1 引言

半个多世纪以来的高超声速研究极大地推动了可压缩空气动力学的发展，在此期间，超声速与高超声速气动设计方法也取得了较大的进展。特征线法、流函数坐标、残差修正、伴随方法、优化算法等在不同的时期引领了超声速气动设计的发展方向^[1~3]。其中特征线法是一种数值求解超声速无粘流动的重要方法。二维特征线法在超声速流场设计中得到广泛应用^[4~6]，相比之下，三维特征线法的推广遇到了较大的困难。

三维特征线主要包括参考平面法^[7~9]和双特征线方法^[10~13]。参考平面法较双特征线方法具有一些优势：在参考平面上的有限差分网格较为简单，可以非常方便地使用特征线网格; 并且可以引入二维特征线中发展较为成熟的激波处理等技术; 与密切理论相类似，将复杂的型面进行降维设计，但参考平面法可以考虑垂直于平面的流动梯度。

国内外对三维超声速复杂气动外形的设计研究有着悠久的历史。Thompson等^[14]采用变分法对三维喷管的设计问题展开了研究。Wang等^{[15, 16]}采用CFD结合模拟退火的优化方法对三维圆转椭圆喷管的过渡型面进行了一些优化分析。此外还有一些三维喷管的直接设计方法，如等长度等截面积转换法，超椭圆面积转换法等^{[17, 18]}。在冲压发动机进排气系统的设计领域，能够实现三维超声速气动设计且应用较为成熟的方法主要包括流线追踪技术和密切理论。

Smart等^{[19, 20]}基于流线追踪技术和Busemann基准流场^[21]提出REST进气道的设计方法，是进气道设计方法的创新。尤延铖等提出内乘波式进气道的设计概念^[22]，且具有较高的三维压缩效果。莫建伟等采用双向流线追踪技术，进行进出口形状均可控的三维变截面(圆转方)非对称喷管的设计^[23]。

密切理论、密切轴对称理论自Jones等^[24]提出以来，广泛应用于内外流压缩流场的组织。在此基础上，密切锥乘波体设计方法^[25]、双密切轴对称近似设计方法^[26]、三维变截面尾喷管设计方法^[27]等不断出现。但在一定几何约束条件下，密切理论难以直接应用于强三维流场的精细组织。

为了进一步探索三维超声速的气动设计方法，本文从降维的三维特征线方法出发，提出三维曲面激波反问题的求解方法。为了验证设计方法的可靠性，采用参考平面法设计了锥形流和平面激波，并用泰勒麦科尔流动和斜激波关系验证了数值方法的精度; 设计了在3°来流攻角下产生圆锥激波的三维型面，并与CFD数值模拟结果进行对比。

2 参考平面解法 2.1 特征方程与相容方程

柱坐标系下三维定常、无粘、可压、等熵流动的控制方程为

$ \begin{array}{*{20}{l}} {u\frac{{\partial \rho }}{{\partial x}} + v\frac{{\partial \rho }}{{\partial r}} + \rho \frac{{\partial u}}{{\partial x}} + \rho \frac{{\partial v}}{{\partial r}} = }\\ { - \frac{w}{r}\frac{{\partial \rho }}{{\partial \phi }} - \frac{\rho }{r}\frac{{\partial w}}{{\partial \phi }} - \frac{{\rho v}}{r}} \end{array} $

(1)

$ u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial r}} = - \frac{w}{r}\frac{{\partial u}}{{\partial \phi }} - \frac{1}{\rho }\frac{{\partial p}}{{\partial x}} $

(2)

$ u\frac{{\partial w}}{{\partial x}} + v\frac{{\partial w}}{{\partial r}} = - \frac{w}{r}\frac{{\partial w}}{{\partial \phi }} - \frac{{vw}}{r} - \frac{1}{{\rho r}}\frac{{\partial p}}{{\partial \phi }} $

(3)

$ u\frac{{\partial w}}{{\partial x}} + v\frac{{\partial w}}{{\partial r}} = - \frac{w}{r}\frac{{\partial w}}{{\partial \phi }} - \frac{{vw}}{r} - \frac{1}{{\rho r}}\frac{{\partial p}}{{\partial \phi }} $

(4)

