在枪炮武器、火箭和航空航天等工程技术领域都涉及燃气射流冲击的问题，射流波系结构的影响因素有很多^[1]。为了减少燃气射流对近场设备和工程设施的冲击危害，需要了解射流流场波系结构和流场的流动状态。由于相关的射流实验费用昂贵，且通过实验不可能详尽地研究各种因素在不同水平对射流波系结构的影响。然而，采用高精度的数值方法能够有效地捕捉燃气射流场的激波波系结构，并且与相同工况下的实验结果能很好地吻合^[2-3]。

近十几年来，基于非结构/结构网格的高精度算法发展迅速，主要有TVD、DGM、ENO、k-exact有限体积方法等^[4-6]。从A.Harten等提出ENO格式以后，许多学者从不同的思路出发，提出了多种形式的ENO格式，如WENO、CENO^[7-8]。徐文灿等^[9]采用高分辨率的ENO格式对喷管流动问题进行了数值计算，得到了推力随摆角变化的规律，与实验结果基本一致，证明了采用ENO格式能够很好地捕获复杂波系，反映激波和边界层之间的相互干扰。陆霄露等^[10]也采用ENO格式计算了一维非定常进排气流动问题，计算结果表明发动机的主要性能参数都和实验结果吻合很好。

为了清晰地捕捉到流场中的间断区域(激波结构)，研究一种有效求解Lagrange坐标系中积分形式欧拉方程的三阶ENO有限体积法，并且采用具有精确解的激波管算例验证该方法的有效性。最后，采用该方法求解典型的欠膨胀射流的流动问题。

1.   控制方程

在Lagrange坐标系中，二维轴对称非稳态可压缩气体流动的积分形式欧拉方程^[11]为：

式中:Ω(t)为由边界Γ(t)包围的运动控制体，守恒变量矢量U=[ρ, M, E]^T，通量矢量F=[0, pn, pu·n]^T。其中，ρ为密度，u为速度矢量，M=(M_x, M_y)=ρu为动量，E为总能量，n=(n_x, n_y)为边界Γ(t)的外法线向量。理想气体的状态方程压力p=(γ-1)ρe, e为气体的质量内能。

2.   数值计算方法

二维空间域划分为m×n个计算单元。Ω_{i+1/2, j+1/2}为一个四边形的计算单元，它的顶点分别为(x_{i, j}, y_{i, j})、(x_{i+1, j}, y_{i+1, j})、(x_{i+1, j+1}, y_{i+1, j+1})、(x_{i, j+1}, y_{i, j+1})。S_{i+1/2, j+1/2}为计算单元Ω_{i+1/2, j+1/2}的面积。采用格心型有限体积法，所有的物理量都存储在计算单元的中心。因此，采用均值形式的密度ρ_{i+1/2, j+1/2}、动量M_{i+1/2, j+1/2}以及总能量E分别如下式所示：

式中:M_x、M_y分别为x方向和y方向的动量。

采用有限体积法，对控制方程式(1)进行离散，得到轴对称Euler方程组的半离散守恒方程形式:

式中:k为计算单元4个顶点的序号，按逆时针方向排序；p_k，k+1为计算单元k顶点和k+1顶点边界上的压强。由于采用有限体积方法得到的是网格平均值p，由p可以构造高阶插值多项式p(x, y)，但不能保证p(x, y)在网格边界上是连续的，因此p_k，k+1不能直接从构造的插值多项式p(x, y)求得。这里需要解一个Riemann问题，通过求精确Riemann解，可以得到p_k，k+1的值。而法向速度u_n和网格侧面积的乘积取为:

式中:u_k为计算单元Ω_{i+1/2，j+1/2}中序号为k的顶点在x方向的速度。同理，v_k为y方向的速度，(x_k, y_k)为该顶点的坐标值。Δt时间后网格中心新速度$\mathit{\boldsymbol{\hat u = }}\left( {\hat u, \hat v} \right)$和能量$\hat E$可由式(4)、(5)联立求得。计算单元顶点的新速度采用格点型格式直接计算，构造顶点处速度的控制体为图 1中虚线边界Γ′_k包围的控制体Ω′_{i+1/2，j+1/2}，在二维空间上进行三阶ENO重构，则有插值多项式：

图  1  速度的控制体

Figure  1.  Control volume of velocity

下载: 全尺寸图片 幻灯片

式中:q₀、q_x、q_y、q_xx、q_xy、q_yy为待定系数。将Ω′_{i+1/2, j+1/2}边界的线段用参数方程x=x₁+(x₂-x₁)t, y=y₁+(y₂-y₁)t，根据文献[10]，计算单元顶点运动的位置$\left( {{{\hat x}_{ij}}, {{\hat y}_{ij}}} \right), {\hat x_{ij}} = {x_{ij}} + \Delta t{\hat u_{ij}}, {\hat y_{ij}} = {y_{ij}} + \Delta t{\hat v_{ij}}$。根据每个顶点的新位置求得计算单元新的面积${\hat S_{i + 1/2, j + 1/2}}$，进而得到新的密度${\hat \rho _{i + 1/2, j + 1/2}} = {w_{i + 1/2, j + 1/2}}/{\hat S_{i + 1, j + 1/2}}$，其中w_{i+1/2, j+1/2}为计算单元内流体质量。

联立式(2)~(4)，在结构网格下采用高阶ENO-FVM格式，方程(1)能够表示^[4]为：

