格子ボルツマン法と粒子法

なぜ「速度場を直接解かない」流体解析があるのか

流体シミュレーションの主流は、ナビエ・ストークス方程式を有限差分・有限体積・有限要素で離散化し、圧力と速度を連立して解く手法です。しかし、この直接的なアプローチには弱点があります。非圧縮性を保証するための圧力補正（poisson方程式の求解）は大域的な連立一次方程式を伴い、並列化すると通信コストが跳ね上がります。また、複雑形状や自由表面、大変形を扱うにはメッシュの生成・追従が重荷になります。

格子ボルツマン法(LBM)とSPH(smoothed particle hydrodynamics、平滑化粒子流体力学)は、どちらも「巨視的な流速場を直接解かない」という発想でこの弱点を回避します。LBMは分布関数という中間的な量を格子上で動かし、SPHはメッシュそのものを捨てて粒子だけで場を表現します。着眼点は正反対ですが、どちらも「局所的な情報のやり取りだけで大域的な流れを再現する」ことを狙っており、並列計算との相性を強く意識して設計された手法だと言えます。

格子ボルツマン法(LBM) ── 分布関数の衝突と移流

LBMは流体を粒子分布関数 f_i(x, t) で表します。これは位置 x、時刻 t において、離散化された有限個の速度方向 i（典型的には2次元9方向のD2Q9、3次元19方向のD3Q19など）へ運動する仮想粒子の密度です。巨視的な密度 ρ と速度 u は、この分布関数のモーメント（速度方向についての重み付き総和）として事後的に計算されます。

巨視量の復元（D2Q9の例、e_i は離散速度ベクトル）

  密度   ρ(x,t)   = Σ_i f_i(x,t)
  運動量 ρu(x,t)  = Σ_i f_i(x,t)・e_i

LBMの時間発展は、各格子点で完結する2段階の演算に分解されます。

LBMの1タイムステップ

  1. 衝突(collision):
     各格子点で、局所平衡分布 f_i^eq へ緩和させる
     f_i*(x,t) = f_i(x,t) - (1/τ)・[f_i(x,t) - f_i^eq(x,t)]
     ※ BGK近似（単一緩和時間モデル）。τ は緩和時間で動粘性係数 ν に対応

  2. 移流(streaming):
     衝突後の分布を、離散速度方向にそのまま1格子分だけ動かす
     f_i(x + e_i・Δt, t + Δt) = f_i*(x, t)

衝突ステップは各格子点が自分の値だけで完結する完全な局所演算です。移流ステップも、隣接格子点への単純なデータ移動に過ぎません。この2段階を繰り返すと、チャップマン・エンスコグ展開という漸近解析により、連続極限でナビエ・ストークス方程式（圧縮性がごく小さい極限での近似）が回復されることが示されます。

LBMの並列親和性

LBMが並列計算と相性が良い理由は、演算が「同一格子点内で完結する衝突」と「最近傍格子点への移流」だけに分解されることに尽きます。領域分割(domain decomposition)を行うと、各ランクが担当する格子ブロックの内部はローカルに衝突・移流を計算でき、通信が必要なのは境界に接するゴースト層（隣接ランクの最外殻1〜2層分の分布関数値）の交換だけです。この通信パターンは有限差分法の隣接通信と同質ですが、LBMは大域的な連立方程式求解を伴わない分、通信と計算の重ね合わせ（オーバーラップ）がしやすく、GPU実装とも親和性が高いという特徴があります。

複雑形状の扱いやすさも実務上重要です。障害物の表面をbounce-back（衝突後の分布を反射させ壁で流速ゼロを模す境界条件）という単純なルールで表現できるため、多孔質媒体や骨梁構造のような複雑形状でも境界適合メッシュを作らずに済みます。ただし、緩和時間 τ が1/2に近づくと数値不安定になりやすく、高レイノルズ数（乱流）領域ではMRT(multiple-relaxation-time、複数緩和時間)モデルなど安定化の工夫が必要になる点は限界として押さえておくべきです。

