求解如下方程 u = [ ρ ρ u E ] , F = [ ρ u ρ u 2 + p u ( E + p ) ] , E = p γ − 1 + 1 2 ρ u 2   ∂ u ∂ t + ∂ F ∂ x = 0 u=\left[\begin{matrix} \rho \\ \rho u\\ E \end{matrix}\right], F=\left[\begin{matrix} \rho u\\ \rho u^2+p\\ u(E+p) \end{matrix}\right], E=\frac{p}{\gamma-1}+\frac{1}{2}\rho u^2\\ \ \\ \frac{\partial u}{\partial t}+\frac{\partial F}{\partial x}=0 u=⎣⎡​ρρuE​⎦⎤​,F=⎣⎡​ρuρu2+pu(E+p)​⎦⎤​,E=γ−1p​+21​ρu2 ∂t∂u​+∂x∂F​=0 的间断初值问题称为激波管问题也叫黎曼问题。 该问题描述的是无粘理想气体在中一根管道中由压力或者速度驱动在突变的初始条件下的变化过程。 通过求解通量的雅各比矩阵 A = ∂ F ⃗ ∂ u ⃗ A=\frac{\partial {\vec F}}{\partial {\vec u}} A=∂u ∂F ​可以将前面的控制方程写成 ∂ u ∂ t + A ∂ u ∂ x = 0 \frac{\partial u}{\partial t}+A\frac{\partial u}{\partial x}=0 ∂t∂u​+A∂x∂u​=0 其中 u = [ ρ ρ u E ] , A = [ 0 1 0 − 3 − γ 2 u 2 ( 3 − γ ) u γ − 1 − u ( 1 γ − 1 c 2 + 2 − γ 2 u 2 ) 1 γ − 1 c 2 + 3 − 2 γ 2 u 2 γ u ] u=\left[\begin{matrix} \rho \\ \rho u\\ E \end{matrix}\right], A=\left[\begin{matrix} 0&1&0 \\ -\frac{3-\gamma}{2}u^2&(3-\gamma)u&\gamma-1\\ -u(\frac{1}{\gamma -1}c^2+\frac{2-\gamma}{2}u^2)&\frac{1}{\gamma -1}c^2+\frac{3-2\gamma}{2}u^2&\gamma u \end{matrix}\right] u=⎣⎡​ρρuE​⎦⎤​,A=⎣⎡​0−23−γ​u2−u(γ−11​c2+22−γ​u2)​1(3−γ)uγ−11​c2+23−2γ​u2​0γ−1γu​⎦⎤​ 其中 c = γ p ρ c=\sqrt{\gamma\frac{p}{\rho}} c=γρp​ ​是声速 对于该方程存在一个特殊关系 以守恒变量表示的守恒形式和非守恒形式的方程之间满足雅各比矩阵 A A A有齐次性 即 d F = A d u dF=Adu dF=Adu且 F = A u F=Au F=Au

任何与 A A A相似的矩阵 B B B以及相应的可逆变换 P P P满足 A = P − 1 B P A=P^{-1}BP A=P−1BP，都可以得到一个新的等价形式 ∂ ( P u ) ∂ t + B ∂ ( P u ) ∂ x = 0 \frac{\partial (Pu)}{\partial t}+B\frac{\partial (Pu)}{\partial x}=0 ∂t∂(Pu)​+B∂x∂(Pu)​=0 于是根据原变量 ρ , u , p \rho,u,p ρ,u,p可以得到相应的等价形式 ∂ u ∂ t + A ∂ u ∂ x = 0 u = [ ρ u p ] , A = [ u ρ 0 0 u 1 ρ 0 ρ c 2 u ] \frac{\partial u}{\partial t}+A\frac{\partial u}{\partial x}=0\\ u=\left[\begin{matrix} \rho \\ u\\ p \end{matrix}\right], A=\left[\begin{matrix} u&\rho&0 \\ 0&u&\frac{1}{\rho}\\ 0&\rho c^2&u \end{matrix}\right] ∂t∂u​+A∂x∂u​=0u=⎣⎡​ρup​⎦⎤​,A=⎣⎡​u00​ρuρc2​0ρ1​u​⎦⎤​ 通过原变量形式的方程可以很容易计算出矩阵 A A A的特征值，从而得到与 A A A相似的对角阵，将方程解耦，得到三个方程。分别求出三个方程的特征线并对应三个变量，沿三条特征线三个变量的值不发生改变，这三个变量称为黎曼不变量。