$ \begin{array}{*{20}{l}} {u\frac{{\partial p}}{{\partial x}} + v\frac{{\partial p}}{{\partial r}} - {a^2}\left( {u\frac{{\partial \rho }}{{\partial x}} + v\frac{{\partial \rho }}{{\partial r}}} \right) = }\\ { - \frac{w}{r}\frac{{\partial p}}{{\partial \varphi }} + {a^2}\frac{w}{r}\frac{{\partial \rho }}{{\partial \varphi }}} \end{array} $

(5)

此外

$ \frac{{{u^2} + {v^2} + {w^2}}}{2} + \frac{{{a^2}}}{{\gamma - 1}} = \frac{{Ma_0^2}}{2} + \frac{1}{{\gamma - 1}} $

(6)

$ \frac{p}{{{\rho ^\gamma }}}{{\rm{e}}^{ - \left( {\gamma - 1} \right)S}} = {\rm{const}} $

(7)

$ \mathit{\boldsymbol{u}} \cdot \nabla S = 0 $

(8)

式(1)为连续方程，式(2)~式(4)为动量方程，式(5)为声速方程。其中，速度$u, {\rm{}}v, {\rm{}}w$由来流声速进行无量纲化，$x, {\rm{}}r$由物理长度进行无量纲化，密度$\rho$，压力$p $，熵$ S$分别由来流密度，$ \gamma {p_0}, \left( {S - {S_0}} \right)/R$进行无量纲化，$\gamma $为气体比热比，${p_0} $为来流总压，$ R$为气体常数。

在参考平面方法中，特征线被定义在一个参考平面内，其物理意义为马赫锥及空间流线在该平面上的投影。如果认为式(1)~式(5)的右边项都是已知的，则这些方程组仅依赖于确定参考平面的两个坐标，即为(x，r)。对这两个独立变量进行特征分析，可以推导出适用于参考平面法的特征方程和相容方程。参考面内的特征曲线也称为近特征线。沿近特征线成立的相容方程足以确定解点上除了垂直于参考平面方向的速度分量以外的所有流动参数。剩下的这一个速度分量由垂直于参考平面方向的动量方程来确定，如式(4)所示。

柱坐标系$\left( {x,r,\phi } \right)$参考平面法的特征方程和相容方程如下：

特征方程

$ {\left( {\frac{{{\rm{d}}r}}{{{\rm{d}}x}}} \right)_0} = {\lambda _0} = \frac{v}{u} $

(9)

$ {\left( {\frac{{{\rm{d}}r}}{{{\rm{d}}x}}} \right)_ \pm } = {\lambda _ \pm } = {\rm{tan}}\left( {\theta \pm \mu } \right) $

(10)

相容方程

$ {{\rm{d}}_ \pm }\theta \frac{{\sqrt {Ma_*^2 - 1} }}{{\rho {q^2}}}{{\rm{d}}_ \pm }p = \frac{{{F_1}c{\rm{os}}\mu \mp {F_2}{\rm{sin}}\mu }}{{{\rm{cos}}(\theta \pm \mu )}}{{\rm{d}}_ \pm }x $

(11)

$\rho q{{\rm{d}}_0}q + {{\rm{d}}_0}p = \left( { - \frac{w}{r}\frac{{\partial q}}{{\partial \phi }} + \frac{{{w^2}}}{r}{\rm{sin}}\theta } \right)\rho {{\rm{d}}_0}s$

(12)

$q{{\rm{d}}_0}w = \left( { - \frac{w}{r}\frac{{\partial w}}{{\partial \phi }} - \frac{{vw}}{r} - \frac{1}{{\rho r}}\frac{{\partial p}}{{\partial \phi }}} \right){{\rm{d}}_0}s$

(13)

$q({{\rm{d}}_0}p - {a^2}{{\rm{d}}_0}\rho ) = ( - \frac{w}{r}\frac{{\partial p}}{{\partial \phi }} + {a^2}\frac{w}{r}\frac{{\partial \rho }}{{\partial \phi }}){\rm{d}}{s_0}$

(14)

式(11)中的${F_1} $，$ {F_2}$为

${F_1} = \frac{{{w^2}}}{{{q^2}r}}{\rm{cos}}\theta - \frac{w}{{qr}}\frac{{\partial \theta }}{{\partial \phi }}$

(15)

