三维超声速流动的压力反问题

蓝庆生, 赵玉新,*, 赵一龙, 刘红阳

(1.国防科技大学 高超声速冲压发动机技术重点实验室, 湖南 长沙 410073;2.中国空气动力研究与发展中心 计算空气动力研究所, 四川 绵阳 621000)



三维超声速流动的压力反问题

蓝庆生1, 赵玉新1,*, 赵一龙1, 刘红阳2

(1.国防科技大学 高超声速冲压发动机技术重点实验室, 湖南 长沙 410073;2.中国空气动力研究与发展中心 计算空气动力研究所, 四川 绵阳 621000)

为了进一步探索三维超声速流道的设计方法，采用一种预设壁面压力分布计算壁面型线的思想，并结合双特征线方法提出一种全三维超声速流动压力反问题的求解方法。在三维超声速流场设计中，可直接根据来流条件和壁面压力分布求解壁面的三维坐标，通过空间步进的方式，使得解在一系列解平面上推进，从而使得所设计的型面与预设的壁面压力分布相容。通过Prandtl-Meyer膨胀波的理论解验证了该格式的设计精度。根据预设的压力分布，设计了圆形和椭圆形入口的三维超声速喷管，并将设计方法与数值模拟进行对比验证。验证结果表明：所设计的流场与CFD计算得到的等值线符合得较好，因此基于双特征线的压力反问题求解方法具备三维超声速气动设计的能力，并具有纯三维、高精度、壁面压力分布可控的优势，对未来高超声速气动设计应用将起到重要的支撑作用。

压力反问题；双特征线方法；三维超声速气动设计；壁面压力分布可控

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)

[u2-(V2-a2)](dx)2+[v2-(V2-a2)](dy)2+[w2-(V2-a2)](dz)2+2uvdxdy

(9)

相容方程：

(10)

(11)

ρanxut+ρanyvt+ρanzwt-pt

(12)

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

1.2 双特征线方法单元过程

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

图1 五面体双特征曲线网格Fig.1 Pentahedral line network

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

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

(14)

其中a表示当地声速。

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

(15)

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

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

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

(16)

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

(17)

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

2 压力反问题求解器

2.1 特征线求解技术

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

(18)

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

图2 压力反问题单元过程Fig.2 Unit process of pressure inverse problem

2.2 方法精度验证

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

图3 精度验证算例模型Fig.3 A case for accuracy validation

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

图4 特征线解与精确解对比Fig.4 Comparison between analytical and bicharacteristics solutions for PM flow表1 设计方法与Prandtl-Meyer精确解对比的最大相对误差Table 1 Maximal relative error of the design method comparing with Prandtl-Meyer flow

p/Paρ/(kg·m-3)u/(m·s-1)v/(m·s-1)MaP0/Pa87.441.38.60.29

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

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

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

(19)

其中：z0=y0=0.1+0.4x2。

图5 预设的壁面压力分布Fig.5 Preassigned pressure distribution

3.1 圆形入口的喷管设计

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

图6 圆形进口喷管沿-x方向的视图Fig.6 View of nozzle with circle inlet along -x direction

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

图7 设计喷管型线与原始喷管对比Fig.7 Comparison between designed and original nozzle shape

图8 设计马赫数与原始马赫数沿流向的对比Fig.8 Comparison between designed and original Mach number

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

图9 出口截面沿45°线的 马赫数分布对比Fig.9 Comparison for Mach number along 45° direction in outlet section

图10 出口截面沿45°线的 压力分布对比Fig.10 Comparison for pressure along 45° direction in outlet section

图11 出口截面沿45°线的x方向 速度分布对比Fig.11 Comparison for x-component of velocity along 45° direction in outlet section

3.2 椭圆形入口的喷管设计

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

图12 椭圆形进口喷管沿-x方向的视图Fig.12 View of nozzle with ellipse inlet along -x direction

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

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

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

(a) y=0截面

(b) z=0截面图13 截面的马赫数等值线图Fig.13 Mach number in different sections at y=0 and z=0

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

(a) 设计流场的三维压力云图

(b) 数值模拟流场的三维压力云图

(c) 压力等值线对比 (d) 马赫数等值线对比图14 三维压力云图及出口截面压力、马赫数的等值线图Fig.14 Contours of three-dimensional pressure, Mach number and pressure in outlet section

(a) 马赫数

(b) 压力图15 出口截面沿45°的马赫数分布、压力分布曲线Fig.15 Comparison for pressure and Mach number along 45° direction in outlet section

