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

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

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

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


今回は、HLL法(Harten, Lax, van Leer法)を用いてダム崩壊問題(dam break problem)を解いてみました。HLL法とは、波速が最大のものと最小のものだけを考え、3つの状態(左側、右側、star region)を仮定して解く手法です。解いている方程式は、摩擦、勾配なしの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^{-7} で、水面は静止しているものとします。このように、dry bedの場合は非常に小さい水深を与えることとします。時刻  t=0 に真ん中にあるしきり(ダム)を取り去った時に、水がどのように動くかを計算します。境界条件としては右端、左端で水深、流量ともに斉次のノイマン条件を課しています。


参考文献としてはToroの"Shock-Capturing Methods for Free-Surface Shallow Flows"のpp.179-181です。バイブルですね。波速を推定するための手法としては、Exact Depth Positivity Riemann Solver (EDPRS)、Two-Rarefaction Riemann Solver (TRRS)、Two-Shock Riemann Solver (TSRS) を選べるようになっています。

Shock-Capturing Methods for Free-Surface Shallow Flows

Shock-Capturing Methods for Free-Surface Shallow Flows


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

f:id:mutsumunemitsutan:20190527211046g:plain:w500

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

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

f:id:mutsumunemitsutan:20190530181953g:plain:w500

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


一応Youtubeのチャンネルも作ってみました!まあ、上のGif動画と同じものをアップロードしただけですが…


では計算コードです。

#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);
}

//function for TRRS//
inline double gk(double h0, double hlr)
{
	return sqrt(g/2.0*(h0+hlr)/(h0*hlr));
}

//Riemann Solver Based on Exact Depth Positivity//
inline void EDPRS(double &hst, double &ust, double hL, double hR, double qL, double qR)
{
	double aL=a(hL);
	double aR=a(hR);
	double uL=qL/hL;
	double uR=qR/hR;
	
	hst=(hL+hR)/2.0-(uR-uL)*(hL+hR)/4.0/(aL+aR);
	ust=(uL+uR)/2.0-(hR-hL)*(aL+aR)/(hL+hR);
}

//Two-Rarefaction Riemann Solver//
inline void TRRS(double &hst, double &ust, double hL, double hR, double qL, double qR)
{
	double aL=a(hL);
	double aR=a(hR);
	double uL=qL/hL;
	double uR=qR/hR;
	
	hst=1.0/g*((aL+aR)/2.0+(uL-uR)/4.0)*((aL+aR)/2.0+(uL-uR)/4.0);
	ust=(uL+uR)/2.0+aL-aR;
}

//Two-Shock Riemann Solver//
inline void TSRS(double &hst, double &ust, double hL, double hR, double qL, double qR)
{
	double uL=qL/hL;
	double uR=qR/hR;
	
	TRRS(hst,ust,hL,hR,qL,qR);
	//EDPRS(hst,ust,hL,hR,qL,qR);
	double hcon=hst;
	hst=(gk(hst,hL)*hL+gk(hst,hR)*hR+uL-uR)/(gk(hst,hL)+gk(hst,hR));
	ust=(uL+uR)/2.0+((hcon-hR)*gk(hcon,hR)-(hcon-hL)*gk(hcon,hL))/2.0;
}

//function for HLL//
inline double qk(double hst, double hlr)
{
	double qk=1.0;
	if(hst>hlr) qk=sqrt((hst+hlr)*hst/(hlr*hlr)/2.0);
	return qk;
}


//left wave estimation//
inline double SLeft(double hst, double hL, double qL)
{
	double uL=qL/hL;
	return uL-a(hL)*qk(hst,hL);
}

//right wave estimation//
inline double SRight(double hst, double hR, double qR)
{
	double uR=qR/hR;
	return uR+a(hR)*qk(hst,hR);
}

//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;

}

//hll flux//
inline void HLLFlux(double FLR[], double hL, double hR, double qL, double qR)
{
	double hst,ust;

	//estimastion of hst and ust//
	TRRS(hst,ust,hL,hR,qL,qR);
	//TSRS(hst,ust,hL,hR,qL,qR);
	//EDPRS(hst,ust,hL,hR,qL,qR);
	
	//wave speed estimation//
	double SL=SLeft(hst,hL,qL);
	double SR=SRight(hst,hR,qR);
	
	if(SL>=0)
	{
		FLR[0]=Fluxcon(hL,qL);
		FLR[1]=Fluxmom(hL,qL);
	}
	else if(SR<=0)
	{
		FLR[0]=Fluxcon(hR,qR);
		FLR[1]=Fluxmom(hR,qR);
	}
	else
	{
		FLR[0]=(SR*Fluxcon(hL,qL)-SL*Fluxcon(hR,qR)+SR*SL*(hR-hL))/(SR-SL);
		FLR[1]=(SR*Fluxmom(hL,qL)-SL*Fluxmom(hR,qR)+SR*SL*(qR-qL))/(SR-SL);
	}
}

inline void output(double h[], 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]<<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=10;//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;
	
	//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;
	}
	
	//IC output//
	output(h,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//
			HLLFlux(FL,h[i-1],h[i],q[i-1],q[i]);
			
			//right flux, between i and i+1//
			HLLFlux(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]);
		}

		//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,Node,dx);
			cout<<"t="<<double(k)*dt<<endl;
		}
	}
	
	return 0;
}


こちらもどうぞ