格子ボルツマン法（LBM）を使った解析

1. 概要

格子ボルツマン法（Lattice Boltzmann Method, LBM）は、 流体解析（CFD）の世界において、 従来のナビエ・ストークス方程式を直接解く手法に代わる手法として注目されている。

これは、「流体を連続体としてではなく、仮想的な粒子の集まりとして捉え、 その動きをシミュレーションする」手法であり、 有限体積法（FVM）、有限要素法（FEM）などの従来の手法に比較して解析プログラムを簡単に実装できる点に特徴がある。

ここでは一般の流体解析用のベンチマーク計算に使われているテイラーグリーン渦流れを対象として LBM法を使った計算を行ってみた。

2. 格子ボルツマン法とは？

一般的な流体シミュレーションでは、 流体の速度や圧力を変化量として扱う「マクロ」な視点（ナビエ・ストークス方程式）を用いる。

このスケール領域では、FDM、FVM、FEMなどの解析手法が使われる。

一方、極めて小さいミクロなスケールに対しては 分子の運動を直接シミュレートする分子動力学（MD法）が使われる。

LBMはそれよりもミクロとマクロの中間に位置する、 **「メゾスコピック（中間的）」**な視点に基づき解析を行う。

図1 — 流体シミュレーションのスケール領域

LBMでは、空間を格子（グリッド）に区切り、各格子点に存在する粒子の分布関数を計算する。

粒子が隣の格子へ移動し（並進）、粒子同士がぶつかって速度を変える（衝突）という シンプルなプロセスを繰り返すことで、結果として流体全体の動きを再現する。

LBM法は気体分子運動論のアナロジーから導かれているが、 漸近理論（S展開）に基づけば、ナビエ・ストークス方程式を近似することが示されている。

LBM法の長所、短所は以下の通りである。

• 長所 (Advantages)

  • 並列計算との親和性が極めて高い:
    LBMの基本演算は「衝突（Collision）」と「移動（Streaming）」の2ステップで構成されている。 各格子点での計算が局所的であるため、 MPIやGPUを用いた大規模な並列計算において非常に高いスケーラビリティを発揮する。

  • 複雑な境界形状への対応が容易:
    格子（セル）がデカルト格子（正方形・立方体）ベースであるため、 複雑な形状の壁面を「Bounce-back」という単純なルールで処理できる。 多孔質媒体や複雑な幾何形状内の流れにおいて、メッシュ生成の手間を大幅に削減できる。

  • アルゴリズムが単純:
    NS方程式でボトルネックとなる圧力に関するポアソン方程式を解く必要がない。 圧力が密度と状態方程式を介して直接求まるため、プログラムの実装が比較的シンプルである。

  • 多相流・微細流動のシミュレーションに強い:
    粒子間相互作用をモデル化しやすいため、 気液二相流、界面の追跡、懸濁液（サスペンション）などの複雑な現象を扱うのに適している。

• 短所 (Disadvantages)

  • 格子依存性とメモリ消費量:
    標準的なLBMは等間隔のデカルト格子を使用するため、 壁面近傍など特定の領域だけ細かくする「局所細分化」の柔軟性が NS方程式ベースの手法（有限体積法など）に比べて低い。 その結果、計算領域全体で格子数が膨大になり、メモリ消費量が増える傾向がある。

  • 高マッハ数（圧縮性流れ）への制限:
    LBMは基本的に低マッハ数（非圧縮性近似）を前提としている。 音速に近い流れや衝撃波を伴う超音速流の計算には、 高度なモデルの拡張が必要であり、安定性の確保が難しくなる。

  • 数値的不安定性:
    高いレイノルズ数の流れを計算する際、緩和時間（Relaxation time）が極端に短くなると、数値計算が発散しやすくなる。 これを回避するためにMRT（Multiple Relaxation Time）モデルなどが使われるが、パラメータ調整の難易度が上がる。

  • 熱流体計算の複雑化:
    標準的なモデルは等温（Isothermal）を前提としている。 温度変化を伴う熱流体を扱う場合、追加の分布関数を導入する必要があり、 LBMの持つ「シンプルさ」という利点が薄れることがある。

