この文書では、WSL2のUbuntu24.04上でLBM(格子ボルツマン法)流体解析ツールXLBを活用する方法を紹介します。
この書き出して、Autodeskが開発している「XLB(Accelerated Lattice Boltzmann)」について、以下10個のQiita記事をまとめてきました。GPUによる高速計算と機械学習との連携を重視したPython製LBM(格子ボルツマン法)ライブラリとして、2次元・3次元の流体解析に対応し、標準的なBGKモデルや高レイノルズ数流れに適したKBCモデルなどを利用できます。
LBM流体解析ツールXLBの活用【01:JAXcpu導入方法と動作確認】
LBM流体解析ツールXLBの活用【02:例題lid_driven_cavity_2d.pyの分析】
LBM流体解析ツールXLBの活用【03:JAXgpu導入方法と動作確認】
LBM流体解析ツールXLBの活用【04:例題flow_past_sphere_3d.pyの分析】
LBM流体解析ツールXLBの活用【05:Neon導入方法とマルチ解像度例題の分析】
LBM流体解析ツールXLBの活用【06:例題multires_windtunnel_3d.pyの分析】
LBM流体解析ツールXLBの活用【07:Warp導入方法と動作確認】
LBM流体解析ツールXLBの活用【08:例題rotating_sphere_3d.pyの分析】
LBM流体解析ツールXLBの活用【09:XLB付属例題概要と分析】
LBM流体解析ツールXLBの活用【10:XLB独自例題の作成と分析】
以上の記事で、それなりに、XLBの活用ノウハウは獲得できたと思います。ただそれでも、そもそもの「格子ボルツマン法(LBM)の基本」については、十分に解説してなかったと思い、ここでは改めて 【00:LBMの基礎の直観的理解】 を目的とします。
この記事でLBMを学ぶ12のステップ
【前半:00A】
STEP 0 流体を知る
STEP 1 拡散方程式を知る
STEP 2 手で計算する
STEP 3 Excelで計算する
STEP 4 Pythonで計算する
STEP 5 格子と時間発展を理解する
STEP 6 分布関数 f を知る
【後半:00B】
STEP 7 Collisionを知る
STEP 8 Streamingを知る
STEP 9 D1Q3でLBMを体験する
STEP 10 D2Q9で2次元化する
STEP 11 境界条件を入れる
STEP 12 実際の流体を解析する
この構成のポイントは、STEP 1~5で「数値計算そのもの」を理解してから、STEP 6以降でLBM特有の考え方に入ることです。これなら、LBMを初めて学ぶ学生にもかなり説明しやすい構成になります。つまり、この構成なら、**「数学 → 手計算 → Excel → Python → LBM → 実際の流体」**という流れが明確で、初心者向けの教材として非常に分かりやすいです。
LBM初心者向け学習ロードマップ
もう少し詳しい【00:LBMの基礎の直観的理解】に関する目次を、以下にまとめます。ただし、解説に図解を入れて、詳しく解説するのですが、1記事のファイル制限100MBのために、この【00:LBMの基礎の直観的理解】は、前半【00A:LBMの基礎の直観的理解】+後半【00B:LBMの基礎の直観的理解】に分割します。
前半【00A:LBMの基礎の直観的理解】 次の記事です
STEP 0 流体の基礎
- 速度・圧力・密度・粘性など、流体を表す基本的な物理量を理解する。
- 流体が「流れる」「拡散する」「壁にぶつかる」といった基本的な現象を理解する。
- Reynolds数を学び、流れの状態と速度・粘性・代表長さの関係を理解する。
STEP 1 1次元拡散方程式
- まず流体ではなく、熱や物質が広がる拡散現象を題材にする。
- 1次元拡散方程式 $\frac{\partial T}{\partial t} = \alpha\frac{\partial^2T}{\partial x^2}$ の意味を理解する。
- 「ある場所の値が、時間とともに左右へ広がっていく」という時間発展を理解する。
STEP 2 手計算
- 空間をいくつかの格子点に分割し、それぞれの値を計算する。
- 1ステップずつ計算して、現在の値 → 次の時間の値という流れを確認する。
- 実際に数値を計算して、拡散によって値が「平らになっていく」ことを確認する。
STEP 3 Excel
- 手計算した拡散計算をExcelの表計算に置き換える。
- 時間ステップを増やし、値が時間とともに変化する様子を確認する。
- グラフを作り、拡散する様子を目で確認する。
STEP 4 Pythonで拡散方程式
- Excelで行った計算を、Pythonで同じように再現する。
-
for文や配列を使って、多数の格子点・時間ステップを計算する。 - Matplotlibなどを使い、温度分布の変化をグラフや動画で表示する。
STEP 5 「格子上で時間発展」の理解
- 解析領域を格子に分割し、格子点ごとの値を更新する考え方を理解する。
- 「現在の状態 → 計算 → 次の状態」という時間発展計算を理解する。
- 拡散方程式とLBMに共通する、格子上で局所的な計算を繰り返す考え方を理解する。
STEP 6 LBMの分布関数 f
- LBMでは温度や速度を直接計算するのではなく、**分布関数
f**を計算することを理解する。 - 分布関数が「各方向へ進もうとしている流体の情報」を表していることを理解する。
- 分布関数から、密度や速度を計算できることを学ぶ。
後半【00B:LBMの基礎の直観的理解】 この記事です
STEP 7 Collision
- 分布関数が局所的に変化する処理としてCollision(衝突)を理解する。
- 分布関数を平衡分布関数
f_eqに近づけるという考え方を学ぶ。 - BGKモデルなどを使って、Collisionの計算を手計算・Pythonで確認する。
STEP 8 Streaming
- 分布関数が隣の格子点へ移動する処理を理解する。
- 「Collision → Streaming」というLBMの基本的な計算順序を理解する。
- 格子上で分布関数が移動する様子を、簡単な図や数値例で確認する。
STEP 9 D1Q3
- **1次元・3速度方向(D1Q3)**の非常に単純なLBMを学ぶ。
-
←、停止、→の3方向の分布関数を使って計算する。 - 手計算やExcelを使って、分布関数 → Collision → Streamingを実際に確認する。
STEP 10 D2Q9
- D1Q3を2次元に拡張し、D2Q9を理解する。
- 1つの格子点にある9方向の分布関数を理解する。
- 2次元の流れを計算し、速度・密度・流れ場を可視化する。
STEP 11 境界条件
- 壁・入口・出口など、流体解析で必要となる境界条件を理解する。
- LBMで代表的なBounce Backによる壁境界を学ぶ。
- 境界条件によって、流れの速度分布や圧力などがどのように変化するか確認する。
STEP 12 実際の流体
- Poiseuille流、円柱周りの流れ、障害物周りの流れなどを計算する。
- Reynolds数を変化させ、層流から複雑な流れへの変化を確認する。
- 最終的に、揚力・抗力・渦放出など、実際の流体解析で扱う現象へ発展させる。
■ LBMの基礎の直観的理解【STEP 7 Collisionを知る】
STEP 7の Collision(衝突) は、STEP 6で学んだ「方向ごとの分布関数 $f$」が、時間とともにどのように変化するのかを説明する段階です。ここでは、まだ「速度モデルD1Q3やD2Q9」の具体的な計算には深く入りません。まず、Collision = 1つの格子点の中で、分布関数の偏りを調整する処理 と理解するのが重要です。
1. STEP 6からのつながり
STEP 6では、「1つの格子点 → 方向ごとの分布関数:f₀, f₁, f₂, ... → 密度 ρ・速度 u」と考えました。しかし、時間が進めば分布関数 $f$ も変化します。そこで、「fをどうやって次の状態に変えるのか?」 という問題が出てきます。その最初の処理が、Collision(衝突) です。
2. Collisionは「衝突」という名前
LBMでは、流体の微視的な運動を直接計算する代わりに、流体粒子同士が衝突して、進む方向の分布が変化するという効果をモデル化します。実際の分子を1個ずつ追跡するわけではありません。イメージとしては、方向ごとの分布が変化すると考えます。
3. Collisionの重要な特徴
Collisionには重要な特徴があります。それは、「その格子点の中だけで計算する」 ことです。例えば、Collisionを計算するとき、「AのCollision → Aの分布関数だけを見る・BのCollision → Bの分布関数だけを見る」となります。つまり、隣の格子点を見るのはCollisionではありません。この点は、後のStreamingとの大きな違いです。
4. Collisionでは何を目指すのか
Collisionでは、現在の分布関数を、平衡分布関数 $f_i^{eq}$ に近づけます。「平衡」というのは、その格子点の密度と速度に対応した、バランスの取れた分布 と考えてください。ただし、全部を一気に平衡にするわけではありません。少しずつ近づけます。
5. Collisionの基本式
代表的なBGKモデルでは、$f_i^{new}=f_i-\omega\left(f_i-f_i^{eq}\right)$を使います。この式は、初心者には次のように読むとよいです。「新しいf = 現在のf + 「平衡状態との差」を少し修正」あるいは、「現在の分布 $f$ を、平衡分布 $f^{eq}$ に近づける」 という式です。
LBMのCollisionモデル:BGKモデルとその他のモデル
LBM(格子ボルツマン法)では、各格子点にある分布関数 $f_i$ を、Collision(衝突)によって更新します。このとき、分布関数をどのような方法で変化させるかを決めるのが、Collisionモデルです。
【1. BGKモデルとは?】
BGKモデルは、LBMで広く使われる、最も基本的なCollisionモデルです。現在の分布関数 $f_i$ を、平衡分布関数 $f_i^{eq}$ に向けて、少しずつ近づけます。$ f_i^{*}=f_i-\frac{1}{\tau}(f_i-f_i^{eq}) $
式の意味としては、以下の通りです。
- $f_i$:衝突前の現在の分布関数
- $f_i^{eq}$:その格子点の密度・速度に対応する平衡分布関数
- $\tau$:緩和時間。分布が平衡に近づく速さを決めるパラメータ
- $f_i^*$:Collision後の分布関数
直観的には、方向ごとの分布の偏りを、一定のルールで少しずつ整える方法です。Collisionの後にはStreaming(隣の格子点への移動)が行われます。なお、BGKは一般に、すべての方向の分布を同じ緩和時間で扱います。この単純さが利点ですが、条件によっては数値的な安定性に課題があります。
ここで、LBM(格子ボルツマン法)の 『BGK」 は、次の人名に由来する略称です。「Bhatnagar–Gross–Krook(バートナーガル・グロス・クルック)」という、3人の研究者の姓の頭文字を取っています。(B:Bhatnagar・G:Gross・K:Krook)
BGKモデルは、もともと気体分子の分布を記述するボルツマン方程式を簡略化するために提案されたモデルです。LBMでは、衝突過程を簡単に計算する方法として広く使われています。ここで重要なのは、すべての方向の分布関数を、共通の緩和時間 $\tau$ を使って平衡分布 $f_i^{eq}$ に近づけるという点です。
【2. BGK以外の代表的なCollisionモデル】
LBMには、BGKの弱点を補い、数値計算の安定性などを改善するためのモデルがあります。代表的なものを比較します。
■ BGK ■ 単一緩和時間:基本モデル・最初に学ぶモデル
すべての分布成分を、基本的に1つの緩和時間で平衡分布に近づけます。計算が簡単で、教材や基本的な流れの計算に向いています。
注意点:条件によっては計算が不安定になることがあります。
■ TRT ■ 2つの緩和時間:BGKを少し発展させたモデル
分布を対称成分と反対称成分に分け、それぞれに異なる緩和時間を設定します。壁付近の計算誤差などを調整しやすくなります。
特徴:BGKより調整の自由度が増えます。
■ MRT ■ 複数緩和時間:各物理量に応じて緩和を調整
分布関数から求めた密度・運動量などの「モーメント」に変換して、各成分の緩和速度を個別に設定します。
特徴:BGKより安定性を改善できる場合がありますが、実装は複雑になります。
■ Cumulant ■ キュムラント:高い安定性を目指す発展的なモデル
分布関数をキュムラントという量に変換して衝突を計算します。高Reynolds数など、計算の安定性が重要な流れで用いられます。
特徴:実装が複雑で、変換や緩和の設定にも注意が必要です。
■ Entropic ■ エントロピー型:エントロピーの条件を利用するモデル
エントロピーに関する条件を使い、分布の更新を制御して数値的な安定性を高めることを目指します。
特徴:通常のBGKより計算や実装が複雑になります。
【3. どれを選べばよいか?】
| 目的 | 検討するモデル |
|---|---|
| LBMの基礎学習、簡単な流れ | BGK |
| 壁境界や誤差の調整 | TRT |
| 数値安定性を改善したい | MRT |
| 高Re流れなどで、より高い安定性を目指す | Cumulant、MRTなど |
| エントロピー条件を利用した安定化 | Entropic |
これらは単純な優劣の順番ではありません。格子モデル、境界条件、緩和パラメータ、流れの条件によって適したモデルは変わります。Autodesk XLBで学ぶ場合は、まずBGKに相当する基本的な衝突処理を理解し、次にTRT・MRT・Cumulantの違いを調べると理解しやすいです。なお、これらのモデルがXLBの特定バージョンで標準実装されているか、追加実装が必要かは別途ソースコードで確認する必要があります。
6. 式をもっと簡単に見て直観的に理解する
例えば、$f_i=10$で、$f_i^{eq}=6$だったとします。すると、$f_i-f_i^{eq}=4$です。つまり、「現在 = 10・平衡 = 6・差= 4」となっています。ここで、$\omega=0.5$なら、$f_i^{new}=10-0.5(10-6)=10-2=8$です。つまり、「10 → Collision → 8 → 6:平衡値」となり、と、10が6に近づきました。
式や値としては、上記のようになります。でもそもそも「平衡分布関数」の 「平衡」 とはなんでしょうか?「平衡」とは、方向ごとの分布が均一になることではなく、その場所の密度や流速に対して、分布が適切な形になっていることです。例えば、右向きに流れている流体では、平衡分布も右向きの流れを反映した形になります。LBMのCollisionは、現在の分布をこの平衡分布に近づけます。ただし、密度や運動量など、保存されるべき量は基本的に保たれます。
7. Collisionのイメージ
【直観的な理解に向けて】
気体や液体の中では、分子がさまざまな方向に運動し、互いに衝突しています。分子一つひとつの動きを見ると、進む方向や速さは絶えず変化しています。一方、流体全体をマクロに見ると、右方向に流れている場合があります。このとき、分子はランダムに運動していますが、その速度分布には、マクロな流れに対応した右向きの偏りがあります。つまり、分子がさまざまな方向に動いていても、全体としては右向きの流れが生じています。
では、分子同士が衝突すると、どうなるのでしょうか。衝突によって分子の運動方向や速さが変わり、右向きや左向きなど、各方向に分布する分子の割合も変化します。このような衝突が繰り返されることで、分子の速度分布は、その場所の密度や流速などに対応した平衡分布へと近づきます。ここでいう平衡分布とは、すべての方向に同じ数の分子が動いている状態ではありません。右向きに流れている流体であれば、その流れを反映した分布になります。また、平衡状態になっても、分子の運動そのものが止まるわけではありません。分子は運動を続けながら、速度分布としては平衡状態を保つことができます。
LBM(格子ボルツマン法)の Collision(衝突) は、この分子の衝突による速度分布の変化と、平衡分布へ近づく効果をモデル化した計算処理です。実際の分子一つひとつの衝突を計算するのではなく、格子点ごとに方向別の分布関数を更新することで、その効果を表現します。この後の Streaming(移動) では、方向別の分布関数が隣の格子点へ移動します。Collisionによる分布の調整とStreamingによる移動を繰り返すことで、LBMは流体の流れや粘性などのマクロな振る舞いを再現します。
つまりCollisionは、分布関数を平衡状態へ「緩和」させる処理 です。「緩和」という言葉は後でよく出てくるので、ここで覚えておくと便利です。LBMにおける「緩和」とは、方向別の分布関数が、その場所の流速や密度に対応する「平衡分布」に近づくことです。
身近なイメージで考えてみます。例えば、流体が右向きに流れているとします。このとき、分子はさまざまな方向に運動しています。方向ごとの分布が、その場所の流れに対応した平衡分布からずれていることがあります。
そこで、衝突による変化がおこり、分子同士の衝突で運動状態が変わり、方向別の分布が変化します。LBMでは、この効果をCollisionの計算でモデル化します。
【平衡分布に近づく緩和】
平衡分布に近づく(緩和)について、説明します。方向別の分布が、その場所の密度や流速に対応する平衡分布へ近づきます。ただし、分子の運動が止まるわけではありません。実際のLBMは、分子一つひとつを追跡するのではなく、格子点ごとの分布関数を計算します。
BGKモデルの衝突計算は、次の式で表されます。
$ f_i^*=f_i-\frac{1}{\tau}(f_i-f_i^{eq}) $
ここで、$f_i$:現在の方向別の分布関数、$f_i^{eq}$:平衡分布関数、$\tau$:緩和時間、$f_i^*$:衝突計算後の分布関数 この式では、現在の分布 $f_i$ と平衡分布 $f_i^{eq}$ の差を、一定の割合で小さくすることで緩和を表現します。
最後に、緩和時間 $\tau$ の意味を考えます。緩和時間は、平衡分布との差をどの程度の割合で修正するかを決めるパラメータです。「$\tau$ が小さい:1回の衝突計算で差を大きく修正する・$\tau$ が大きい:1回の衝突計算で差を小さく修正する」です。ただし、標準的なBGKモデルでは、安定性などの条件があり、$\tau$ を自由に小さくしたり大きくしたりできるわけではありません。また、粘性係数とも関係します。
まとめると、緩和とは「分布のずれを修正し、平衡分布に近づけること」です。 LBMでは、この緩和をCollisionで計算し、その後のStreamingで分布関数を隣の格子点へ移動させます。
8. ω は何を決めるのか
式の中の$\omega$は、どのくらい強く平衡状態へ近づけるかを決めます。例えば、$\omega$ が小さいと、「現在:10 → 9 → 8 → 7 → 6」ゆっくり近づく。逆に$\omega$ が大きいと、「現在:10 → 7 → 6」速く近づく。つまり、$\omega$ は分布関数の変化の強さを決めると理解できます。
9. Collisionでは「空間移動」はしていない
ここは非常に重要です。実は分布関数を考えるとき、Collision前とCollision後も、ある格子点の平衡分布に近づけますが、それは格子点のなかで配分されます。まだ別の格子点へ移動していません。
つまり、「Collision → その場所で分布を変更」です。
これに対して、次のSTEP 8のStreamingは、「● ─────→ ●」のように分布関数を隣の格子へ移します。
10. CollisionとStreamingの違い
ここはSTEP 7で必ず整理しておきます。
| Collision | Streaming | |
|---|---|---|
| 何をする? | 分布関数を変える | 分布関数を移動する |
| 場所 | 同じ格子点 | 隣の格子点へ |
| 隣の格子点を見る? | 基本的に見ない | 移動先が隣接格子 |
| 主な意味 | 緩和・衝突 | 移流・移動 |
図示すると、「Collision: ● fを変更」「Streaming:● ─────→ ● 隣の格子へ」となります。
11. なぜCollisionが必要なのか
ここで、「分布関数をただ移動させるだけではダメなの?」という疑問が出ます。ただ移動するだけだと、「流れ:→ → → → 」という情報がそのまま運ばれるだけになります。
そこでCollisionによって、「方向別の分布 → 相互作用 → 分布の再調整」を行います。これによって、粘性を持つ流体らしい振る舞いを表現できるようになります。
12. 拡散方程式とのつながり
ここでSTEP 1~5へ戻ります。拡散方程式では、「高い温度 → 周囲へ広がる → 分布が平らになる」という現象を学びました。Collisionでも、「偏った分布 → Collision → よりバランスした分布」という変化があります。したがって、「分布を調整することで、最終的にマクロな流れの性質を表現する」 という考え方につながります。ただし、Collisionそのものが拡散方程式の計算というわけではありません。
13. 「保存されるもの」も重要
Collisionでは、勝手に流体を作ったり消したりしてはいけません。そのため、適切なCollisionモデルでは、「質量・運動量」などの保存が重要になります。つまり、「Collision前 → Collision → 分布の形は変わる → でも全体としての密度・運動量は整合する」という考え方です。ここは今の段階では、「分布の形は変わるが、物理的に重要な保存則を守る」 とだけ理解すれば十分です。
STEP 7からSTEP 8への接続
最後は、次の問いでStreamingにつなげるのが自然です。「Collisionによって新しい分布関数ができました。では、その分布関数は次にどこへ行くのでしょうか?」答えは、「それぞれの方向に沿って、隣の格子へ移動する」です。これがSTEP 8のStreamingです。
STEP 7の到達目標
STEP 7では、次の3つが分かれば十分です。
① Collisionとは何か: 分布関数を、その格子点の平衡状態へ近づける処理。
② 何が起こるのか: $f \rightarrow f^{eq}$ の方向に分布が変化する。
③ Streamingとの違い: Collision=その場所で分布を変える・Streaming=分布を隣へ運
■ LBMの基礎の直観的理解【STEP 8 Streamingを知る】
STEP 8では、STEP 7のCollisionの後に起こる**「分布関数 $f$ が格子から格子へ移動する」**ことを理解します。ここでは、まだLBMの速度モデル(D1Q3やD2Q9)の細かな式に入らず、Streaming = Collision後の分布関数を、それぞれの方向に従って隣の格子へ運ぶ処理 と理解することを中心にします。
1. STEP 7から接続する
STEP 7では、「分布関数 f → Collision → その格子点の中で f を変化させる」と学びました。ここで、「新しくなった分布関数 $f$ は、そのあとどうなるのか?」という疑問が出てきます。そこで登場するのがStreamingです。つまり、「[1] Collision → 分布関数 f が変化 → [2] Streaming → 隣の格子へ移動」この2段階がLBMの基本です。
2. Streamingは「運ぶ」処理
例えば、1次元の格子を考えます。
●────●────●────●────● ●────●────●────●────●
0 1 2 3 4 0 1 2 [→]3 4
格子点2から、右方向→の分布関数があれば、右隣の格子点3へ移動します。
3. 左向きなら反対に移動する
同じように、左方向←の分布関数があれば、右隣の格子点3へ移動します。
●────●────●────●────●
0 1[←] 2 3 4
4. Streamingでは「方向」が重要
STEP 6で、分布関数には「方向」があると学びました。その方向が、Streamingでそのまま使われます。つまり、「fがどちらへ進むか」を決めているのが方向情報です。
f↑:上の格子へ
左の格子へ:f← ● f→:右の格子へ
f↓:下の格子へ
5. CollisionとStreamingを比較する
ここはSTEP 8で必ず整理します。
【Collision】 その格子点の中で分布を変える:Collision = 変える
● → fを変更:新しいf
【Streaming】 分布を隣の格子へ運ぶ:Streaming = 運ぶ
● ─────────→ ●
6. 1次元で考えると非常に分かりやすい
D1Q3につながるように、以下のように考えます。
← 2 →
●──●──●──●──●
各格子点には、「左へ進む分布・停止している分布・右へ進む分布」があります。
例えば格子点2で、「f← = 2・f₀ = 5・f→ = 8」だったとします。Collision後に、「f←* = 3・f₀* = 5・f→* = 7」になったとします。この * は、Collision後の分布 を表すことにします。するとStreamingでは、「f←* = 3」は左へ、「f→* = 7」は右へ移動します。以下のイメージです。
●────●────● → ●────●────● collision後
2← 5 8→ 3← 5 7→
7. Streamingの後は、隣の格子の状態が変わる
ここが重要です。Streamingによって、例えば、「格子2 → 格子3」へ分布関数が移動すると、「格子3」に新しい分布が入ってきます。つまり、「格子1の情報 → 格子2へ・格子2の情報 → 格子3へ・格子3の情報 → 格子4へ」というように、情報が格子全体を移動していく ことになります。
8. これが「流れ」の基本になる
例えば、「→ → → → →」という分布があれば、「● → ● → ● → ● → ●」と情報が右へ進みます。これを時間ステップごとに繰り返すことで、流体の流れを表現できます。 つまり、Streamingは単なるデータ移動ではありますが、LBMではこの方向性を持った情報の移動が、流体の移流・輸送につながります。
9. CollisionとStreamingを連続して考える
ここまで来たら、LBMの基本ループが見えてきます。これが非常に重要です。
現在の分布関数 f → 【Collision】 → Collision後の分布関数 f
→ 【Streaming】 → 次の時刻の分布関数 f → 繰り返す
10. 「格子と時間発展」ともつながる
STEP 5で、「現在 → 計算 → 次の状態」という時間発展を学びました。LBMでは、「現在の f → Collision → Streaming → 次の f」となります。したがって、「STEP 5:格子上で時間発展 → STEP 6:分布関数 f → STEP 7:Collision → STEP 8:Streaming」が一本につながります。
11. STEP 8では「1ステップ」を理解する
ここでは、LBMの1時間ステップを次のようにまとめるとよいです。「時刻 n → 各格子点に f がある → Collision → 各格子点で f を修正 → Streaming → f を隣の格子へ移動 → 時刻 n+1」つまり、CollisionとStreamingを1組として、1ステップの計算を行う ということです。
12. 拡散方程式との違いもここで整理する
STEP 1~5では、拡散方程式で、「温度 T → 周囲の値を使う → 次の T」と計算していました。LBMでは、「分布関数 f → Collision → Streaming → 次の f」です。したがって、両者には、現在の状態から次の時間の状態を作るという共通点があります。
ただし、LBMでは単なる温度ではなく、方向を持った分布関数を時間発展させています。
13. Streamingで「何が移動しているのか」
初心者が間違えやすいところです。Streamingで移動しているのは、「格子そのもの」ではありません。 また、「流体粒子1個1個」でもありません。 移動しているのは、分布関数 $f_i$ が表している流体の情報です。この理解が大切です。
14. D1Q3への接続
ここで次のSTEP 9へつなげます。STEP 8では、「右へ・左へ」といった方向別のStreamingを理解しました。では、「1次元で、具体的に何方向の分布を持たせればよいでしょうか?」 という問いにつなげます。そこで、「← ● →」の3つを使うD1Q3が登場します。
「D1 = 1次元・Q3 = 3つの速度方向」です。つまり、「STEP 6:分布関数 f → STEP 7:Collision → STEP 8:Streaming → STEP 9:D1Q3」となります。
STEP 8の到達目標
この段階では、次の3つが分かれば十分です。
① Streamingとは何か:Collision後の分布関数を、決められた方向へ隣の格子へ移動させる処理。
② Collisionとの違い:Collision → 分布を変える + Streaming → 分布を運ぶ
③ 1ステップの流れ:分布関数 f → Collision → Streamin → 次の分布関数 f
この3点が理解できれば、次のSTEP 9では、いよいよD1Q3を使ってCollisionとStreamingを実際の数字で手計算することができます。そして、ここで初めて「分布関数 (f) が実際にどう動くのか」を紙の上で追跡すると、LBMの仕組みがかなり明確になります。
もう一度「Collision:緩和」「Streaming:移動」を直感的に理解する
これまでの説明だと「Collision=緩和」「Streaming=移動」という言葉が抽象的すぎます。ここでは数式をいったん忘れて、「3つの箱」と「矢印のついたカード」 で説明すると分かりやすいです。
【まず、LBMを「カードの移動」として考える】
1次元の簡単な世界を考えます。
格子0 ● 格子1 ● 格子2 ● 格子3 ● 格子4 ●
それぞれの格子点には、「← / ○ / →」という3種類のカードがあるとします。つまり、「←カード:左へ進む情報・○カード:その場に残る情報・→カード:右へ進む情報」です。ここで大切なのは、カード1枚=流体分子1個ではありません。 「その方向へ進む流体の情報」を表しています。
【1. まず、ある格子点だけを見る】
例えば格子2に、「← / → / ○」のカードが何枚か入っているとします。数字(枚数)で書くと、「←左方向:2・停止: 6・→右方向:6」です。つまり、合計は、$2+6+6=14$です。
【Collisionは「同じ格子の中でカードの配分を変える」】
ここがCollisionです。格子2の中だけで、Collision前「← = 2・○ = 6・→ = 6」だったものを、Collision後「← = 3・○ = 4・→ = 7」に変えたとします。格子2から外にはまだ何も出ていません。 これが非常に重要です。
【では、なぜカードの枚数を変えるのか?】
「なぜ勝手に2枚を3枚にするの?」と思いますよね。LBMでは、流体が衝突した結果、進む方向の分布が変化する という現象を数式で表しています。
実際には、「2 → 3・6 → 4・6 → 7」という単純な整数操作をしているわけではなく、$f$ という数値を計算します。でも初心者の段階では、Collision = 「その場所で、方向別の流体の割合を組み替える」 と考えてください。これならかなり具体的です。
【Collisionの前後を並べる】
格子2の中だけを見るて、Collisionによって、「← = 2・○ = 6・→ = 6」だったものを、「← = 3・○ = 4・→ = 7」にして、Collisionは「その場で配分を変える」 だけです。
【次がStreaming】
ここで初めて、格子と格子の間を考えます。「格子1 ● 格子2 ● 格子3 ●」Collision後の格子2は、「← = 3・○ = 4・→ = 7」でした。それぞれのカードには進む方向があります。したがって、「← 3枚は左へ行きます」「→ 7枚は右へ行きます」「○ 4枚はその場に残ります。」
【Streamingするとこうなる】
格子1 ● 格子2 ● 格子3 ●
←3枚 ○4枚 →7枚
つまり最終的に、「格子1 ← 3・格子2 ○ 4・格子3 → 7」という形になります。これがStreamingです。
【つまり、ものすごく簡単に言うと】
Collision = 格子の中で「向きを組み替える」
Streaming = 組み替えたものを「隣の格子へ運ぶ」 このような機能です。
【「人」で考えるとさらに分かりやすい】
例えば駅が3つあります。「駅A ─── 駅B ─── 駅C」駅Bに14人いるとします。その人たちが、「左へ行く人:2人・駅Bに残る人:6人・右へ行く人:6人」だったとします。
=Collision:駅Bの中で、人の行動が変化して、「左へ行く人:3人・駅Bに残る人:4人・右へ行く人:7人」になった。これがCollisionです。
=Streaming:次に、「左へ3人 → 駅A・その場4人 → 駅B・右へ7人 → 駅C」と移動する。これがStreamingです。
【実際のLBMでは何が違うのか】
ここで重要な注意があります。LBMで本当に計算しているのは、「人」ではなく、「分布関数 f」です。例えば、「格子2:f← = 2・f0 = 6・f→ = 6」という数値を持っています。
Collisionによって、「f← = 3・f0 = 4・f→ = 7」のように変化する。
Streamingによって、「格子1:f←・格子2:f0・ 格子3:f3」と隣の格子へ移ります。
【これを時間ごとに繰り返す】
1回だけで終わりません。時刻 n において、「各格子に f がある → Collision → 各格子で f を組み替える → Streaming → f が隣の格子へ移動 → 時刻 n+1 → また Collision → また Streaming → 繰り返す」これがLBMの基本的な時間発展です。
【ここで「流れ」が見えてくる】
例えば、全体として「→方向の f が多い」なら、「→ → → →」と情報が右へ進みます。反対に、「←方向の f が多い」なら、「← ← ← ←」となります。つまり、Streamingによって方向を持った情報が格子から格子へ移動することが、流れを表現する わけです。
【拡散方程式とのつながり】
ここでSTEP 1~5に戻ります。拡散方程式では、「格子点の温度 → 周囲との関係から新しい温度 → 次の時間」でした。LBMでは、「格子点の分布関数 f → Collision → Streaming → 次の時間」です。つまり、「拡散方程式:値を更新する」、「LBM:分布を更新して運ぶ」という違いがあります。
【CollisionとStreamingを一枚で見る】
初心者には、この図を覚えてもらうのが一番よいと思います。もっと簡単にすると、「Collision = 混ぜ直す・Streaming = 運ぶ」です。
① Collision:「その場所で組み替える」 ② Streaming:「方向に沿って運ぶ」
格子 ● ● ←──── ● ────→ ●
←2 ○6 →6 格子1 格子2 格子3
←3 ○4 →7 3 4 7
【なぜ2つに分けるのか】
これも直感的に考えられます。流体には、「① その場所で他の流体と影響し合う・② その結果、別の場所へ移動する」という2つの側面があります。LBMでは、それを計算上、「① Collision → ② Streaming」に分けています。
【ここまで理解したら、次はD1Q3】
ここまでの説明で、「← / ○ / →」という3種類だけを使えば、かなり具体的にLBMが見えてきます。そこで次のSTEP 9では、D1Q3を使って本当に数字を入れて計算します。ここまで来ると、「Collision」と「Streaming」が単なる言葉ではなく、実際に数字がどこへ行くのかとして理解できるようになります。
今回の学習ルートでは、私は**STEP 9のD1Q3を「LBMを初めて実際に計算する回」**にして、5個程度の格子だけで手計算する構成をおすすめします。
■ LBMの基礎の直観的理解【STEP 9 D1Q3でLBMを体験する】
STEP 9では、「D1Q3を実際に数字で計算する」 ところまで進めます。ここでは、これまでのSTEP 6~8を全部使って、「分布関数 f → Collision → Streaming」を手計算します。
1. まずD1Q3とは何か
D1Q3は、「D1:1次元・Q3:3方向」という意味です。1つの格子点に、「左 ←・停止 ○ ・→ 右」の3種類の分布関数があります。数式では、$f_{-}, f_0, f_{+}$と書きます。
2. まず1個の格子点だけを見る
例えば格子点1を考えます。「格子点1 ● : ← ○ → 」この格子点に、$f_- = 4$・$f_0 = 4$・$f_+ = 4$という分布関数があるとします。つまり、「左方向:4・停止:4・右方向:4」です。
3. まず密度を計算する
分布関数を全部足します。$\rho=f_-+f_0+f_+$ したがって、$\rho=4+4+4=12$ です。
改めて物理的な意味を解説します。格子ボルツマン法(LBM)における質量、速度、運動エネルギーの関係は、マクロ(巨視的)な流体としての物理量と、ミクロ・メゾ(微視的)な仮想粒子の統計的な挙動を結ぶ「統計のモーメント(平均操作)」によって定義されます。
「基本的な定義(速度分布関数との関係)」
LBMでは、格子点ごとに特定の離散速度 $c_i$ を持った仮想粒子の「粒子速度分布関数 $f_{i}$」を計算します。この $f_{i}$ を用いて、マクロな質量密度 $\rho $、流速 $u$ は次のように結びついています。
• 質量密度(マクロ質量) $\rho $
すべての方向の分布関数を足し合わせたものです(0次モーメント)。
$\rho=\sum_{i}f_{i}$
• 流動量(マクロ質量 × マクロ速度) $\rho u$
分布関数に粒子の離散速度をかけて足し合わせたものです(1次モーメント)。
$\rho u=\sum_{i}f_{i} c_{i}$
4. 次に速度を計算する
D1Q3では速度方向(離散速度)は、$c_-=-1, c_0=0, c_+=1$ です。したがって、流動量 $\rho u=(-1)f_-+0f_0+(1)f_+$ となります。今回、$\rho u=-4+0+4=0$ なので、質量密度(マクロ質量)$\rho =12$ですから、マクロ速度 $u=0$ です。つまり、「← 4 ○ 4 →」で左右が同じなので、全体としては流れがない状態です。
5. 平衡分布 feq を求める
Collisionでは、現在の分布を、平衡分布に近づけるのでした。D1Q3の標準的な重みは、$w_0=\frac23$、$w_- =w_+=\frac16$です。今回、$\rho=12, u=0$なので、平衡分布は非常に簡単になります。
$f_-^{eq}=\frac16\times12=2$ $f_0^{eq}=\frac23\times12=8$ $f_+^{eq}=\frac16\times12=2$ したがって、$f_-^{eq}=2, f_0^{eq}=8, f_+^{eq}=2$ です。
6. Collisionをする
今回、分布関数は、「Collision前:← 4・○ 4・→ 4」でした。平衡分布は、「平衡状態:← 2・○ 8・→ 2」です。つまり、「Collision前:4 4 4 →Collision:平衡に近づく→ 2 8 2」となります。
7. Collisionの式を使う
BGKモデルでは、$f_i^{*}=f_i-\frac{1}{\tau}(f_i-f_i^{eq})$ です。今回は計算を簡単にするため、$\tau=1$とします。すると、$f_i^{*}=f_i-(f_i-f_i^{eq})$ なので、$f_i^*=f_i^{eq}$ となります。つまり、「Collision前:←4 ○4 →4 →Collision後→ ←2 ○8 →2」です。ここで、格子点の中だけで分布を変えています。 まだどこにも移動していません。
8. ここからStreaming
ここがSTEP 8で学んだStreamingです。Collision後は、「← 2・○ 8・→ 2」です。それぞれの進む方向に従って移動します。「左向き2 ← ←・停止8 ○○○○○○○○・右向き2 → →」
9. 3個の格子を用意する
ここで、「格子0 ● ─────格子1 ● ─────格子2 ●」とします。Collision後の格子1は、「格子1:← 2・○ 8・→ 2」でした。
10. Streamingするとどうなるか
左向きの2は、左隣へ移動します。「格子0 ● ── ←2 ──格子1 ● ─────格子2 ●」
右向きの2は、右隣へ移動します。「格子0 ● ─────格子1 ● ── 2→ ──格子2 ●」
停止している8は動きません。 「格子0 ● ─────格子1 ● 8 ─────格子2 ●」
11. ここで「Collision」と「Streaming」がはっきり分かれる
今回の計算を並べると、以下のようになります。
【最初】格子1: ←4 ○4 →4 →【Collision】→ 格子1:←2 ○8 →2
【Streaming】格子0 [←2] 格子1[8] 格子2 [2→]
つまり、「Collision = その格子点の中で数字を変える」「Streaming = 変わった数字を隣の格子へ運ぶ」です。ここまでなら、かなり具体的にイメージできます。
12. 「1ステップ」とは、この一連の動作
LBMの1ステップは、「現在の f → ① 密度・速度を計算 → ② 平衡分布 feqを計算 → ③ Collision → ④ Streaming → 次の時刻の f」です。D1Q3では、この動作を毎回繰り返します。
13. D1Q3を時間発展として見る
例えば中央の格子を考えると、以下のようになります。
時刻 n → Collision → Streaming → 時刻 n+1
格子1 格子0 格子1 格子2
←4 ○4 →4 ←2 ○8 →2 ←2 8 2→
そして次の時刻では、以下のように繰り返します。
時刻 n+1 → Collision → Streaming → 時刻 n+2
14. STEP 5とのつながり
ここで、これまでの学習が一本につながります。
【拡散方程式】:現在の温度 → 計算 → 次の温度
【D1Q3】:現在の分布関数 f → Collision → Streaming → 次の分布関数 f
どちらも、「現在の状態 → 計算 → 次の状態」という時間発展です。
15. ただし、ここで大きな違いがある
拡散方程式では、「格子点 → 温度 T」だけでした。D1Q3では、「格子点 → f₋ f₀ f₊」の3つの値を持っています。つまり、「拡散方程式:● → T」「D1Q3:● → f₋, f₀, f₊」です。これが、LBMの非常に重要な特徴です。
16. 「なぜ3方向なのか」が見えてくる
D1Q3では、「← ○ →」の3つだけです。したがって、「f₋ → 左へ進む情報・f₀ → その場に残る情報・f₊ → 右へ進む情報」と考えられます。つまりD1Q3は、「方向を持った分布関数を最も簡単な形で学ぶためのLBM」 です。
17. STEP 9では、まずこの計算を紙に書く
教材としては、次の問題を実際に手計算すると非常に良いです。
問題:ある格子点で、$f_- = 4,\quad f_0=4,\quad f_+=4$とする。$\tau=1,\quad \rho=12,\quad u=0$とする。
計算;
① 密度:ρ = 4 + 4 + 4 = 12
② 速度:u = (-4 + 4)/12 = 0
③ 平衡分布:feq = (2, 8, 2)
④ Collision:f* = (2, 8, 2)
⑤ Streaming:左2 → 左の格子・8 → その場・右2 → 右の格子
これを自分の手で紙に書くことがSTEP 9の中心です。
18. ここで初めて「LBMの計算」が見える
これまでのSTEP 6~8では、「分布関数・Collision・Streaming」を概念として説明してきました。STEP 9では、それが、「4, 4, 4 → 2, 8, 2 → 左へ2・中央に8・右へ2」という具体的な数字の動きになります。ここで初めて、**「ああ、これがLBMの計算なのか」**と理解できるようになります。
STEP 9の到達目標
STEP 9では、次の3つができれば十分です。
① 1つの格子点に3つの (f) があることを説明できる
← f₋ ○ f₀ → f₊
② Collisionを手計算できる
現在の f → 平衡分布 feq → Collision
③ Streamingを図に描ける
← 左へ・○ その場・→ 右へ
これができれば、STEP 10のD2Q9では、この「←・○・→」という考え方を2次元に広げて、9方向の分布関数へ自然に進めます。
■ LBMの基礎の直観的理解【STEP 10 D2Q9で2次元化する】
STEP 10は、STEP 9のD1Q3で理解した**「左・停止・右の3方向の分布関数を、CollisionしてからStreamingする」**という考え方を、2次元へ拡張する段階です。
ここでは、単に「3方向が9方向になる」と覚えるのではなく、「なぜ9方向になるのか → 9個の分布関数は何を意味するのか → CollisionとStreamingはどうなるのか → 3次元ではどう拡張するのか」という順番で説明するのがよいと思います。それでは、「D2Q9 ― 1次元LBMを2次元へ拡張する」の説明に進みます。
1. D1Q3を思い出す
STEP 9では、1次元なので方向は3つでした。「左 ← 停止 ○ 右 →」したがって、1つの格子点に$f_-,\quad f_0,\quad f_+$の3つの分布関数があります。そして、「f → Collision → Streaming → 次の時間」という計算をしました。
2. 2次元になると「上下」が必要になる
2次元の流体では、「右・左・上・下」だけでなく、「右上・左上・右下・左下」という斜め方向も考える必要があります。したがって、1つの格子点を中心にすると、以下のような9方向になります。これがD2Q9です。「D2 = 2次元・Q9 = 9個の離散速度」です。
↖ ↑ ↗
← ● →
↙ ↓ ↘
3. 9方向を具体的に見る
D2Q9では、例えば次のように番号を付けます。
5 ↖ ↑2 ↗ 3
← 1 ● 6 → 0 停止
8 ↙ ↓7 ↘ 4
重要なのは番号そのものではなく、1つの格子点に「9方向の流体情報」が存在することです。
4. 1次元から2次元への拡張
D1Q3では、「← ○ →」でした。2次元にすると、説明図のように「8方向+停止」となります。つまり、「D1Q3 → 左右+停止 → D2Q9 → 上下+斜め方向を追加」です。この「方向を増やす」という考え方が、D2Q9を理解するポイントです。
5. 1つの格子点に9個のfがある
例えば格子点 $(x,y)$ について、$f_0,f_1,f_2,\ldots,f_8$があります。具体的には、説明図の項目5となります。ここで、$f_0$は停止成分、$f_1,\ldots,f_8$は移動方向を持つ成分です。
6. D2Q9では「方向ベクトル」を持つ
D1Q3では、$c_-=-1,\quad c_0=0,\quad c_+=1$でした。D2Q9では2次元ベクトルになります。例えば、
$\mathbf c_0=(0,0)$ $\mathbf c_1=(-1,0)$ $\mathbf c_2=(0,1)$ $\mathbf c_3=(1,1)$などです。つまり、各分布関数 $f_i$ に「どちらへ進むか」という速度ベクトル $\mathbf c_i$ が対応しているわけです。この考え方がStreamingの中心になります。
7. 例えば右方向の分布
右方向の分布を$f_6$とします。するとStreamingでは、以下のように移動します。つまり、$f_6(x,y)\rightarrow f_6(x+1,y)$です。
現在 →
● ●──────●隣の格子
8. 上方向ならどうなるか
上方向の分布なら、説明図の項目8に、示すようになります。つまり、$f_2(x,y)\rightarrow f_2(x,y+1)$です。
9. 斜め方向も同じ
例えば右上方向なら、説明図の項目9に、示すようになります。したがって、$f_3(x,y)\rightarrow f_3(x+1,y+1)$となります。
ここがD2Q9の面白いところです。Streamingは9方向それぞれについて独立に行われます。
これまでの1から9までの、図解を以下に示します。
10. CollisionはD1Q3と基本的に同じ
ここが非常に重要です。D1Q3からD2Q9になっても、Collisionの考え方は変わりません。D1Q3では、$f_i^*=f_i-\frac{1}{\tau}(f_i-f_i^{eq})$でした。
D2Q9でも同じです。$f_i^*=f_i-\frac{1}{\tau}\left(f_i-f_i^{eq}\right)$ ただし、$i=0,1,\ldots,8$について計算します。
つまり、「D1Q3:3個のfをCollision・D2Q9:9個のfをCollision」です。
11. では「平衡分布」はどうなるのか
D2Q9では、平衡分布関数は一般に、$f_i^{eq}=w_i\rho\left[1+\frac{\mathbf c_i\cdot\mathbf u}{c_s^2}+\frac{(\mathbf c_i\cdot\mathbf u)^2}{2c_s^4}-\frac{\mathbf u\cdot\mathbf u}{2c_s^2}\right]$と書けます。
初心者の段階では、この式を暗記する必要はありません。重要なのは、「密度 ρ → 速度 u → 平衡分布 feq → 各方向の f の基準を決める」という関係です。
12. D2Q9の「重み」が重要
D2Q9では方向ごとに重み $w_i$ が違います。代表的には、$w_0=\frac49$、軸方向は、$w_{1,2,6,7}=\frac19$ 斜め方向は、$w_{3,4,5,8}=\frac1{36}$です。
図にすると、以下のようになり、これらは、2次元の等方的な流れをうまく表現するための重みと考えればよいです。
1/36 1/9 1/36
↖ ↑ ↗
1/9 ← ● → 1/9 4/9(停止)
↙ ↓ ↘
1/36 1/9 1/36
13. 密度と速度も9個のfから求める
D1Q3では、$\rho=f_-+f_0+f_+$でした。D2Q9では、$\rho=\sum_{i=0}^{8}f_i$です。
つまり、「9方向のf → 全部足す → 密度ρ」です。
14. 速度はベクトルになる
D1Q3では速度は1方向だけでした。$u$ D2Q9では2次元なので、$\mathbf u=(u_x,u_y)$になります。式は、$\rho\mathbf u=\sum_i f_i\mathbf c_i$です。つまり、「各方向のf × 各方向のc → 全部足す → ux, uy」です。
15. 例えば右向きの流れ
右方向の分布が多いとします。数値的には、「右向き f が大きい・左向き f が小さい」なら、$u_x>0$となります。上向きと下向きも比較して、$u_y$が決まります。そのため、「fの分布 → 平均すると → 速度ベクトル」が現れます。
16. D2Q9の1時間ステップ
ここまでをつなげると、「時刻 n → 9方向のf → ρ , u を計算 → feqを計算 → Collision → 9方向のf*ができる → Streaming → それぞれの方向へ移動 → 時刻 n+1」です。D1Q3と構造は同じです。ただし、「3方向 → 9方向」になっています。
17. 実際の格子全体ではこうなる
2次元格子を考えます。各格子点 ● に、「9個の f(分布関数)」があります。Collisionでは、それぞれの格子点の中で処理します。
2次元格子 その後Streamingによって情報が周囲へ移動します。
●──●──●──● ● → ● → ● → ● → ●
│ │ │ │ 格子点での ↑ ↗ ↓ ↘ ↘
●──●──●──● Collision ● → ● → ● → ● → ●
│ │ │ │ ↓ ↘ ↑ ↗ ↗
●──●──●──● ● → ● → ● → ● → ●
18. ここで「流れ」が2次元になる
例えば左から右へ流体を流すと、以下の主流ができます。そこに障害物を置くと、流体は、のように上下に回り込めます。これがD2Q9で初めて表現できる2次元的な流れです。
→ → → → → → → → → → → → → → ↗ ↘
→ → → → → → → → → → ● → → → → → → ● → → →
→ → → → → → → → → → → → → → ↘ ↗
19. D1Q3からD2Q9への本質的な変化
| D1Q3 | D2Q9 |
|---|---|
| 1次元 | 2次元 |
| 3方向 | 9方向 |
| $f_-,f_0,f_+$ | $f_0\sim f_8$ |
| 速度 $u$ | 速度 $(u_x,u_y)$ |
| 左右へ移動 | 上下・左右・斜めへ移動 |
| 直線的なモデル | 2次元流れを表現 |
20. ここで非常に重要な「D2Q9のイメージ」
私は教材では、D2Q9をまず次のように覚えてもらうのがよいと思います。つまり、**「1つの格子点から、8方向へ移動する情報+その場に残る情報」**です。これを覚えれば、D2Q9の意味がかなり分かります。
↖ ↑上↗
左 ← ● → 右 + その場に残る
↙下↓ ↘
これまでの10から20までの、図解を以下に示します。
21. では、3次元ではどうなるか?
ここから3次元に拡張します。D2Q9では、「x方向・y方向」の2次元でした。3次元では、「x方向・y方向・z方向」が必要です。したがって、流体は、3次元空間を動きます。
22. 3次元LBMの考え方
3次元でも、基本的な考え方は変わりません。ここで、「D3 = 3次元・Q19 = 19方向・Q27 = 27方向」です。「【1次元】D1Q3:3方向 → 【1次元】D2Q9:9方向 → 【1次元】D3Q19:19方向・D3Q27:27方向」
23. D3Q19
D3Q19では、代表的に、「1つの停止 + 6つの軸方向 + 12個の面対角方向」を使います。したがって、1+6+12=19 です。イメージとしては、「停止 → 6方向 ±x, ±y, ±z → 面の対角方向 → 合計19」です。
24. D3Q27
D3Q27では、3次元の格子で考えられます。「停止 + 軸方向 + 面対角方向 + 立体対角方向」を、全部含めます。数えると、「1:停止 + 6:軸方向 + 12:面対角 + 8:立体対角 = 27方向」です。
25. D2Q9からD3Q27を見る
これは非常に分かりやすい関係です。
D2Q9 ↖ ↑ ↗ D3Q27
← ● → これをさらにz方向へ立体的に拡張します。
↙ ↓ ↘ D2Q9の「平面上の方向」を、3次元空間へ拡張したもの
26. D3Q19とD3Q27の違い
初心者向けには、まず次の理解で十分です。D3Q27は立体対角方向まで含めるのが大きな違いです。
一般に、D3Q27はより多くの離散方向を持つため、方向性の表現をより完全にしたモデルとして使えます。一方、D3Q19は方向数が少ないため、計算量・メモリ使用量とのバランスを取りやすく、3次元LBMで広く使われています。
| D3Q19 | D3Q27 | |
|---|---|---|
| 次元 | 3D | 3D |
| 速度方向 | 19 | 27 |
| 軸方向 | 6 | 6 |
| 面対角 | 12 | 12 |
| 立体対角 | なし | 8 |
| 停止 | 1 | 1 |
27. 3次元でもCollision → Streaming
3次元になっても計算の基本は変わりません。「各格子点:f₀~f₁₈ → Collision → f₀~f₁₈ → Streaming → 次の格子へ」D3Q27なら、「f₀~f₂₆」を扱います。
つまり、D1Q3 → D2Q9 → D3Q19/D3Q27と進んでも、「Collision → Streaming → 次の時間」という基本構造は変わりません。
28. STEP 10で一番重要なのは「数字が増えた」ことではない
ここが最も重要です。以下のように拡張されています。
D1Q3: → D2Q9: → D3Q19 / D3Q27:
← ○ → ↖ ↑ ↗ 3次元のあらゆる方向へ
← ● →
↙ ↓ ↘
つまり、LBMは、格子点に「方向を持った分布関数」を置き、その方向に従って情報を移動させることで、流体の運動を表現しているということです。
これまでの21から28までの、図解を以下に示します。
29. STEP 10の学習順
ここまでを一度に全部やるのではなく、次の順番が良いと思います。
① D1Q3を思い出す ② 1次元 → 2次元 ③ 3方向 → 9方向
④ D2Q9の9つのfを見る ⑤ D2Q9のCollision ⑥ D2Q9のStreaming
⑦ 速度・密度を求める ⑧ 2次元流れを考える ⑨ 2次元 → 3次元
⑩ D3Q19 / D3Q27
STEP 10の到達目標
最終的に、次の3つを説明できれば十分です。
① D2Q9とは何か: 2次元空間で、8方向への移動+停止の計9種類の分布関数を使うLBMモデル。
② D1Q3から何が変わったか: 方向が増えて、2次元の流れを表現できるようになった。
③ 3次元ではどうなるか: 方向を3次元へ拡張するが、Collision → Streamingという基本構造は変わらない。
そして、このSTEP 10を終えると、次のSTEP 11「境界条件」 が非常に自然になります。例えば2次元で、以下の流れを考えたとき、「障害物の格子にStreamingしてきた分布関数を、どう処理するのか?」 という問題が出てきます。ここで初めて、Bounce Back(跳ね返し) などの境界条件が登場します。これがSTEP 11へのきれいな接続になります。
■ LBMの基礎の直観的理解【STEP 11 境界条件を入れる】
STEP 11は、これまで学んだLBMを**「実際の流体解析」にするための重要な段階です。STEP 10までで、「D2Q9 → 9方向の分布関数 f → Collision → Streaming」まで理解しました。
しかし、ここで障害物や壁を置くと、「Streamingしてきた分布関数が、壁にぶつかったらどうするのか?」** という問題が出てきます。これが境界条件です。
1. そもそも境界条件とは何か
例えば、2次元の流体を考えます。ここには、「左:流体が入る・右:流体が出る・上下:壁・中央:障害物」があります。
流入口 → 流出口 →
┌─────────────────────────────┐ 壁
│ → → → → → → → → → → → → → → │
│ → → → → ■ → → → → → → → │
│ → → → → ■ → → → → → → → │
│ → → → → ■ → → → → → → → │
│ → → → → → → → → → → → → → → │
└─────────────────────────────┘ 壁
LBMの計算領域の中だけでは、「Collision → Streaming」を計算できます。しかし、「計算領域の外」 にStreamingしてしまった分布関数は、どうするのでしょうか?そのため、境界で分布関数 $f$ をどう処理するか を決める必要があります。
2. 境界条件は「壁で何が起こるか」を決める
例えば固体壁なら、流体は壁を突き抜けない という条件が必要です。さらに通常の粘性流体では、壁面では流速が0になるというNo-slip条件(すべりなし条件) を使います。これをLBMで実現する代表的な方法が、Bounce Back(反射) です。
壁 ██████████████████
流体 → → → → → → → → → 壁面速度:u = 0
→ → → → → → → → →
3. Bounce Backを直感的に考える
1次元D1Q3で考えると、とても簡単です。格子1から右へ向かう分布「→」があったとします。そのままでは壁の中へ入ってしまいます。そこで、「→」を「←」に反転させます。つまり、壁に向かって進む分布関数を、反対方向へ跳ね返す という処理です。
格子0 格子1 壁 反転:Bounce Back
● ● ────→ █ ● ←──── █
4. 数値で見ると分かりやすい
例えば壁に向かう分布が、$f_{\rightarrow}=5$ だったとします。壁でBounce Backすると、壁から戻ることになり、$f_{\leftarrow}^{new}=5$とします。つまり、値5が消えるのではなく、方向が反転すると考えると分かりやすいです。
5. D2Q9では8方向が対象になる
2次元では、9個の方向があります。
壁 █████████ 例えば上壁にぶつかる分布 これらを反対方向へ跳ね返します。
流体 ↖ ↑ ↗ ↖ ↑ ↗ ↙ ↓ ↘ Bounce Back
← ● →
↙ ↓ ↘
6. これで「壁」が作れる
D2Q9では、障害物を単純に、「■ = 固体・○ = 流体」と判定します。例えば、以下の流体空間を考えます。Streamingで、「○ → ■」となった分布関数に対して、「← ■」と反射させます。これで、流体が障害物を通り抜けなくなります。
○ ○ ○ ○ ○ ○
○ ○ ○ ■ ■ ○
○ ○ ○ ■ ■ ○
○ ○ ○ ○ ○ ○
7. 「Halfway Bounce Back」が重要
実際のLBMでは、Halfway Bounce Backがよく使われます。考え方は、実際の壁を格子点そのものではなく、格子点と格子点の中間に置く というものです。右図のように、実際のLBMでは、Halfway Bounce Backがよく使われます。考え方は、実際の壁を格子点そのものではなく、格子点と格子点の中間に置く というものです。右図のように、
流体格子 ● ───────── ■ 壁格子 ●──────│──────■
この間に壁がある ↑ より正確には 壁
8. FullwayとHalfway
ここでは名前だけ整理しておきます。初心者の段階では、「Halfway Bounce Backは、壁を格子点の間に置くことで、壁位置をより自然に扱う方法」 と理解しておけば十分です。
| 方法 | 考え方 |
|---|---|
| Fullway Bounce Back | 格子点上で反射させる考え方 |
| Halfway Bounce Back | 格子点間の中間に壁を置いて反射させる |
9. 流入口には別の境界条件が必要
壁だけではありません。例えば左から流体を入れたい場合、「流入部 → → 」としたいわけです。そこで入口では、流速 $u$ を指定する という条件を与えます。例えば、$u_x=0.1,\quad u_y=0$とします。つまり、「入口: → 0.1:流速を指定」です。
10. 流入口では「fの一部が分からない」
ここはLBM独特です。入口では、「→」方向の分布だけが分かっていても、「←・→・↑・↓・↗・↘」など、Streamingとの関係で一部の $f_i$ が未知になる場合があります。そこで、「指定した速度・密度と整合するように、未知の分布関数を計算する」 という処理が必要になります。
代表的な方法の一つが、Zou–He境界条件です。これは後で具体的な式を学べばよく、STEP 11では、> 入口では、指定した流速に合うように未知の $f$ を決めると理解すれば十分です。
11. 出口はさらに難しい
出口では、「→ → → → │ 出口」となります。ここで、「流体を自然に外へ出したい」 という条件が必要です。単純に、「出口のf = その一つ手前」などとする簡単な方法もあります。
ただし、流れによっては、下図のような*逆流(backflow)**が発生することがあります。そのため、出口条件は入口より少し難しくなります。
→ → → │
← ← ↙ │ 逆流
12. 圧力を指定する方法もある
入口で速度ではなく、$p = p_{in}$を指定する方法もあります。例えば、「入口:p = 101325 Pa」出口で、「p = 101300 Pa」とすると、「圧力差 → 流体が流れる」という設定ができます。LBMでは圧力と密度の関係が重要で、$p=c_s^2\rho$という関係を使います。したがって、圧力条件を密度条件として扱うという考え方も出てきます。
13. 周期境界条件
もう一つ非常に便利な境界条件があります。Periodic Boundaryです。例えば右端から出たものが、左端に戻ってくると考えます。右端から出たら、つまり、「右端から出る → 左端に入る」です。これは、「無限に続く流れ・繰り返し構造・チャンネル流れ」などで便利です。
左端 右端
● → → → → → → → → ●
↑ ↓
└─────────────────┘
14. 境界条件を4種類に整理
STEP 11では、まず次の4つを覚えるとよいです。
| 境界 | 代表的な条件 | 目的 |
|---|---|---|
| 壁 | Bounce Back | 流体を壁から反射 |
| 入口 | Velocity BC | 流速を指定 |
| 出口 | Pressure / Outlet BC | 流体を外へ出す |
| 周期境界 | Periodic BC | 反対側とつなぐ |
15. ここから「数値計算上の工夫」が重要になる
ここまでは物理的な意味です。しかし実際のLBMでは、境界条件を正しく計算機上に実装することが難しくなります。代表的な工夫を順番に見ます。
【工夫① 固体を「マスク」で管理する】
例えば、「fluid 流体 = 0・solid 固体 = 1」として、以下のような配列を作ります。とすれば、Streamingのとき、「流体 → 流体」なら普通に移動。「流体 → 固体」ならBounce Back、と分岐できます。これは実装上かなり重要です。
0 0 0 0 0
0 0 1 1 0
0 0 1 1 0
0 0 0 0 0
【工夫② 角では特別な処理が必要】
例えば四角い障害物では、4つのという角があります。そのため、どの分布をどの方向へ反射させるかを間違えると、角付近に不自然な速度が出ます。初心者の段階では、「壁より角の方が境界条件の実装が難しい」 と覚えておけば十分です。
┌─────┐ 角では、
│ │ ↗ → ↑
│ │ このような複数方向が同時に固体へ向かいます。
└─────┘
【工夫③ 曲面は「階段状」になってしまう】
例えば円柱を格子上に表現すると、本当は、円弧なのに、格子では、階段状の壁になりやすいです。これを、staircase boundaryなどと考えます。粗い格子では、本当の円柱と計算上の円柱が違ってしまうという誤差が生じます。
【工夫④ 格子を細かくする】
最も単純な対策は、「粗い格子から細かい格子」にすることです。すると曲面をより正確に表現できます。ただし、「格子数 → 計算量 [増加↑] → GPUメモリ使用量 [増加↑]」となります。これは、これまでXLBで扱ってきた大規模計算にも直結する問題です。
【工夫⑤ 曲面境界を補間する】
格子を細かくするだけでなく、「格子点から本当の壁までの距離を考えて反射させる」 という方法もあります。例えば、壁までの距離に応じて分布関数を補間する 方法があります。これにより、円柱や複雑形状をより正確に表現できます。ただし、計算はかなり複雑になります。
【工夫⑥ 入口・出口で「未知のf」を注意して扱う】
LBMでは、「境界の位置 → 一部のfが不足 → 他の物理量からfを復元」ということが起こります。したがって、「速度・密度・圧力」と、「$f_i$」の関係を使って補います。これが、「Zou–He・非平衡外挿・convective outlet」などの境界条件につながります。
【工夫⑦ 境界条件が「質量保存」を壊さないようにする】
例えば、「入口から100入った・出口から90しか出なかった」となれば、「内部で10増えてしまう」ことになります。その結果、$\rho$が徐々に増えてしまいます。したがって実際の計算では、質量保存が成立しているかを常に確認することが重要です。これは非常に重要なチェック項目です。例えば計算中に、平均密度が大きく変化したら、境界条件や計算条件に問題がある可能性があります。
【工夫⑧ 速度を大きくしすぎない】
LBMでは、一般に低Mach数で計算することが重要です。つまり、流体速度を格子音速に対して十分小さくします。非常に簡単に言えば、**「LBMでは、あまり大きな速度を設定しない」**ということです。初心者向けの例題なら、「u = 0.01」程度の小さな速度から始めることが多いです。
【工夫⑨ 粘性と $\tau$ の関係】
BGK-LBMでは、$\nu=c_s^2\left(\tau-\frac12\right)\Delta t$ という関係があります。つまり、「τ → 粘性 ν → Reynolds数」につながります。そして、$\omega=\frac{1}{\tau}$ です。したがって、「Collision → τ → 粘性 → Re」という関係があります。これはSTEP 7で説明したCollisionと、STEP 11の実際の流体解析をつなぐ重要な部分です。
16. ここで「壁+入口+出口」を全部組み合わせる
最終的には、例えば、以下のような問題になります。そしてLBMでは、① Collision → ② Streaming → ③ 壁ならBounce Back → ④ 入口なら速度条件 → ⑤ 出口なら圧力・流出条件 → ⑥ 次の時間、という処理を繰り返します。
┌────────────────────────┐壁
│ → → → → → → → → → → → →│
入口 │ → → → ■■■ → → → → →│ 出口
速度指定 │ → → → ■■■ → → → → →│ 圧力指定
│ → → → ■■■ → → → → →│
│ → → → → → → → → → → → →│
└────────────────────────┘壁
17. STEP 11の全体像
ここまでを一枚にまとめると、以下のようになります。
LBM
┌─────────┴─────────┐
流体領域 境界
│ ┌────────┼────────┐
Collision 壁 入口 出口
↓ │ │ │
Streaming Bounce 速度 圧力
↓ Back 指定 指定
境界で処理
↓
次の時間
STEP 11で特に重要な「3段階」
この段階では、全部を一度に覚えない方がよいです。
【レベル1:物理的な意味】 壁 → 流体を通さない・入口 → 流体を入れる・出口 → 流体を出す
【レベル2:LBMでの処理】 壁 → Bounce Back・入口 → fを速度に合わせて決める・出口 → fを流出条件に合わせる
【レベル3:数値計算上の工夫】 角の処理・曲面の処理・格子解像度・質量保存・低Mach数・安定性
STEP 10からSTEP 11へのつながり
STEP 10では、「D2Q9では、9方向へ分布関数が移動する」ことを学びました。するとSTEP 11では、「では、その9方向の分布が壁にぶつかったらどうする?」 となります。これがBounce Backです。つまり、「D2Q9 → 方向別の分布 → 壁との衝突 → Bounce Back」という流れで理解できます。
STEP 11からSTEP 12へ
最後に、次のSTEP 12「実際の流体」へ、「ここまでで、D2Q9・Collision・Streaming・境界条件がそろいました。では、この仕組みを使って、実際に流体を流してみましょう」とつなぎます。そこで初めて、「Poiseuille流 → 円柱周りの流れ → 抗力 CD → 揚力 CL → 渦放出 → Reynolds数による変化 」へ進むことができます。
この順番なら、「なぜXLBの境界条件にFullwayBounceBackBCやHalfwayBounceBackBCがあるのか」も、STEP 11まで来れば意味を持って理解できるようになります。
■ LBMの基礎の直観的理解【STEP 12 実際の流体を解析する】
これまでのSTEP 0~11をすべてつなげると、STEP 12「実際の流体」 は、単なる「LBMを動かしてみる」段階ではなく、「これまで学んだLBMの仕組みを使って、実際の流体現象を解析し、結果が物理的に正しいか確認する」 という最終段階にすると、きれいにまとまります。
1. STEP 0~11を振り返る
ここまでの学習を一本につなぐと、以下のようになり、STEP 12では、このすべてを使います。
STEP 0 流体の基礎 → STEP 1 1次元拡散方程式 → STEP 2 手計算 → STEP 3 Excel
→ STEP 4 Python → STEP 5 格子と時間発展 → STEP 6 分布関数 f → STEP 7 Collision
→ STEP 8 Streaming → STEP 9 D1Q3 → STEP 10 D2Q9 → STEP 11 境界条件 →
STEP 12 実際の流体
2. 最初は「流路を流れる流体」
いきなり複雑な形状にせず、まずは単純な流路から始めます。ここで、「入口で流体を入れる・壁では流体を止める・出口から流体を出す」という、STEP 11の境界条件を実際に使います。
┌─────────────────────────┐ 壁 \ 流速分布
│ → → → → → → → → → → → → │ \
入口 │ → → → → → → → → → → → → │ 出口 │ 中央
│ → → → → → → → → → → → → │ /
└─────────────────────────┘ 壁 /
3. 最初の実用例:Poiseuille流
最初の実際の流体問題として、Poiseuille流を使うのがおすすめです。例えば、上図のような流路です。壁では、$u=0$とします。すると中央付近の流速が大きく、壁付近では小さくなります。この速度分布が理論解と一致するかを確認します。ここがSTEP 12の最初の重要なポイントです。「LBMで計算できた」ではなく、「理論解と一致した」ことを確認する。
Poiseuille流 :ハーゲン・ポアズイユ流れ(ポアズイユ流れ)とは、一定の半径を持つ円管の中を、非圧縮性のニュートン流体が「定常層流」で流れるときの現象およびその解析解のことです。
4. 次に障害物を入れる
流れが正しく計算できることを確認したら、障害物を置きます。すると流体は障害物を避けて流れます。ここで、「D2Q9・Collision・Streaming・Bounce Back」が全部登場します。
┌─────────────────────────┐
│ → → → → → → → → → → → → │ → → → ↗ ↘ → → →
│ → → → → ■ → → → → → │ → → → → ■ → → → →
│ → → → → → → → → → → → → │ → → → ↘ ↗ → → →
└─────────────────────────┘
5. 円柱周りの流れ
次に、四角い障害物から円柱へ進みます。これが、LBMの代表的な教材問題です。なぜ円柱がよいかというと、流れの変化を目で確認しやすい からです。
→ → → → → → → → →
→ → → ○ → → →
→ → → → → → → → →
6. Reynolds数を変えてみる
ここでSTEP 0のReynolds数が生きてきます。$Re=\frac{UL}{\nu}$ 例えば、Re = 10, 50, 100, 200 と変化させます。すると、「低Re → 比較的穏やかな流れ・高Re → 後流が複雑になる」という変化を確認できます。
7. 渦放出を見る
Reynolds数がある程度高くなると、円柱の後ろに、以下のような交互の渦が現れます。これは、カルマン渦列です。ここで、「LBMで実際の非定常流れが再現できている」 ことを実感できます。
○ ↗ ○
→ → → ○ → → → → →
↘ ○ ↘ ○
8. 速度分布を見る
STEP 12では、単に流線を見るだけではなく、速度場を可視化します。例えば、「低速 ──→ 高速」を色で表示します。概念的には、以下のようにどこで流速が速く、どこで遅いか
→ → → → → → →
→ → ↗ ○ ↘ → →
→ ↗ ○ ↘ → →
→ → ↘ ○ ↗ → →
→ → → → → → →
9. 圧力分布を見る
次に、$p=c_s^2\rho$を使って圧力を調べます。例えば円柱の前面では圧力が高く、後ろでは低くなる、といった分布が現れます。これによって、**「なぜ物体に抗力が生じるのか」**という構造解析・CAEにも近い考え方へ進めます。
→ → → 高圧 ○ 低圧
10. 抗力 CD と揚力 CL
ここから本格的な工学解析になります。物体に働く力を、$C_D=\frac{F_D}{\frac12\rho U^2A}$ $C_L=\frac{F_L}{\frac12\rho U^2A}$のような無次元量で評価します。つまり、「流体解析 → 物体表面の力 → 抗力 FD揚力 FL → CD・CL」です。これによって、「流れがどうなっているか」 だけでなく、「その流れによって物体にどんな力が働くか」 まで評価できます。
11. 実験値・理論値と比較する
STEP 12で非常に重要なのが検証です。例えば、以下のように、「理論解・実験値・LBM計算値」を比較します。ここで初めて、「計算結果は正しいのか?」 を確認できます。
CL CD │ 理論解* ● 実験値
│ /
│ ○─── LBM
└────────────→ Re
12. 格子解像度の影響を見る
同じ問題を、「粗い格子:50 × 20」 → 「普通の格子:100 × 40」 → 「細かい格子:200 × 80」のように計算します。すると、「格子が粗い → 形状表現が粗い → 誤差が大きい・格子が細かい → 形状が正確 → 結果が改善」ということが分かります。これはSTEP 11で学んだ、「曲面・境界の表現」 につながります。
13. 時間刻みの影響も見る
同様に、$\Delta t$を変えることで、「安定性・精度・計算速度」の関係を学びます。つまり、数値解析で非常に重要な考え方に進みます。
精度 ↑
│ ●
│ ●
│ ●
└────────→ 計算コスト
14. STEP 12では「正しく計算できているか」を確認する
最終段階では、次のようなチェックをします。
【物理的なチェック】
・流体は壁を通り抜けていないか ・入口から流体が入っているか
・出口から流体が出ているか ・速度分布は妥当か ・圧力分布は妥当か
【数値的なチェック】
・密度が大きく変化していないか ・計算が発散していないか
・格子を細かくすると結果が収束するか ・時間を進めると適切な状態になるか
【工学的なチェック】
・理論値と一致するか ・実験値と一致するか ・CD・CLが妥当か
15. STEP 12のおすすめ例題
初心者向けの最終章としては、次の順番がよいと思います。この⑤までできれば、「LBMを学んだ」という実感がかなり得られると思います。
| 例題 | 学ぶこと |
|---|---|
| ① Poiseuille流 | 速度分布・理論解との比較 |
| ② 正方形障害物 | 境界条件・流れの回り込み |
| ③ 円柱周りの流れ | 圧力・速度・抗力 |
| ④ 円柱+Re数変更 | 流れの変化・渦 |
| ⑤ 円柱+渦放出 | 非定常流れ・(C_D,C_L) |
16. 最終的な学習体系
今回の学習ルートを最終的に整理すると、非常にきれいです。
STEP 0 【基礎】流体の基礎 STEP 1 【数学】1次元拡散方程式
STEP 2 【手で理解】手計算 STEP 3 【表計算】Excel
STEP 4 【プログラム】Python STEP 5 【数値解析の考え方】格子と時間発展
STEP 6 【LBM】分布関数 f STEP 7 Collision
STEP 8 Streaming STEP 9 D1Q3
STEP 10 D2Q9 STEP 11 境界条件
STEP 12 【実問題】実際の流体
┌────────────────────────────────────────────────┐
│ Poiseuille流 障害物周りの流れ 円柱周りの流れ │
│ Reynolds数 渦放出 CD・CL 理論・実験との比較 │
└────────────────────────────────────────────────┘
この学習ルートの「ゴール」
最終的には、「LBMのコードを実行できる」 ではなく、「なぜその計算をしているのかを説明でき、計算結果が物理的に妥当か判断できる」 ところまで持っていくのが理想です。
特に今回の構成では、**STEP 1~5が「数値計算の基礎」、STEP 6~10が「LBMの仕組み」、STEP 11が「実際の計算方法」、STEP 12が「工学的な流体解析」**という4つのまとまりになります。これは、初心者向けのLBM教材としてかなりきれいな構成になります。