(a) 马赫数

(b) 压力图16 沿中心线及长短轴壁面的马赫数分布、压力分布曲线Fig.16 Comparison for pressure and Mach number along central axis, minor axis, and major axis

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

表2 出口截面各项参数的最大相对误差Table 2 Maximal relative error for parameters in outlet section

4 结 论

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

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

[1]Yi S H, Zhao Y X, He L, et al.Supersonic and hypersonic nozzle design[M].Beijing: National Defense Industry Press, 2013.(in Chinese)易仕和, 赵玉新, 何霖, 等.超声速与高超声速喷管设计[M].北京: 国防工业出版社, 2013.

[2]Zucrow M J, Hoffman J D.Gas dynamics[M].Beijing: National Defense Industry Press, 1984.(in Chinese)左克罗M J, 霍夫曼J D.气体动力学[M].北京: 国防工业出版社, 1984.

[3]Tong B G, Kong X Y, Deng G H.Gas dynamics[M].Beijing: Higher Education Press, 1900.(in Chinese)童秉纲, 孔祥言, 邓国华.气体动力学[M].北京: 高等教育出版社, 1990.

[4]Zhang M L, Yi S H, Zhao Y X.The design and experimental investigation of supersonic length-shorted nozzle[J].Acta Aerodynamica Sinica, 2007, 25(4): 500-503.(in Chinese)张敏莉, 易仕和, 赵玉新.超声速短化喷管的设计与实验研究[J].空气动力学学报, 2007, 25(4): 500-503.

[5]Ransom V H, Hoffman J D, Thompson H D.A second-order bicharacteristics method for three-dimensional, steady, supersonic flow[J].AIAA Journal, 1972, 10(12): 1573-1581.

[6]Beckel S A, Garrett J L, Gettinger C G.Technologies for robust and affordable scramjet propulsion[R].AIAA 2006-7980, 2006.

[7]Bulman M J, Siebenhaar A.The rebirth of round hypersonic propulsion[R].AIAA 2006-5035, 2006.

[8]Yates L A, Chapman G T.Streamlines, vorticity lines, and vortices around three-dimensional bodies[J].AIAA Journal, 1992, 30(7): 1819-1826.

[9]Billig F S, Kothari A P.Streamline tracing: technique for designing hypersonic vehicles[J].Journal of Propulsion and Power, 2000, 16(3): 465-471.

[10]Lu X, Yue L J, Xiao Y B, et al.Design of scramjet nozzle

based on streamline tracing technique[J].Journal of Propulsion Technology, 2011, 32(1): 91-96.(in Chinese)卢鑫, 岳连捷, 肖雅彬, 等.超燃冲压发动机尾喷管流线追踪设计[J].推进技术, 2011, 32(1): 91-96.

[11]You Y C, Liang D W.Design concept of three-dimensional section controllable internal waverider hypersonic inlet[J].Sci China Ser E-Tech Sci, 2009, 52(7): 2017-2028.(in Chinese)尤延铖, 梁德旺.基于内乘波概念的三维变截面高超声速进气道[J].中国科学: E 辑, 2009, 39(8): 1483-1494.doi: 10.1007/s11431-009-0125-1.

[12]Lu X, Yue L J, Xiao Y B, et al.Design of three-dimensional section controllable scramjet nozzle[C]//The 3th National conference on Hypersonic Technology.2010, Wuxi.(in Chinese)卢鑫, 岳连捷, 肖雅彬, 等.超燃冲压发动机三维变截面尾喷管设计[C]//第三届高超声速科技学术会议, 2010, 无锡.

[13]Smart M K.Design of three-dimensional hypersonic inlets with rectangular-to-elliptical shape transition[J].Journal of Propulsion and Power, 1999, 15(3): 408-416.

[14]Quan Z B, Xu J L, Mo J W.Design and experiment validation of single expansion ramp nozzle based on the control of wall pressure distribution[J].Journal of Propulsion Technology, 2013, 34(3): 307-310.(in Chinese)全志斌, 徐惊雷, 莫建伟.基于控制壁面压力分布的分级单边膨胀喷管设计及试验验证[J].推进技术, 2013, 34(3): 307-310.

[15]Zhao Y X, Liu H Y, Zhao Y L.The inverse problem of the three-dimensional supersonic flow[C]//The 8th National conference on Hypersonic Technology, 2015, Harbin.(in Chinese)赵玉新, 刘红阳, 赵一龙.三维超声速流动的反问题[C]//第八届全国高超声速科技学术会议, 2015, 哈尔滨.

