0 引言

特征线法是一种数值求解无粘流动的重要方法。关于二维特征线的方法和应用，国内外都有了很多相关成果^[1-4]，相比之下三维特征线的应用较少^[5]。而现阶段随着机体/发动机一体化设计概念的提出，对超声速流道高效转弯变形的需求也日渐增长，圆形和椭圆形等流道设计再次成为国际上的研究热点^[6-7]，这类流道通常需要实现三维空间复杂曲面的均匀过渡，因此流动的三维特性进一步增加了对高精度、纯三维数值方法的需求。

在三维气动设计的研究过程中出现过许多精妙的构思，如Busemann进气道，流线追踪技术及密切理论、三维乘波体设计等。但是流线追踪方法主要限制于基准流场的设计^[8-10]；在一定强几何约束条件下，密切理论难以直接应用于三维流场的精确组织^[11-12]；流线几何融合能够实现变截面压缩，但其融合主要是从几何的角度出发，并未考虑到流场的特征，使得压缩不好控制而出现流向涡^[13]。因此寻找一种新的三维超声速流场的设计方法至关重要。

2013年南京航空航天大学的全志斌等人设计了壁面压力分布可控的单边膨胀喷管^[14]。2015年，国防科学技术大学的赵玉新等人提出一种三维超声速流动的压力反问题，并构造了边界反Riemann求解器^[15]，实现了入口为矩形、椭圆形、扇环形流管的三维反设计。2016年厦门大学的李怡庆等人设计了壁面压力分布可控的高超声速进气道^[16]。此外文献[17]中提出了全三维超声速流场的求解方法，以及基于特征线追踪的气动反设计方法，现已成功运用于圆转方超声速流道的设计中^[18-19]。

为了进一步探索纯三维超声速流场压力反问题的求解方法，本文通过分析三维特征线方法^[20-21]，采用预设壁面压力分布计算壁面型线的思想，提出了基于双特征线的压力反问题求解方法。为了验证设计方法的可靠性，本文对设计的喷管构型进行CFD数值模拟，并将其与设计流场进行分析比较。结果表明：基于双特征线方法的压力反问题求解器具备三维超声速与高超声速气动设计的能力，同时具有高精度、纯三维、易控制的特性，对未来高超声速气动设计应用将起到重要的支撑作用。

1 双特征线算法 1.1 特征方程及相容方程

三维定常、无粘、可压、等熵流动的控制方程为：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

其中，u、v、w为x、y、z方向的速度，p为静压，ρ为密度，a为声速，P为总压，H为总焓，下标x、y、z表示对x、y、z坐标方向的方向导数。这些方程可以写成如下所示的定常三维等熵超声速流动的特征方程和相容性方程^[2]。

特征方程：

(8)

(9)

相容方程：

(10)

(11)

(12)

其中：n=in_x+jn_y+kn_z，为特征法线，a为当地声速，下标t表示在特征方向的方向导数。

1.2 双特征线方法单元过程

在进行数值求解时，特征方程决定离散网格点，相容方程决定解点流动参数。本文采用的是Butler提出的五面体双特征线网格^[20-21]，如图 1所示。该格式采用空间步进的方法，使得解在一系列平行的解平面上推进。

图中点5为微团轨迹与初值面的交点，点6为解点，其位置由点5向前延伸微团轨迹来确定。随后四条间隔相同的双特征线从点6向后延伸，与初值面在点1、2、3和点4处相交。点1~点5上的流动参数由拟合了点5和八个相邻点的双变量内插多项式来确定。

此外，Butler还引进了双特征线的参数化，如下式所示。

(13)

其中dx_i是双特征射线，t是描述射线长度的参数，θ是规定的双特征线的参数角，q表示速度大小，在点6上的u_i/q，α_i，β_i组成一个正交单位向量。c为马赫锥在与微团轨迹垂直的平面上的扩张速度，它由下式确定：

(14)

其中a表示当地声速。

Butler还从原先偏微分方程的线性组合导出一个在微团轨迹上能够成立的非特征线关系式：

(15)

其中d_up定义为压力沿流线切向的导数。