3. LBMの基本的なアルゴリズム

LBMの計算は、非常にシンプルな2つのステップ（衝突ステップと並進ステップ）の繰り返しで構成されている。

図2 — LBMの計算フローチャート

💡 LBMの計算は「衝突」と「並進」の2ステップが基本であり、NS方程式ベースの手法で必要となる圧力ポアソン方程式の求解が不要であるため、アルゴリズムが非常にシンプルである。

3.1. 衝突ステップ (Collision)

格子点に集まった粒子が衝突し、統計的に安定な状態（平衡状態）に近づくプロセスである。 最も一般的なモデルはLBGKモデルと呼ばれ、以下の式で表される。 添字の$i$は、粒子の速度方向を表すインデックスである。

$$f_i(x, t + \Delta t) = f_i(x, t) - \frac{1}{\tau} (f_i(x, t) - f_i^{eq}(x, t))$$

• $f_i$: 粒子の分布関数
• $f_i^{eq}$: 平衡分布関数
• $\tau$: 緩和時間（粘性に関係するパラメータ）

3.2. 並進ステップ (Streaming)

衝突した後の粒子が、それぞれの速度方向に沿って隣接する格子点へと移動するプロセスである。 $$f_i(x + e_i \Delta t, t + \Delta t) = f_i(x, t + \Delta t)$$

3.3. 平衡分布関数

平衡分布関数の定義式

$$f_i^{eq}(\rho, \mathbf{u}) = w_i \rho \left[ 1 + \frac{\mathbf{e}_i \cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{e}_i \cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u} \cdot \mathbf{u}}{2c_s^2} \right]$$

各変数の意味は以下の通り。

3.4. 離散速度モデル

LBMでは、粒子の移動方向を離散化するために、D2Q9モデルなどの離散速度モデルが用いられる。 D2Q9モデルは、2次元空間において9つの離散速度を持つモデルで、以下のように定義される。

図3 — D2Q9離散速度モデル（2次元・9速度）

3.5. 密度、流速の計算

LBMでは、分布関数からマクロな物理量である密度と流速を計算することができる。 密度 $\rho$ と流速 $\mathbf{u}$ は以下の式で求められる。 $$ \begin{align} \rho(x, t) &= \sum_{i} f_i(x, t) \notag \\ \mathbf{u}(x, t) &= \frac{1}{\rho(x, t)} \sum_{i} f_i(x, t) \mathbf{e}_i \notag \end{align} $$

4. Taylor-Green渦流れの解析

Taylor–Green vortex の「理論解」は、 粘性ありのナビエ・ストークス方程式の厳密解（特定条件下）として与えられる。 特に有名なのは2次元・非圧縮・周期境界の場合である。

4.1. 理論解

２次元のナビエ・ストークス方程式は以下のように表される。

• 連続の式（発散ゼロ）

$$ \begin{align} \nabla \cdot \mathbf{u} &= 0 \notag \end{align} $$

• 運動方程式

$$ \begin{align} \frac{\partial \mathbf{u}}{\partial t} +(\mathbf{u} \cdot \nabla)\mathbf{u} &= -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} \notag \end{align} $$

Taylor–Green vortexの理論解は次式で表される。

これは上記のナビエ・ストークス方程式を満たす特定の解になっている。

• 速度場

$$ \begin{align} u(x,y,t) &= U_0 \sin(x)\cos(y) e^{-2\nu t} \notag \\ v(x,y,t) &= -U_0 \cos(x)\sin(y) e^{-2\nu t} \notag \end{align} $$

• 圧力場

$$ \begin{align} p(x,y,t) &= -\frac{\rho U_0^2}{4} \left(\cos(2x) + \cos(2y)\right) e^{-4\nu t} \notag \end{align} $$

4.2. LBMでの計算

流速場の初期条件は以下のように指定した。

