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

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

【浅水流方程式】Roe法でダム崩壊問題を解いてみました(wet bedとdry bed) 動画とC++コード付き

【浅水流方程式】サイトマップ(ここから関連記事が探せます)
http://mathlang.hatenablog.com/entry/2018/12/18/213650


今回は、Roe法を用いてダム崩壊問題(dam break problem)を解いてみました。Roe法とは、ヤコビアン固有値の正負によって分割し、Roe平均(Roe average)と呼ばれる特別な平均を用いてそれらを評価し、正負に合わせて風上化をかける手法です。解いている方程式は、摩擦、勾配なしの1次元浅水流方程式(Shallow Water Equations)です。

 \displaystyle
\begin{eqnarray} 
&\frac{\partial h}{\partial t} + \frac{\partial q}{\partial x} = 0 \\
&\frac{\partial q}{\partial t} + \frac{\partial}{\partial x} \left( \frac{q^2}{h} + \frac{1}{2}gh^2 \right) = 0
\end{eqnarray}

初期状態としては、wet bedの場合は、領域のちょうど左側で水深が1.0、領域のちょうど右側で水深が0.1で、水面は静止しているものとします。dry bedの場合は、領域のちょうど左側で水深が1.0、領域のちょうど右側で水深が  10^{-4} で、水面は静止しているものとします。このように、dry bedの場合は非常に小さい水深を与えることとします。時刻  t=0 に真ん中にあるしきり(ダム)を取り去った時に、水がどのように動くかを計算します。境界条件としては右端、左端で水深、流量ともに斉次のノイマン条件を課しています。


参考文献としてはToroの"Shock-Capturing Methods for Free-Surface Shallow Flows"のpp.187-191です。

Shock-Capturing Methods for Free-Surface Shallow Flows

Shock-Capturing Methods for Free-Surface Shallow Flows


こちらがwet bedの場合の計算結果です。水深が青い線で表示されています。

f:id:mutsumunemitsutan:20190605231210g:plain:w500

右へと進行する衝撃波(不連続な波、Shock wave)と左へと進行する膨張波(連続な波、Rarefaction wave)が観察できます。よい感じです。

こちらがdry bedの場合の計算結果です。

f:id:mutsumunemitsutan:20190605232945g:plain:w500

右へと進行する膨張波(連続な波、Rarefaction wave)が観察できます。よい感じです。しかしHLL法と比べて小さい水深がとれません。また、やはり膨張波は苦手のようです。

Youtube


では計算コードです。

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

using namespace std;

const double g=9.8;

//celerity//
inline double a(double h)
{
	return sqrt(g*h);
}

//flux for continuity equation//
inline double Fluxcon(double h, double q)
{
	return q;
}

//flux for momentum equation//
inline double Fluxmom(double h, double q)
{
	return q*q/h+g*h*h/2.0;

}

//Roe average for h
double Roeh(double hL, double hR)
{
	return sqrt(hL*hR);
} 

//Roe average for q
double Roeq(double hL, double hR, double qL, double qR)
{
	return (sqrt(hL)*qL+sqrt(hR)*qR)/(sqrt(hL)+sqrt(hR));
}

//Roe flux//
inline void RoeFlux(double FLR[], double hL, double hR, double qL, double qR)
{
	double rh=Roeh(hL,hR);
	double rq=Roeq(hL,hR,qL,qR);
	double lamp=rq/rh+a(rh);//plus
	double lamm=rq/rh-a(rh);//minus
	double lampR=qR/hR+a(hR);
	double lampL=qL/hL+a(hL);
	double lammR=qR/hR-a(hR);
	double lammL=qL/hL-a(hL);
	
	//entropy fix, Chakravarthy//
	if((lampL<0.0)&&(0.0<lampR)) lamp=lamp/abs(lamp)*(abs(lamp)+(lampR-lampL)/2.0);
	if((lammL<0.0)&&(0.0<lammR)) lamm=lamm/abs(lamm)*(abs(lamm)+(lammR-lammL)/2.0);

	double M11=(lamp*abs(lamm)-abs(lamp)*lamm)/(2.0*sqrt(g*rh));
	double M12=(abs(lamp)-abs(lamm))/(2.0*sqrt(g*rh));
	double M21=(lamp*lamm*(abs(lamm)-abs(lamp)))/(2.0*sqrt(g*rh));
	double M22=(abs(lamp)*lamp-abs(lamm)*lamm)/(2.0*sqrt(g*rh));
	
	double dh=hR-hL;
	double dq=qR-qL;
	
	FLR[0]=(Fluxcon(hR,qR)+Fluxcon(hL,qL)-(M11*dh+M12*dq))/2.0;
	FLR[1]=(Fluxmom(hR,qR)+Fluxmom(hL,qL)-(M21*dh+M22*dq))/2.0;
}

