1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@ShibataRyoichi(Ryoichi Shibata)

LBM流体解析ツールXLBの活用【00A:LBMの基礎の直観的理解(STEP 0-6】

1
Last updated at Posted at 2026-10-12

この文書では、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の基礎の直観的理解】に関する目次を、以下にまとめます。ただし、解説に図解を入れて、詳しく解説するのですが、１記事のファイル制限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**を計算することを理解する。
  • 分布関数が「各方向へ進もうとしている流体の情報」を表していることを理解する。
  • 分布関数から、密度や速度を計算できることを学ぶ。

image.png

後半【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 0　流体の基礎】

STEP 0では、次の3つを理解します。ここでは、LBMの数式を覚える前に、「流体をどんな量で表すのか」 を理解する段階です

① 速度・圧力・密度・粘性　　② 流体で起こる基本現象　　③ Reynolds数

LBMでは最終的に、分布関数 $f_i$ から密度・速度などの物理量を取り出します。そのため、STEP 0でこれらの意味をしっかり理解しておくと、後のLBMがかなり分かりやすくなります。

1. まず「流体」とは何か

流体は、簡単に言えば、形を保ったままではなく、流れることができる物質 です。代表例は、「液体：水、油・気体：空気」です。例えば水を容器に入れると、液体は容器の形に合わせて水の形が変わります。一方、固体は基本的には自分自身の形をある程度保ちます。LBMは、この流体の運動を計算する方法です。

2. 速度 u

最初に理解するのは速度です。例えば、水が左から右へ流れているとします。この「どの方向へ、どのくらい速く動いているか」を表すのが速度です。1次元なら、$u$ 、2次元なら、$\mathbf{u}=(u_x,u_y)$、3次元なら、$\mathbf{u}=(u_x,u_y,u_z)$ です。なお、「速度」は大きさだけでなく方向も持つことが重要です。単位は、$\mathrm{m/s}$ です。

3. 圧力 p

圧力は、単位面積あたりに働く力 です。基本的には、$p=\frac{F}{A}$ と考えます。単位は、$\mathrm{Pa}$ （パスカル）です。流体では、圧力差があると流れが生じます。例えば、「高圧 → → 低圧」のように、圧力の高い側から低い側へ流れやすくなります。パイプなら、「高圧：p₁ → → 低圧：p₂」です。LBMでも圧力は重要で、後で $p=c_s^2\rho$ という関係が出てきます。

4. 密度 ρ

密度は、単位体積あたりにどれだけ質量があるか です。$\rho=\frac{m}{V}$ と考えます。単位は、$\mathrm{kg/m^3}$ です。例えば、同じ体積の中に、物質がぎっしり詰まっていれば密度が高くなります。LBMでは、密度は後で $\rho=\sum_i f_i$ として、分布関数 (f_i) を足し合わせて求めることになります。
したがって、「STEP 0：密度とは何か　→ STEP 6：f を足すと密度になる」というつながりがあります。

5. 粘性とは何か

ここはLBMを学ぶ上で特に重要です。粘性は簡単に言えば、「流体の流れにくさ」 です。例えば、「水 → 比較的さらさら・はちみつ → どろっとしている」という違いがあります。ただし、厳密には「粘性＝流れにくさ」だけではなく、流体内部で速度差があるとき、その差をならそうとする性質 と考えると、LBMにつながりやすくなります。

6. 粘性を速度分布で考える

例えば、以下のような、上下に壁がある流路を考えます。壁の近くでは流体と壁面との摩擦により速度が小さく、中央では速度が大きくなります。このような速度差をならそうとする作用に、粘性が関係しています。

【STEP 0 から 6 のまとめ】

image.png

7. 動粘性係数 ν

LBMでは、「粘性」という言葉だけでなく、動粘性係数 $\nu$ が非常によく登場します。動粘性係数は、$\nu=\frac{\mu}{\rho}$ です。ここで、$\mu$：動粘度ではなく粘性係数（動粘性係数ではない）＝粘性係数・dynamic viscosity、$\rho$：密度、$\nu$：動粘性係数・kinematic viscosity です。
単位は、$\mathrm{m^2/s}$ です。LBMでは後で、$\nu=c_s^2\left(\tau-\frac12\right)\Delta t$ が出てくるので、粘性と $\tau$ の関係につながります。

8. 「流れる」という現象

流体の基本現象として、まず「流れる」を見ます。これは速度 $u$ を持つ流体が移動している状態です。LBMでは、後でこの流れを、「分布関数　→　Streaming　→　隣の格子へ移動」という形で表現します。つまり、「流体が流れる」→「分布関数が方向に沿って運ばれる」 というつながりがあります。

9. 「拡散する」という現象

次に拡散です。例えば温度が中央だけ高いと、例えば、最初は中央が「温度100度 ● 」ならば、時間が経つと「左横20度・中央60度・右側20度」のように周囲へ広がります。これはSTEP 1で学んだ1次元拡散方程式につながります。重要なのは、「高いところから低いところへ広がり、分布がなだらかになる」 ということです。また流体力学では、温度や物質だけでなく、速度の拡散＝運動量の拡散という考え方もあり、これが粘性と深く関係します。

10. 「壁にぶつかる」という現象

流体が壁に近づくと、「流体→ → → ＊█████壁面」のように、流体は壁面に＊の位置で衝突します。物理的には、流体は壁を通り抜けません。さらに粘性流体では、通常、壁面において$u=0$というすべりなし条件を考えます。これが後のSTEP 11で、Bounce Back として登場します。
つまり、「STEP 0：壁では流体が通り抜けない　→　STEP 11：Bounce Backで実現する」という関係です。

11. ここまでの4つを整理する

この4つがSTEP 0の中心です。

物理量 意味 単位 LBMとの関係
速度 $u$ 流体の速さ・方向 m/s $f_i$から求める
圧力 $p$ 面積あたりの力 Pa 密度と関係
密度 $\rho$ 単位体積あたりの質量 kg/m³ $f_i$の総和
粘性 $\nu,\mu$ 流れの変化をならす性質 m²/s, Pa·s $\tau$と関係

12. ここで「流れの状態」を考える

次に、「同じ速度でも、流れ方は違うのでは？」という話に進みます。例えば、以下のような「ゆっくりした流れ」と「複雑な流れ」では、流れ方が違います。この流れの特徴を表す代表的な指標が、Reynolds数 です。

ゆっくりした流れ　　　複雑な流れ
→ → → → → → →　　　→ → ↗ → ↘ →
→ → → → → → →　　　→ ↘ → ↗ → →
→ → → → → → →　　　→ → ↘ → ↗ →

【STEP 7 から 12 のまとめ】
image.png

13. Reynolds数

Reynolds数は、$Re=\frac{UL}{\nu}$ です。ここで、$U$：代表速度、$L$：代表長さ、$\nu$：動粘性係数です。Reynolds数は無次元量です。

14. Reynolds数の意味

Reynolds数は簡単に言うと、「流れの慣性の影響」と「粘性の影響」の強さの比 です。概念的には、$Re\sim\frac{\text{慣性の影響}}{\text{粘性の影響}}$ と考えられます。したがって、「Reが小さい：粘性の影響 ＞ 慣性の影響」で比較的穏やかな流れになります。一方、「Reが大きい：慣性の影響 ＞ 粘性の影響」では、流れが複雑になりやすくなります。

15. 例えば円柱の後ろを見る

円柱を流れの中に置きます。低いReでは比較的穏やかです。Reが高くなると渦放出が現れるようになります。STEP 12で扱う、カルマン渦列につながります。

　　　　　　　　　　　低いReでは　　　　  Reが高くなると
→ →  ○  → →　　　　 → →  ○  → → 　　　→ →  ○
　　　　　　　　　       ↘   ↗　　　　　　　↘  → ↘
                       ↘ ↗　　　　　　　　　↖ ↙︎

16. Reynolds数は何を変えるのか

式を見ると、$Re=\frac{UL}{\nu}$なので、「速度 $U$ を大きくすると：Re ↑大きく」「代表長さ $L$ を大きくすると：Re ↑大きく」「動粘性係数 $\nu$ を小さくすると：Re ↑大きく」、のようになります。つまり、「速度を変える → Reが変わる・粘性を変える → Reが変わる・物体サイズを変える → Reが変わる」ということです。

17. LBM学習とのつながり

STEP 0で学んだものが、その後どこに使われるかを整理すると非常に分かりやすいです。「速度 u → STEP 6　分布関数から速度を求める」「密度 ρ → STEP 6　fの総和から求める」「粘性 ν → STEP 7　Collision → 緩和時間 τ」「壁 → STEP 11　Bounce Back」「Reynolds数 → STEP 12　実際の流れ」
つまりSTEP 0は、後のすべての基礎になっています。

18. STEP 0 → STEP 1への接続

ここで「流体の基礎」から「拡散方程式」へつなぎます。STEP 0で、「流体には、速度・圧力・密度・粘性がある」 ことを学びました。しかし、「物理量が時間とともに変化する」 という現象を、まだ数式では扱っていません。そこでSTEP 1で、$\frac{\partial T}{\partial t}=\alpha\frac{\partial^2T}{\partial x^2}$ という最も簡単な拡散方程式を学びます。
つまり、「STEP 0：流体にはどんな物理量があるか → STEP 1：物理量が時間・空間でどう変化するか」という接続です。

【STEP 13 から 18 のまとめ】
image.png

STEP 0の到達目標

STEP 0では、細かい流体力学を覚える必要はありません。次の3つが説明できれば十分です。

【① 流体を表す量】$u,\quad p,\quad \rho,\quad \nu$：それぞれが何を意味するか。
【② 基本的な現象】「流れる・拡散する・壁にぶつかる」：それぞれ何が起こっているか。
【③ Reynolds数】$Re=\frac{UL}{\nu}$ ：「速度・大きさ・粘性によって、流れ方が変わる」 と理解する。
今回の学習ルートでは、このSTEP 0が「物理の言葉」を学ぶ段階です。その後、「STEP 0：流体の言葉を知る → STEP 1：拡散方程式で数学的に表す → STEP 2：手で計算する → STEP 3：Excelで計算する → STEP 4：Pythonで計算する」と進めると、「物理 → 数学 → 数値計算 → プログラム」 という流れがきれいにできます。

■ LBMの基礎の直観的理解【STEP 1 拡散方程式を知る】

STEP 1　1次元拡散方程式を知る — 15項目の詳しい解説

STEP 1では、手計算・Excel・Pythonに進む前に、拡散という物理現象と、それを表す数式の意味を理解します。この段階では、微分方程式を解く必要はありません。まず、細い棒の中央を熱したとき、なぜ中央の温度が下がり、左右に熱が広がるのかを考えます。

1. 「拡散」とは何か

拡散とは、温度や物質の濃度などの分布が、時間とともに空間的に広がり、なだらかになる現象です。例えば、細い棒の中央だけを熱すると、中央の熱が周囲へ伝わります。

初期状態：中央だけ高温　→  時間経過：熱が左右へ広がる
              100
               ●
               │
───●────●────●────●────●───
   低温   低温   高温   低温   低温

このように、ある場所だけに集中していた温度が周囲に広がる現象が熱拡散です。拡散は熱だけでなく、インクが水中に広がる現象や、物質の濃度が均一になっていく現象にも見られます。このSTEPでは、熱の拡散を例にして学びます。

2. なぜ「1次元」で考えるのか

実際の熱の伝わり方は、上下・左右など複数の方向に広がる3次元現象です。しかし、最初から3次元を考えると、計算も説明も複雑になります。そこで、細く長い棒を考え、棒の長さ方向だけに温度が変化すると仮定します。

 左端                           右端
 x=0                           x=L
  │                             │
  ●────●────●────●────●────●────●
              棒の長さ方向 x →

このように、変化する空間方向を1つに限定したモデルを1次元モデルと呼びます。1次元モデルでは、温度は位置( $x$) と時間 ($t$) の関数として表せます。$T=T(x,t)$ これは、温度 $T$ が場所と時間によって変わることを意味します。

3. 最初の温度分布を考える

最初の状態を、中央だけ温度が高い状態にします。例えば、次のように考えます。

位置 $x$ 0 1 2 3 4
初期温度 0 0 100 0 0

ここでの数値は、分かりやすくするための温度の例です。単位は任意の温度単位とします。グラフにすると、中央だけが高い山のような形になります。この初期状態から、時間とともに温度分布がどのように変わるかを考えます。

温度 T
100          ●
             │
             │
  0  ●──●────┴────●──● → 位置 x
       0    1   2   3   4

4. 時間が経つと温度はどう変化するか

中央が高温で、隣接する場所が低温なら、熱は中央から周囲へ伝わります。その結果、中央の温度は下がり、左右の温度は上がります。例えば、温度分布の変化を概念的に表すと、次のようになります。

初期状態
0    0   100    0    0
       ↓ 拡散
少し後
0   20    60   20    0
       ↓ さらに時間経過
さらに後
10  20    40   20   10

これは現象を理解するための模式例であり、特定の拡散方程式を数値的に解いた結果ではありません。重要なのは、中央の高い温度が下がり、その周囲の温度が上がるという変化です。

5. 1次元拡散方程式を知る

このような温度変化を表す代表的な式が、1次元の拡散方程式です。$\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}$ この式は、温度 $T$ が、位置 $x$ と時間 $t$ に応じて、どのように変化するかを表しています。各記号の意味は次のとおりです。