$$ \begin{align} u(x,y,0) &= U_0 \sin(x)\cos(y) \notag \\ v(x,y,0) &= -U_0 \cos(x)\sin(y) \notag \end{align} $$

解析領域は $[0, 100] \times [0, 100]$ であり、$x,y$領域の両端には周期境界条件を設定する。 計算条件を以下の表に示す。

計算領域の長さLを 周期 $2\pi$ に相当する格子数である NX = 100 とすると、レイノルズ数は以下のようになる。 $$ Re = \frac{U_0 L}{\nu} = \frac{0.1 \times 100}{0.0333…} = 300 $$

4.3. 解析結果と理論解の比較

時刻 $t=5000$ まで解析を行い、流れ場の最大流速を理論解と比較した。 理論解は次式で表される。

$$ U_{max}(t) = U_0 e^{-2\nu\left(k_x^2 + k_y^2 \right) t} $$

ここで、$k_x$ と $k_y$ は空間周波数であり、今回のケースでは $k_x = k_y = 2\pi/NX$ となる。

流速場の渦構造はそのまま保たれ、振幅だけが指数関数的に減衰していく解になっている。

図4 — 最大流速の時間履歴（LBM計算結果 vs 理論解）

💡 LBMの計算結果と理論解はよく一致しており、LBMがナビエ・ストークス方程式を正しく近似できていることが確認できる。

時刻 $t=3000$ での流速ベクトル図は以下のようになった。

図5 — 時刻 t=3000 における流速ベクトル図

4.4. 解析に使用したコード

解析は以下のC++コードを用いて行った。

#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>

// D2Q9 モデルの定数
const int NX = 100;
const int NY = 100;
const int Q = 9;
const double TAU = 0.6;
const int STEPS = 5000;
const double U0 = 0.1;

const double w[Q] = {4.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0};
const int cx[Q] = {0, 1, 0, -1, 0, 1, -1, -1, 1};
const int cy[Q] = {0, 0, 1, 0, -1, 1, 1, -1, -1};

struct Lattice {
    double f[Q];
    double ftmp[Q];
    double rho, ux, uy;
};

// 動粘性係数の計算 (LBMの定義)
double get_viscosity() {
    return (TAU - 0.5) / 3.0;
}

// Taylor-Green Vortex 初期化 (t=0)
void init_tgv(std::vector<std::vector<Lattice>>& grid) {
    double kx = 2.0 * M_PI / NX;
    double ky = 2.0 * M_PI / NY;

    for (int i = 0; i < NX; ++i) {
        for (int j = 0; j < NY; ++j) {
            // 解析解に基づく初期流速
            grid[i][j].ux = U0 * sin(kx * i) * cos(ky * j);
            grid[i][j].uy = -U0 * cos(kx * i) * sin(ky * j);
            grid[i][j].rho = 1.0;

            for (int k = 0; k < Q; ++k) {
                double cu = cx[k] * grid[i][j].ux + cy[k] * grid[i][j].uy;
                double u2 = grid[i][j].ux * grid[i][j].ux + grid[i][j].uy * grid[i][j].uy;
                grid[i][j].f[k] = w[k] * grid[i][j].rho * (1.0 + 3.0*cu + 4.5*cu*cu - 1.5*u2);
            }
        }
    }
}

