はじめに
共同研究を行っている大学において、しばしば孔隙内の油・ガス・水の流動を解析するために格子ボルツマン法(Lattice Boltzaman Mehtod, LBM)が使われているので、勉強もかねてC++でLBMシミュレータを作ってみました。
使用しているライブラリは以下の通りです。
- 行列計算：Eigen
- JSONファイルの読み込み：JSON for Modern C++
- コマンドライン解析：CLI11
- フォーマット：{fmt}
- 結果の可視化：Jupyter Notebook + matplotlib.pyplot
参考文献は以下の2つです。
稲室・吉野・鈴木の「格子ボルツマン法入門」は、基礎から丁寧に解説されています。一方、瀬田の「格子ボルツマン法」は応用面について詳しいです。お互いが補完しあっているので、両方共読んでおくといいでしょう。
C++での実装
ここでは瀬田の「格子ボルツマン法」の第2章で解説されている二次元Poiseuille流れをC++で実装しました。著者もC言語で書かれたプログラムを配布していますが、Windowsでコンパイルするには色々と修正が必要なのでご注意ください。
以下では、衝突過程、並行分布関数、外力項の実装方法のみを解説しています。並進過程、境界条件については簡単なので除外しています。
衝突過程
単一緩和時間(Single Relaxation Time, SRT)衝突則は次式で表されます。
f_k(\boldsymbol{x} + \boldsymbol{c} _ k \delta_t, t + \delta _ t ) =
f_k (\boldsymbol{x}, t) - \frac{1}{\tau} \left( f_k (\boldsymbol{x}, t) - f_k^\mathrm{eq} (\boldsymbol{x}, t) \right)
ここで$\boldsymbol{f} _ k$は分布関数(distribution function)、$\boldsymbol{f} _ k ^\mathrm{eq}$は平衡分布関数(equilibrium distribution function)、$\tau$は緩和時間(relaxation time)、$\boldsymbol{x}$は位置ベクトル、$\boldsymbol{c} _ k$は離散速度ベクトル、$\delta _ t$は時間ステップ幅、下付き文字の$k$は格子ベクトル番号を表します。
右辺の衝突過程の結果の分布関数を$f _ k ^*$とすると
$$
f_k^* (\boldsymbol{x}, t) =
f_k (\boldsymbol{x}, t) - \frac{1}{\tau} \left( f_k (\boldsymbol{x}, t) - f_k^\mathrm{eq} (\boldsymbol{x}, t) \right)
$$
となります。これをEigenライブラリを用いて以下のようにC++で実装します。
まずEigen::Matrixクラスを用いて分布関数$f _ k$と平衡分布関数$f _ k^\mathrm{eq}$を$Q \times N$行列として定義します。ここで$Q$は格子ベクトル数、$N$はグリッド数です。
using Eigen::Matrix, Eigen::Dynamic;
Matrix<double, Q, Dynamic> f(Q, N);
Matrix<double, Q, Dynamic> feq(Q, N);
すると衝突過程は行列演算として次のように表すことができます。
Matrix<double, Q, Dynamic> fstar = f - (f - feq)/tau;
ここでtauはSRTの緩和時間$\tau$です。
平衡分布関数
次に平衡分布関数の計算に用いられる式
$$
f _ k^\mathrm{eq} = \omega _ k \rho \left[ 1 + 3 \boldsymbol{c} _ k \cdot \boldsymbol{u} + \frac{9}{2} \left( \boldsymbol{c} _ k \cdot \boldsymbol{u} \right)^2 - \frac{3}{2} \boldsymbol{u} \cdot \boldsymbol{u} \right]
$$
について、重み係数ベクトル$\omega _ k$を$Q \times 1$行列、離散速度ベクトル$\boldsymbol{c} _ k$を$D \times Q$行列、速度ベクトルを$D \times N$行列、密度を$N \times 1$行列として表します。ここで$D$は次元です。固定長行列はEigen::Matrixクラスを使って表すことができます。
using Eigen::VectorXd;
Matrix<double, Q, 1> w;
Matrix<double, D, Q> c;
Matrix<double, D, Dynamic> u(D, N);
VectorXd rho(N);
これらを用いると平衡分布関数の計算を以下の様に表すことができます。
const Matrix<double, Q, Dynamic> cu = c.transpose() * u;
const Matrix<double, 1, Dynamic> u2 = u.colwise().squaredNorm();
feq = ((w * rho.transpose()).array() * (
1.0 + 3.0 * cu.array() +
4.5 * cu.array().square() -
1.5 * u2.replicate<Q, 1>().array())
).matrix();
ただこのように表すとforループでべた書きするよりもパフォーマンスが低かったです。
おそらく最初の直積w * rho.transpose()で一時変数が生成されているのではないかと推察していますが、アセンブリまで確認したわけではありません。今後の課題ですね。
外力項
外力項$F _ k$が
$$
F _ k = 3 \omega _ k \rho \boldsymbol{G} \cdot \boldsymbol{c}_k
$$
で定式化されているとします。ここで$\boldsymbol{G}$は体積力ベクトルです。体積力ベクトル$\boldsymbol{G}$を$D \times 1$行列、外力項$F _ k$を$Q \times N$行列として実装します。
Matrix<double, D, 1> g; // 体積力ベクトル
Matrix<double, Q, Dynamic> fg(Q, N); // 外力項
これを用いると外力項を次式のように計算できます。
Matrix<double, Q, 1> tmp = 3.0 * w.cwiseProduct(c.transpose() * g);
fg = tmp * rho.transpose();
コード
実際のコードはレポジトリをご確認ください。パフォーマンス向上のために若干修正している場合があります。
結果
$Y$軸に沿った$x$方向の流速$u_x$を流速の最大値$u_c$
$$
u_c = \frac{G_x h^2}{8 \nu}
$$
で正規化したプロットを作りました。ここで$G_x$は$x$方向の外力、$\nu$は動粘性係数です。Y軸は縦の長さ$h = n_y - 1$（$n_y$はY方向のセル数）で正規化してあります。
瀬田の「格子ボルツマン法」の第2章と同じ結果を得られました。
余談
本論とは関係ありませんが、プログラムをGitHub等で公開せずにダウンロード形式で配布するのはやめて欲しいと思っています。