数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。noteにも書いています。https://note.mu/kumag_min

【格子ボルツマン法】格子ボルツマン法で1次元移流拡散方程式を解きます(D1Q2) C++コード付き

LBMで1次元移流拡散方程式を解くC++コードを作りました。D1Q2でやってます。

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

Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes

Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes


今回解く1次元移流拡散方程式は

 \displaystyle \frac{\partial T}{\partial t} = -u\frac{\partial T}{\partial x} + \alpha \frac{\partial^2 T}{\partial x^2}

とあらわされます。境界条件は簡単のため両端Dirichlet条件とします。初期条件はsin波で与えます。

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

f:id:mutsumunemitsutan:20190511013958p:plain:w500

移流速度  u = 0.1 で一定とし、両端斉次の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 V[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=500;
	double DL=0.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));
		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*V[i]/cs/cs);
			feq2[i]=w2*x[i]*(1.0+c2*V[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%100==0)
		{
			for(int i=0;i<Node;i++)
			{
				fk<<double(i)*dx<<" "<<x[i]<<endl;
			}
			fk<<endl;
			fk<<endl;
		}
	}
	
	return 0;
}


次書くのはBurgers方程式でしょうか。

『格子ボルツマン法・差分格子ボルツマン法』を読み始めました。べらぼうに面白いですね!Streaming processの移流方程式をTVDとか使って解くと面白いんじゃないかと考えています。

格子ボルツマン法・差分格子ボルツマン法- DVD付き -

格子ボルツマン法・差分格子ボルツマン法- DVD付き -