算例收集5：双马赫反射
1.问题描述
双马赫反射（Double Mach Reflection，DMR）是计算流体力学（CFD）和高超声速空气动力学中一个经典的验证算例，主要用于测试数值方法（特别是激波捕捉格式）处理强间断、复杂波系相互作用和多维效应的能力。它模拟了一个特定的激波反射现象，具有清晰的物理背景和独特的数值挑战性。
该问题模拟了一个强平面斜激波（通常马赫数很高，如 Mach=10）撞击刚性倾斜壁面时发生的现象。当入射激波与壁面的夹角（入射角）大于某个临界值时，会发生马赫反射（Mach Reflection），而不是规则反射（Regular Reflection）。在 DMR 算例的特定设置下，这种马赫反射过程会演化出一种对称的双马赫杆结构。
流场特点有：
1、强间断性： 包含极强的激波（高压比、高密度比、高温度比），对数值格式的激波捕捉能力、鲁棒性和非振荡性（无伪振荡） 要求极高。
2、多波系复杂相互作用： 两个强马赫杆（强间断）、入射激波（部分区域）反射激波（可能不止一条）滑移线（接触间断，两侧速度相等但密度、温度、熵不等）
3、典型结构清晰：成功的模拟应该清晰地分辨出：两个竖直（或接近竖直）的密集等值线/色块带（马赫杆）、连接两个马赫杆顶部的、向上游倾斜的反射激波、连接两个马赫杆底部（靠近壁面）的波（可能是反射激波的延伸或另一条波）、从两个三波点发出的滑移线（通常表现为密度或压力的间断线）、滑移线下方（靠近壁面）的复杂涡结构（由滑移线不稳定性及波系相互作用引起）
4、高数值挑战性：激波分辨率、接触间断分辨率、多尺度问题、网格收敛性（特别是滑移线的发展、涡结构的细节、三波点的位置对网格分辨率非常敏感）、边界条件处理、鲁棒性。
双马赫反射算例通过模拟一个强斜激波在倾斜壁面上诱发的复杂马赫反射及其特殊边界条件下形成的对称双马赫杆结构，为 CFD 数值方法（尤其是欧拉方程求解器）提供了一个极其重要的基准测试。
1.1控制方程
控制方程采用二维Euler方程：
\(\frac{ \partial q}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y}= 0\)
其中：
\(q=\begin{pmatrix}
\rho \\\rho u \\\rho v
\\\rho e
\end{pmatrix},F=\begin{pmatrix}
\rho u \\\rho u^2+p \\\rho uv
\\u(\rho e+p)
\end{pmatrix},G=\begin{pmatrix}
\rho v \\\rho uv \\\rho v^2+p
\\v(\rho e+p)
\end{pmatrix}\)
1.2 初始条件及边界条件
双马赫反射问题的初场如图所示：
平板位于x=1/6, y=0处；t=0时，M=10的右行正激波置于x=1/6, y=0处,与平板之间成60°夹角。该算例中气体比热比为1.4。
算例中波前和波后参数分别取为：
1、左边界处设为波后参数；
2、对于下边界，0<x<1/6设为波后参数，1/6<x<4处为壁面，即反射边界条件,具体来说，这里还是把\(\rho,u,\rho e\)延拓为关于y的偶函数，把\(v\)延拓为关于y的奇函数；
3、右边界为出口条件；
4、上边界为准确描述马赫数为10的激波运动的条件，即以激波速度移动。
当\(x<1/6+(1+20t_{RK})/\sqrt{3}\)时，上边界设定为波后值，反之则设定为波前值。这里\(t_{RK}\)表示3阶Runge-Kutta格式中每个Stage里的时间，具体写出为：
\(\begin{matrix}
stage1:t_{RK}=t^n\\stage2:t_{RK}=t^n+dt
\\stage3:t_{RK}=t^n+dt/2
\end{matrix}\)
计算域为[0,4]x[0,1],计算至\(t_{end}==0.2s\)。
初始条件：
gama=1.4
nx=804
ny=661
tf=0.2
dt=1e-4
num_iter=round(Int,tf/dt)
xb=0
xe=3.2
yb=0
ye=1.0
dx=(xe-xb)/(nx-1)
dy=(ye-yb)/(ny-1)
pr = 116.5
rho = 8.0
u = 8.25*cos(pi/6)
v = -8.25*sin(pi/6)
q1=turn_var(gama,rho,u,v,pr)
pr = 1.0
rho = 1.4
u = 0.0
v = 0.0
q2=turn_var(gama,rho,u,v,pr)
x=Array{Float64}(undef,nx)
y=Array{Float64}(undef,ny)
qn=Array{Float64}(undef,4,ny,nx)
for i=1:nx
x[i]=(i-1)*dx
end
for j=1:ny
y[j]=(j-1)*dy
end
for i=1:nx
for j=1:ny
if(x[i] < 1/6.0+y[j]/sqrt(3.0))
qn[:,j,i] = q1
else
qn[:,j,i] = q2
end
end
end
2.数值方法
采用格点格式有限差分法。
2.1 插值格式
采用了WENO格式。
WENO（Weighted Essentially Non-Oscillatory，加权本质无振荡）格式是一种高精度、高分辨率且能稳定捕捉强间断的数值方法，广泛应用于计算流体力学（CFD）中求解含有激波、接触间断等复杂结构的偏微分方程（如欧拉方程）。
WENO格式继承自 ENO（Essentially Non-Oscillatory）格式的思想：在间断附近自适应选择最光滑的子模板重构数值通量，避免跨过间断时产生非物理振荡。计算时将所有候选模板的插值结果通过非线性权重融合，形成最终重构，其权重设计要求为：对光滑区域的模板赋予大权重（接近理想线性组合）；对包含间断的模板赋予接近零的权重，抑制振荡；权重由模板的光滑度量因子（Smoothness Indicator）动态计算，实现自适应。
本算例采用5阶精度的WENO-JS格式进行插值，计算单元的左右界面的守恒变量，具体格式如下：
\(\hat{q}_{i+1/2}=\sum_{k=0}^{2} \omega _{k} q_{k}\)
其中qk是包含三个格点的子模板\(S_{k}^3=(i+k-2,i+k-1,i+k)\) , qk由下式给出：
\(\left\{\begin{matrix}
q_0=\frac{1}{3}q_{i-2}-\frac{7}{6}q_{i-1}+\frac{11}{6}q_i\\q_1=-\frac{1}{6}q_{i-1}+\frac{5}{6}q_{i}+\frac{1}{3}q_{i+1}
\\q_2=\frac{1}{3}q_{i}+\frac{5}{6}q_{i+1}-\frac{1}{6}q_{i+2}
\end{matrix}\right.\)
权重的计算公式为：
\(\omega _{k}=\frac{\alpha _k}{\alpha _0+ \alpha _1 + \alpha _2},\alpha _k=\frac{c_k}{(IS_k+\epsilon )^2},k=0,1,2\)
其中\(IS_k\)就是光滑度量因子（Smoothness Indicator），c0=0.1,c1=0.6，c2=0.3分别为模板最佳权重， \(\epsilon=10^{-6}\)引入一个很小的正数使分母不为零。
对于5阶格式，\(IS_k\)由下式给出：
\(\left\{\begin{matrix}
IS_0=\frac{13}{12}(q_{i-2}-2q_{i-1}+q_{i})^2+\frac{1}{4}(q_{i-2}-4q_{i-1}+3q_{i})^2 \\IS_1=\frac{13}{12}(q_{i_1}-2q_{i}+q_{i+1})^2+\frac{1}{4}(q_{i_1}-q_{i+1})^2
\\IS_2=\frac{13}{12}(q_{i}-2q_{i+1}+q_{i+2})^2+\frac{1}{4}(3q_{i}-4q_{i+1}+q_{i+2})^2
\end{matrix}\right.\)
插值示意图如下：
以正向为例（L方向）计算过程如下：
function ML(q1,q2,q3,q4,q5)
eps=1.0e-6
s1=(13.0/12.0)*(q1-2.0*q2+q3)^2 + 0.25*(q1-4.0*q2+3.0*q3)^2
s2=(13.0/12.0)*(q2-2.0*q3+q4)^2 + 0.25*(q2-q4)^2
s3=(13.0/12.0)*(q3-2.0*q4+q5)^2 + 0.25*(3.0*q3-4.0*q4+q5)^2
a1=0.1/(s1+eps)^2
a2=0.6/(s2+eps)^2
a3=0.3/(s3+eps)^2
sum=a1+a2+a3
w1=a1/sum
w2=a2/sum
w3=a3/sum
f1=1.0/3.0*q1 -7.0/6.0*q2 +11.0/6.0*q3
f2=-1.0/6.0*q2 +5.0/6.0*q3 + 1.0/3.0*q4
f3=1.0/3.0*q3 + 5.0/6.0*q4 -1.0/6.0*q5
qL=w1*f1+w2*f2+w3*f3
return qL
end
2.2 重构格式
采用AUSMPW+格式，具体可以见算例收集4。
2.3 时间推进
采用3步Runge-Kutta格式。
for iter=1:num_iter
global tn
t_Rk = tn
rhs1=rhs(nx,ny,dx,dy,qn,q1,q2,t_Rk,x)
for k=1:4
for j=1:ny
for i=1:nx
qt[k,j,i]=qn[k,j,i]+dt*rhs1[k,j,i]
end
end
end
t_Rk=tn+dt
rhs2=rhs(nx,ny,dx,dy,qt,q1,q2,t_Rk,x)
for k=1:4
for j=1:ny
for i=1:nx
qt[k,j,i]=0.75*qn[k,j,i]+0.25*qt[k,j,i]+0.25*dt*rhs2[k,j,i]
end
end
end
t_Rk=tn+dt/2.0
rhs3=rhs(nx,ny,dx,dy,qt,q1,q2,t_Rk,x)
for k=1:4
for j=1:ny
for i=1:nx
qt[k,j,i]=(1.0/3.0)*qn[k,j,i]+(2.0/3.0)*qt[k,j,i]+(2.0/3.0)*dt*rhs3[k,j,i]
end
end
end
tn = tn+dt
for k=1:4
for j=1:ny
for i=1:nx
qn[k,j,i]=qt[k,j,i]
end
end
end
end
3.计算结果
计算时间等于0.2s时的流场密度分布云图如下：
参考文献
J. Shi et al.Resolution of high order WENO schemes for complicated flow structures[J].Journal of Computational Physics 186 (2003) 690–696.
Y. Shen, S. Li, S. Liu et al.A robust common-weights WENO scheme based on the flux vector splitting for Euler equations[J].Communications in Nonlinear Science and Numerical Simulation 119 (2023) 107112
保姆级间断Galerkin方法(10)：二维Euler方程
浙公网安备 33010602011771号