四个双特征线相容性方程可以用θ=0，π/2，π，3π/2写出，再将这四个方程与式(15) 进行有限差分。这五个有限差分方程可进行线性组合，从而获得三个相互独立的方程，同时也能消除解点上的交叉倒数。这三个独立的方程加上在微团轨迹上成立的两个相容性方程，用改进的欧拉预估-校正法进行求解。为达到二阶精度，在校正步时对相容性方程系数取平均。

对于无粘情况下的固体边界，流体必须与边界相切。用流线相切的条件代替四条双特征线之一的相容性方程，内点单元即可修改为边界单元。流动相切条件为：

(16)

其中n_i是点6固体边界的单位外法线。

在沿流动方向进行空间推进时，逐个截面间的距离必须调整到使网格上所有点均满足CFL稳定性准则。可允许的步长大小是当地流动参数和点间距离的函数，由下式给出：

(17)

其中R_min是流线与初值面的交点和用来构成有限差分网格凸包的点中离得最近的那个点之间的距离。

2 压力反问题求解器 2.1 特征线求解技术

如图 2所示，在点5建立基准面，该基准面由向量OP₁和向量OP₂共同构成，其中OP₁与点5左右两点的连线相互平行，OP₂为向量(1, 0, 0)。直线L56为流线微元。密切面由矢量V和上游的β_i共同构成，V为点5和点6的平均速度矢量，β_i为流动相切条件中的n_i。流线L56为由密切面和旋转θ角的基准面相贯获得，旋转角度的正方向满足旋转轴为OP₁的右手定则。因此存在一个单调关系：

(18)

因此，P_6i会随着θ_i的增大而减小，反之亦然。通过给定点6压力的设计值，反求对应的θ值，从而获得点6坐标，这就是压力反问题。

2.2 方法精度验证

本小节对该求解器的精度进行验证。设计的算例模型如图 3所示，外折角为20°，来流马赫数为1.4，总压为101 325Pa，总温为288K，气体常数为287J/(kg·K)，比热比为1.4。假设图中矩形区域内未受到上壁面反射膨胀波的影响，故其精确流场可用普朗特-迈耶函数求得。图 4中黑色等值线为马赫线的精确解，从顶点处(0.2449，0.25) 追踪一条流线，将其上的压力分布作为反问题求解器的输入项。

初始截面所在位置已在图 3中标示，设z方向垂直于纸面向外，则初始截面在yz平面内，分别沿y、z方向设定40、20个计算节点。x为流动方向，沿流向步进80步，获得如图 4所示的马赫数等值线，其中红色等值线为求解器求得，黑色等值线为精确解。从图 4中可以看出，该方法所求得的马赫线和理论几乎完全重合。壁面上各流动参数与精确解的最大相对误差见表 1。从表 1中可以看出，求解器继承了特征线高精度的特性。

3 三维超声速喷管设计算例

为了验证压力反问题求解方法的有效性，针对进口为圆形、椭圆形的喷管进行设计，入口条件为马赫数为1.05的平行均匀流，总压为101 325Pa，总温为288K，比热比为1.4，气体常数为287J/(kg·K)。在壁面上调用压力反问题求解器，内部流场采用双特征线方法进行计算。

若预设的壁面压力变化过于剧烈，容易出现强压缩或者膨胀波，因此本文计算了一个超椭圆喷管，将其壁面压力分布作为一个较为合理的壁面压力分布，如图 5所示。其中超椭圆的壁面函数为：

(19)

其中：z₀=y₀=0.1+0.4x²。

3.1 圆形入口的喷管设计

喷管的进口截面设为R=0.1m的圆，考虑空间的对称性，为了减少计算量，仅计算1/4圆。初始截面的结构网格采用放射状的方式划分。x方向为流动方向，喷管沿-x向的视图如图 6所示。沿流动方向空间步进400次，步进距离为0.347m，网格量为21×21×400，在主频为2.6GHZ的PC上运行时间需要120s。图 6中可以看出，每一条沿j向的网格线均表示一个流面切片，每一个流面均通过轴线，说明气体在喷管内部严格地沿着径向流动，符合轴对称规律。

