【格子ボルツマン法】格子ボルツマン法で2次元拡散方程式を解きます(D2Q4) C++コード付き
今回から2次元に入っていきます。2次元拡散方程式を解くC++コードを作りました。D2Q4です。どうやら拡散方程式を解く分にはこれで十分みたいです。
参考にしたのは"Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes"の2-4章とAppendixのコード(まったく洗練されていないが故に読みやすい不思議なコード)です。
Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes
- 作者: A. A. Mohamad
- 出版社/メーカー: Springer
- 発売日: 2019/05/16
- メディア: ハードカバー
- この商品を含むブログを見る
今回解く2次元拡散方程式は
とあらわされます。領域を長方形(正方形)として、境界条件は左側にディリクレで1、右側にディリクレで0、上側にディリクレで0、下側に斉次のノイマンです。
計算結果です。
うん、良さそうです。
以下にコードを示します。久しぶりに3次元配列使いました。あと本のコードは端の境界条件が怪しいですね。きちんと場合分けすべきだと私は思いました。コードではそうしてあります。
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <string> #include <sstream> using namespace std; int main() { int xNode=101; int yNode=101; double f[4][xNode][yNode];//f_k at (i,j) double phi[xNode][yNode]; double feq; double dt=1.0; double dx=1.0; double dy=dx; double csq=dx*dx/(dt*dt); double alpha=0.25; double omega=1.0/(2.0*alpha/(dt*csq)+0.5); int step=500; double DL=1.0;//dirichlet, left double DR=0.0;//dirichlet, right double DT=0.0;//dirichlet, top double DB=0.0;//dirichlet, bottom double w[4]; w[0]=0.25; w[1]=0.25; w[2]=0.25; w[3]=0.25; double c1=dx/dt; double c2=-dx/dt; double cs=1.0/sqrt(2.0)*dx/dt; ofstream fk; fk.open("phi.txt"); //initialization// for(int i=0;i<xNode;i++) { for(int j=0;j<yNode;j++) { phi[i][j]=0.0; for(int k=0;k<4;k++) { f[k][i][j]=w[k]*phi[i][j]; } } } //main loop// for(int k=0;k<step;k++) { cout<<k<<endl; //collision process// for(int i=0;i<xNode;i++) { for(int j=0;j<yNode;j++) { for(int k=0;k<4;k++) { feq=w[k]*phi[i][j]; f[k][i][j]=(1.0-omega)*f[k][i][j]+omega*feq; } } } //Streaming process// //x// for(int i=1;i<xNode;i++) { for(int j=0;j<yNode;j++) { f[0][xNode-i][j]=f[0][xNode-i-1][j]; f[1][i-1][j]=f[1][i][j]; } } //y// for(int i=0;i<xNode;i++) { for(int j=1;j<yNode;j++) { f[2][i][yNode-j]=f[2][i][yNode-j-1]; f[3][i][j-1]=f[3][i][j]; } } //BC// for(int j=1;j<yNode-1;j++) { //left// f[0][0][j]=w[0]*DL+w[1]*DL-f[1][0][j]; //right// f[1][xNode-1][j]=w[0]*DR+w[1]*DR-f[0][xNode-1][j]; } //left upper corner// f[0][0][yNode-1]=w[0]*DL+w[1]*DL-f[1][0][yNode-1]; f[3][0][yNode-1]=w[2]*DL+w[3]*DL-f[2][0][yNode-1]; //right upper corner// f[1][xNode-1][yNode-1]=w[0]*DR+w[1]*DR-f[0][xNode-1][yNode-1]; f[3][xNode-1][yNode-1]=w[2]*DR+w[3]*DR-f[2][xNode-1][yNode-1]; for(int i=1;i<xNode-1;i++) { //top// f[3][i][yNode-1]=w[2]*DT+w[3]*DT-f[2][i][yNode-1]; //bottom// for(int k=0;k<4;k++) f[k][i][0]=f[k][i][1]; } //left lower corner// f[0][0][0]=w[0]*DL+w[1]*DL-f[1][0][0]; f[2][0][0]=w[2]*DL+w[3]*DL-f[3][0][0]; //right lower corner// f[1][xNode-1][0]=w[0]*DR+w[1]*DR-f[0][xNode-1][0]; f[2][xNode-1][0]=w[2]*DR+w[3]*DR-f[2][xNode-1][0]; //summation// for(int i=0;i<xNode;i++) { for(int j=0;j<yNode;j++) { double sum=0.0; for(int k=0;k<4;k++) { sum+=f[k][i][j]; } phi[i][j]=sum; } } } //output// for(int i=0;i<xNode;i++) { for(int j=0;j<yNode;j++) { fk<<double(i)*dx<<" "<<double(j)*dy<<" "<<phi[i][j]<<endl; } fk<<endl; fk<<endl; } return 0; }
次は『格子ボルツマン法・差分格子ボルツマン法』を読みながらD2Q9作ります。
- 作者: 蔦原道久
- 出版社/メーカー: コロナ社
- 発売日: 2018/12/01
- メディア: 単行本
- この商品を含むブログを見る