【LBM】ナビエ・ストークス方程式の復元（Chapman-Enskog展開）【格子ボルツマン法講座 #6】

はじめて格子ボルツマン法を学ぶ初心者のための講座、第６回になります。

格子ボルツマン法は歴史が浅いこともあり、まだ情報が少ないのが正直なところです。そこで、本サイトでは、ゼロからでも理解できるように丁寧に進めていくようにします。数学的な部分も極力簡単にするので、ぜひ最後までついてきていただけると幸いです。

全く知識のない人でも一応理解できるように説明しますが、可能であればナビエストークス方程式やほかのCFD（流体解析）の手法を知っておくとベターです。

前回は格子ボルツマン方程式の計算フローについて説明しました。下記からどうぞ。

【LBM】計算フロー【格子ボルツマン法講座 #5】

重要なポイントは下記のとおりです。

計算は、「初期化」「衝突」「並進」「マクロ値の更新」で構成される
衝突と並進を分けるため、一時変数を導入する
境界条件は並進の後に入れる

前回は計算フローに関する説明をしました。今回までの講座で格子ボルツマン方程式の計算について説明しましたが、今回は一味変わって証明です。格子ボルツマン方程式がなぜあの形であるのかを理解しましょう。

ナビエ・ストークス方程式の復元（Chapman-Enskog展開）

さて、格子ボルツマン法で流体解析をするとナビエ・ストークス方程式と同等の答えが得られますが、これはなぜでしょうか？

その理由は、格子ボルツマン法は「ナビエ・ストークス方程式と同じ答えが出るような式として作られているから」です。

格子ボルツマン法は元々ボルツマン方程式という気体分子運動論の式からできていますが、このボルツマン方程式に対して「ナビエ・ストークス方程式と一致するような変数設定を行う」ことによってナビエ・ストークス方程式と同等の結果が得られます。

例えば、粘性$\nu$や平衡分布$f^{eq}$などは、ナビエ・ストークス方程式と一致させるために一見複雑な式となっています。

$$ \nu = c_s^2 (\tau – \frac{\Delta t}{2}) $$

$$ f_i^{eq} (x,t) = w_i \rho ( 1 + \frac{u \cdot c_i}{c_s^2} + \frac{ (u \cdot c_i)^2 }{ 2 c_s^4} – \frac{u \cdot u}{ 2 c_s^2}) $$

Chapman-Enskog理論

格子ボルツマン法は、ナビエ・ストークス方程式とほぼ同等の結果が得られますが、厳密に一致するわけではありません。

流れの主要成分が一致するように式を組み立てています。

摂動展開

流れの主要成分を得るための手法として「摂動展開」が使われます。

摂動とは「主要な力による運動が副次的な力（摂動解）によって乱される現象」のことです。これだとわかりにくいので、格子ボルツマン法に置き換えてみましょう。

格子ボルツマン法では、主要な力とは平衡分布$f^{eq}$のことです。そして、副次的な力（摂動解）とは非平衡分布関数$f^{neq} ( = f -f^{eq})$のことです。

つまり、格子ボルツマン法とは平衡分布関数$f^{eq}$というマクロな流れが非平衡分布関数$f^{neq}$というミクロな成分によって乱される現象ということになります。

摂動展開のイメージは掴めたでしょうか？

摂動展開とは、副次的な力（摂動解）を展開したもの、つまり非平衡分布関数$f^{neq}$を展開したものになります。

イメージ的には、テイラー展開のように高次の項に分けることを想像してもらうと良いです。ただ、摂動展開で「クヌーセン数」という値を使うことに注意してください。

クヌーセン数

クヌーセン数は下記で計算されます。元々は気体分子運動論の値で使われる値なので、分子運動を想定しています。

$$ クヌーセン数 = \frac{\lambda ( 平均自由行程 )}{ L ( 代表長さ ) } $$

平均自由行程$\lambda$とは、分子が邪魔されることなく進める距離のことです。代表長さ$L$は流れ場の大きさです。

つまり、クヌーセン数とは「分子がどれだけ自由に動けるかを表す値」です。

クヌーセン数が大きいとスカスカで、分子同士の衝突がほとんどなく、分子と壁の衝突が多い状態です。一方でクヌーセン数が小さいと、分子同士の衝突が多い状態となります。

言うまでもなく、格子ボルツマン法はクヌーセン数が小さい状態を対象としています。

すべての計算店で毎ステップ衝突計算をしていることからも、分子の衝突がたくさん起こると想定しているのは明らかでしょう。

また、大前提として格子ボルツマン法は「クヌーセン数が十分小さい」という仮定のもと成り立っています。

分子が多く、衝突が起こることで、エネルギーが平均化されます。その結果、計算対象を連続体（流体）としてみなすことができるようになります。

