#概要
格子ボルツマン法(LBM)の基礎について，物理の視点からお話します．
最初に，よくある誤解と論文では触れられない基礎的な前提について述べます．
最も基本的なD2Q9モデルにおけるNavier-Stokes方程式の復元の一部を具体的に追いかけ，LBMの基礎的な部分の理解を目指します．高度な内容には踏み込みませんが，その足掛かりとなる内容だと思います．
###この記事に関する注意
1:Guoらの論文[2]を読破している方にとっては退屈な記事です．
2:全ての計算はフォローしません．
3:僕の思想が所々に強く入っています．
4:2つのベクトル$\vec{V}$を$\vec{V}\vec{V}$のようにくっつけて2階のテンソルと表記する場合があります．
5:テンソル同士を潰すような縮約はコロンを用います$T_{ij} U_{ij} = \boldsymbol{T}:\boldsymbol{U}$
6:僕がMarkdownに不慣れなため，ボールド体をテンソルとするなど表記が揺れる部分もありますが，文脈から推察してください．
#よくある誤解に関する注意
LBMがトリッキーな手法であるために誤解されやすいので，簡単に触れておきます．
###注意1：LBMは粒子法ではない
仮想「粒子」の分布関数を計算するために粒子法と誤解されることがあります．
LBMは「lattice」とあるように(相空間中に)格子を生成してから話が始まりますし，離散化された場を考えるので無論Lagrange描像ではないです．格子ボルツマン方程式(LBE)はボルツマン方程式を離散化したものであることも思い出しましょう．
詳細は後述しますが，そもそも仮想粒子は物理としては何の意味も持ちません．
###注意2：D2Q9モデルは厳密には等温モデルではない
「D2Q9モデルは自動的に非圧縮流体になる」という表現を時折見かけますが，「D2Q9は運動量空間の自由度が9つしかないためそもそも圧縮性のコントロールも温度の表現もできないが，(自由度が足りないが故に)実効的に非圧縮性流体として振る舞う」という表現の方が正確です．
論文のイントロとかでは「実在流体の僅かな圧縮性も再現できる！」みたいな書き方がされることがありますが，個人的にはこの表現は嫌いです．できないことをできると言うな
蔦原氏も[1]で「非熱モデルと呼ぶべき」と仰っています．
面倒な言い回しをしましたが，D2Q9の計算を全部やればわかります．（本末転倒）
#LBMの思想
具体的な計算に入る前に，皆さんに僕の思想を擦り込みます．
長々と書いてしまいましたが，
「流体は高い普遍性をもつ．LBMではその性質と統計物理の枠組みを悪用する」
ということを伝えたいだけです．
##統計物理に関して
粗視化は，統計力学を勉強したことがある人なら馴染みのある概念だと思います．
これを抑えられるとメインの計算の理解が捗るので触れておきます．
物理には階層性があります．というのは，それぞれの長さや時間スケールによって支配的な法則が異なります．(ミクロなら量子力学，人間のスケールなら古典力学や熱力学，宇宙スケールなら一般相対論，など)
より大きなスケールの系は，それを構成する小さいスケールの系の集団として原理的には記述できるはずです．しかし，実際に観測されるマクロな振る舞いはそれ自体が一つの実体を持つように見え，ミクロな世界は多くの場合に顔を出しません．また，マクロな系の記述に必要な自由度は，それを構成するミクロな系のそれよりも少なくなります．
ですので，マクロな系にミクロな意味づけを与える際は，ミクロな情報を程よく「塗り潰し」て自由度を大幅に削る必要があります．
量子力学と熱力学を橋渡しする統計力学の場合，$10^{23}$オーダーの自由度を持つミクロな分子の集団を，僅か6つの量(温度，内部エネルギー，圧力，体積，化学ポテンシャル，粒子数)で記述される熱力学と無矛盾になるように自由度を塗りつぶします．
粗視化を行うと，当然ですが元のミクロな情報の大部分は失われます．つまり，一般に粗視化の前後で全然違う系に変身してしまいます．上の例では，統計力学を通じて量子力学から熱力学っぽいものが捻り出されています．これは現代科学では広く受け入れられているものの，よくよく考えてみると全く非自明なプロセスです．粗視化の適当な手法は問題設定によって様々です．
我々の身の回りに流体が溢れているのは，粒子の集団を粗視化すると容易に流体になってしまうためです．
つまり流体は物理学的なオブジェクトの中でも非常に高い普遍性を持ち，このことは**「粒子の集団が巨視的に流体として振る舞う条件は結構ゆるい」**とも表現できます．
##LBMって何をしてるの？
統計物理の思想と対比しつつLBMの思想について述べます．
###Boltzmann方程式と連続体方程式の関係
気体分子運動論において，Boltzmann方程式が希薄な気体を記述します．これは分子の分布関数に関する偏微分方程式です．
Chapman-Enskog理論によると，系のKnudsen数が十分小さいならばBoltzmann方程式の摂動解から連続体方程式を導くことができます．
摂動は局所平衡分布(Maxwell分布)からのずれです．
巨視的な物理量は分布関数のモーメントとして求められます．このモーメントの計算(速度に関する積分)は相空間の自由度を落とす，ある種の粗視化です．
粗視化によって得られるストレステンソルの形は衝突項の具体形に強く依存します．
(衝突項の具体形が何であれ，Boltzmann方程式がGalilei不変なNewton力学を参照しているので連続体方程式の慣性項はちゃんと$(\partial_t +\vec{v}\cdot\nabla)\vec{v}$になります)
一旦まとめると，ロジックの流れとしては，「Boltzmann方程式→摂動展開&粗視化→連続体方程式」です．
###逆転の発想
Boltzmann方程式は希薄な気体の方程式だから液体には利用できないのでしょうか？
いいえ，どこかのクレイジーな天才は，上の話を逆走するような逆転の発想をします．
「ええ感じのNavier-Stokes方程式が出てくるような衝突項をもつ集団を作ってしまえ」
「運動量空間の自由度が9つだけだと数値計算が楽だよね」
粗視化の性質が悪用されています．と言うのは，ありもしない希薄気体，仮想粒子の集団を自分勝手に作っているのですが，どうせ粗視化してしまうので正当化されます．
現実とズレるならそれを打ち消すような項を好きな形で付け足せばいいんです．
運動量空間の自由度も，流体の等方性を守り，なおかつ必要な巨視的な量を表現できる最低限の数で良いのです．
このように，統計物理の手法の応用ではあるものの，ロジックが全然違います．
しつこくおさらいをしますが，統計物理では，
「ミクロとマクロな系が先にあって，それらの橋渡しをうまいことやる」
のですが，LBMでは
「欲しいマクロな系と橋渡しの手法が先にあって，それらに辻褄が合うような「LBEに従うミクロな系」を強引に用意し，そのミクロな系を数値的に解く」
となっています．
多分ここが初学者にとってのネックなのではないかと思います．(1年前は僕も初学者でしたけどよく分かんないです)
無理矢理作ったミクロな系ではアルゴリズムも単純で境界条件も相対的にラクです．(というかそうなるように作ってあります)
仮想粒子の集団は，粗視化という「橋渡し」ありきの嘘みたいなものです．
局所平衡分布関数の形も，欲しい答えが出てくるように決めてしまいます．
最終的に辻褄が合えばなんでも良いので，再現したい系や求める精度に応じて適切な手法は変わります．「橋渡し」を少し変えると分布関数もデザインしなおす必要がありますし，希薄気体効果の抑え方や外力の加え方も何通りもあります．
今日もLBMの手法の開発に関する論文が出続けるのは，上記のようなLBMのヘンテコな性質が理由だと思います．
#Navier-Stokes方程式の復元
実際の手続きでは，摂動の1次のゼロ次と1次のモーメントで密度保存とEuler方程式，摂動の2次のゼロ次と1次のモーメントで(より遅い時間スケールでの)密度保存とNS方程式の拡散項が出ます．
この記事ではEuler方程式の計算のみを扱います．
##LBE
D2Q9モデルを導入します．BGK近似(単衝突)でのLBEは次式である．
$$
f_i(\vec{x} +\vec{e}_i \Delta t, t+\Delta t) -f_i(\vec{x}, t)
= -\frac{\Delta t}{\tau}\left[f_i(\vec{x}, t) -f_i^{eq}\right] +\Delta tF_i
$$
時間と空間はそれぞれ$\Delta t, \Delta x$で離散化され，音速は$c = \Delta x/\Delta t,c_s = c/\sqrt{3}$．
$\vec{e}_0 = (0, 0)$, $i=1,...,4$の時は最近接格子点への移動$\vec{e}_i = c(\cos[\pi(i-1)/2], \sin[\pi(i-1)/2])$,
$i=5,...,8$の時は次近接格子点への移動$\vec{e}_i = \sqrt{2}c(\cos[\pi(i-9/2)/2],\sin[\pi(i-9/2)/2])$，
$f_i(\vec{x}, t)$は，仮想粒子が時刻$t$，座標$\vec{x}$において速度$\vec{e}_i$である確率分布を表す．
巨視的な物理量は，分布関数のモーメントとして求められる．
$$
\rho = \sum_i f_i, \ \ \rho \vec{u}^* = \sum_i\vec{e}_if_i +m\Delta t \vec{\rm F},\ \
\frac{1}{2}\rho\vec{u}^2 +e = \sum_i\frac{1}{2}f_i\vec{e}_i^2
$$
$\rho\vec{u} = \sum_i f_i\vec{e}_i $であり，$e$は内部エネルギー密度，$\vec{\rm F}$は外力，定数$m$は後のマルチスケール解析で決める．
$\vec{u} $の「*」の有無は外力の印加の有無を表す．
$f_i^{eq}$は局所平衡分布関数で，Maxwell分布を$\mathcal{O}(u^2)$までで打ち切ったものである．
$$
f_i^{eq} = w_i\rho\left[1 +\frac{\vec{e}_i\cdot\vec{u}}{c_s^2} +\frac{\vec{u}\vec{u}:(\vec{e}_i\vec{e}_i-c_s^2\boldsymbol{I})}{2c_s^4}\right]
$$
$w_i$は$w_0 = 4/9$,$i=1,..,4$の時$w_i=1/9$, $i=5,...,8$の時$w_i = 1/36$である．
$F_i$はforcing termと呼ばれ，外力の印加に起因するspurious flowを打ち消すためにあらかじめ仕込む項である．外力の$\vec{\rm F}$と混同しないように注意．(spurious flowは2次の摂動で粘性応力を計算するときに出るもので，$\vec{u}\vec{\rm F} +\vec{\rm F}\vec{u}$という格好をしてます．)
ボルツマン方程式における$\vec{\rm F}\cdot \frac{\partial f}{\partial \vec{p} }$という項に対応するこの$F_i$は，スカラー$A$，ベクトル$\vec{B}$，テンソル$\boldsymbol{C}$を用いて次の$\vec{e}_i$の級数で表される．$A,\vec{B},\boldsymbol{C}$は$\vec{\rm F}$の関数として後で決定する．
$$
F_i = w_i\left[A +\frac{\vec{e}_i\cdot\vec{B}}{c_s^2} +\frac{\boldsymbol{C}:(\vec{e}_i\vec{e}_i-c_s^2\boldsymbol{I})}{2c_s^4}\right]
$$
##Chapman-Enskog理論，マルチスケール解析
マルチスケール解析によってLBEの摂動解を2次まで求めます．1，2次それぞれの摂動解のモーメントから移流と拡散が出ます．減衰振動での振動と減衰に対応します．
摂動の計算では，分布関数以外は下付き添字の1,2で摂動の次数を表します．分布関数の添字との混同に注意．
###準備
LBEの左辺を$\Delta t $に関してTaylor展開し，2次で打ち切る．
$$f_i(\vec{x}+\Delta t\vec{e}_i, t+\Delta t) -f_i = \Delta t \left(\frac{\partial }{\partial t} + \vec{e}\cdot\nabla\right) f_i +\frac{\Delta t^2}{2}\left(\frac{\partial }{\partial t} + \vec{e}\cdot\nabla\right)^2f_i$$
Knudsen数を摂動パラメータ$\epsilon$とし，分布関数と演算子を展開する．外力は1次の摂動として扱う．$\vec{\rm F} = \epsilon \vec{\rm F}_1 $分布関数は2次で打ち切る．
$$
f_i = f_i^{(0)} +\epsilon f_i^{(1)} +\epsilon^2 f_i^{(2)} +\mathcal{O}(\epsilon^3)
$$
巨視的な時間スケールは，仮想粒子の時間スケールより遅い$t_1$と，さらに遅い時間$t_2$の2つを用いる．この2つの時間は独立変数として扱うことに注意．空間は1段階だけ粗視化を行う． $F_i $は外力と同様に摂動の1次であると仮定する．
$$
\frac{\partial}{\partial t} = \epsilon\frac{\partial}{\partial t_1} +\epsilon^2 \frac{\partial}{\partial t_2}, \ \nabla = \epsilon\nabla_1
$$
以上の$\epsilon $で展開した分布関数と演算子をLBEに代入し，$\epsilon $の各次の項を比較する．すると，分布関数の摂動解を逐次求めることができる．ゼロ次は次式．
$$\mathcal{O}(\epsilon^0): \ f_i^{(0)} = f_i^{eq} $$
この$f_i^{0}$のモーメントを具体的に計算してみると，巨視的な変数(密度と運動量，エネルギー密度)がきちんと出てくることがわかる．
###摂動の1次
2次まで計算しないとNS方程式は出ないのですが，今回は割愛して1次での計算の一部だけを追いかけます．
摂動の1次の演算子を導入する．
$$D_{1i} = \frac{\partial }{\partial t_1} + \vec{e}\cdot\nabla_1$$
LBEと$A=\epsilon A_1,\vec{B} = \epsilon \vec{B}_1, \boldsymbol{C} = \epsilon \boldsymbol{C}_1$より，分布関数の1次の摂動解は次式より求めることができる．
$$\mathcal{O}(\epsilon^1): \ D_{1i}f_i^{(0)} = -\frac{1}{\tau}f_i^{(1)} +F_{1i}$$
この式に$\vec{e}_i$を掛けて$i$に関する和を取る．
外力は摂動の1次であること，$\vec{u}$は摂動のゼロ次であることより，$\sum_i \vec{e}_i f_i^{(1)} = -m\Delta t\vec{\rm F} $及び$\sum_i \vec{e}_i f_i^{(r>1)} = 0$が要請される．
更に係数$n $を用いて$\vec{B}_1 = n\vec{\rm F}_1 $と書けると仮定する．$\sum_i\vec{e}_i\vec{e}_if_i^{eq} $の計算に注意すると，次式が得られる．
$$
\sum_i D_{1i}f_i^{(0)} = \frac{\partial }{\partial t_1 }(\rho \vec{u}) +\nabla_1\cdot\left(c_s^2\rho\boldsymbol{I} +\rho\vec{u}\vec{u} \right) = \left(n +\frac{m\Delta t}{\tau }\right)\vec{\rm F}_1
$$
この式がEuler方程式になるには，$n +m\Delta t/\tau = 1$が満足されなければならない．
2次の摂動まで計算すると$m = 1/2$となるため，Euler方程式を得るには$n = 1 -\Delta t/(2\tau) $と選ぶ必要がある．
#終わりに
Twitterで「LBMの記事書こうかな〜(チラッチラッ」と言ったら予想以上に反応があったので書きました．僕もLBMの基礎でかなり苦しんだので，誰かの助けになればいいかなあと思います．
全ての計算は追えませんでしたが，誰でも自力で計算できるようなるための取っ掛かりににはなると思います．
回収されずに放り投げられた計算が気になる人は[2]を読んで計算全部フォローしようね❤️
質問や間違いがあればご指摘ください．
###参考文献
[1]蔦原道久， 格子ボルツマン法・差分格子ボルツマン法.
FDLBMについても書いてあります．おすすめの一冊です．
[2]Z. Guo, et.al. Phys. Rev. E, 65, 046308 (2002).
これの計算を全部追うことができれば基礎はばっちりです．