【格子ボルツマン法】格子ボルツマン法で1次元バーガース方程式を解きます（D1Q2）　C++コード付き

LBMで1次元バーガース方程式を解くC++コードを作りました。D1Q2でやってます。

参考にしたのは"Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes"の2-5章、Appendixのコード（まったく洗練されていないが故に読みやすい不思議なコード）です。

Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes

作者: A. A. Mohamad
出版社/メーカー: Springer
発売日: 2019/05/16
メディア: ハードカバー
この商品を含むブログを見る

今回解く1次元バーガース方程式は

$\displaystyle \frac{\partial u}{\partial t} = -u\frac{\partial u}{\partial x} + \alpha \frac{\partial^2 u}{\partial x^2}$
とあらわされます。境界条件は簡単のため両端Dirichlet条件とします。初期条件はsin波で与えます。

まず計算結果を示します。

f:id:mutsumunemitsutan:20190511124841p:plain:w500

両端斉次のDirichlet条件にしています。初期に与えたsin波が右に流れながら拡散していく様子が確認できます。その際解の値が大きいほど移流速度が大きいこともわかります。

以下にコードを示します。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>

using namespace std;

int main()
{
	int Node=1001;
	double f1[Node];
	double f2[Node];
	double x[Node];
	double feq1[Node];
	double feq2[Node];
	double V[Node];
	double dt=0.01;
	double dx=0.1;
	double csq=dx*dx/(dt*dt);
	double alpha=0.5;
	double omega=1.0/(alpha/(dt*csq)+0.5);
	int step=5000;
	double DL=1.0;//dirichlet, left
	double DR=0.0;//dirichlet, right
	double w1=0.5;
	double w2=0.5;
	double c1=dx/dt;
	double c2=-dx/dt;
	double cs=1.0/sqrt(2.0)*dx/dt;
	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));
		//if(i<Node/3) x[i]=1.0;
		//else x[i]=0.0;
		f1[i]=w1*x[i];
		f2[i]=w2*x[i];
		V[i]=0.1;
		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]*(1.0+c1*x[i]/cs/cs);
			feq2[i]=w2*x[i]*(1.0+c2*x[i]/cs/cs);
			
			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%500==0)
		{
			for(int i=0;i<Node;i++)
			{
				fk<<double(i)*dx<<" "<<x[i]<<endl;
			}
			fk<<endl;
			fk<<endl;
		}
	}
	
	return 0;
}