算例收集4：二维高马赫数喷流
1.问题描述
高马赫数喷流问题是一个在流体力学、航空航天工程、声学和推进系统设计中至关重要的复杂物理现象。它指的是高速（马赫数显著大于1，通常指马赫数3以上）的气体从喷管（如火箭发动机、超燃冲压发动机、喷气发动机加力燃烧室或实验设施）排出进入静止或低速环境（通常是大气）时，所发生的一系列高度非线性、非稳态的物理过程。
下图是拉瓦尔喷管形成激波的过程
1、收缩段（亚声速加速）
2、喉部（声速临界点）在最小截面（喉部），气流达到马赫数 Ma=1（声速）
3、扩张段（超声速加速）：进入扩张段后，截面积增大，气流继续膨胀加速，马赫数提升至设计值（如 Ma>3）。
以下是该问题的主要特征和关键方面：
1）高马赫数： 喷流出口速度远高于当地声速（Ma >> 1）。这导致喷流内部和与环境气体相互作用时，惯性力远大于粘性力（高雷诺数），并且压缩性效应（气体密度变化）极为显著。
2）强欠膨胀： 在火箭和高超音速推进系统中，为了获得高推力，燃烧室压力通常远高于环境大气压。当高压气体通过喷管膨胀加速到超声速时，如果喷管出口处的静压（p_e）远高于环境背压（p_b），即 p_e > p_b，则喷流处于“欠膨胀”状态。高马赫数喷流通常伴随着极强的欠膨胀。
高欠膨胀超声速喷流最显著的特征之一。在喷流中心线上，膨胀波汇聚导致一道近乎垂直的正激波，称为马赫盘（或终端激波）。马赫盘下游，超声速流突然转变为亚声速流，压力和密度骤升，温度也急剧升高（气动加热）。
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 初始条件及边界条件
计算域大小为 \(\Omega =[0,1]\times [=0.25,0.25]\), 内部填充 \(（\rho,u,v,p）=(0.5,0,0,0.4172)\)的气体，边界条件为：上边、右边、下边为出流边界条件，左边 $y\in [-0.05,0.05] $的部分设置为 \(（\rho,u,v,p）=(5,800,0,0.4172)\)的喷流，其余部分与内部气体设置一致。另外，该算例中气体比热比为5/3。
这是一个很明显的高马赫数喷流+强欠膨胀过程问题，喷流马赫数很高，$ Ma=\frac{\left | U\right |}{\sqrt{\gamma p/\rho }}\approx 2156.92$。
计算过程中容易出现负密度和负压强，可以用来检验程序的保正性。
初始条件：
gama=5/3;
nx=320;
ny=180;
xb=0;
xe=1.0;
yb=-0.25;
ye=0.25;
dx=(xe-xb)/nx;
dy=(ye-yb)/ny;
dt=5e-7;
tm=0.001;
nt=tm/dt;
rho1=0.5;
u1=0.0;
v1=0.0;
p1=0.4127;
rho2=5;
u2=800;
v2=0.0;
p2=0.4127;
2.数值方法
总体采用格点格式有限差分法。
2.1 插值格式
采用了MUSCL格式，该格式是由Bram van Leer于1979年提出的高阶数值格式，主要用于求解双曲型守恒律方程（如欧拉方程）。其核心思想是通过高阶插值重构和限制器技术，在保持数值解单调性的同时提高空间离散精度。
该格式的特性是能够在激波和光滑区域间自适应切换，计算效率高于高阶WENO/TENO，但是达到的空间精度一般在2阶左右。
插值示意图：
MUSCL格式表达式为：
上式中W即为控制方程中的守恒变量q
插值精度及两侧的信息取舍选择取决于κ的值
采用MUSCL格式的插值方法在流场光滑区可以得到很好的结果，但当流场中存在大的梯度或间断时，这种方法就会引起振荡，出现上冲或下冲，即插值过大或过小。为了避免这一不利现象，引入了各种限制器来限制插值时的梯度。
将\(\Delta^- , \Delta^+\)替换为限制器函数，这里使用Minmod限制器，其特性为最保守，强烈抑制振荡。
\(\Delta^+ = minmod(r_L,1)\)，\(\Delta^- = minmod(r_R,1)\)
以正向（从左流向右 或 从下流向上）为例：
function qL=ML(q1,q2,q3)
qL=zeros(4,1);
K=1/3;
dp=delta(q3(:,1,1),q2(:,1,1));
dm=delta(q2(:,1,1),q1(:,1,1));
rl=zeros(4,1);
rl1=zeros(4,1);
for i=1:4
rl(i)=dp(i)/(dm(i)+1e-8);
rl1(i)=dm(i)/(dp(i)+1e-8);
end
v1=limiter(rl);
v2=limiter(rl1);
for i=1:4
qL(i)=q2(i,1,1) + 1/4 * ((1-K)*v1(i,1) + (1+K)*rl(i,1)*v2(i,1)) * dm(i,1);
end
end
function v=limiter(R) %Minmod limiter
v=zeros(4,1);
for i=1:4
v(i)=(sign(R(i))+1.0)/2.0*min(abs(R(i)),1.0);
end
end
2.2 重构格式
采用AUSMPW+格式，AUSM格式在理论上将流动对流特征中的线性场（与特征速度u有关）和非线性场（与特征速度u±a有关）相区别，并且将压力项与对流通量项分别开来分裂。从格式构造来讲，AUSM格式是Van Leer格式的一种发展改进，但从其耗散项来分析，它是一种FVS（矢量通量分裂格式）与FDS（通量差分格式）的复合格式。
这种格式构造简单，无矩阵运算，激波分辨率高而且稳定性好，既具有FDS格式在边界层中解的精确性，也具有FVS格式在捕捉强间断时的健壮性，具有非常优秀的性能。经过十多年的发展，已经应用于从低速到高超声速的各种气体模型的流动计算研究，形成了一系列在CFD研究中应用广泛的AUSM类格式。
自AUSM格式提出以来，有了两个分支，其一是Liou等通过修改马赫数分裂函数和压力项分裂项函数，发展了AUSMDV、AUSM+；另外一支是Kim等引入压力权函数修正等技术，发展了AUSMPW、AUSMPW+格式。
AUSMPW+格式界面上的通量重构表达式为：
\(F_{1/2,AUSM+}=c_{1/2}(M_{L}^{+}|_{\beta =1/8}\Phi _L+M_{R}^{-}|_{\beta =1/8}\Phi _R)+(P_{L}^{+}|_{\alpha =3/16} P_L+P_{R}^{-}|_{\alpha =3/16} P_R)\)
\(\Phi^x_{L/R}=\Phi ^y_{L/R}=\begin{pmatrix}
\rho \\\rho u \\\rho v
\\\rho e
\end{pmatrix},P^x_{L/R}=\begin{pmatrix}
0 \\p \\0
\\0
\end{pmatrix},P^y_{L/R}=\begin{pmatrix}
0 \\0 \\p
\\0
\end{pmatrix}\)
马赫数和压力项分裂函数分别为：
\(\tilde{c}=(c^*)^2 / max(\left | U\right |,c^*)\)
$ c^*=\sqrt{\frac{2(\gamma -1)}{\gamma +1}H} , H=\rho h=\rho (e+p/\rho) $
以x方向为例：
function residual_x=AUSM_plus_x(nx,ny,gama,qL_x,qR_x,dx)
Fx=zeros(4,nx+1);
residual_x=zeros(4,nx,ny);
vector(1)=1.0;
vector(2)=0.0;
for j=1:ny
for num_faces=1:nx+1
rhoL=qL_x(1,num_faces,j);
uL=qL_x(2,num_faces,j)/qL_x(1,num_faces,j);
vL=qL_x(3,num_faces,j)/qL_x(1,num_faces,j);
unL=uL*vector(1)+vL*vector(2); %velocity component
eL=qL_x(4,num_faces,j)/qL_x(1,num_faces,j);
pL=(gama-1)*rhoL*(eL-0.5*(uL^2+vL^2));
aL=sqrt(2*(gama-1)/(gama+1)*(eL+pL/rhoL));
rhoR=qR_x(1,num_faces,j);
uR=qR_x(2,num_faces,j)/qR_x(1,num_faces,j);
vR=qR_x(3,num_faces,j)/qR_x(1,num_faces,j);
unR=uR*vector(1)+vR*vector(2); %velocity component
eR=qR_x(4,num_faces,j)/qR_x(1,num_faces,j);
pR=(gama-1)*rhoR*(eR-0.5*(uR^2+vR^2));
aR=sqrt(2*(gama-1)/(gama+1)*(eR+pR/rhoR));
aL_ave=(aL^2)/max(abs(unL),aL);
aR_ave=(aR^2)/max(abs(unR),aR);
a_mid=min(aL_ave,aR_ave);
MaL=unL/a_mid;
MaR=unR/a_mid;
[phiL,PL_V]=AUSM_turnvar_x(gama,qL_x(:,num_faces,j));
[phiR,PR_V]=AUSM_turnvar_x(gama,qR_x(:,num_faces,j));
% Mach number and pressure splitting functions
M0L=ausmsplit_Ma(MaL);
M0R=ausmsplit_Ma(MaR);
P0L=ausmsplit_P(MaL);
P0R=ausmsplit_P(MaR);
for k=1:4
Fx(k,num_faces)=a_mid*(M0L(1)*phiL(k)+M0R(2)*phiR(k))+...
P0L(1)*PL_V(k)+P0R(2)*PR_V(k);
end
end
for s=1:nx
residual_x(:,s,j)=-(Fx(:,s+1)-Fx(:,s))/dx;
end
end
end
function M0=ausmsplit_Ma(M)
M0=zeros(1,2);
if(abs(M)>1.0)
Mp=1/2*(M+abs(M));
Mm=1/2*(M-abs(M));
else
Mp=1/4*((M+1)^2)+1/8*(M^2-1)^2;
Mm=-1/4*((M-1)^2)-1/8*(M^2-1)^2;
end
M0(1)=Mp;
M0(2)=Mm;
end
function P0=ausmsplit_P(M)
P0=zeros(1,2);
if(abs(M)>1.0)
Pp=1/2*(1+sign(M));
Pm=1/2*(1-sign(M));
else
Pp=1/4*(M+1)^2*(2-M)+3/16*M*(M^2-1)^2;
Pm=1/4*(M-1)^2*(2+M)-3/16*M*(M^2-1)^2;
end
P0(1)=Pp;
P0(2)=Pm;
end
function [phi,P0]=AUSM_turnvar_x(gama,q)
phi=zeros(1,4);
P0=zeros(1,4);
rho=q(1);
u=q(2)/q(1);
v=q(3)/q(1);
e=q(4)/q(1);
p=rho*(gama-1.0)*(e-0.5*(u^2+v^2));
phi(1)=rho;
phi(2)=rho*u;
phi(3)=rho*v;
phi(4)=q(4)+p;
P0(1)=0;
P0(2)=p;
P0(3)=0;
P0(4)=0;
end
2.3 时间推进
采用显式的四步Runge-Kutta格式来完成：
3.计算结果
密度云图
时间=7.5^-5s:
参考文献：潘沙.高超声速气动热数值模拟方法及大规模并行计算研究[D].长沙:国防科学技术大学,2010,4.
Kim K H,Kim C,Rho O H. Methods for the Accurate Computations of Hypersonic Flows[J].Journal of Computational Physics 174, 38–80 (2001) doi:10.1006/jcph.2001.6873
浙公网安备 33010602011771号