计算流体力学简介(八)——激波管问题
问题描述
求解如下方程
u=[ρρuE],F=[ρuρu2+pu(E+p)],E=pγ−1+12ρu2∂u∂t+∂F∂x=0u=\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}=0u=⎣⎡ρρ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=[ρρuE],A=[010−3−γ2u2(3−γ)uγ−1−u(1γ−1c2+2−γ2u2)1γ−1c2+3−2γ2u2γ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(γ−11c2+22−γu2)1(3−γ)uγ−11c2+23−2γu20γ−1γu⎦⎤
其中c=γpρc=\sqrt{\gamma\frac{p}{\rho}}c=γρp是声速
对于该方程存在一个特殊关系
以守恒变量表示的守恒形式和非守恒形式的方程之间满足雅各比矩阵AAA有齐次性
即dF=AdudF=AdudF=Adu且F=AuF=AuF=Au
任何与AAA相似的矩阵BBB以及相应的可逆变换PPP满足A=P−1BPA=P^{-1}BPA=P−1BP，都可以得到一个新的等价形式
∂(Pu)∂t+B∂(Pu)∂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=0u=[ρup],A=[uρ00u1ρ0ρc2u]\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ρc20ρ1u⎦⎤
通过原变量形式的方程可以很容易计算出矩阵AAA的特征值，从而得到与AAA相似的对角阵，将方程解耦，得到三个方程。分别求出三个方程的特征线并对应三个变量，沿三条特征线三个变量的值不发生改变，这三个变量称为黎曼不变量。
问题求解
一根两端封闭的管道内左右两侧有不同的气体，中间以一物体将两侧气体分隔开，左侧气体满足ρL=1,uL=0,pL=1\rho_L=1,u_L=0,p_L=1ρL=1,uL=0,pL=1右侧气体满足ρL=2,uL=0,pL=2\rho_L=2,u_L=0,p_L=2ρL=2,uL=0,pL=2现在突然撤去中间分隔物，两侧气体将在压力作用下相互混合，最终达到平衡状态。
这里以二阶迎风的通量差分格式计算。
由原变量表达的方程中空间导数项的系数矩阵AAA的特征值为u+c,u−c,uu+c,u-c,uu+c,u−c,u分别表示三个黎曼不变量的传播速度，其中本问题中显然一定有亚音速的区域，因此必然有u+c>0,u−c<0u+c>0,u-c<0u+c>0,u−c<0说明方程本身既有向左传播的信息又有向右传播的信息。这时就需要根据信息传播方向决定使用哪一侧的迎风格式。
通量分裂的方法有很多，这里使用Lax-Friedrichs分裂
前面提到了雅各比矩阵的特征值是u,u+c,u−cu,u+c,u-cu,u+c,u−c那么根据当地的ρ,u,p\rho,u,pρ,u,p可以很容易计算出三个特征值的值，找到三者中模最大的特征值∣λ∣max|\lambda|_{max}∣λ∣max然后构造新的雅各比矩阵A+=A+∣λ∣maxI2,A−=A−∣λ∣maxI2A^+=\frac{A+|\lambda|_{max}I}{2},A^-=\frac{A-|\lambda|_{max}I}{2}A+=2A+∣λ∣maxI,A−=2A−∣λ∣maxI这样必然有A=A++A−A=A^++A^-A=A++A−且A+A^+A+的全部特征值大于0，A−A^-A−的全部特征值小于0，利用守恒型变量雅各比矩阵的齐次性有
F+=A+u=F+∣λ∣maxu2,F−=A−u=F−∣λ∣maxu2F^+=A^+u=\frac{F+|\lambda|_{max}u}{2},F^-=A^-u=\frac{F-|\lambda|_{max}u}{2}F+=A+u=2F+∣λ∣maxu,F−=A−u=2F−∣λ∣maxu
于是得到分裂的守恒方程
∂u∂t+∂F+∂x+∂F−∂x=0\frac{\partial u}{\partial t}+\frac{\partial F^+}{\partial x}+\frac{\partial F^-}{\partial x}=0∂t∂u+∂x∂F++∂x∂F−=0
式中第二项表示向右传播的信息，第三项表示向左传播的信息，分别使用3I−4E−1+E−22Δx\frac{3I-4E^{-1}+E^{-2}}{2\Delta x}2Δx3I−4E−1+E−2和−E2−4E1+3I2Δx-\frac{E^2-4E^1+3I}{2\Delta x}−2ΔxE2−4E1+3I进行离散。
时间使用单步推进。
Δt=0.005,Δx=0.1\Delta t=0.005,\Delta x=0.1Δt=0.005,Δx=0.1
具体代码如下
计算结果如下
蓝线是压力，黑线是密度，红线是速度
由于我也没有激波管解析解的结果。所以也没办法判断具体结果差多少，不过从网上看到的一些激波管问题的解来看基本规律应该没有太大的问题，不过由于二阶迎风格式的耗散较小，因此间断处有一些震荡。
要想消除震荡可以利用限制器在间断附近使用耗散更大的低阶通量或者使用ENO和WENO格式来抑制震荡。
现在先测试一下一阶格式，使用如下advance函数
void advance(vector<double>& w1,vector<double>& w2,vector<double>& w3,vector<double>& F_p,vector<double>& F_m) {vector<double> tF(3,0);double l=0;vector<double>::iterator f_p=F_p.begin(),f_m=F_m.begin();for(int i=0;i<w1.size();i++){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.size();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<w1.size()-2;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]; }结果如下
可以看到使用高耗散的低阶格式间断处的震荡就可以被抑制住。但是耗散非常明显。
这里使用限制器来抑制间断震荡，首先相邻确定差分的比值根据信息的传播方向分成两个方向的比值
ri+12−=ui+2−ui+1ui+1−uiri+12+=ui−ui−1ui−1−ui−2r^-_{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−uiui+2−ui+1ri+21+=ui−1−ui−2ui−ui−1
限制器函数φ(r)=min(1,∣r∣)\varphi(r)=min(1,|r|)φ(r)=min(1,∣r∣)
格式如下
uin+1=uin−Δt2Δx(φ(ri+12+)Fi+12+(1)−φ(ri−12+)Fi−12+(1)+(1−φ(ri+12−))Fi+12−(2)−(1−φ(ri+12−))Fi+12−(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))
这里做了一个简化，直接令
ri+=ui−ui−1ui−1−ui−2ri−=ui+2−ui+1ui+1−uir^+_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−2ui−ui−1ri−=ui+1−uiui+2−ui+1
用φ(ri)\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函数如下
Δt=0.001,Δx=0.1\Delta t=0.001,\Delta x=0.1Δt=0.001,Δx=0.1
计算结果如下
对比前面一阶格式和二阶格式可以看到，使用限制器之后二阶格式的间断震荡被明显抑制住了，而耗散也没有一阶格式那么明显。
总结
以上是生活随笔为你收集整理的计算流体力学简介(八)——激波管问题的全部内容，希望文章能够帮你解决所遇到的问题。
- 上一篇： 计算机游戏与动漫设计大赛,西安综合职业中
- 下一篇： 2021-09-27 网安实验-文件恢复