$\begin{array}{*{20}{l}} {{F_2} = \frac{w}{{\gamma qpr}}\frac{{\partial p}}{{\partial \phi }} - \frac{w}{{{q^2}r}}\frac{{\partial q}}{{\partial \phi }} + }\\ {\frac{1}{{qr}}\frac{{\partial w}}{{\partial \phi }} + \frac{{{\rm{sin}}\theta }}{r}\left( {1 + \frac{{{w^2}}}{{{q^2}}}} \right)} \end{array}$

(16)

式中θ，q分别为参考平面内的流动方向角、平面速度大小，$M{a_*}$为平面马赫数，定义为

$Ma_{\rm{*}}^2 = \left( {{u^2} + {v^2}} \right)/{a^2}$

(17)

因此，从式(11)中可以看出，只有当平面马赫数大于1时，平面内的特征锥才是实的。故该方法不适用于横向速度过大的流动，因为尽管来流为超声速，但垂直于参考平面速度太大时，$M{a_*}$仍可能小于1。

从式(9)~式(16)中可以看出：参考平面法的特征方程和相容方程与二维(或者轴对称)特征线方法具有相同形式，当方程右手端的周向偏导数($ \frac{\partial }{{\partial \phi }}$项)均为0时，方程退化为二维/轴对称特征线方程。

参考平面法的主要优势是简单，仅用了三条近特征线，而且流动参数在初值点上的内插是简单的线性内插，解可以在参平面内灵活地采用正步或者逆步进法。这个方法已被Moretti^[28]，Katskova等^[29]，Rakich^{[7, 8]}成功应用于定常三维超声速流动中。

2.2 参考平面解法

在参考平面内，已知图 1中直线L12为激波微元，点1和点2为已知值，其流动参数可根据斜激波关系式求得。点2发出的左行特征线${C_ + } $与1点发出的流线${C_ 0 } $相交，交点即为待求点，标记为红色点4。从点4返回右行特征线${C_ - } $与L12相交，交点为点3，点3的流动参数可于初始线段L12上进行线性内插。预设点4中垂直于参考面的流动梯度均等于点1，而点1的周向偏导数($ \frac{\partial }{{\partial \phi }}$项)均可求，故式(9)~式(16)的右手端均为常数，采用欧拉预估校正法求解，可以获得待求点4的位置和流动参数。

在参考平面内重复图 1所示的单元过程，可以追踪出从激波点发的左行特征线，直至参考平面内的壁面点，即为图 2中激波点至壁面流线点的红色实、虚线，记为${C_ + }\left( 1 \right)、\left( i \right){C_ + }\left( 2 \right)\left( i \right)、{C_ + }\left( 3 \right)\left( i \right)、{C_ + }\left( 4 \right)\left( i \right) $，i表示第i个参考面。

n个参考平面内的左行特征线拼合成特征曲面，可用下式表达

${F_ + }\left( 1 \right) = \{ {C_ + }\left( 1 \right)\left( i \right)|i \le n, {\rm{}}i \in {\rm{N}}\} $

(18)

即${C_ + }\left( 1 \right)\left( i \right) $拼合成曲面，记为$ {F_ + }\left( 1 \right)\left( i \right)$，同理，${C_ + }\left( 2 \right) $对应${F_ + }\left( 2 \right) $，${C_ + }\left( 3 \right) $对应${F_ + }\left( 3 \right) $，${C_ + }\left( 4 \right) $对应$ {F_ + }\left( 4 \right)$。直观地说，轴对称流动中，${C_ + }\left( 1 \right) $绕轴线旋转一周形成的曲面即为${F_ + }\left( 1 \right) $。通过特征曲面的空间步进逐渐生成三维流场，即从${F_ + }\left( 1 \right), {F_ + }\left( 2 \right), {F_ + }\left( 3 \right) $至${F_ + }\left( 4 \right) $逐渐向下游步进，可用下式表达

${\rm{Flow~field}} = \left\{ {{F_ + }\left( i \right)|i \le k, {\rm{}}i \in {\rm{N}}} \right\}$

(19)

式中k表示步进次数。

以下详细介绍三维流场的生成过程，以锥形激波的波后流场为例，流场如图 3所示，这里旨在描述算法中相关概念的三维空间结构，流场的相关参数详见2.3.1节。