記号 意味
$T$ 温度
$t$ 時間
$x$ 位置
$\alpha$ 熱拡散率
$\partial T/\partial t$ 時間による温度変化
$\partial^2T/\partial x^2$ 空間方向の温度分布の曲がり具合

最初は、この式を解く必要はありません。左右の温度差などによって生じる温度変化を、数式で表したものと理解してください。なお、この式は、熱拡散率が一定で、内部発熱がない場合などの基本的なモデルです。

6. 左辺は「時間による温度変化」

方程式の左辺は、$ \frac{\partial T}{\partial t} $です。これは、同じ場所において、温度が時間とともにどれくらい変化するかを表します。例えば、ある位置の温度が、「時間 t=0：温度100　→ 時間 t=1：温度80」となった場合、温度は時間とともに低下しています。

一定の短い時間の変化を考えると、温度変化の速さは概念的に、$\frac{\Delta T}{\Delta t}$ で表せます。ここで $\Delta T$ は温度の変化量、$\Delta t$ は時間の変化量です。微分の記号 $\partial T/\partial t$ は、この考え方を、時間を限りなく小さくした場合へ拡張したものです。

7. 右辺は「温度分布の曲がり具合」

方程式の右辺は、$\alpha\frac{\partial^2T}{\partial x^2}$ です。このうち、特に重要なのが、$\frac{\partial^2T}{\partial x^2}$ です。これは、温度分布が空間方向にどのように曲がっているかを表します。例えば、中央だけ温度が高い場合、グラフは山の形になります。

山の頂上付近では、温度分布の曲がり具合は負になります。そのため、拡散方程式ではその場所の温度が下がる方向に変化します。反対に、周囲より低温の谷では、曲がり具合が正になり、温度が上がる方向に変化する場合があります。このように、温度分布の曲がり具合が、温度の時間変化に関係しているのです。

この例は、「温度分布の曲がり具合が、温度の時間変化を決める」 ということを理解するのに適しています。この状態を、具体的に$T(x)=100-x^2$ を使って、順番に確認すると、以下のという関係になります。
①の山型の温度分布は、②の一階微分で変化率（傾き）を見ると一様に下がっており温度が下がる方向（-2x）に変化し、③の傾きの変化率はこの場合は一定（-2）です。

