【有限体積法】1次元移流方程式を風上差分で解く　C++コード付き

今回は1次元移流方程式を有限体積法（風上差分）で解いていきます。C++コード付きです。

1次元移流方程式とは

$\displaystyle \frac{\partial u}{\partial t} + V\frac{\partial u}{\partial x} =0, \quad x \in [0, L]$

のような一階の偏微分方程式のことをいいます。ここで、 $u$ は溶質の濃度、 $V$ は風などによる移流速度です。

移流方程式はその名前の通り、初期の波形（プロファイル）が $V$ の速度で波形を変えずに移動していきます。移流方程式はこのように非常に直観的な解を持つのですが、これを数値計算で扱うのは非常に難しいです。いかに移流方程式をきとんと解くか、これがCFD（Computational Fluid Dynamics, 数値流体力学）における究極的な問いです。移流方程式をきちんと解くために今まで多くの労力がさかれてきました。なので、CFDを学ぼうと考えている方は移流方程式の数値解法をきちんと勉強することをおすすめします。今までにたくさんのスキームが移流方程式を解くために作られてきたのです。私も今後そうしたスキーム達を紹介していきます。

さて、風上差分というのは、特性曲線の向き（情報の向き）を考えることによって移流方程式を離散化しようという試みです。移流速度が正の場合は波形は左からくるので左のステンシルの値を用いるようにします。移流速度が負の場合はその逆です。詳しくは別記事で説明します。お待ちください。

境界条件は2パターンあります。すなわち、

Case 1, $\boldsymbol{x=0}$ でDirichlet条件を課す
（ $u(0)$ の値を指定）
Case 2, $\boldsymbol{x=L}$ でDirichlet条件を課す
（ $u(L)$ の値を指定）

です。それぞれコード上ではbc=1、bc=2に対応しています。特性曲線の向きを考えるとわかるのですが、 $\boldsymbol{x=0}$ でDirichlet条件を課すのは移流速度が正の場合です。直観的には波が境界の外から内部へ入って来る様子を表しています。なので、 $\boldsymbol{x=L}$ では境界条件は必要ありません。これが移流方程式の境界条件の独特な部分です。逆に移流速度が負の場合は $x=L$ でDirichlet条件を課します。

今回の計算条件は、移流速度が正 $V=0.1$ で、分割数は400、時間刻みは0.001、領域の広さは $L=1$ 、初期条件は領域の左端の矩形波（長方形の波）です。

以下では2つコードを示します。絶対値を使うものとC++の組み込み関数maxを使うコードです。どちらを使っても大丈夫です。絶対値を使う方法は離散化方程式を理論的に調べる際によく使われます。一方maxを使うとプログラムがとても簡単に書け、読みやすくなります。絶対値はこちらの記事を見てください。

まずは絶対値のコードです。

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

using namespace std;


inline void solve(double x[], int &Node, double &dt, double &dx, double &V)
{
	int i;
	
	for(i=1;i<Node-1;i++)
	{
		x[i]=x[i]-dt/dx*(-(V+abs(V))/2.0*x[i-1]+abs(V)*x[i]+(V-abs(V))/2.0*x[i+1]);
	}

}

inline void boundary(int &bc, double &D0, double x[], double &V, int &Node, double &dt, double &dx)
{
	if(bc==1)
	{
		x[0]=D0;
		x[Node-1]=x[Node-1]-dt/dx*V*(x[Node-1]-x[Node-2]);
	}
	if(bc==2)
	{
		x[0]=x[0]-dt/dx*V*(x[1]-x[0]);
		x[Node-1]=D0;
	}
}

inline void text(int &i,double x[], double &dx, int &Node)
{
	int j;
	
	stringstream ss;
	string name;
	ofstream fo;
	
	ss<<i;
	name=ss.str();
	name="Answer_" + name + ".txt";
	fo.open(name.c_str ());
	
	for(j=0;j<Node;j++)
	{
		fo<<dx*float(j)<<" "<<x[j]<<endl;
	}
}

