今日から「DSMC ソルバーをイチから作ろう！」と題して、DSMC ソルバーを作ることをやっていこうと思います。
僕も勉強しながらなので、今回は DSMC とはなんぞや、というところから解説していこうと思います。

一応、大学の研究で dsmcFoam+ や SPARTA と呼ばれる OSS の DSMC ソルバーを使って研究していたのですが、イチから作ったことは無いので、卒業した今、作ろうとしている次第です。

ちなみに、大学の卒業論文は Github に置いてありますので、興味のある方は眺めてみてください。
(130 ページ、 150 式以上というボリュームなので、ちゃんと読もうとしないほうがいいです😅)

なんとなく以下のようなスケジュールで進めようと思っていますが、都度変更しながら全体としてなるべく分かりやすいように作っていきます。

DSMC って？ ← 今ここ
DSMC の支配方程式
DSMC の導入
分子モデル
分子間衝突
多次元化してみる
ソルバーをコーディングして作る (多分複数回になる)

また、対象読者としてはなんとなく数値流体力学を知ってる方を対象としています。

気圧がとてもとても低い気体、希薄流

気体 (空気等) の圧力が低下する (≒ 分子の密度が減る) と、気体は連続体とみなすことができなくなります。
これを「希薄流 (きはくりゅう)」と言います。

どんなシチュエーションで希薄流になるかというと、

高高度大気
ミクロスケールな半導体製造装置内の流れ
薄膜
真空

等です。まとめると、

宇宙のように、そもそも気体分子が少ない
大気圧化でも、対象のスケールがとても小さい

のようなときです。

Navier–Stokes 方程式を数値的に解く手法は様々ありますが、どれも「連続体近似」が成り立っているという前提の元解かれます。
しかし、希薄流では連続体とみなせず、連続体近似が破綻してしまうので、その前提が成り立たず、同様の方法では Navier–Stokes 方程式を解くことができません。

今回取り扱う、DSMC とは、このような Navier–Stokes 解析が上手く行かない希薄流の解析に効果を発揮する手法になります。

希薄流になる条件は？

希薄流という概念を導入しました。とにかく薄い気体です。
では、基準はあるのか？というと、あります。
基準になる量が Knudsen 数 (クヌッセン数、クヌーセン数) という量です。

ほんとにざっくりな説明だと「気体のスカスカ具合」を表した量で、Knudsen 数が大きいほどスカスカ = 希薄ということになります。
Knudsen 数は $Kn$ と書いたりしますが、

$$ Kn > 0.01 $$

程度で希薄流とみなす場合が多いです。
逆に、 $ Kn \ll 1 $ の場合は連続体とみなせます。

この Knudsen 数については、論文をいくつか読んでまとめた記事を以前自身のブログで書いたので、詳しくはそちらを参照してください。
→ クヌッセン数について日本一詳しく解説します！

(このクヌッセン数にも種類があり、しっかり理解するのは大変かもしれませんが、そんな方は、以下の論文をまず読むことをおすすめします。: Boyd, I. D., Chen, G., and Candler, G. V., “Predicting Failure of the Continuum Fluid Equations in Transitional Hypersonic Flows”, Phys. Fluids, 7, pp. 210–219, 1995.)

Knudsen 数の定義

この Knudsen 数 $Kn$ は平均自由行程 $\lambda$ と流れの代表長さ $L$ の比で、

$$ Kn = \dfrac{\lambda}{L} $$

と表されます。
いきなり平均自由行程が出てきましたが、これは「分子同士が衝突して、次に別の分子と衝突するまでの距離」を全分子で平均化したものです。
自由に 進む距離 (行程) を平均したものなので平均自由行程と呼びます。
なぜ、流れの代表長さとの比を出すかというと、相対的に見てスカスカかどうかということが重要だからです。

砂場を例に挙げると、人間から見ればギチギチに砂が詰まっていますが、微生物にとって見れば間を通れるほどスカスカと見ることができます。
なので、代表長さとの比で表す必要があるわけです。

平均自由行程は、様々な形で表せますが、分子直径 $ d $, 気体の数密度 $n$ を用いて、

$$ \lambda = \dfrac{1}{\sqrt{2}\pi d^2n} $$

と表されます。
気体の数密度 $n$ は圧力 $p$ と温度 $T$ が与えられれば、状態方程式 (懐かしい) によって、

$$ p = nk_{\mathrm{B}}T $$

から求めることができます。
$k_\mathrm{B}$ はボルツマン方程式で、$k_\mathrm{B} = 1.38 \times 10^{-23}~\mathrm{J/K}$ です。

DSMC って？

DSMC とは、 Direct Simulation Monte Carlo 法の略で、日本語ではモンテカルロ直接法などと呼ばれます。
モンテカルロとあるように、確率統計的な手法を用いて希薄流の流体場を解析する手法です。

どの部分を確立的に解くかというと、分子同士の衝突です。
細かくはこの記事で書きませんが、サンプル粒子というものを導入します。
これはいくつかの分子を代表した粒子になっていて、例えば 1000 個の粒子を 1 つの粒子で代表させて、そのサンプル粒子だけ追跡し計算します。
代表した粒子のみを解くので、全分子を追跡して解析するよりも計算コストが低く抑えることができます。

計算の流れ

基本的には以下の流れになっています。

流れ場のセルを分割する
セルに分子を配置する (初期条件設定)
分子を $\Delta t$ の間動かす (+ 境界処理)
分子間衝突を模擬する
各セルで物理量を計算する
流れが定常か確認: Yes → 次ステップへ, No → 3. へ戻る
物理量の時間平均を計算する
ゆらぎを判定: Yes → 次ステップへ, No → 3. へ戻る
終了

分子の動かし方、すなわち $(\boldsymbol{x}, \boldsymbol{c})$ の与え方は、簡単なものでよく、マクスウェル分布によって与えたりします。

参考文献

日本機械学会, 『原子・分子モデルを用いる数値シミュレーション』, コロナ社, 1996.
林大地, 東北大学卒業論文, 2020.

全体的に、こちらの文献を参考にしています。