$\underbrace{T(x)=100-x^2}_{\text{温度分布}}
\quad\Rightarrow\quad
\underbrace{T'(x)=-2x}_{\text{温度が下がる方向}}
\quad\Rightarrow\quad
\underbrace{T''(x)=-2}_{\text{下向きの曲がり}}
$

image.png

以上の、１次元拡散方程式の説明「STEP 1 から 7」について、以下に図示します。
【STEP 1 から 7 のまとめ】
image.png

8. 「山が崩れる」と考える

拡散方程式の理解では、温度分布を山や谷に例えると分かりやすくなります。

初期     　→拡散→          少し後  　 →拡散→          さらに後
           /\ 
          /  \                    /────\                /──────────\
─────────/────\───────── ────────/      \──────── ─────/            \─────

時間が経つと、山は低くなり、幅が広がります。数学的には、例えば中央の高温部分で、$\frac{\partial^2T}{\partial x^2}<0$ となることが、温度低下の方向を説明します。逆に、温度分布の谷では、$\frac{\partial^2T}{\partial x^2}>0$ となることがあり、温度は上昇する方向に変化します。
ただし、温度分布のすべての場所で同じ変化が起こるわけではありません。温度の上昇・下降は、その場所の温度分布の形によって決まります。

9. 熱拡散率 α とは何か

式の中の $\alpha$ は、熱拡散率です。熱拡散率は、材料の中で熱がどれくらい速く広がるかに関係する物性値です。単位は、$[\alpha]=\mathrm{m^2/s}$ です。一般に、同じ温度分布、同じ時間、同じ長さの条件で比較すると、熱拡散率が大きい材料ほど、温度分布が速く変化します。
例えば、「熱拡散率が小さい → 温度分布がゆっくり変わる・熱拡散率が大きい → 温度分布が速く変わる」と理解できます。

10. α を変えると温度分布はどうなるか

同じ初期温度分布を考えます。$T(x,0)=T_0(x)$ この初期状態を変えずに、熱拡散率 $\alpha$ だけを変えます。
熱拡散率が小さい場合は、ある時間が経過しても山が比較的残ります。熱拡散率が大きい場合は、同じ時間で山がより低く、広くなります。

α が小さい                α が大きい
          /\
         /  \                  /────────\
────────/────\────────  ──────/          \──────

このように、熱拡散率は、温度分布が時間とともに変化する速さを決めます。ここで「熱が広がる速さ」を、単純に熱が一定の速度で移動することと混同しないようにしましょう。拡散は、温度分布が時間とともに変化する現象です。

11. 空間 x の意味

拡散方程式では、位置 $x$ が重要です。$T=T(x,t)$ の $x$ は、棒のどの位置かを表します。

空間 x →
0────1────2────3────4────5
●    ●    ●    ●    ●    ●

例えば、$x=0$ は棒の左端、$x=5$ はその右側の位置を表すものとします。ここでの位置の数値は説明用です。STEP 1では、棒は連続した空間として考えます。まだ $x=0,1,2,\ldots$ のような格子点で計算する必要はありません。

12. 時間 t と「時間発展」の意味

時間 $t$ は、現象が始まってからの経過時間を表します。

t = 0  →  t = Δt  →  t = 2Δt  →  t = 3Δt
 初期       少し後       さらに後      さらに後

例えば、最初に温度100だった位置が、時間が経つにつれて80、60と変化することがあります。このように、ある時刻の状態から、その後の状態へ変化していくことを時間発展と呼びます。拡散方程式は、位置 \(x\) における温度の時間変化を計算するための式です。

13. STEP 1では連続した空間として考える

この段階では、棒を細かい格子に分けません。実際の棒は連続した空間であり、温度も位置に応じて連続的に変化すると考えます。

連続した棒
──────────────────────────

温度分布
       ／＼
      ／  ＼
─────／────＼──────────→ x

そのため、拡散方程式では、位置についての微分を使います。$ \frac{\partial T}{\partial t} = \alpha\frac{\partial^2T}{\partial x^2} $ ここで使われる微分は、温度の変化を連続的に表現するための数学的な道具です。格子点に分割して数値を計算するのはSTEP 2以降です。

14. 初期条件・境界条件と、差分式との違い

拡散方程式を実際に解くには、方程式に加えて、最初の温度分布や棒の端の条件を決める必要があります。初期条件の例は、$ T(x,0)=T_0(x) $です。これは「時間 $t=0$ において、各位置の温度がどうなっているか」を指定します。
境界条件の例としては、棒の両端の温度を一定に保つ条件などがあります。一方、STEP 2以降では、微分方程式を格子点で計算できる式に変換します。例えば、$ T_i^{n+1} = T_i^n+ r\left(T_{i-1}^n-2T_i^n+T_{i+1}^n\right) $です。
ここで、$i$：格子点の位置番号、$n$：時間ステップ番号、$r=\alpha\Delta t/\Delta x^2$：計算条件から決まる係数、を表します。この式は、拡散方程式を数値計算するための差分式です。

ただし、この差分式はSTEP 1では使いません。STEP 1では連続体としての拡散方程式を理解し、STEP 2で初めて差分式を使って手計算します。

【STEP 8 から 14 のまとめ】
image.png

15. STEP 1から手計算、Excel、Python、LBMへ

STEP 1の目的は、拡散方程式を使って物理現象を表す考え方を理解することです。その後、計算の具体的な方法へ進みます。

STEP 1：拡散方程式を知る$ \frac{\partial T}{\partial t} = \alpha\frac{\partial^2T}{\partial x^2} $：温度分布が時間とともにどう変化するか理解します。
STEP 2：手計算：格子点を使い、中央と左右の温度を1ステップずつ計算します。
STEP 3：Excel：手計算の式をセルに入力し、時間発展を自動計算してグラフ化します。
STEP 4：Python：配列と繰り返し処理で、Excelと同じ計算をプログラムにします。
STEP 5：格子と時間発展を理解する：空間を格子に分け、現在の状態から次の状態を計算するという数値解析の基本を整理します。
STEP 6以降：LBMを学ぶ

温度 $T$ を直接更新する考え方から、方向を持つ分布関数 $f_i$ を計算する考え方へ進みます。CollisionとStreamingを理解し、D1Q3、D2Q9へと発展させます。

STEP 1の到達目標

STEP 1を終えた段階で、次の3点を自分の言葉で説明できれば十分です。

  1. 拡散とは何か。 温度分布などが時間とともになだらかになる現象である。
  2. 拡散方程式は何を表すか。 温度の時間変化と、空間方向の温度分布の曲がり具合を関係づける式である。
  3. なぜ次に手計算するのか。 連続的な方程式を格子点で計算できる形に変換し、時間ステップごとに温度が変わることを確かめるためである。

以上を理解すると、STEP 2の手計算は単なる数値計算ではなく、拡散方程式で表した物理現象を、格子点の数値として確かめる作業だと理解できるようになります。

■ LBMの基礎の直観的理解【STEP 2 手で計算する】

STEP 2は、ExcelやPythonを使わず、自分の手で「1回の時間ステップ」を計算する段階です。ここで大切なのは、計算量を増やすことではなく、「ある格子点の値が、左右の格子点の値からどう変化するのか」 を実際に計算して理解することです。

1. まず1次元の格子を考える

位置を0～20の21個の格子点に分けます。今回は、格子間隔： $Delta x=1$、時間刻み：$\Delta t=0.1$、拡散係数：$\alpha=1$、とします。

位置
 0   1   2   3   4   5   6   7   8   9  10  11  12  13 ... 20
 ●───●───●───●───●───●───●───●───●───●───●───●───●───・・・●

2. 最初の状態を決める

最初の時刻 $t=0$ では、位置10だけ100 とします。

位置   0  1  2  3  4  5  6  7  8  9  10 11 12 13 ... 20
温度   0  0  0  0  0  0  0  0  0  0 100  0  0  0  ...  0

イメージすると、以下の状態です。

温度
100 │                    ●
    │                   / \
    │                  /   \
  0 └──────────────────────────→ 位置
                         10

3. 拡散方程式を使う

1次元拡散方程式は、$\frac{\partial T}{\partial t}=\alpha\frac{\partial^2T}{\partial x^2}$ です。しかし、このままでは手計算しにくいので、格子点で計算できる形にします。$T_i^{n+1}=T_i^n+r\left(T_{i-1}^n-2T_i^n+T_{i+1}^n\right)$
ここで、$r=\frac{\alpha\Delta t}{\Delta x^2}$ です。今回は、$r=\frac{1.0\times0.1}{1.0^2}=0.1$ です。
したがって、今回の計算式は、$T_i^{n+1}=T_i^n+0.1(T_{i-1}^n-2T_i^n+T_{i+1}^n)$ になります。

【差分式の直感的な理解】
image.png

4. まず「位置10」を計算する

中央の位置10を考えます。

位置      9       10       11
         ↓        ↓        ↓
温度      0       100       0

式に入れます。$T_{10}^{1}=100+0.1(0-2\times100+0)$ まず、$0-200+0=-200$ なので、$T_{10}^{1}=100+0.1\times(-200)=100-20=80$ となります。つまり、中央の100が80に下がりました。

5. 次に「位置9」を計算する

位置9を見ます。

位置      8       9       10
         ↓       ↓        ↓
温度      0       0       100

式は、$T_9^1=T_9^0+0.1(T_8^0-2T_9^0+T_{10}^0)$ です。値を入れると、$=0+0.1(0-2\times0+100)=0+10=10$ です。つまり位置9には、10の温度が現れました。

6. 位置11も計算する

位置11は左右対称なので、以下になります。

位置      10      11      12
          ↓       ↓       ↓
温度     100       0       0

$T_{11}^1=0+0.1(100-2\times0+0)$ したがって、$10$ です。

7. 1回目の計算結果

したがって、最初の時間ステップが終わると、以下の状態になります。

位置       8   9   10  11  12
初期       0   0  100   0   0
             ↓ 計算
1ステップ   0  10   80  10   0

全体では、以下の状態です。これが**「拡散した最初の状態」**です。

位置   0 1 2 ... 8  9  10  11 12 ... 20
温度   0 0 0 ... 0 10  80  10  0  ... 0

8. ここで重要なこと

この計算から、最初は「位置10●のみ値100」が、1ステップ後は「位置9●値10・位置10●値80・位置11●値10」となりました。つまり、高温の場所から、隣の場所へ熱が広がった ことが分かります。これが「拡散」です。

9. もう1ステップ計算してみる

次は、現在の状態から、さらに1ステップ進めます。

位置      8   9   10  11  12
温度      0  10   80  10   0

【位置10】：$T_{10}^{2}=80+0.1(10-2\times80+10)=80+0.1(10-160+10)=66$
【位置09】：$T_9^2=10+0.1(0-2\times10+80)=10+0.1(60)=16$
【位置08】：$T_8^2=0+0.1(0-2\times0+10)=1$
同様に対称となる右側も入れると、以下のようになります。

位置      8   9   10  11  12
温度      1  16   66  16   1

10. さらに進める

同じ計算を繰り返します。このように、中央の温度が下がり、左右に広がっていくことが確認できます。

時間0：       位置10 = 100
時間1：    0     10    80     10     0
時間2：    1     16    66     16     1
時間3：  2.4   19.5    56   19.5   2.4
時間4： 3.88  21.44  48.7  21.44  3.88

11. 手計算では全部を計算しなくてよい

ここは初心者向け教材として重要です。21個全部を毎回計算する必要はありません。最初は、「位置 9・位置10・位置11」の3点だけを計算すれば十分です。なぜなら、最初の状態では他の場所はすべて0だからです。つまり、「 位置9：0・位置10：100・位置11：０」です。その後、「位置8～12」まで広がってきます。したがって、手計算ではまず中央付近だけを見ると理解しやすいです。

12. 「1回の計算」の意味

このSTEP 2で一番理解してほしいのは、現在の「左・中央・右」から、次の時間の「中央」を計算する だけです。

ここまでの手計算の手順の全体を、図示します。
【差分式による計算の具体例】
image.png

14. さらにLBMにつながる部分

そして、この「手計算」はLBM学習にかなり重要です。今回、以下の手順を進めました。

現在の状態　→ 周囲の値を見る　→ 次の状態を計算　→ 時間を1つ進める

対応して、LBMでも、以下のように時間を1ステップずつ進めます。

現在の分布関数 f　→ Collision　→ Streaming　→ 次の時間の分布関数

したがって、STEP 2ではLBMそのものをまだ勉強しませんが、「格子を使って、時間を1ステップずつ進める」 という、LBMの基本的な数値計算の考え方を体験しています。

STEP 2の到達目標

この段階では、次の3つが分かれば十分です。

  1. 空間を格子に分ける
  2. 左右の値を使って、次の時間の値を計算する
  3. これを繰り返すと、温度が広がっていく

ここまで理解できれば、次のSTEP 3のExcelはかなり自然に理解できます。さらにそのExcelをPythonに置き換えたSTEP 4が、「数値計算をプログラムにする」練習になります。

■ LBMの基礎の直観的理解【STEP 3 Excelで計算する】

STEP 0～2で学んだ内容をもとに、今回のSTEP 3では、完成したExcelファイルを使って、手計算で理解した拡散現象を、Excelで自動的に繰り返し計算する方法を学びます。

完成版Excelファイル：　STEP3_1D_diffusion_Excel_final.xlsx

このファイルをPCにダウンロードして、Excelなどの表計算ソフトで開きます。「1D Diffusion」と「学習ポイント」の2つのシートがあります。最初は「1D Diffusion」シートを開いて、表の数値とグラフを確認してください。

1. まず、STEP 0～3のつながりを理解する

ここまでの学習を振り返ると、それぞれのSTEPには役割があります。

【STEP 0　物理現象を知る】　流体の速度・圧力・密度・粘性などを学びます。また、熱や物質が広がる「拡散」という現象を知ります。
【STEP 1　拡散現象を数式で表す】　$\frac{\partial T}{\partial t}=\alpha\frac{\partial^2T}{\partial x^2}$ 温度が位置や時間によってどう変化するかを学びます。
【STEP 2　差分式を手で計算する】　中央の温度100が80になり、左右に温度10が現れることを、実際に数値を代入して確かめます。
【### STEP 3　Excelで自動計算する】　STEP 2で手計算した式をExcelのセルに入力し、時間を1ステップずつ進めます。これを何度も繰り返して、温度分布の変化を確認します。

つまり、STEP 3では新しい物理法則を学ぶのではありません。STEP 2の手計算をExcelで自動化することが中心です。

2. Excelで何を計算するのか

最初に、完成版Excelの計算条件を確認します。この条件では、20ステップ計算すると、最終時刻は $t=20\times0.1=2.0$ となります。今回の値は、数値計算の仕組みを理解するための教材用設定です。実際の材料の熱物性値を再現することが目的ではありません。

計算条件 設定値 意味
格子点数 (N) 21 位置0～20の21点
時間ステップ数 20 20回計算する
格子間隔 (\Delta x) 1.0 隣り合う格子点の間隔
時間刻み (\Delta t) 0.1 1回で進む時間
拡散係数 (\alpha) 1.0 今回使用する拡散の係数
係数 (r) 0.1 1回の計算での変化の大きさ

3. なぜ格子点を21個にするのか

完成版Excelでは、空間を位置0～20の21個の格子点で表します。横方向は空間の位置です。各格子点に温度 $T$ を与え、その温度が時間とともに変わるようにします。ここで、格子点は、連続した棒の温度を数値で表すために選んだ位置です。

STEP 1では連続した棒を考えましたが、STEP 2では手計算できるように格子点に分け、STEP 3ではその格子点をExcelの列に対応させています。

位置 i                         ↓中央位置10
0──1──2──3──4──5──6──7──8──9──10──11──12──…──20
●  ●  ●  ●  ●  ●  ●  ●  ●  ●   ●   ●   ●       ●

4. 最初の温度をExcelに入力する

今回の初期条件は、中央の位置10だけ高温で温度100とし、その他はすべて0です。

位置 i    0   1   2   3  …   9   10   11   12  …  20
初期温度   0   0   0   0  …   0  100    0    0  …   0

Excelでは「1D Diffusion」シートの15行目が時間ステップ0の状態です。この行では、「B列～V列：位置0～20の温度、L列：中央の位置10、L15：初期温度100」です。したがって、初期状態では L15が100、その他の位置が0 になっています。この「中央だけ100」という条件が、これから広がる温度分布の出発点になります。

5. STEP 2で学んだ差分式を使う

STEP 2では、次の差分式を手計算しました。$T_i^{n+1}=T_i^n+r\left(T_{i-1}^n-2T_i^n+T_{i+1}^n\right)$ この式は、次のように読みます。「次の温度 ＝ 現在の温度 ＋ 周囲との温度差による変化」記号の意味を整理します。この式の重要な点は、ある格子点の次の温度を、現在の左・中央・右の3点から計算することです。

記号 意味
$T_i^n$ 時間ステップ (n) の位置 (i) の温度
$T_i^{n+1}$ 次の時間ステップの温度
$T_{i-1}^n$ 現在の左隣の温度
$T_{i+1}^n$ 現在の右隣の温度
$r$ 温度変化の大きさを決める係数

6. まず中央の位置10を計算する

時間ステップ0では、中央の温度が100、左右が0です。ここで (r=0.1) を差分式に代入します。

位置         9        10        11
温度         0        100        0

ExcelのL行において、$T_{10}^{1}=100+0.1(0-2\times100+0)$ かっこの中は、$0-200+0=-200$ なので、$T_{10}^{1}=100+0.1(-200)=80$ となります。中央の温度は、100から80に下がりました。完成版Excelでは、これは セルL16 に対応します。実際の数式は、「=L15+$B$11*(K15-2*L15+M15)」です。
この数式では、L15が中央の現在温度、K15が左隣、M15が右隣、B11が係数 (r) を表します。

7. 次に位置9と位置11を計算する

中央の温度が下がるだけではありません。左右の隣接位置にも温度が現れます。
【位置9】：現在の左・中央・右の温度は、次のとおりです。

位置          8        9        10
温度          0        0        100

差分式に代入します。$T_9^1=0+0.1(0-2\times0+100)=10$
【位置11】：同様に、

位置          10       11       12
温度          100       0        0

なので、$T_{11}^1=0+0.1(100-2\times0+0)=10$ です。

1ステップ後の中央付近は、次のようになります。

位置 8 9 10 11 12
時間ステップ0 0 0 100 0 0
時間ステップ1 0 10 80 10 0

これがSTEP 2で手計算した内容です。Excelでは、この計算を各格子点に対して一括して行います。

8. Excelが自動計算する仕組み

Excelで最も注目してほしいのは、数式を次の時間ステップの行へ適用していることです。「1D Diffusion」シートでは、次のように表を構成しています。

Excelの行 時間ステップ 時刻
15行目 0 0.0
16行目 1 0.1
17行目 2 0.2
18行目 3 0.3
19行目 4 0.4
… … …
35行目 20 2.0

横方向は空間、縦方向は時間です。セルL16にある数式は、L15・K15・M15を参照して中央の次の温度を計算します。次の行のL17はL16・K16・M16を参照します。

つまり、「時間ステップ0：15行目　→ 時間ステップ1：16行目　→ …　→　時間ステップ20：35行目」というように、直前の時間ステップの温度を使って、次の時間ステップを計算する仕組みです。また、位置0と位置20は、各時間ステップで温度0に固定しています。これは両端の温度を一定にする境界条件です。

9. もう1ステップ進めるとどうなるか

時間ステップ1の状態は、$[0,\ 10,\ 80,\ 10,\ 0]$ です。ここから時間ステップ2を計算します。
【中央の位置10】：$T_{10}^{2}=80+0.1(10-2\times80+10)=80+0.1(-140)=66$
【位置9】：$T_9^2=10+0.1(0-2\times10+80)=16$
【位置8】：$T_8^2=0+0.1(0-2\times0+10)=1$
左右対称なので、右側も同じ値になります。

位置 8 9 10 11 12
時間ステップ1 0 10 80 10 0
時間ステップ2 1 16 66 16 1

時間ステップ2では、熱の影響が位置8と位置12にも広がっています。

10. Excelの表から分かる温度の変化

完成版Excelでは、時間ステップ0～20までの計算結果が表になっています。中央付近を抜き出すと、次のようになります。

時間ステップ 位置8 位置9 位置10 位置11 位置12
0 0 0 100 0 0
1 0 10 80 10 0
2 1 16 66 16 1
3 2.4 19.5 56 19.5 2.4
4 3.88 21.44 48.7 21.44 3.88

image.png

この表から、次の2つの変化が分かります。
中央の温度は下がる：$100\rightarrow80\rightarrow66\rightarrow56\rightarrow48.7$
左右の温度は上がり、より外側へ広がる：例えば位置9では、$0\rightarrow10\rightarrow16\rightarrow19.5\rightarrow21.44$ と上昇します。

この結果はSTEP 1で学んだ「温度分布の山が低くなり、幅が広がる」という現象に対応しています。

11. 最後にグラフを確認する

完成版Excelのグラフでは、次の時刻の温度分布を重ねて表示しています。$t=0,\quad0.5,\quad1.0,\quad1.5,\quad2.0$ 時間刻みが $\Delta t=0.1$ なので、それぞれ時間ステップ0、5、10、15、20に対応します。グラフの横軸は位置0～20、縦軸は温度 $T$ です。

最初は中央の位置10に高さ100の山があります。時間が進むと、中央の山が低くなり、左右へ広がっていきます。ここで大切なのは、時間の違う温度分布を、同じ位置軸で比較できることです。表だけでは分かりにくい変化も、グラフなら直感的に理解できます。
image.png
なお、Excelのグラフでは、指定点を曲線でつないで表現するため、位置8-9の部分で負にグラフが進みますが、計算の値が負になるわけではないです。

12. STEP 3で実際に確認してほしいこと

Excelを開いたら、次の順番で確認するのがおすすめです。

  1. 計算条件を見る。 格子点21個、時間刻み0.1、$r=0.1$ などの設定を確認します。
  2. 15行目を見る。 初期状態で、位置10だけが100になっていることを確かめます。
  3. 16行目を見る。 位置9・10・11の温度が10・80・10に変化したことを確認します。
  4. 17行目を見る。 温度がさらに左右へ広がり、中央が66に下がることを確認します。
  5. L16の数式を見る。 左・中央・右の値を使う差分式と一致することを確かめます。
  6. グラフを見る。 時間とともに中央の山が低く、幅広くなることを確認します。

補足： 現在の完成版Excelでは、係数 $r=0.1$ はB11セルに値として入力されています。今回の計算条件では正しく計算されますが、将来 $\alpha$、$\Delta t$、$\Delta x$ を変更する場合には、B11の値も$r=\frac{\alpha\Delta t}{\Delta x^2}$ に従って更新する必要があります。

Excelシートの、パラメータ「$r=\frac{\alpha\Delta t}{\Delta x^2}$」を、「0.3」に変化させてみます。以下のように、自動的に再計算されて、グラフも時間ステップ２の橙色線の最大が「24.7」に低くなっています。
image.png

image.png

13. STEP 0～5、そしてLBMへつなげる

今回のExcelで行っていることは、次の学習につながります。

STEP 0　流体・物理現象の基礎：「何が変化するのか」
STEP 1　拡散方程式：「変化を数式でどう表すか」
STEP 2　手計算：「実際に数値をどう計算するか」
STEP 3　Excel：「その計算を自動で繰り返す」
STEP 4　Python：「プログラムで同じ計算を行う」
STEP 5　格子と時間発展：「数値計算の基本を整理する」
STEP 6以降　LBM：「分布関数 f を使って流体を計算する」

STEP 3では、格子点の温度 $T$ を直接計算しています。一方、LBMでは、各格子点に方向別の分布関数 $f_i$ を持たせて計算します。

計算する量や更新方法は異なりますが、現在の状態から次の状態を計算し、それを時間方向に繰り返すという基本的な考え方は共通しています。まずは完成版Excelで、この計算の流れを理解してください。STEP 4ではExcelの各セルをPythonの配列に置き換え、同じ計算式をプログラムで再現します。

image.png

■ LBMの基礎の直観的理解【STEP 4 Pythonで計算する】

STEP 4では、STEP 1で学んだ拡散方程式、STEP 2で行った手計算、STEP 3のExcelによる自動計算を踏まえて、同じ計算をPythonで実行し、温度分布が時間とともに変化する様子をPNG画像で確認することを目標にします。

ここでは、LBM解析システム（Autodesk XLB、Google JAX、NVIDIA Warp、Autodesk Neon) は使用しません。WSL Ubuntu 24.04上で、Python・NumPy・Matplotlibだけの、専用計算環境を使用します。

