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

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

【格子ボルツマン法】格子ボルツマン法で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のコード(まったく洗練されていないが故に読みやすい不思議なコード)です。

混相流の数値シミュレーション

混相流の数値シミュレーション

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} = \alpha \frac{\partial^2 T}{\partial x^2}

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

LBMでは、速度分布関数  f_i i は方向をあらわす添え字)の時間発展を以下の格子ボルツマン方程式によって計算します。


 \displaystyle f_i (x+\Delta x, t+\Delta t) = f_i (x,t) (1-\omega) + \omega f_i^{eq} (x,t)

ここで、 \Delta x = c_i \Delta t c_i i 方向の速度、 \omega は緩和時間、 f_i^{eq} は平衡分布関数です。どうやらこの平衡分布関数を変えるだけで、色々な問題(拡散、移流拡散、ナヴィエ・ストークス方程式)などが計算できるようです。上記の格子ボルツマン方程式ではBGK(Bhatnagar-Gross-Krook)近似にもとづく衝突則を採用しています。これがスタンダードみたいですね。

この格子ボルツマン方程式によって速度分布関数の時間発展を計算し、速度分布関数を用いてマクロな物理量(密度や運動量)を計算します。


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

f:id:mutsumunemitsutan:20190510000841p:plain:w500

両端斉次の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;
}


取りあえず合ってそうな動くものが作れました。ただ理論はまだ全然です。統計力学の闇が垣間見えます…理論に関しても『格子気体法・格子ボルツマン法』や『格子ボルツマン法・差分格子ボルツマン法』を読んで勉強していくことにしましょう。

格子気体法・格子ボルツマン法―新しい数値流体力学の手法

格子気体法・格子ボルツマン法―新しい数値流体力学の手法

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

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