そのため、分子の少ない希薄な気体などの高クヌーセン数流れには、基本的な格子ボルツマン法は適用できないので注意しましょう。

摂動展開

クヌーセン数を理解したところで、摂動展開の説明に入っていきましょう。

摂動展開ではクヌーセン数の次数に応じて、分布関数を展開します。

$$ f_i = f_i^{(0)} + \epsilon f_i^{(1)} + \epsilon^2 f_i^{(2)} + O (\epsilon^3) $$

クヌーセン数の次数が高い成分は誤差として、計算上は無視します。

基本的には、高次項は微小なので無視しても問題ありません。しかし、希薄流体のような特殊な流れだとクヌーセン数が大きくなり、高次項の影響が無視できず、問題が起こることもあります。

各項のイメージとしては、分子運動のスケールが異なる＝流れのスケールが異なる、くらいに理解しておくと良いと思います。

$f^{(0)}$という主な流れに対して、副次的な流れである $f^{(1)}$ があり、さらに副次的な流れである $f^{(2)}$ があるといったイメージです。

時間に関しても摂動展開を行います。

$$ \frac{\partial}{\partial t} = \epsilon \frac{\partial}{\partial t_1} + \epsilon^2 \frac{\partial}{\partial t_2} $$

こちらも時間スケールの違いとしてイメージを持つと良いです。

格子ボルツマン方程式のテイラー展開

ここからは、格子ボルツマン方程式をテイラー展開と摂動展開することによって、ナビエ・ストークス方程式を復元する手順について説明します。

本来やるべきなのは、ナビエ・ストークス方程式から格子ボルツマン方程式を作ることなのですが、これはあまりに難しいです。そのため、一般的には格子ボルツマン方程式からナビエ・ストークス方程式を復元して、正しさを証明することがよく行われます。

実用性で言えばあまり重要ではないので、あくまで格子ボルツマン方程式の理解の助けになると思って気楽に見てください。

格子ボルツマン方程式の左辺

格子ボルツマン方程式の左辺をテイラー展開すると下記のようになります。左辺は並進に関する部分です。

$$ f_i (x+e_i \Delta t, t + \Delta t) – f_i (x,t) = \Delta t (e \cdot \nabla + \frac{\partial }{\partial t}) f_i + \frac{\Delta t^2}{2} (e \cdot \nabla + \frac{\partial }{\partial t})^2 f_i + O $$

三次以降の項は誤差として無視しています。

上式に対して摂動展開を行います。

$$ \frac{\partial}{\partial t} = \epsilon \frac{\partial}{\partial t_1} + \epsilon^2 \frac{\partial}{\partial t_2} $$

$$ f_i = f_i^{(0)} + \epsilon f_i^{(1)} + \epsilon^2 f_i^{(2)} + O (\epsilon^3) $$

すると、格子ボルツマン方程式の左辺（LHS）は下記のようになります。

$$ LHS = \Delta t ( e \cdot \Delta + \epsilon \frac{\partial }{\partial t_1} + \epsilon^2 \frac{\partial }{\partial t_2} ) ( f_i^{(0)} + f_i^{(1)} + f_i^{(2)} ) + \frac{\Delta t^2}{2} ( e \cdot \Delta + \epsilon \frac{\partial }{\partial t_1} + \epsilon^2 \frac{\partial }{\partial t_2} )^2 ( f_i^{(0)} + f_i^{(1)} + f_i^{(2)} ) $$

上の式を展開すると、$\epsilon^3$以上の係数の項が出てきますが、ナビエ・ストークス方程式の復元には $\epsilon^2$ までしか必要ないので、全て消してOKです。

格子ボルツマン方程式の右辺

格子ボルツマン方程式の右辺も摂動展開します。右辺は衝突に関する部分で、今回はBGKモデルを使用します。

$$ – \frac{1}{\tau} [ f_i – f_i^{eq}] \Delta t F_i = – \frac{1}{\tau} [ f_i^{(0)} + \epsilon f_i^{(1)} + \epsilon f_i^{(2)} – f_i^{0}] \Delta t \epsilon F_i $$

平衡分布関数$f^{eq}$は$f^{(0)}$と同じ意味なので、置き換えられています。

今回は参考にした下記論文と合わせるため、外力項ありで進めます。

https://journals.aps.org/pre/abstract/10.1103/PhysRevE.65.046308

格子ボルツマン方程式の摂動展開

両辺を摂動展開したので、（展開後の左辺）＝（展開後の右辺）の形にして、クヌーセン数の次数で整理します。

左辺の展開が非常に大変ですが、最終的に下記の３つの式が得られます。（$\epsilon^3$以上の項は使用しないので消すと、展開後の項の整理がかなり楽になるのでオススメです。）

