【格子ボルツマン法】格子ボルツマン法で1次元拡散方程式を解きます(D1Q2) C++コード付き
幅を広げようと思い、格子ボルツマン法(Lattice Boltzmann Method, LBM)の勉強を始めました!LBMは流体の方程式に対する離散化手法です。
そこで、取りあえず1次元拡散方程式を解くC++コードを作りました。D1Q2というやつです。Dの後の数字が次元で、Qの後の数字が粒子の自由度です。なのでこの場合流体粒子は右か左に進むだけ、となります。
参考にしたのは『混相流の数値シミュレーション』の5章と"Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes"の2章と3章、Appendixのコード(まったく洗練されていないが故に読みやすい不思議なコード)です。
- 作者: 太田光浩,酒井幹夫,島田直樹,本間俊司,松隈洋介,気泡・液滴・微粒子分散工学分科会
- 出版社/メーカー: 丸善出版
- 発売日: 2015/07/27
- メディア: 単行本(ソフトカバー)
- この商品を含むブログを見る
Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes
- 作者: A. A. Mohamad
- 出版社/メーカー: Springer
- 発売日: 2019/05/16
- メディア: ハードカバー
- この商品を含むブログを見る
今回解く1次元拡散方程式は
とあらわされます。境界条件は簡単のため両端Dirichlet条件とします。初期条件はsin波で与えます。
LBMでは、速度分布関数 ( は方向をあらわす添え字)の時間発展を以下の格子ボルツマン方程式によって計算します。
ここで、 で は 方向の速度、 は緩和時間、 は平衡分布関数です。どうやらこの平衡分布関数を変えるだけで、色々な問題(拡散、移流拡散、ナヴィエ・ストークス方程式)などが計算できるようです。上記の格子ボルツマン方程式ではBGK(Bhatnagar-Gross-Krook)近似にもとづく衝突則を採用しています。これがスタンダードみたいですね。
この格子ボルツマン方程式によって速度分布関数の時間発展を計算し、速度分布関数を用いてマクロな物理量(密度や運動量)を計算します。
まず計算結果を示します。
両端斉次のDirichlet条件にしています。初期に与えたsin波が拡散していく様子が確認できます。
以下にコードを示します。
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <string> #include <sstream> using namespace std; int main() { int Node=101; double f1[Node]; double f2[Node]; double x[Node]; double feq1[Node]; double feq2[Node]; double dt=1.0; double dx=1.0; double csq=dx*dx/(dt*dt); double alpha=0.5; double omega=1.0/(alpha/(dt*csq)+0.5); int step=3000; double DL=0.0;//dirichlet, left double DR=0.0;//dirichlet, right double w1=0.5; double w2=0.5; ofstream fk; fk.open("x.txt"); //initialization// for(int i=0;i<Node;i++) { x[i]=sin(atan(1.0)*4.0*double(i)*dx/(double(Node-1)*dx)); f1[i]=w1*x[i]; f2[i]=w2*x[i]; fk<<double(i)*dx<<" "<<x[i]<<endl; } fk<<endl; fk<<endl; //main loop// for(int k=0;k<step;k++) { //collision process// for(int i=0;i<Node;i++) { x[i]=f1[i]+f2[i]; feq1[i]=w1*x[i]; feq2[i]=w2*x[i]; f1[i]=(1.0-omega)*f1[i]+omega*feq1[i]; f2[i]=(1.0-omega)*f2[i]+omega*feq2[i]; } //Streaming process// for(int i=1;i<Node;i++) { f1[Node-i]=f1[Node-i-1]; f2[i-1]=f2[i]; } //BC// f1[0]=DL-f2[0];//dirichlet, left f2[Node-1]=DR-f1[Node-1];//dirichlet, right //f1[Node-1]=f1[Node-2];//neumann, right //f2[Node-1]=f2[Node-2];//neumann, right //output// if(k%100==0) { for(int i=0;i<Node;i++) { fk<<double(i)*dx<<" "<<x[i]<<endl; } fk<<endl; fk<<endl; } } return 0; }
取りあえず合ってそうな動くものが作れました。ただ理論はまだ全然です。統計力学の闇が垣間見えます…理論に関しても『格子気体法・格子ボルツマン法』や『格子ボルツマン法・差分格子ボルツマン法』を読んで勉強していくことにしましょう。
- 作者: 蔦原道久,片岡武,高田尚樹
- 出版社/メーカー: コロナ社
- 発売日: 1999/08
- メディア: 単行本
- クリック: 14回
- この商品を含むブログ (1件) を見る
- 作者: 蔦原道久
- 出版社/メーカー: コロナ社
- 発売日: 2018/12/01
- メディア: 単行本
- この商品を含むブログを見る