int main()
{
	int i;
	int Partition=400;
	int Node=Partition+1;
	double LL=1.0;
	double dx=LL/(Node-1);
	double dt=0.001;
	double NT=7000;
	double eps=pow(2.0,-50);
	double V=0.1;
	int bc=1;//bc=1 left Dirichlet (V>=0), bc=2 right Dirichlet (V<0)
	double D0=0.0;//for Dirichlet
	double *x=new double[Node];
	
	ofstream fk;
	fk.open("Answer_0.txt");
	
	//initial condition//
	for(i=0;i<Node;i++)
	{
		x[i]=0.0;
		if((i>=40)&&(i<=80)) x[i]=1.0;
		fk<<dx*float(i)<<" "<<x[i]<<endl;
	}
	
	for(i=1;i<=NT;i++)
	{
		solve(x,Node,dt,dx,V);
		boundary(bc,D0,x,V,Node,dt,dx);
		
		if(i%1000==0)
		{
			cout<<i<<endl;
			text(i,x,dx,Node);
		}
	}
	
	delete[] x;
	
	return 0;
}

次にmaxを用いるコードです。

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

using namespace std;


inline void solve(double x[], int &Node, double &dt, double &dx, double &V)
{
	int i;
	
	for(i=1;i<Node-1;i++)
	{
		x[i]=x[i]-dt/dx*(-max(V,0.0)*x[i-1]+(max(V,0.0)+max(-V,0.0))*x[i]-max(-V,0.0)*x[i+1]);
	}

}

inline void boundary(int &bc, double &D0, double x[], double &V, int &Node, double &dt, double &dx)
{
	if(bc==1)
	{
		x[0]=D0;
		x[Node-1]=x[Node-1]-dt/dx*V*(x[Node-1]-x[Node-2]);
	}
	if(bc==2)
	{
		x[0]=x[0]-dt/dx*V*(x[1]-x[0]);
		x[Node-1]=D0;
	}
}

inline void text(int &i,double x[], double &dx, int &Node)
{
	int j;
	
	stringstream ss;
	string name;
	ofstream fo;
	
	ss<<i;
	name=ss.str();
	name="Answer_" + name + ".txt";
	fo.open(name.c_str ());
	
	for(j=0;j<Node;j++)
	{
		fo<<dx*float(j)<<" "<<x[j]<<endl;
	}
}

int main()
{
	int i;
	int Partition=400;
	int Node=Partition+1;
	double LL=1.0;
	double dx=LL/(Node-1);
	double dt=0.001;
	double NT=7000;
	double eps=pow(2.0,-50);
	double V=0.1;
	int bc=1;//bc=1 left Dirichlet (V>=0), bc=2 right Dirichlet (V<0)
	double D0=0.0;//for Dirichlet
	double *x=new double[Node];
	
	ofstream fk;
	fk.open("Answer_0.txt");
	
	//initial condition//
	for(i=0;i<Node;i++)
	{
		x[i]=0.0;
		if((i>=40)&&(i<=80)) x[i]=1.0;
		fk<<dx*float(i)<<" "<<x[i]<<endl;
	}
	
	for(i=1;i<=NT;i++)
	{
		solve(x,Node,dt,dx,V);
		boundary(bc,D0,x,V,Node,dt,dx);
		
		if(i%1000==0)
		{
			cout<<i<<endl;
			text(i,x,dx,Node);
		}
	}
	
	delete[] x;
	
	return 0;
}

計算結果を見てみましょう。

f:id:mutsumunemitsutan:20181111231051p:plain:w500

初期条件である左端の矩形波が時間がすすむにつれて右へと流されていきます。ただし、注意して頂きたいのは、厳密解だと初期波形がそのまま形を保ったまま流されていくはずですが、数値計算ではだんだんと角がとれてなまってきている点です。これは数値拡散と呼ばれています。風上差分は非常に安定して計算を行うことができるのですが、精度が1次しかなく、数値拡散が大きすぎます。そのため、この欠点を克服するために今までたくさんのスキームが作られてきました。

離散化の詳細は後でまとめます。