STEP 4の成果物

Pythonスクリプト：diffusion_step4.py：完成版Excelと同じ計算条件を使用
計算結果のグラフ：diffusion_result.png：時間0～2.0の温度分布を比較

1. STEP 1～4のつながりを整理する

まず、これまでの学習を振り返ります。
【STEP 1】　拡散方程式を知る $ \frac{\partial T}{\partial t} =\alpha\frac{\partial^2T}{\partial x^2} $ 温度分布が時間とともに変化する物理現象を理解します。
【STEP 2】　手計算で確かめる　差分式を使って、中央の温度100が80に下がり、左右の温度が10になることを計算しました。
【STEP 3】　Excelで自動計算する　差分式をセルの数式にし、時間ステップ0～20の温度分布を表とグラフで確認しました。
【STEP 4】　Pythonで自動計算する　Excelと同じ差分式をプログラムで繰り返し計算し、温度分布をPNG画像として保存します。

今回の重要なポイントは、STEP 2・3・4では、基本的に同じ計算を 「手計算・Excel・Python」 という3つの方法で実行していることです。

2. 計算する問題と条件を決める

STEP 3の完成版Excelと同じ条件を使用します。

項目 設定値 意味
格子点数 $N$ 21 位置0～20
時間ステップ数 $n_t$ 20 20回計算する
格子間隔 $\Delta x$ 1.0 隣接する格子点の間隔
時間刻み $\Delta t$ 0.1 1回で進める時間
拡散係数 $\alpha$ 1.0 今回使用する係数
差分係数 $r$ 0.1 1回の計算における変化の大きさ