图 7为给定压力的原始喷管型线与设计型线的对比，其中绿色点表示原始喷管轴对称面上的壁面型线，红色线条为设计喷管的壁面型线。图 8为对比了轴线和壁面的马赫数分布，同样绿色的点表示给定的流场数据，红色线条表示相应的设计值。通过计算，中心马赫数和壁面马赫数的最大相对误差分别为0.07%和0.1%。因此，该方法高精度地还原了给定压力的原始构型和原始流场。

为了进一步证明该方法对原始流场的还原程度，本文比较了出口截面沿45°的参数分布，如图 9~图 11所示，对比了马赫数、压力、以及x轴的速度分量随y的变化曲线。图中可以看出，点和曲线几乎完全重合，其最大相对误差如表 2所示，进一步证明了设计程序的可靠性。

3.2 椭圆形入口的喷管设计

喷管的入口截面为椭圆形，长轴长0.15m，短轴长0.1m，流动参数与轴对称喷管一致，同样是马赫数为1.05的平行均匀流。初始截面采用放射状网格，x方向表示流动方向。图 12为喷管沿-x方向的视图。沿流动方向空间步进400次，步进距离为0.318m，网格量为21×21×400。图中的轴线、边界流线、进口截面沿j向的网格线，出口截面沿j向的网格线，这四条线共同组成的曲面为流动的流面。从流面的三维特性可以看出，椭圆喷管内部存在从y轴向z轴的周向流动。由于初始截面为椭圆，长轴膨胀效率高，膨胀快，为了在周向获得相同的压力，长轴的扩张角就要小于短轴的扩张角，使短轴长度获得更快的增长速率，因此出现喷管横截面从椭圆逐渐向圆形过渡的现象。

图 13为截面马赫数等值线图。从图 13中可以看出，椭圆喷管内部的马赫数等值线均匀平滑，无明显的膨胀波，具有较高的流场品质。

为验证设计方法的优劣，对椭圆喷管进行CFD数值模拟，将CFD计算的流场与设计的结果进行对比分析。图 14(a)~(d)分别为三维压力云图和出口截面的压力等值线图、马赫数等值线图，其中CFD计算得到的等值线为红色，3D-MOC设计流场的等值线为黑色。

为了获得更加精确地对比，本文分别在设计流场和数值模拟流场的出口截面，提取了沿45°方向的流动参数，如图 15所示，分别对比了压力、马赫数。可以发现，在壁面和中心处各项参数均有不同程度的误差，但最大相对误差均为千分之几量级，甚至为万分之几(表 2)。

图 16分别对比了沿中心线及长短轴壁面的马赫数和压力分布。图 16中可以看出，CFD结果与设计结果匹配得好，同时长轴壁面和短轴壁面对应的压力曲线吻合程度高，说明能够实现壁面沿流向的一维压力分布和一维马赫数分布。

表 2中列出了喷管出口截面与CFD数值模拟对比的最大相对误差，由于圆形进口的喷管具有较好的轴对称特性，流动三维特性小，网格的正交性能够一直保持至出口，因此精度相对较高。而对于椭圆形进口，喷管内部存在周向流动，网格也随着流面的扭曲发生相应的扭动，使得网格的正交性略有下降，影响了单元过程中双线性内插的插值精度，同时该误差数值在一定程度上包含了CFD求解的部分误差，因此误差相对较大，但是最大相对误差在千分之几量级，其设计精度仍能够保证要求。

4 结论

综上所述，基于双特征线方法的压力反问题求解技术具备三维超声速与高超声速气动设计的能力，可应用于三维超声速喷管设计。该方法未引入任何平面假设，为纯三维的气动设计方法；最大误差在千分之几量级；可通过设计压力进行壁面构型的设计，使得壁面流场参数控制更加灵活。

在设计过程中仍然存在许多问题还有待深入探究，如：流管内部的三维激波问题较为复杂、控制喷管出口形状需要复杂的壁面压力设计。因此，还需要继续发展新的技术来完善该求解技术。