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

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

【格子ボルツマン法】格子ボルツマン法で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

Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes


今回解く2次元拡散方程式は

 \displaystyle \frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right)

とあらわされます。領域を長方形(正方形)として、境界条件は左側にディリクレで1、右側にディリクレで0、上側にディリクレで0、下側に斉次のノイマンです。


計算結果です。

f:id:mutsumunemitsutan:20190515222141p:plain:w500

うん、良さそうです。


以下にコードを示します。久しぶりに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作ります。

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

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