クヌーセン数ゼロ次　$O(\epsilon^0)$

$$ f_i^{(0)} = f_i^{eq} $$

クヌーセン数一次　$O(\epsilon^1)$

$$ D_1 f_i^{(0)} = – \frac{1}{\tau \Delta t} f_i^{(1)} + F_{1i} $$

クヌーセン数二次　$O(\epsilon^2)$

$$ \frac{\partial f_i^{(0)}}{\partial t_2} + ( 1 – \frac{1}{2 \tau}) D_{1i} f_i^{(1)} = – \frac{1}{\tau \Delta t} f_i^{(2)} – \frac{\Delta t}{2} D_{1i} F_{1i} $$

ゼロ次は平衡分布に当たるので、すぐに分かります。一次と二次に含まれる$D_{1i}$は、下記で定義される物質微分です。

$$ D_{1i} = \frac{\partial }{ \partial t} + e_i \nabla_i $$

物質微分はラグランジュ微分とも呼ばれ、流体が移動する座標系から見たときの時間変化を表します。例えば、粒子法では計算点が移動するので、物質微分で表現されます。

上の三つの式に対して、モーメントを計算することで、ナビエ・ストークス方程式が復元できます。モーメントとは、粒子速度による積分で下記式のように今まで密度を求める場合などに行っていたものです。

$$ \rho = \sum f_i$$

直感的な理解としては、マクロ値への変換と思ってください。

下記の式は、クヌーセン数一次の項に対してモーメントを計算することで得られます。

$$ \frac{\partial \rho}{\partial t_1} + \nabla_1 \cdot (\rho u^* ) = A_1 $$

$$ \frac{ \partial (\rho u^*)}{ \partial t_1} + \nabla_1 \cdot \Pi^{(0)} = (n+\frac{m}{\tau}) F_1$$

一段目の式は、クヌーセン数が一次の式のモーメントを計算したものです。二段目の式は、クヌーセン数が二次の式に粒子速度$e_i$をかけて、モーメントを計算したものです。

上の２つの式は、ナビエ・ストークス方程式の連続の式と輸送式に非常によく似ていることがわかります。

ここで、ナビエ・ストークス方程式に近づけるために、下記を与えます。

$$A = 0$$

$$ n + \frac{m}{\tau} = 1 $$

次に、クヌーセン数二次の式対して、モーメントを計算します。

$$ \frac{\partial \rho}{\partial t_2} = \Delta t (m – \frac{1}{2}) \nabla_1 \cdot F_1$$

$$ \frac{\partial (\rho u^*)}{\partial t_2} = \Delta t (m – \frac{1}{2}) \frac{\partial F_1}{\partial t_1} + \nabla_1 \cdot \sigma_1$$

上の段の式は、クヌーセン数が二次の式のモーメントを計算したものです。下の段の式は、クヌーセン数が二次の式に粒子速度$e_i$をかけてモーメントを計算したものです。

これらの式に対しても、書きを与えてナビエ・ストークス方程式に近づけます。

$$ m = \frac{1}{2}$$

ストレステンソル$\sigma_1$は下記で与えられます。

$$\sigma_1 = – (1 – \frac{1}{2 \tau}) \nabla_1 \cdot \Pi^{(1)} – \frac{\Delta t}{4} (C_{1\beta \alpha}) $$

$$ = (\tau – \frac{1}{2}) c_s^2 \Delta t \rho (\nabla_{\alpha} u^*_{\beta} + \nabla_{\beta} u^*_{\alpha} ) + \Delta t [ (\tau – \frac{1}{2}) ( u^*_{\alpha} F_{1 \beta} + u^*_{\beta} F_{1 \alpha} ) – \frac{\tau}{2} (C_{1 \alpha \beta} + C_{1 \beta \alpha} ) ]$$

ストレステンソルもナビエ・ストークス方程式と合わせる必要があります。そのために、右辺第一項の係数を粘性係数とし、右辺第二項を消去するように変数を設定します。

$$ \nu = (\tau – \frac{1}{2}) c_s^2 \Delta t $$

$$ C = (1 – \frac{1}{2 \tau}) 2 u^* F$$

以上のように変数を設定することで、ナビエ・ストークス方程式が得られます。

$$\frac{\partial \rho}{\partial t} = \nabla \cdot \rho v = 0$$

$$\frac {\partial}{\partial t} (\rho v) + \nabla \cdot (\rho v v)= – \nabla p + \nu \nabla \cdot [\rho (\nabla v + ( \nabla v )^T)] +F$$

圧力と粘性係数は下記で得られます。

$$ p = \frac{\rho c_s^3}{3} $$

$$ \nu = \frac{1}{3}(\tau – \frac{1}{2}) c_s^2 \Delta t $$