求解过程忽略3阶以上高阶项的影响，方程(1)的求解转化为常微分方程的求解问题。常微分方程的求解采用具有TVD性质的三阶Runge-Kutta方法，如下式^[10]所示：

式中:时间步长Δt根据式Δt^(r)=λmin(Δl_{i+1/2, j+1/2}^(r)/c_{i+1/2, j+1/2}^(r))来确定，其中，r表示第r个计算步；Δl_{i+1/2, j+1/2}^(r)为Γ_{i+1/2，j+1/2}边界中最短边的长度；c_{i+1/2, j+1/2}^(r)为计算单元Ω_{i+1/2, j+1/2}内的声速；计算中Courant数λ取为0.5。

3.   计算结果与分析

3.1   Sod激波管算例

计算域为[0, 2]，所采用的网格数为100，初始条件为:

采用Newmann边界条件，取CFL系数为0.5，计算推进到t=0.5 s时终止计算。此时，流场中包含一个激波、一个接触间断和一个膨胀波。图 2给出了采用数值方法得到的密度分布曲线。从图 2可以看出，三阶ENO有限体积法的激波分辨率很高，能准确捕捉到激波结构，将激波被抹平的厚度限制在一个网格单元尺度。表 1给出了不同网格尺度下，三阶ENO有限体积法的密度误差L1范数及其收敛精度^[6]。从表 1可以看出，随着网格尺度的减小，计算方法体现出了网格收敛性，并且收敛的精度大致相当，计算时间以指数倍增长。

图  2  Sod激波管密度分布

Figure  2.  Density profiles for the Sod problem

下载: 全尺寸图片 幻灯片

3.2   喷管出口欠膨胀射流问题

根据文献[1]，计算了喷管出口欠膨胀射流例子，计算域为喷管子午面截得的一半区域，如图 3所示。图中AB表示喷管的出口半径d_e，AC表示喷管的中心轴线，上游边界BD和外边界DE为静止大气条件，下游边界CE采用外推。采用无量纲格式计算，参考长度d_e=10 mm，参考压力p₀=101.325 kPa，参考温度T=300 K。计算区域为[0, 9]×[0, 4]计算区域内布置270×160个网格，整个计算域初始时为静止大气条件，在t=0时刻射流突然喷出。计算了两种工况，工况1轻度欠膨胀，喷管出口压力比为2.06；工况2重度欠膨胀，喷管出口压力比为26.4。

表  1  三阶ENO有限体积法密度误差

Table  1.  Density error for third-order ENO finite volume method

网格尺度	L1误差	计算时间/s	收敛精度
0.02	3.03×10-3	0.072 6
0.01	4.09×10-4	0.296 2	2.89
0.005	3.90×10-4	1.102 2	2.96

下载: 导出CSV 

| 显示表格

图  3  计算区域

Figure  3.  Computational domain

下载: 全尺寸图片 幻灯片

图 4(a)~(c)分别给出了工况1在第1 000个计算时间步时刻的马赫数、压力、密度等值线分布。从图 4可以看出，射流的流场结构类似钻石状，马赫盘略有呈现，入射激波、反射激波、马赫盘等构成的波胞交替产生。在下游边界处，由于射流与静止大气之间的对流作用，流场中产生了间断面。由于Euler方程可以看作是高Reynolds数下的近似，因此边界层区域附近将出现Kelvin-Helmholtz不稳定性，越往下游不稳定性越明显。计算结果与正格式数值计算的结果^[12]和实验纹影照片图 4(d)反映的流场结构吻合较好。

图  4  工况1下计算结果等值线图和纹影照片

Figure  4.  Contour maps and schlieren photograph in case 1

下载: 全尺寸图片 幻灯片

图 5(a)~(c)分别为工况2下马赫数、压力和密度等值线分布。与工况1相比，由于喷管出口压力比增大，图 5中有明显的马赫盘结构，在马赫盘的边缘与入射激波和反射激波交汇形成三叉激波，出现三波交点和滑移线，马赫盘上游射流的结构较稳定，下游流场出现了明显的Kelvin-Helmholtz不稳定性。由于出口压力比增大，计算域内只能形成一个波胞，并且马赫盘的过渡区间很窄，说明本文中采用的方法具有较强的捕捉激波的能力。图 5(d)为相同流动条件下的实验纹影照片^[13]，实验结果中第一个马赫盘距喷管出口平面的距离为4.56d_e。采用数值方法得到的马赫盘的位置为4.51d_e，文献[12]中采用正格式的结果为4.68d_e，两种算法相比较，采用本文方法得到的马赫盘位置更精确。与实验得到的马赫盘位置比较，本文计算结果的误差小于1.1%，且与实验纹影照片所反映的流场结构吻合较好，表明该方法具有较高的精度，能精确捕捉激波位置，并且在激波面上不会产生振荡或抹平间断现象。

图  5  工况2下计算结果等值线图和纹影照片

Figure  5.  Contour maps and schlieren photograph in case 2

下载: 全尺寸图片 幻灯片

4.   结论

以ENO格式为基础，通过在单元交界面处进行三阶ENO格式插值，构造了一类求解Lagrange坐标系中积分形式欧拉方程的三阶ENO有限体积法。数值计算结果表明，该方法具有较高的数值精度，适用于非稳态流场的数值计算，并且具有较强的激波捕捉能力，能够清晰地模拟出复杂流场的射流结构，与相同工况下实验结果吻合较好。