C++による格子ボルツマン法を用いた二次元キャビティ流れのシミュレーション

はじめに

前回の記事の続きで、格子ボルツマン法を用いて二次元キャビティ流れをシミュレーションするプログラムをC++で作成してみました。

使用しているライブラリは前回同様以下の通りです。

• 行列計算：Eigen
• JSONファイルの読み込み：JSON for Modern C++
• コマンドライン解析：CLI11
• フォーマット：{fmt}
• 結果の可視化：Jupyter Notebook + matplotlib.pyplot

参考文献

C++での実装

境界条件以外の部分は、外力項がなくなっただけで、前回とほとんど変わりません。詳細はcavity_flow_simulator.cppのlbm::CavityFlowSimulator::run関数を参照してください。

境界条件については、瀬田の「格子ボルツマン法」の2.7節と同様にハーフウェイバウンスバックスキームを用いています。実際にプログラムを書いていて、異なる境界条件の影響を受ける角のセルの処理方法について本の中で解説されていなかったことが気になりました。色々試したところあまり影響がなかったので、最終的に特別な処理は行わないこととしました。

衝突則については、Single Relaxation Time (SRT) modelとMultiple Relaxation Time (MRT) modelを実装しました。

結果

セル数$n _ x = n _ y = 52$、上面の壁の速度を$u_x=0.1$、レイノルズ数を$\mathit{R\kern-.04em e}=400$として、MRT衝突則を用いて計算したときの速度分布が下図となります。

$x=0.5$および$y=0.5$における$x$・$y$方向速度をプロットしたものが下図となります。

Ghia et al. (1982)¹とほぼ同じ結果が得られていることがわかります。

SRT衝突則を用いたときにVisual Studio 2022のパフォーマンスプロファイラーを使って分析した結果は下図です。SIMDやOpenMPは有効にしていない状態での測定結果です。

平衡分布関数(calc_equilibrium_distribution_function)および並進過程(run_propagation_process)の計算が全体の58%を占めており、格子ボルツマン法ではこれらの2つの部分の最適化がパフォーマンスの観点で重要であることがわかります。

1. V.Ghia, K.N.Ghia and C.T.Shin. 1982. "High-Re solutions for incompressible flow using the Navier-Stokes equations and a multi-grid method", J. Comp. Phys., vol. 48, pp. 387-411 ↩