差分係数は、次の式で求めます。$r=\frac{\alpha\Delta t}{(\Delta x)^2}$ 今回の条件を代入すると、$ r=\frac{1.0\times0.1}{1.0^2} =0.1 $となります。初期条件は、中央の位置10だけ温度100とし、その他はすべて0にします。

0  1  2  3  4  5  6  7  8  9  10 11 12 ... 20  位置 x
●──●──●──●──●──●──●──●──●──●──●──●──●──…──●
0  0  0  0  0  0  0  0  0  0 100  0  0 ...  0  温度
                              ↑ 中央位置10

両端の位置0と20は、計算中も温度0に固定します。これは両端の温度を一定にする境界条件です。

3. WSL UbuntuにPythonの解析環境を作る

今回必要なのは、Pythonと2つのライブラリだけで、本検証の主目的である「Autodesk XLB」や、そのバックエンドとなる「Jax, Warp, Neon」、さらに高速計算用の「GPU, CUDA」などは不要です。単純に、CPUのみで動作するUbuntu環境で検証します。

  • NumPy：格子点の温度を配列として管理し、数値計算を行う。
  • Matplotlib：計算した温度分布をグラフにし、PNG画像として保存する。

なお、XLBのCPU検証環境を、この手順「LBM流体解析ツールXLBの活用【01:JAXcpu導入方法と動作確認】」で、「Ubuntu-24.04-XLBjaxCPU」として用意してある場合には、その上に、ここでのLBMの検証環境を追加することも可能です。ただし、この場合でもPython仮想環境「」は区別にして作ります。

【手順⓪　WSLにUbuntuを用意する】

ここでは、「WSL2上のUbuntu24.04LTS（Ubuntu-24.04-LBM）」を用いることにして、LBM用のPython環境を作ることにします。初めて取り組む方は、以下の手順でWSL2上にUbuntu24.04を、導入してください。
https://qiita.com/ShibataRyoichi/items/b9e1bb2e9923e4e2c9bd

WSL2上に複数のUbuntu24.04を利用する場合には、基盤となるUbuntuを導入して、これをexportしてから、目的の環境としてimportします。具体的な以下の手順で行います。環境名としては「Ubuntu-24.04-LBM」とします。
https://qiita.com/ShibataRyoichi/items/31bd2df1091ee8687775

PS C:\Users\ryos> wsl --import Ubuntu-24.04-LBM .\WSL2VM\Ubuntu-24.04-LBM .\Ubuntu-24.04-Base.tar
[========================  41.6%                     ]
↓
この操作を正しく終了しました。

PS C:\Users\ryos> wsl -l -v
  NAME                         STATE           VERSION
* Ubuntu-24.04                 Stopped         2
  Ubuntu-24.04-LBM             Stopped         2

この導入で、作成したUbuntu-24.04-LBMを起動したとき、プロンプトが＃でroot環境の時には、先のQiitaに示した手順で、一般ユーザー（ここではryos）にしておきます。

Ubuntu環境を起動して、最新「Ubuntu 24.04.5 LTS」に更新しておきます。

PS C:\Users\ryos> wsl -d Ubuntu-24.04-LBM
ryos@dalab2403sv01:/mnt/c/Users/ryos$ cd

ryos@dalab2509dt01:/mnt/c/Users/ryos$ cd
ryos@dalab2509dt01:~$ sudo apt update
[sudo] password for ryos:
：
ryos@dalab2509dt01:~$ sudo apt upgrade
：

ryos@dalab2509dt01:~$ lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 24.04.5 LTS
Release:        24.04
Codename:       noble

【手順①　UbuntuのPython環境の準備】

まず、Pythonのバージョンを確認します。次に、venv が使える状態にして、Pythonの環境を整るために、次のコマンドで準備します。最後に「STEP 4用」で用いるLBMのPython検証用の仮想環境を作成します。

ryos@dalab2403sv01:~$ python3 --version
Python 3.12.3
ryos@dalab2403sv01:~$ sudo apt install -y python3-venv python3-pip python-is-python3
：

ryos@dalab2403sv01:~$ python3 -m venv LBM
ryos@dalab2403sv01:~$ ls
LBM
ryos@dalab2403sv01:~$ source ~/LBM/bin/activate
(LBM) ryos@dalab2403sv01:~$

仮想環境を有効にすると、ターミナルの先頭に (LBM) のように表示されます。これからLBMの演習では、必ずこのPython仮想環境(LBM)の中で行い、各種のツールもこの仮想環境でpip導入します。

【手順②　必要なライブラリをインストールする】