给定图 3中的锥形激波，求解对应的壁面。采用一系列参考平面(等$\phi $平面)划分锥形激波的波后流场，图 3中的流场切片即为其中的一个参考平面，图 4(d)中黑色网格所在流场切片同样表示为参考平面，且黑色网格的求解同图 2。

在求解第一个步进面$ {F_ + }\left( 1 \right)$时，每个参考面调用一次图 1所示的单元过程，从而确定了$ {C_ + }\left( 1 \right)\left( i \right)$，进而确定了${F_ + }\left( 1 \right) $，如图 4(d)中绿色网格所示; 继续步进，分别在第i个参考面中确定${C_ + }\left( 2 \right)\left( i \right) $，并拼合成${F_ + }\left( 2 \right) $; 步进至最后一个激波离散点，完成给定激波求解壁面的过程，红色曲面即表示最终计算求得的壁面。图 4的(a)~(c)分别为图 4(d)的正视图、左视图和俯视图。

采用如图 5所示的预估校正向下游空间步进，具有较好的收敛性，迭代过程为

(1) 预估${F_ + }\left( {k + 1} \right) $上的点的周向偏导数与曲面$ {F_ + }\left( k \right)$相等。

(2) 完成$ {C_ + }\left( {k + 1} \right)\left( i \right)$的求解，确定曲面$ {F_ + }\left( {k + 1} \right)$。

(3) 在曲面$中求解周向偏导数(式(9)~式(16)右手端的$ \frac{\partial }{{\partial \phi }}$项)。

(4) 返回步骤(2)直至收敛。

(5) 进行${F_ + }\left( {k + 2} \right) $的求解。

综上所述，参考平面法通过参考面内的局部迭代和垂直于参考面的周向迭代，以及彼此间的相互嵌套完成曲面的空间步进，从而逐渐生成三维复杂型面。

2.3 精度验证

为了验证设计方法的可靠性和精度，本节设计了两个算例：无攻角的锥形激波、平面激波，并与理论解进行对比验证。

2.3.1 锥形激波

给定均匀来流，马赫数为4，比热比为1.25，气体常数为287J/(kg·K)，流动方向为x方向。再设定锥形激波的半锥角为19.92°，锥顶点位置在(0，0，0)，轴线为x轴。在相同来流马赫数下。锥面激波开始脱体时的半锥角要大于平面激波的二维斜劈角，因此所设定的锥形激波具有物理解。

为了减少计算，考虑空间的对称性，仅取锥形激波的一半进行计算：$ \phi $∈[0，π]，$ \phi $为柱坐标系下的周向坐标分量。选取121个轴对称面作为参考平面。初值线为x=0.1m与圆锥激波相贯的环线，沿x轴空间步进80步，轴向最大距离为0.7m，网格量约为100万，所得流场及特征网络如图 6所示，其三维结构如图 4所示。

无攻角锥形激波的波后流场是无旋的，为锥型流，其控制方程是一个二阶非线性常微分方程(Taylor-Maccoll方程^{[30, 31]})，因此可将其作为评判图 6流场精度的标准。图 7中展示了流场压力随球面角的关系，其中黑色点为数值积分Taylor-Maccoll方程获得的无量纲压力分布，红色的线表示参考平面法解，可以看出，对于轴对称流动，参考平面法具有非常高的计算精度，其中压力的最大相对误差为0.128%，速度误差小至3.2×10^-5。

2.3.2 平面激波

上一小节中采用锥形流验证参考平面法的有效，但在该流动中，垂直于参考平面的流动梯度均为零，参考平面法实则退化成了轴对称的特征线方法。考虑到平面激波有理论解，且通过坐标变换(直角坐标系转化到柱坐标系)使得速度矢量具有周向梯度，因此采用平面激波验证周向偏导数不全为零时参考平面法的有效性。

如图 8所示，来流马赫数为4，流动参数均与上一小节相同，图中红色的平面激波过z轴，并与zx平面成19.92°夹角。初值线为x=0.1m与平面激波相贯的直线。将波后$ \phi $∈[-π/6，π/6]的流场用41个过x轴的参考平面划分，其中第1个参考面和第n个参考面与xy平面的夹角分别为-π/6和π/6。由于坐标系的变换，使得不同参考面具有相同的流场标量、不同的速度矢量，例如垂直于参考面方向的速度：参考面1 > 参考面2 > 参考面3 > …… > 参考面41。空间步进80步后，网格量约为33万，轴向最大距离为0.9m，如图 9所示，平面ABCD表示平面激波，红色面表示对应的壁面，绿色曲面为${F_ + }\left( {80} \right) $，即为