inline void output(double h[], double z[], int Node, double dx)
{
	stringstream ss;
	string name;
	ofstream fo;
	static int count=0;
	
	ss<<count;
	name=ss.str();
	name="h_" + name + ".txt";
	fo.open(name.c_str ());
	
	for(int i=0;i<Node;i++)
	{
		fo<<dx*float(i)<<" "<<h[i]+z[i]<<endl;
	}
	fo<<endl;
	fo<<endl;
	
	for(int i=0;i<Node;i++)
	{
		fo<<dx*float(i)<<" "<<z[i]<<endl;
	}
	
	count+=1;
}

int main()
{
	const int n=200;//number of cells
	const int Node=n+1;
	double L=1.0;
	double dx=L/double(n);
	double dt=0.0001;
	double TotalTime=0.15;
	int div=20;//for output
	double q[Node];
	double h[Node];
	double qn[Node];
	double hn[Node];
	double FL[2],FR[2];
	double InithL=1.0;
	double InithR=0.1;
	double z[Node];
	double I[Node];
	double mn=0.01;
	
	//IC//
	for(int i=0;i<Node;i++)
	{
		q[i]=0.0;
		h[i]=InithL;
	}
	for(int i=Node/2;i<Node;i++)
	{
		q[i]=0.0;
		h[i]=InithR;
	}
	
	//bottom//
	for(int i=0;i<Node;i++)
	{
		z[i]=0.0;
		double a=0.5-20.0*(float(i)*dx-0.5)*(float(i)*dx-0.5);
		if(a>0.0) z[i]=a;
		z[i]=0.0;
	}
	
	//slope//
	for(int i=1;i<Node-1;i++)
	{
		I[i]=-(z[i]-z[i-1])/dx;
	}
	//at boundary// 
	I[0]=-(z[1]-z[0])/dx;
	I[Node-1]=-(z[Node-1]-z[Node-2])/dx;
	
	//IC output//
	output(h,z,Node,dx);
	
	//time step// 
	for(int k=1;double(k)*dt<=TotalTime;k++)
	{	
		//inner cells//
		for(int i=1;i<Node-1;i++)
		{
			//left flux, between i-1 and i//
			RoeFlux(FL,h[i-1],h[i],q[i-1],q[i]);
			
			//right flux, between i and i+1//
			RoeFlux(FR,h[i],h[i+1],q[i],q[i+1]);
			
			hn[i]=h[i]-dt/dx*(FR[0]-FL[0]);
			qn[i]=q[i]-dt/dx*(FR[1]-FL[1]);//+dt*g*h[i]*(I[i]-mn*mn*q[i]*abs(q[i])/pow(h[i],10.0/3.0));
		}

		//BC//
		//left boundary//
		hn[0]=hn[1];
		qn[0]=qn[1];
		
		//right boundary//
		hn[Node-1]=hn[Node-2];
		qn[Node-1]=qn[Node-2];
		
		//update//
		for(int i=0;i<Node;i++)
		{
			q[i]=qn[i];
			h[i]=hn[i];
		}
		
		//output// 
		if(k%div==0)
		{
			output(h,z,Node,dx);
			cout<<"t="<<double(k)*dt<<endl;
		}
	}
	
	return 0;
}


こちらもどうぞ