ラグランジュ的粒子法 ── SPHの発想

LBMが「格子に固定された分布関数」というオイラー的（空間に固定した視点）な表現を使うのに対し、SPHは完全にラグランジュ的（物質点に追従する視点）です。流体を離散的な粒子の集合として表現し、各粒子が質量・運動量・エネルギーを持って移動します。任意の物理量 A(x) は、近傍粒子の値を平滑化カーネル関数 W で重み付けした総和として補間されます。

SPHの基本補間（粒子 j の近傍集合について和を取る）

  A(x_i) ≈ Σ_j  m_j・(A_j / ρ_j)・W(x_i - x_j, h)

  m_j : 粒子jの質量
  ρ_j : 粒子jの密度
  h   : 平滑化長（カーネルの影響半径）
  W   : カーネル関数（|x_i - x_j| > h で0になるコンパクト台を持つ）

密度や圧力勾配、粘性項もこの補間の微分形で計算され、運動方程式（ニュートンの第二法則）を各粒子について時間積分することで流体全体の運動を追跡します。メッシュを持たないため、自由表面（液面の分裂・合体）、大変形、破砕といった現象を境界追跡なしに自然に扱えるのがSPHの最大の強みです。決壊した水の挙動や、天体物理での星形成シミュレーションなど、境界形状が刻々と変わる問題で使われるのはこのためです。

近傍探索という並列化のボトルネック

SPHの並列計算における最大の課題は、格子法のような固定した通信パターンを持たない点にあります。各粒子は毎ステップ移動するため、「誰が自分の近傍か」を毎回探索し直す必要があります。この近傍探索を愚直に全粒子対で行うと計算量は粒子数の2乗に比例してしまうため、実用上は空間分割データ構造で近傍候補を絞り込みます。

近傍探索の高速化に使われる代表的なデータ構造

  1. セルリスト(linked-cell法):
     空間を平滑化長 h 程度の一様なセルに分割し、
     各粒子が属するセルを管理。近傍探索は隣接する
     26セル(3次元の場合)のみ走査すればよい → 計算量は概ね粒子数に比例

  2. k-d木・八分木(オクトツリー):
     粒子密度が空間的に不均一な場合に有効。
     セルリストより探索コストは高いが、疎密差が大きい系に強い

並列環境では、この空間分割そのものを領域分割の単位として使うのが一般的です。しかし格子法と違い、粒子は時間とともに担当ランクをまたいで移動するため、次のような追加コストが発生します。

とりわけ粒子の移動(migration)は格子法には存在しない負荷で、担当領域の境界を越えた粒子をランク間でやり取りするたびにデータ構造の再構築が必要になります。さらに、自由表面や崩壊現象では粒子分布が急激に偏るため、静的な領域分割では負荷不均衡が発生しやすく、実運用では空間充填曲線（ヒルベルト曲線など）に基づく動的な再分割で負荷を均す実装が広く使われています。

まとめ

• 格子ボルツマン法(LBM)は流速場でなく粒子分布関数を扱い、衝突（局所演算）と移流（最近傍への単純移動）の2段階だけでナビエ・ストークス方程式を近似的に回復する。大域的な圧力求解が不要なため通信は最近傍に限定され、並列化・GPU化に強い。
• SPHに代表されるラグランジュ的粒子法はメッシュを持たず、近傍粒子のカーネル補間で物理量を表現するため自由表面や大変形に強いが、近傍関係が時々刻々変わる。
• SPHの並列化ではセルリストやk-d木で近傍探索の計算量を抑えつつ、粒子の担当ランク間移動(migration)とゴースト粒子同期、動的負荷分散が性能を左右する固有の課題となる。
• 格子法(LBM)は通信パターンが静的で予測しやすく、粒子法は動的な近傍管理が必要という対比が、両者の並列化戦略の違いの本質である。