${F_ + }\left( {80} \right) = \{ {C_ + }\left( {80} \right)\left( i \right)|i \le 41, {\rm{}}i \in {\rm{N}}\} $

(20)

通过与斜激波关系式对比，压力、马赫数的最大相对误差分别为6.1×10^-5，4.8×10^-6，在一定程度上验证了参考平面法的有效性。

3 结果与讨论

高超声速飞行器在进行有攻角巡航时，进气道第一道激波的波后流场是三维的，而按照轴对称理论设计的第一道激波通常不能完全封住唇口。因此，可采用参考平面法求解三维激波反问题，使得三维激波的形状和强度得以控制。本节设计了在3°来流攻角情况下，产生圆锥激波的三维型面，即在设计状态进行3°攻角巡航时，高超声速飞行器头部产生圆锥形封口激波。

3.1 流场设计

为了保证解的存在，预设锥形激波形状与2.3.1节相同，即锥形激波的半锥角为19.92°，锥顶点位置在(0，0，0)。来流马赫数为4，流动参数均与上节相同。来流存在攻角，流动方向平行于xy平面，与x轴夹角为3°，如图 10所示。

为了减少计算，考虑空间的对称性，仅取锥形激波的一半进行计算：$\phi $∈[0，π]，$\phi $=0为背风面，$\phi $=π为迎风面。选取121个通过物理轴线的子午面作为参考平面。考虑到尖锥附近流动参数变化剧烈，为了保证参考面内的流动为超声速，暂时不考虑奇点附近的流动，故初值线选定为x=0.1m与圆锥激波相贯的环线，沿x轴向下游步进80步，轴向最大距离为0.7m，网格量约为100万，所得流场及特征网络如图 11所示，三维结构如图 12所示。

3.2 与CFD解对比

采用商用软件Fluent对设计的三维型面进行数值模拟，模拟的控制方程为Euler方程，网格划分如图 13(a)所示，其网格量约为262万; 作为对比，降低沿壁面法向的网格量，约为178万。三维压力场如图 13(b)所示。提取x=0.15，0.2，0.25，0.3，0.35m处的激波形状，并与设计值对比，如图 14所示，图中点表示不同空间位置处预设的激波位置，黑线表示数值模拟所得激波，可以看出，设计激波和计算激波基本一致，即：完成了在3°来流攻角情况下，产生圆锥激波的三维型面设计。

对比x=0.3m截面，y=0m和z=0.01m方向的马赫数分布，图 15所示为z=0.01m方向的马赫数分布，其中红线为参考平面法解，黑点为262万网格量下的CFD解，绿点为178万网格量下的CFD解。从图中可以看出，除了激波由于网格尺度的影响使其厚度有所不同外，结果基本一致，因此CFD解具有一定的可信性。通过对比，可得在y=0m方向压力、马赫数的最大相对误差分别为1.6%，0.27%;在z=0.01m方向压力、马赫数的最大相对误差分别为1.56%，0.35%。

3.3 三维型面分析

图 16为有攻角的锥形激波对应的三维型面，且在x=0.25m处叠加了一个虚线圆，定性地表明，迎风面壁面曲率半径大于背风面。图 17为x=0.25m处壁面轮廓的曲率半径，从图中可以看出，从迎风面至背风面，壁面的曲率半径逐渐减小。

4 结论

通过本文研究，得到如下结论：

(1) 本文提出了曲面激波反问题的参考平面解法，可在每个平面内设计激波曲线，通过考虑平面法向的偏导数保证三维流动计算的准确性。

(2) 采用采用泰勒-麦科尔流动和斜激波关系式验证了该设计方法的精度，其中压力的相对误差分别为1.3×10^-3，6.1×10^-5。

(3) 基于参考平面法设计了在3°来流攻角下产生圆锥激波的型面，并与传统的数值模拟进行对比，流场参数最大误差不超过2%。

(4) 在有攻角超声速自由流条件下，产生圆锥激波的三维型面具有从迎风面至背风面曲率半径逐渐降低的特征。

致谢: 感谢国家自然科学基金资助。