一根两端封闭的管道内左右两侧有不同的气体，中间以一物体将两侧气体分隔开，左侧气体满足 ρ L = 1 , u L = 0 , p L = 1 \rho_L=1,u_L=0,p_L=1 ρL​=1,uL​=0,pL​=1右侧气体满足 ρ L = 2 , u L = 0 , p L = 2 \rho_L=2,u_L=0,p_L=2 ρL​=2,uL​=0,pL​=2现在突然撤去中间分隔物，两侧气体将在压力作用下相互混合，最终达到平衡状态。 这里以二阶迎风的通量差分格式计算。 由原变量表达的方程中空间导数项的系数矩阵 A A A的特征值为 u + c , u − c , u u+c,u-c,u u+c,u−c,u分别表示三个黎曼不变量的传播速度，其中本问题中显然一定有亚音速的区域，因此必然有 u + c > 0 , u − c < 0 u+c>0,u-c0,u−c *f=F[0]; f++; *f=F[1]; f++; *f=F[2]; f++; } double max(double x1,double x2) { if(x1>x2) return x1; else return x2; } void advance(vector& w1,vector& w2,vector& w3,vector& F_p,vector& F_m) { vector tF(3,0); double l=0; vector::iterator f_p=F_p.begin(),f_m=F_m.begin(); for(int i=0;i F_p[3*i]+=l*w1[i]; F_m[3*i]-=l*w1[i]; F_p[3*i+1]+=l*w2[i]; F_m[3*i+1]-=l*w2[i]; F_p[3*i+2]+=l*w3[i]; F_m[3*i+2]-=l*w3[i]; } f_p=F_p.begin()+3,f_m=F_m.begin()+3; w1[1]=w1[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx; f_p++; f_m++; w2[1]=w2[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx; f_p++; f_m++; w3[1]=w3[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx; f_p++; f_m++; for(int i=2;i const vector &w1,&w2,&w3; val(const vector& _w1,const vector& _w2,const vector& _w3):w1(_w1),w2(_w2),w3(_w3){}; }; void init(vector &w1,vector &w2,vector &w3) { int i=1; for(;i w1[i]=1; w2[i]=0; w3[i]=1/(gamma-1); } w1[0]=w1[1]; w2[0]=-w2[1]; w3[0]=w3[1]; w1[w1.size()-1]=w1[w1.size()-2]; w2[w2.size()-1]=-w2[w2.size()-2]; w3[w3.size()-1]=w3[w3.size()-2]; } ostream& operator double rho=Q.w1[i],u=Q.w2[i]/Q.w1[i],p=(gamma-1)*(Q.w3[i]-Q.w2[i]*Q.w2[i]/Q.w1[i]/2); out advance(w1,w2,w3,F_p,F_m); if(i%SKIP_STEP==0)cout F(tF,w1[i],w2[i],w3[i]); double u=w2[i]/w1[i],p=(gamma-1)*(w3[i]-w2[i]*w2[i]/w1[i]/2),c=sqrt(gamma*p/w1[i]); l=max(max(abs(u+c),abs(u-c)),l); F_div(f_p,tF,w1[i],w2[i],w3[i]); F_div(f_m,tF,w1[i],w2[i],w3[i]); } for(int i=0;i w1[i]=w1[i]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx; f_p++; f_m++; w2[i]=w2[i]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx; f_p++; f_m++; w3[i]=w3[i]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx; f_p++; f_m++; } w1[w1.size()-2]=w1[w1.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx; f_p++; f_m++; w2[w2.size()-2]=w2[w2.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx; f_p++; f_m++; w3[w3.size()-2]=w3[w3.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx; f_p++; f_m++; w1[0]=w1[1]; w2[0]=-w2[1]; w3[0]=w3[1]; w1[w1.size()-1]=w1[w1.size()-2]; w2[w2.size()-1]=-w2[w2.size()-2]; w3[w3.size()-1]=w3[w3.size()-2]; }

结果如下 可以看到使用高耗散的低阶格式间断处的震荡就可以被抑制住。但是耗散非常明显。

这里使用限制器来抑制间断震荡，首先相邻确定差分的比值根据信息的传播方向分成两个方向的比值 r i + 1 2 − = u i + 2 − u i + 1 u i + 1 − u i r i + 1 2 + = u i − u i − 1 u i − 1 − u i − 2 r^-_{i+\frac{1}{2}}=\frac{u_{i+2}-u_{i+1}}{u_{i+1}-u_{i}}\\ r^+_{i+\frac{1}{2}}=\frac{u_{i}-u_{i-1}}{u_{i-1}-u_{i-2}} ri+21​−​=ui+1​−ui​ui+2​−ui+1​​ri+21​+​=ui−1​−ui−2​ui​−ui−1​​ 限制器函数 φ ( r ) = m i n ( 1 , ∣ r ∣ ) \varphi(r)=min(1,|r|) φ(r)=min(1,∣r∣) 格式如下 u i n + 1 = u i n − Δ t 2 Δ x ( φ ( r i + 1 2 + ) F i + 1 2 + ( 1 ) − φ ( r i − 1 2 + ) F i − 1 2 + ( 1 ) + ( 1 − φ ( r i + 1 2 − ) ) F i + 1 2 − ( 2 ) − ( 1 − φ ( r i + 1 2 − ) ) F i + 1 2 − ( 2 ) ) u_i^{n+1}=u_i^n-\frac{\Delta t}{2\Delta x} (\varphi(r^+_{i+\frac{1}{2}})F^{+(1)}_{i+\frac{1}{2}}-\varphi(r^+_{i-\frac{1}{2}})F^{+(1)}_{i-\frac{1}{2}}+ (1-\varphi(r^-_{i+\frac{1}{2}}))F^{-(2)}_{i+\frac{1}{2}}-(1-\varphi(r^-_{i+\frac{1}{2}}))F^{-(2)}_{i+\frac{1}{2}}) uin+1​=uin​−2ΔxΔt​(φ(ri+21​+​)Fi+21​+(1)​−φ(ri−21​+​)Fi−21​+(1)​+(1−φ(ri+21​−​))Fi+21​−(2)​−(1−φ(ri+21​−​))Fi+21​−(2)​) 这里做了一个简化，直接令 r i + = u i − u i − 1 u i − 1 − u i − 2 r i − = u i + 2 − u i + 1 u i + 1 − u i r^+_i=\frac{u_{i}-u_{i-1}}{u_{i-1}-u_{i-2}}\\ r^-_i=\frac{u_{i+2}-u_{i+1}}{u_{i+1}-u_{i}} ri+​=ui−1​−ui−2​ui​−ui−1​​ri−​=ui+1​−ui​ui+2​−ui+1​​ 用 φ ( r i ) \varphi(r_i) φ(ri​)来决定使用一阶还是二阶格式 即 ∂ u ∂ t + φ ( r + ) ∂ F + ( 2 ) ∂ x + φ ( r − ) ∂ F − ( 2 ) ∂ x + ( 1 − φ ( r + ) ) ∂ F + ( 1 ) ∂ x + ( 1 − φ ( r − ) ) ∂ F − ( 1 ) ∂ x = 0 \frac{\partial u}{\partial t}+\varphi(r^+)\frac{\partial F^{+(2)}}{\partial x}+\varphi(r^-)\frac{\partial F^{-(2)}}{\partial x}+(1-\varphi(r^+))\frac{\partial F^{+(1)}}{\partial x}+(1-\varphi(r^-))\frac{\partial F^{-(1)}}{\partial x}=0 ∂t∂u​+φ(r+)∂x∂F+(2)​+φ(r−)∂x∂F−(2)​+(1−φ(r+))∂x∂F+(1)​+(1−φ(r−))∂x∂F−(1)​=0 具体的advance函数如下