以下の手順で、pipを更新してから、必要となるPythonライブラリ「numpy, matplotlib」を導入して、インストールを確認します。OK と表示されれば、準備完了です。

(LBM) ryos@dalab2403sv01:~$ python -m pip install --upgrade pip
：
Successfully installed pip-26.2.1

(LBM) ryos@dalab2403sv01:~$ python -m pip install numpy matplotlib
：
Successfully installed contourpy-1.4.0 cycler-0.12.1 fonttools-4.66.1 kiwisolver-1.5.1 matplotlib-3.11.2 numpy-2.5.4 packaging-26.3 pillow-12.3.0 pyparsing-3.3.3 python-dateutil-2.9.0.post0 six-1.17.0

(LBM) ryos@dalab2403sv01:~$ python -c "import numpy, matplotlib; print('OK')"
OK

4. Pythonスクリプトを作成する

LBM作業用フォルダー「~/LBMwork」を作り、この中に、この段階での演習「Work1」フォルダを作り、その中にスクリプト「diffusion1d-1.py」を作成します。

mkdir -p ~/diffusion_work
cd ~/diffusion_work
nano diffusion_step4.py

次のコードをすべて貼り付けてください。

diffusion1d-1.py
from pathlib import Path

import numpy as np
import matplotlib
matplotlib.use("Agg")  # WSL上で画面表示せず、PNGファイルに保存する
import matplotlib.pyplot as plt

# --------------------------------------------------
# 1. 計算条件（完成版Excelと同じ設定）
# --------------------------------------------------
nx = 21          # 格子点数：位置0～20
nt = 20          # 時間ステップ数：20回計算
dx = 1.0         # 格子間隔
dt = 0.1         # 時間刻み
alpha = 1.0      # 拡散係数（教材用の値）

# rは、1回の計算で温度をどのくらい変化させるかを決める係数
r = alpha * dt / dx**2
if not (0 < r <= 0.5):
    raise ValueError(f"この陽解法の安定条件を満たしていません: r={r:.4f}（0 < r <= 0.5）")