[16]Li Y Q, Han W Q, You Y C, et al.Integration waverider design of hypersonic inlet and forebody with preassigned pressure distribution[J].Acta Aeronautica et Astronautica Sinica, 2016, 37(9): 2711-2720.(in Chinese)李怡庆, 韩伟强, 尤延铖, 等.压力分布可控的高超声速进气道/前体一体化乘波设计[J].航空学报, 2016, 37(9): 2711-2720.

[17]Zhao Y X, Liu H Y.The characteristic tracing method for supersonic flow field design[C]//The 8th National Conference on Hypersonic Technology, 2015, Harbin.(in Chinese)赵玉新, 刘红阳.基于特征线追踪的气动反设计[C]//第八届全国高超声速科技学术会议, 2015, 哈尔滨.

[18]Liu H Y.The characteristics tracing method and its applications[D].National University of Defense Technology, 2015.(in Chinese)刘红阳.特征线追踪方法和应用[D].国防科技大学, 2015.

[19]Liu H Y, Zhao Y X.Design of three-dimensional supersonic flow based on the characteristics tracing method with circle -to-rectangular shape transition.[C]//The Chinese Congress of Theoretical and Applied Mechanics, 2015, Shanghai.(in Chinese)刘红阳, 赵玉新, 基于特征线追踪的三维圆转方超声速流道设计[C]//中国力学大会, 2015, 上海.

[20]Ransom V H, Hoffman J D, Thompson H D.A second-order numerical method of characteristics for three-dimensional supersonic flow[D].Purdue University, 1970.

[21]Cline M C, Hoffman J D.Comparison of characteristic schemes for three-dimensional, steady, isentropic flow[J].AIAA Journal, 1972, 10(11): 1452-1458.

Pressure inverse problem of three-dimensional supersonic flow

LAN Qingsheng1, ZHAO Yuxin1,*, ZHAO Yilong1, LIU Hongyang2

(1.ScienceandTechnologyonScramjetLaboratory,NationalUniversityofDefenseTechnology,Changsha410073,China;2.ComputationalAerodynamicsInstitutesofChinaAerodynamicsResearchandDevelopmentCenter,Sichuan621000,China)

In order to explore valid methods for three-dimensional supersonic flow design, a technology was proposed to solve the pressure inverse problem by finding inviscid contour on the basis of a preassigned pressure distribution along the wall.A bicharacteristic algorithm is used to solve the inverse problem.Three-dimensional wall coordinates can be calculated directly according to the flow condition and predetermined boundary pressure.The designed three-dimensional contour can be compatible with preassigned pressure distribution by moving forward the solution plane along thexcoordinate.The accuracy order of the solver was tested by comparing numerical solutions with analytical results of Prandtl-Meyer expansion wave.Three-dimensional supersonic nozzles with round or oval inlet were designed in this paper on the basis of one-dimensional pressure distributions.The reliability of the solver was validated by corresponding CFD numerical simulations.It can be found that the numerical solutions of the present solver are in good agreement with those of CFD simulations.The solver has the ability for three-dimensional supersonic aerodynamic design by solving the inverse problem based on the bicharacteristics method.Moreover, this solver has advantages of high accuracy, pure three-dimensional solutions, and controllable wall pressure distribution.The method has promising potential of playing a supporting role in the applications of further hypersonic aerodynamic design.

pressure inverse problem;bicharacteristics method;three-dimensional supersonic aerodynamic design;controllable wall pressure distribution

0258-1825(2017)03-0429-07

2016-12-02;

2017-01-11

国家自然科学基金(11472304)

蓝庆生(1992-),男(畲族),福建龙岩人,硕士研究生,研究方向:高超声速气动设计.E-mail:lanqingsheng2015@163.com

赵玉新*(1980-),男(满族),吉林舒兰人,教授,研究方向:高/超声速气动设计、超声速边界层、超声速流动成像测量等.E-mail:zyx_nudt@163.com

蓝庆生, 赵玉新, 赵一龙, 等.三维超声速流动的压力反问题[J].空气动力学学报, 2017, 35(3): 429-435.

10.7638/kqdlxxb-2016.0156 LAN Q S, ZHAO Y X, ZHAO Y L, et al.Pressure inverse problem of three-dimensional supersonic flow[J].Acta Aerodynamica Sinica, 2017, 35(3): 429-435.

V434+.1

A doi: 10.7638/kqdlxxb-2016.0156