int main() {
    std::vector<std::vector<Lattice>> grid(NX, std::vector<Lattice>(NY));
    init_tgv(grid);

    double kx = 2.0 * M_PI / NX;
    double ky = 2.0 * M_PI / NY;
    double nu = get_viscosity();

    for (int t = 0; t <= STEPS; ++t) {
        // 1. 衝突ステップ & マクロ量計算
        for (int i = 0; i < NX; ++i) {
            for (int j = 0; j < NY; ++j) {
                double r = 0, usx = 0, usy = 0;
                for (int k = 0; k < Q; ++k) {
                    r += grid[i][j].f[k];
                    usx += grid[i][j].f[k] * cx[k];
                    usy += grid[i][j].f[k] * cy[k];
                }
                grid[i][j].rho = r;
                grid[i][j].ux = usx / r;
                grid[i][j].uy = usy / r;

                for (int k = 0; k < Q; ++k) {
                    double cu = cx[k] * grid[i][j].ux + cy[k] * grid[i][j].uy;
                    double u2 = grid[i][j].ux * grid[i][j].ux + grid[i][j].uy * grid[i][j].uy;
                    double feq = w[k] * r * (1.0 + 3.0*cu + 4.5*cu*cu - 1.5*u2);
                    grid[i][j].ftmp[k] = grid[i][j].f[k] - (grid[i][j].f[k] - feq) / TAU;
                }
            }
        }

        // 2. ストリーミングステップ
        for (int i = 0; i < NX; ++i) {
            for (int j = 0; j < NY; ++j) {
                for (int k = 0; k < Q; ++k) {
                    int ni = (i + cx[k] + NX) % NX;
                    int nj = (j + cy[k] + NY) % NY;
                    grid[ni][nj].f[k] = grid[i][j].ftmp[k];
                }
            }
        }

        // 減衰のモニタリング (中心付近の流速を表示)
        if (t % 100 == 0) {
            double u_theory = U0 * exp(-nu * (kx*kx + ky*ky) * t);
            std::cout << "Step: " << t
                      << " | Max Ux (Sim): " << grid[NX/4][0].ux
                      << " | Theory: " << u_theory << std::endl;
        }

        // 計算結果のファイル出力
        if (t % 1000 == 0) {
            std::string filename = "result" + std::to_string(t) + ".dat";
            std::ofstream ofs(filename);
            for (int i = 0; i < NX; ++i) {
                for (int j = 0; j < NY; ++j) {
                    ofs << i << " " << j << " " << grid[i][j].ux << " " << grid[i][j].uy << "\n";
                }
            }
        }

    }

    return 0;
}

流速ベクトル図は以下のPythonコードを用いて描画した。

import numpy as np
import matplotlib.pyplot as plt

# 1. データの読み込み (data.txt は 4列のテキストファイルを想定)
# 形式: x1, y1, u1, v1
data = np.loadtxt('result3000.dat')

nx = 100  # x方向の格子数（適宜変更してください）
ny = 100  # y方向の格子数（適宜変更してください）

# 2. 各成分を抽出して2次元配列に変換
X = data[:, 0].reshape((ny, nx))
Y = data[:, 1].reshape((ny, nx))
U = data[:, 2].reshape((ny, nx))
V = data[:, 3].reshape((ny, nx))

# 流速の絶対値を計算（色付け用）
speed = np.sqrt(U**2 + V**2)

plt.figure(figsize=(8, 8))

# データを間引いて表示
skip = 1
plt.quiver(X[::skip, ::skip], Y[::skip, ::skip],
           U[::skip, ::skip], V[::skip, ::skip],
           speed[::skip, ::skip], cmap='jet', scale=1)

plt.colorbar(label='Velocity Magnitude')
plt.title("Velocity Vectors (Quiver)")
plt.xlabel("X")
plt.ylabel("Y")
plt.axis('equal')
plt.show()

5. まとめ

今回、格子ボルツマン法（LBM）を用いて、テイラーグリーン渦流れの解析を行った。

LBMは、流体を仮想的な粒子の集まりとして捉えることで、 従来のナビエ・ストークス方程式を直接解く手法に代わる新しいアプローチを提供する。

今回の解析では、LBMを用いてテイラーグリーン渦流れの理論解を再現することができ、 流速の減衰挙動も理論解と良く一致することが確認された。

LBMは、並列計算との親和性が高く、複雑な境界形状への対応が容易であるなどの長所がある一方で、 格子依存性や高マッハ数への制限などの短所も存在する。

今後は、より高レイノルズ数の流れや、熱流体計算など、LBMの適用範囲を広げるための研究が期待される。

弊社の解析事例

弊社の流体解析事例については、下記のリンクからご覧ください。

👉 解析事例のご紹介 OpenFOAM

👉 流体解析 実績及び事例

• 格子ボルツマン法
• LBM
• Taylor-Green渦