# --------------------------------------------------
# 2. 初期条件：位置10だけ温度100、ほかは0
# --------------------------------------------------
x = np.arange(nx) * dx            # 位置：0, 1, 2, ..., 20
T = np.zeros(nx, dtype=float)     # 全格子点の温度を0にする
T[nx // 2] = 100.0                # 中央位置10の温度だけ100にする

# 各時間ステップの温度を保存する（あとでグラフにするため）
history = [T.copy()]

# --------------------------------------------------
# 3. 時間発展：差分式を20回繰り返す
# --------------------------------------------------
for n in range(nt):
    T_new = T.copy()  # 現在の温度をコピー。両端の温度0もそのまま保つ

    # 位置1～19について、左・中央・右の温度から次の温度を計算する
    for i in range(1, nx - 1):
        T_new[i] = T[i] + r * (
            T[i - 1] - 2 * T[i] + T[i + 1]
        )

    T = T_new                  # 次の時刻の温度に更新
    history.append(T.copy())   # 各時刻の温度を保存

# --------------------------------------------------
# 4. グラフ：Excelと同じく、ステップ0, 5, 10, 15, 20を比較
# --------------------------------------------------
fig, ax = plt.subplots(figsize=(9, 5.5))
plot_steps = [0, 5, 10, 15, 20]

for n in plot_steps:
    ax.plot(
        x,
        history[n],
        marker="o",
        markersize=3,
        linewidth=1.8,
        label=f"t = {n * dt:.1f} (step {n})",
    )

ax.set_title("1D Diffusion: Temperature Profiles Over Time")
ax.set_xlabel("Position x (grid point)")
ax.set_ylabel("Temperature T")
ax.set_xlim(0, 20)
ax.set_ylim(0, 110)
ax.set_xticks(np.arange(0, 21, 2))
ax.set_yticks([0, 20, 40, 60, 80, 100])
ax.grid(True, alpha=0.3)
ax.legend()
fig.tight_layout()

# スクリプトと同じフォルダーに画像を保存する
output_path = Path(__file__).resolve().parent / "diffusion_result.png"
fig.savefig(output_path, dpi=180, bbox_inches="tight")
plt.close(fig)

# --------------------------------------------------
# 5. 計算条件と、中央付近（位置8～12）の数値を表示
# --------------------------------------------------
print("=== STEP 4: 1次元拡散方程式 ===")
print(f"格子点数 nx = {nx}（位置0～{nx - 1}）")
print(f"時間ステップ数 nt = {nt}")
print(f"dx = {dx}, dt = {dt}, alpha = {alpha}")
print(f"r = alpha * dt / dx^2 = {r:.3f}")
print(f"最終時刻 = {nt * dt:.1f}")
print("\n位置8～12の温度（Excelとの比較用）")
print("step | time |       T[8:13]")
for n in range(5):
    values = ", ".join(f"{v:.2f}" for v in history[n][8:13])
    print(f"{n:>4} | {n * dt:>4.1f} | [{values}]")
print(f"\nグラフを保存しました: {output_path}")

5. スクリプトの処理を詳しく理解する

ここでは、上記で準備した、１次元拡散方程式を、手計算やExcelで試したように、単純に数値計算で解く手順をPythonスクリプトで、演習します。ここでは、格子ボルツマン法のような分布関数を用いた計算ではなく、拡散方程式を差分式に変換し、その式を繰り返し計算する「有限差分法」で解く方法になります。
今回のスクリプトは、各格子点に温度 $T$ を持たせ、差分式で直接更新します。一方、LBMでは、各格子点に方向別の分布関数 $f_i$ を持たせます。この2つは、計算方法が異なります。ただし、どちらも空間を格子で表し、時間を少しずつ進めていくという共通点があります。
スクリプト全体は、５つの処理に分かれています。それぞれの部分を、スクリプトを提示しながら、解説します。

【① 計算条件を設定する】
ここは、Excel上部の「計算条件」に相当します。Excelでは数値をセルに入れましたが、Pythonでは変数に代入しています。例えば、nx = 21 は21個の格子点を使うことを意味します。

# --------------------------------------------------
# 1. 計算条件（完成版Excelと同じ設定）
# --------------------------------------------------
nx = 21          # 格子点数：位置0～20
nt = 20          # 時間ステップ数：20回計算
dx = 1.0         # 格子間隔
dt = 0.1         # 時間刻み
alpha = 1.0      # 拡散係数（教材用の値）

【② 初期状態を作る】
ここは、位置Xに対する温度Tの初期状態を設定します。

# --------------------------------------------------
# 2. 初期条件：位置10だけ温度100、ほかは0
# --------------------------------------------------
x = np.arange(nx) * dx            # 位置：0, 1, 2, ..., 20
T = np.zeros(nx, dtype=float)     # 全格子点の温度を0にする
T[nx // 2] = 100.0                # 中央位置10の温度だけ100にする

それぞれの役割は次のとおりです。

コード 意味
np.arange(nx) 0～20の位置番号を作る
np.zeros(nx) 21点すべての温度を0にする
T[nx // 2] = 100.0 中央の位置10だけを100にする

なお、グラフを描くために初期状態を保存しています。計算の途中で温度配列 T は何度も更新されるので、時間ごとの状態を残しておく必要があります。

# 各時間ステップの温度を保存する（あとでグラフにするため）
history = [T.copy()]

【③ 差分式を繰り返し計算する】

ここがプログラムの最も重要な部分です。

for n in range(nt):
    T_new = T.copy()  # 現在の温度をコピー。両端の温度0もそのまま保つ

    # 位置1～19について、左・中央・右の温度から次の温度を計算する
    for i in range(1, nx - 1):
        T_new[i] = T[i] + r * (
            T[i - 1] - 2 * T[i] + T[i + 1]
        )
    T = T_new                  # 次の時刻の温度に更新
    history.append(T.copy())   # 各時刻の温度を保存

2つの for 文の役割は違います。

  • 外側の for n in range(nt)：時間を1ステップずつ進める。
  • 内側の for i in range(1, nx - 1)：各格子点の次の温度を計算する。

つまり、以下の処理を行います。

時間ステップ0：内部の19点を計算　→　時間ステップ1：内部の19点を計算　→　…　→　時間ステップ20

なお、最初に T_new = T.copy() とするのは重要です。新しい温度を計算している途中で、計算済みの値がほかの格子点の計算に混ざらないようにするためです。これは、Excelで「前の行の温度を参照し、次の行の温度を計算する」仕組みに対応します。

【④ 温度変化のグラフを作り画像に保存する】

計算した21時刻分のデータから、0、5、10、15、20ステップの温度分布を選んで描きます。

plot_steps = [0, 5, 10, 15, 20]

for n in plot_steps:
    ax.plot(
        x,
        history[n],
        marker="o",
        markersize=3,
        linewidth=1.8,
        label=f"t = {n * dt:.1f} (step {n})",
    )

時間に換算すると、次の5つです。これは、完成版Excelのグラフと比較するための設定です。具体的には、時間ステップ：時刻 $t$」の対応だと「0:0.0 5:0.5 …　20:2.0」となります。

グラフでは画面表示を使わないための設定を行い、PNG画像として保存します。

output_path = Path(__file__).resolve().parent / "diffusion_result.png"
fig.savefig(output_path, dpi=180, bbox_inches="tight")
plt.close(fig)

【⑤ 計算条件中央付近の数値を表示する】

print("=== STEP 4: 1次元拡散方程式 ===")
print(f"格子点数 nx = {nx}（位置0～{nx - 1}）")
print(f"時間ステップ数 nt = {nt}")
print(f"dx = {dx}, dt = {dt}, alpha = {alpha}")
print(f"r = alpha * dt / dx^2 = {r:.3f}")
print(f"最終時刻 = {nt * dt:.1f}")
print("\n位置8～12の温度（Excelとの比較用）")
print("step | time |       T[8:13]")
for n in range(5):
    values = ", ".join(f"{v:.2f}" for v in history[n][8:13])
    print(f"{n:>4} | {n * dt:>4.1f} | [{values}]")
print(f"\nグラフを保存しました: {output_path}")

6. Pythonスクリプトを実行する

仮想環境(LBM)が有効な状態で、作業フォルダー「~/LBMwork/Work1/」へ移動します。pythonでスクリプトを実行します。正常に実行されると、ターミナルに計算条件と中央付近の温度が表示されます。

(LBM) ryos@dalab2403sv01:~$ cd ~/LBMwork/Work1/
(LBM) ryos@dalab2403sv01:~/LBMwork/Work1$ ls
diffusion1d-1.py
(LBM) ryos@dalab2403sv01:~/LBMwork/Work1$ python diffusion1d-1.py
=== STEP 4: 1次元拡散方程式 ===
格子点数 nx = 21（位置0～20）
時間ステップ数 nt = 20
dx = 1.0, dt = 0.1, alpha = 1.0
r = alpha * dt / dx^2 = 0.100
最終時刻 = 2.0

位置8～12の温度（Excelとの比較用）
step | time |       T[8:13]
   0 |  0.0 | [0.00, 0.00, 100.00, 0.00, 0.00]
   1 |  0.1 | [0.00, 10.00, 80.00, 10.00, 0.00]
   2 |  0.2 | [1.00, 16.00, 66.00, 16.00, 1.00]
   3 |  0.3 | [2.40, 19.50, 56.00, 19.50, 2.40]
   4 |  0.4 | [3.88, 21.44, 48.70, 21.44, 3.88]

グラフを保存しました: /home/ryos/LBMwork/Work1/diffusion_result.png
(LBM) ryos@dalab2403sv01:~/LBMwork/Work1$ tree
├── diffusion1d-1.py
└── diffusion_result.png

最後の表示のtreeコマンドは、「sudo apt install tree」で導入できます。

7. 計算結果のグラフを確認する

作成された画像はこちらです。「計算結果：diffusion_result.png」
diffusion_result.png

グラフで確認できることは、横軸は位置0～20、縦軸は温度 $T$ です。5本の線は、時間ステップ0、5、10、15、20における温度分布を表します。初期状態では、位置10の温度が100で、周囲は0です。時間が進むと、中央の山が低くなり、左右の温度が上がって広がります。これはSTEP 1で学んだ拡散現象そのものです。
なおExcelのグラフでは、指定点を曲線でつないで表現するため、位置8-9の部分で負にグラフが進みますが、このグラフでは指定点を単純に直線でつないでいます。

ここで、Excelの確認と同じように、拡散係数の値を「1→3」に変更した結果が、以下の通りです。
diffusion_result.png
同じように、急速に拡散して、全体的に値が低くなっています。

8. このPythonプログラムの核心

プログラム全体の中で、最も重要なのは次の計算です。

T_new[i] = T[i] + r * ( T[i - 1] - 2 * T[i] + T[i + 1] )

これは、以下のSTEP 2の差分式と同じものです。
$ T_i^{n+1} = T_i^n+ r\left( T_{i-1}^n-2T_i^n+T_{i+1}^n \right) $

記号とPythonの変数を対応させると、次のようになります。

数学の記号 Python 意味
$T_i^n$ T[i] 現在の中央の温度
$T_{i-1}^n$ T[i-1] 現在の左隣の温度
$T_{i+1}^n$ T[i+1] 現在の右隣の温度
$T_i^{n+1}$ T_new[i] 次の時間の中央の温度
$r$ r 温度変化の大きさを決める係数

つまり、「Excel：セルの数式で次の温度を計算する　→　Python：配列とfor文で次の温度を計算する」という違いがありますが、計算の中身は同じです。

9. STEP 4からSTEP 5へ

STEP 4まで終えたら、次に理解したいのは「格子」と「時間発展」です。今回のPythonプログラムは、次の構造になっています。この考え方は、これまでのSTEP 0～4で繰り返し確認してきたものです。

■ 空間を21個の格子点に分ける
↓ 各格子点に温度 T を持たせる
↓ 左・中央・右の温度を使って計算する
↓ 次の時間ステップの温度を得る
これを20回繰り返す

ただし、LBMでは温度 $T$ を直接更新するのではなく、各格子点に方向別の分布関数 $f_i$ を持たせます。

STEP 4：拡散方程式　温度 T　→　差分計算　→　次の温度 T
STEP 5：格子と時間発展　→　LBMへ進む
STEP 6：分布関数 f
STEP 7：Collision
STEP 8：Streaming

差分式がPythonの計算コードになること、時間ステップを繰り返すことで温度分布が変化すること、そして結果をグラフで確認できること。この3点を理解できれば、次のSTEP 5へ進む準備ができています。

■ LBMの基礎の直観的理解【STEP 5 格子と時間発展を理解する】

ここは「LBMの基礎の直観的理解」のルートの中でも非常に重要な接続部分です。STEP 1～4では「拡散方程式を有限差分法で解く」ことを学びましたが、STEP 5では、拡散方程式そのものから一歩離れて、「コンピュータで物理現象を計算するとき、なぜ格子を使いなぜ時間を少しずつ進めるのか」 を整理します。この説明を入れておくと、STEP 6以降のLBMが自然につながります。

STEP 1～4からSTEP 5へのつなぎ方

これまで学習してきた内容の全体を、まず1本につなげます。

STEP 1：拡散方程式を知る　→　「温度が時間とともに変化する」
STEP 2：手で計算する　→　「連続した空間を、格子点に置き換える」
STEP 3：Excelで計算する　→　「格子点の計算を、時間ごとに表にする」
STEP 4：Pythonで計算する　→　「同じ計算をプログラムで繰り返す」
STEP 5：格子と時間発展を理解する　→　「物理現象を、格子上で時間方向に進めて計算する」
→　STEP 6 LBMの分布関数 f：LBMの基本概念への展開

つまり、STEP 5は新しい計算方法を学ぶステップではありません。STEP 1～4でやってきたことを、「数値解析の基本的な考え方」 として整理するステップです。

1. STEP 1では「現象」を見ていた

STEP 1では、$\frac{\partial T}{\partial t}=\alpha\frac{\partial^2 T}{\partial x^2}$ という拡散方程式を学びました。ここでは、「温度　→　時間とともに変化　→　熱が広がる」という物理現象が中心でした。この段階では、まだコンピュータでどう計算するかは重要ではありません。

2. STEP 2で「連続空間」を格子にした

次に手計算をするため、ここで初めて、以下のように空間を離散化する という考え方が登場します。
これがSTEP 5の「格子」の原点です。

連続した棒　　　　　　　→　　格子に分ける
───────────────────　　　　●──●──●──●──●──●──●

3. STEP 3のExcelでは「時間」も並べた

Excelにすると、以下のような数値の変化がわかりました。

          空間 [→]
        0  1  2  3  4  5 ...
時間0   0  0  0 100  0  0 ...
時間1   0  0 10  80 10  0 ...
時間2   0  1 16  66 16  1 ...
時間3   0  2 19  56 19  2 ...
時間[↓]　

ここで非常に重要なことが見えてきます。「横方向＝空間・縦方向＝時間」です。つまり、Excelの表そのものが、**「格子上で時間が進んでいる様子」**を表しています。

4. STEP 4のPythonでは「繰り返し」を見た

Pythonでは、以下のスクリプトで計算しました。

for n in range(nt):
    for i in range(1, nx - 1):
        T_new[i] = ...

ここで、「外側のfor　→　時間を進める」「内側のfor　→　各格子点を計算する」となっています。つまりPythonコードを日本語にすると、**「各格子点を計算し、それが終わったら時間を1ステップ進め、またすべての格子点を計算する」**です。

5. ここでSTEP 5に入る

ここまで来たところで、一度Pythonのコードから離れます。そして、「今まで何をしていたのでしょうか？」 と振り返ります。答えは、以下の手順で、これがSTEP 5の中心です。

① 空間を格子に分ける
② 各格子点に値を持たせる
③ 周囲の値を使って新しい値を計算
④ 時間を１ステップ進める
⑤ また計算する
⑥ これを繰り返す

6. 「格子」とは何か

STEP 5では、格子を単なる「計算のための点」として整理します。例えば、以下のような格子があったら、

●──●──●──●──●
0  1  2  3  4

以下のように、各格子点に物理量を持たせると考えます。今回なら「温度」です。

位置0 → 温度 ・ 位置1 → 温度　...　位置4 → 温度

7. 「時間発展」とは何か

時間発展は難しい言葉ですが、簡単です。現在の状態から、次の時間の状態を計算することです。例えば、以下の変化です。

■ t = 0：　0  0  0 100  0  0  0
↓ t = 1：　0  0 10  80 10  0  0
↓ t = 2：　0  1 16  66 16  1  0

つまり、「現在　→　計算　→　次の状態　→　計算　→　次の状態」という繰り返しです。

8. ここで「格子」と「時間」を組み合わせる

STEP 5では、次の図を中心に説明すると非常に分かりやすいです。この図を見ると、「格子」と「時間発展」が同時に見えるようになります。

    空間 x →  0   1   2   3   4   5
              ●   ●   ●   ●   ●   ●
■ 時間 t=0    0   0   0 100   0   0
↓ 時間 t=1    0   0  10  80  10   0
↓ 時間 t=2    0   1  16  66  16   1
↓ 時間 t=3    0   2 19.5 56 19.5  2
↓ ...

9. ここで一つ大事な言葉を導入する

STEP 5では、「状態」 という言葉を導入すると、LBMにつながりやすくなります。例えば拡散方程式では、「t = 0：[0, 0, 0, 100, 0, 0, 0]」が「現在の状態」です。計算すると、「t = 1：[0, 0, 10, 80, 10, 0, 0]」という「次の状態」になります。つまり、「現在の状態 → 計算 → 次の状態」という考え方です。

10. そしてLBMへの橋を架ける

ここで初めて、「では、LBMでは何を状態として持つのでしょうか？」 と問いかけます。拡散方程式では、「各格子点：温度 T」を持っていました。LBMでは、「各格子点：分布関数 f」を持ちます。これがSTEP 6への接続です。
拡散方程式と格子ボルツマン法（LBM）を比較すると、

  • 拡散方程式の差分法：　格子点が温度情報を持つ
  • 格子ボルツマン法：　格子点には分布関数$f$を持つ
    しかし、「どちらも、格子の上で状態を時間発展させる方法」 と説明できます。

STEP 5からSTEP 6へ

最後は、この一言でつなぐのがよいと思います。「これまで私たちは、各格子点に温度を持たせて時間発展させてきました。ではLBMでは、各格子点に何を持たせるのでしょうか？」答えは 「温度 T → 分布関数 f」 と切り替えることです。そしてSTEP 6で初めて、「1つの格子点に、実は1つの値ではなく、方向ごとの複数の分布関数が存在する」 というLBM独特の考え方に入ります。この流れなら、STEP 1～5が「LBMを理解するための土台」になり、STEP 6から本格的なLBMの学習に入るという構成になります。

以上のSTEP 5の内容を説明図にまとめました。
image.png

■ LBMの基礎の直観的理解【STEP 6 分布関数 f を知る】

STEP 6は、今回の学習ルートの中で最も重要な転換点です。STEP 1～5では、「格子点に1つの値を持たせて、時間とともに更新する」という考え方を学びました。STEP 6では、そこから一歩進んで、「LBMでは、1つの格子点に1つの値ではなく、“方向を持った複数の分布関数 f”を持たせる」 という考え方を導入します。

まず、「STEP 6」では、「分布関数 $f$ を理解」 することを目指します。STEP 5までの拡散方程式では、「格子点　→　温度 T　→　時間が進む　→　新しい温度 T」でした。つまり、1つの格子点に、1つの温度 (T) があるという世界でした。

1. LBMでは「1つの値」ではない

ここで、「流体は単に温度のような1つの値だけ持っているのでしょうか？」と考えます。流体には自由に移動する 「方向」 があります。

例えば、２次元の1つの格子点に対して、「右へ進む・左へ進む・上へ進む・下へ進む・斜めに進む」など、いろいろな方向の流体の情報があります。そこでLBMでは、それぞれを$f_0,\ f_1,\ f_2,\ldots$ のような分布関数で表します。

2. 「分布関数 f」とは何か

初心者には、「その格子点にある、ある方向へ進もうとしている流体の量」 と説明するのが分かりやすいです。例えば、以下のようなイメージです。ここでの $f$ は、実際の分子を1個ずつ追跡しているわけではありません。「その方向へ進む流体の分布（割合）」を数学的に表したものです。

       上へ進む　　　　　　　　 １つの格子点
           ↑ f₂　　　　　      ↓
左へ f₁ ←  ●  → f₃ 右へ　　　　方向ごとの流体情報
           ↓ f₄　　　　　      ↓
       下へ進む　　　　　      f₀, f₁, f₂, ...

3. 「分布」という言葉を理解する

例えば、ある格子点で、「右向き：多い（→ → → → ）・左向き：少ない（←）・上向き：少ない（↑）・下向き：少ない（↓）」という分布ならは、その結果、全体として右向きの流れがある と考えられます。
逆に、「右向き ≒ 左向き」なら、左右の流れが打ち消し合って、平均速度は小さい となります。

以上の分布関数$f$について、改めて図解で説明します。
image.png

4. ここで「平均」という考え方が出てくる

STEP 6では、まだ難しい式をたくさん使わない方がよいですが、分布関数を全部合わせると、流体全体の状態が分かる という考え方を導入します。まず密度は、$\rho=\sum_i f_i$です。つまり、この場合は$f$は質量として、「すべての方向の f　→　合計する　→　密度 ρ」です。
もう少し簡単に言うと、全ての方向に分配する関数値$f$（質量）の合計は、そこにある質量の総量なので、単位となる格子１つにおける質量として、密度になります。

5. 速度も分布関数から求める

速度については、$\rho\mathbf{u}=\sum_i f_i\mathbf{c}_i$ と考えます。ここで、$f_i$：その方向の分布（質量）、$\mathbf{c}_i$：その方向（速度ベクトル）、$\rho$は全体の密度、$\mathbf{u}$：流体の平均的な速度、です。初心者向けには、「各方向の流体の量に、その方向を掛けて全部足すと、流体全体の流れの方向と強さが分かる」 と説明するとよいです。
もう少し別の説明をすると、ポイントは、「量 × 速度」を全部足すと、全体の「量 × 平均速度」になるという考え方です。実際のLBMでは、右・左だけでなく、上下や斜めなどの複数の方向を扱います。そのため、速度はベクトルで表します。覚えておきたいことは、以下の３つの物理量です。この関係は、LBMで流体の密度と速度を求める基本原理です。

  • $f_i$ をそのまま足す → 流体の総量（質量の格子毎なので密度）
  • $f_i$ に $\mathbf{c}_i$ を掛けて足す → 方向を考慮した流体の運動量
  • 運動量を密度で割る → 平均速度$\mathbf{u}$

速度を分布関$f$から求めることは 「運動量」 を考えると良いです。運動量とは、物体や流体が「どれくらいの勢いで動いているか」を表す量です。例えば、同じ速度で動いていても、量が多い流体ほど、全体として大きな運動量を持ちます。つまり、「運動量 = 質量$\times$速度」です。
LBMでは、各方向の分布 $f_i$ に、その方向の速度 $\mathbf{c}_i$ を掛けて足し合わせて、「$\rho\mathbf{u}=\sum_i f_i\mathbf{c}_i$」のように求めます。LBMでは、このように各方向の運動量を合計することで、格子点における流体全体の運動量を求めます。そしてこの運動量を、密度 $\rho$ で割ると、平均速度 $\mathbf{u}$ が得られます。覚え方は、密度は「どれだけあるか」、運動量は「どちら向きに、どれだけの勢いで動いているか」 です。

6. STEP 5の差分法とSTEP 6のLBMとの違い

ここは表にすると非常に分かりやすいです。

STEP 5：差分法の拡散方程式 STEP 6：格子ボルツマン法
格子点に温度 $T$ 格子点に分布関数 $f_i$
1つの値 複数の方向の値
温度を直接更新 分布関数を更新
拡散現象 さまざまな流体の運動
$T$を見る 分布関数$f_i$から 密度$\rho,平均速度\mathbf{u}$ を求める

7. なぜ「方向」を持たせるのか

ここが非常に重要です。温度の場合、「温度100」という情報だけでも、その場所の状態を表せます。しかし流体では、「右へ流れている・左へ流れている・上へ流れている・下へ流れている」という方向の情報が必要です。
そのためLBMでは、「どの方向に、どれくらい流体があるか」を、分布関数 $f_i$として持たせます。これによって、流体の流れそのものを表現できるようになります。これにより、「LBMは、各格子点の分布関数を計算し、その結果から密度や速度を得る方法」 と説明します。

8. そしてSTEP 7へ

STEP 6の最後に、「では、この方向別の分布関数 (f) は、時間が進むとどう変化するのでしょうか？」 と問いかけます。そこでSTEP 7のCollisionにつなげます。「STEP 6：方向別の分布関数 f　→ 「fは時間とともに変化する」 → では、何がfを変える？ → STEP 7：Collision」となり、そしてさらに、「Collision → 分布関数を変える → Streaming → 隣の格子へ移動する」という流れに進めます。

STEP 6の到達目標

この段階では、次の3つが分かれば十分です。
① $f$ は何か： 各方向へ進む流体の分布を表す量。
② なぜ複数あるのか：　流体には「方向」があるから。
③ $f$ から何が分かるのか：　分布関数を集計すると、密度や速度を求められる。

この3点が理解できれば、次の 「STEP 7 Collision＝分布関数をどう変化させるのか」 へ非常に自然につながります。

以上のSTEP 6の内容を説明図にまとめました。
image.png

以上で、前半【00A:LBMの基礎の直観的理解(STEP 0-6】 が終わりです。
後半【00B:LBMの基礎の直観的理解(STEP 7-12】は、以下のリンクから進んでください。

1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?