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

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

今回は、HLLC法（Harten, Lax, van Leer, Contact surface法）を用いてダム崩壊問題（dam break problem）を解いてみました。HLLC法とは、HLL法を発展させたもので、波速が最大のものと最小のものに加えて、contact discontinuity（接触不連続、shear wave）を考慮し、全体を4つの状態に分けて解く手法です。Toro教授が考案しました。2×2の普通の1次元浅水流方程式を解く上では、contact discontinuityはあらわれないのですが、1次元浅水流方程式と溶質の濃度を考える場合、1次元浅水流方程式と移動床を考える場合、2次元浅水流方程式などの3×3のシステムを解く際には、contact discontinuityがあらわれます。HLLC法を用いると、このcontact discontinuityを捉えることができるようになります。解いている方程式は、摩擦、勾配なしの1次元浅水流方程式（Shallow Water Equations）に溶質の濃度を考慮した3×3の以下のシステムです。

$\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 \\ &\frac{\partial (h \phi) }{\partial t} + \frac{\partial (uh \phi)}{\partial x} = 0 \end{eqnarray}$

ここで $\phi$ はpassiveな溶質の濃度です。参考文献はLeVequeの"Finite Volume Methods for Hyperbolic Problems"のpp.283-285です。

Finite Volume Methods for Hyperbolic Problems (Cambridge Texts in Applied Mathematics)

作者: Randall J. Leveque
出版社/メーカー: Cambridge University Press
発売日: 2002/08/29
メディア: ペーパーバック
この商品を含むブログを見る

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

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

Shock-Capturing Methods for Free-Surface Shallow Flows

作者: Eleuterio F. Toro
出版社/メーカー: Wiley
発売日: 2001/03/23
メディア: ハードカバー
この商品を含むブログを見る

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

f:id:mutsumunemitsutan:20190630231223g:plain:w500

右へと進行する衝撃波（不連続な波、Shock wave）と左へと進行する膨張波（連続な波、Rarefaction wave）が観察できます。よい感じです。これはHLLとほぼ同じ結果です。

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

f:id:mutsumunemitsutan:20190706152517g:plain:w500

初期時刻にステップ上に分布していた溶質の濃度が、ダムが崩壊すると同時に右方向へと流れていきます。厳密解はステップの形状を保つものです。

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

f:id:mutsumunemitsutan:20190706153010g:plain:w500

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

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

f:id:mutsumunemitsutan:20190706152905g:plain:w500

では計算コードです。

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

//shear wave estimation//
inline double SShear(double hst, double ust, double hL, double hR, double qL, double qR, double SL, double SR)
{
	double uR=qR/hR;
	double uL=qL/hL;
	return (SL*hR*(uR-SR)-SR*hL*(uL-SL))/(hR*(uR-SR)-hL*(uL-SL));
}

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

}

//flux for shear wave//
inline double Fluxshear(double h, double q, double p)
{
	return p*q/h;
}

//hllc flux//
inline void HLLCFlux(double FLR[], double hL, double hR, double qL, double qR, double pL, double pR)
{
	double hst,ust;
	double uR=qR/hR;
	double uL=qL/hL;

	//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);
	
	//estimation of shear wave//
	double Sst=SShear(hst,ust,hL,hR,qL,qR,SL,SR);
	//double Sst=ust;
	
	if(SL>=0)
	{
		FLR[0]=Fluxcon(hL,qL);
		FLR[1]=Fluxmom(hL,qL);
		FLR[2]=Fluxshear(hL,qL,pL);
	}
	else if((SL<=0)&&(0<=Sst))
	{
		double UstL[3];
		UstL[0]=hL*(SL-uL)/(SL-Sst)*1.0;
		UstL[1]=hL*(SL-uL)/(SL-Sst)*Sst;
		UstL[2]=hL*(SL-uL)/(SL-Sst)*pL/hL;
		
		FLR[0]=Fluxcon(hL,qL)+SL*(UstL[0]-hL);
		FLR[1]=Fluxmom(hL,qL)+SL*(UstL[1]-qL);
		FLR[2]=Fluxshear(hL,qL,pL)+SL*(UstL[2]-pL);
	}
	else if((Sst<=0)&&(0<=SR))
	{
		double UstR[3];
		UstR[0]=hR*(SR-uR)/(SR-Sst)*1.0;
		UstR[1]=hR*(SR-uR)/(SR-Sst)*Sst;
		UstR[2]=hR*(SR-uR)/(SR-Sst)*pR/hR;
		
		FLR[0]=Fluxcon(hR,qR)+SR*(UstR[0]-hR);
		FLR[1]=Fluxmom(hR,qR)+SR*(UstR[1]-qR);
		FLR[2]=Fluxshear(hR,qR,pR)+SR*(UstR[2]-pR);
	}
	else
	{
		FLR[0]=Fluxcon(hR,qR);
		FLR[1]=Fluxmom(hR,qR);
		FLR[2]=Fluxshear(hR,qR,pR);
	}
}

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

int main()
{
	const int n=1000;//number of cells
	const int Node=n+1;
	double L=1.0;
	double dx=L/double(n);
	double dt=0.0001;
	double TotalTime=0.2;
	int div=20;//for output
	double q[Node];
	double h[Node];
	double p[Node];//solutants
	double qn[Node];
	double hn[Node];
	double pn[Node];
	double FL[3],FR[3];
	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;
		p[i]=InithL;
	}
	for(int i=Node/2;i<Node;i++)
	{
		q[i]=0.0;
		h[i]=InithR;
		p[i]=0.0;
	}
	
	//bottom//
	for(int i=0;i<Node;i++)
	{
		z[i]=0.0;
		double a=0.1-3.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,p,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//
			HLLCFlux(FL,h[i-1],h[i],q[i-1],q[i],p[i-1],p[i]);
			
			//right flux, between i and i+1//
			HLLCFlux(FR,h[i],h[i+1],q[i],q[i+1],p[i],p[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));
			pn[i]=p[i]-dt/dx*(FR[2]-FL[2]);
		}

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