张涛 郑晓刚 汤祎麒 李怡庆 尤延铖

摘要：基于传统吻切理论，本文提出了一种非轴对称吻切技术，并在此基础上，完善了传统二维特征线技术，可对复杂三维曲面激波进行逆向求解。通过事先指定三维激波曲面形状，根据气流方向与激波曲面当地曲率方向，能够求出各离散激波点对应的离散微吻切平面。以该微吻切平面与激波曲面的交线作为各微吻切平面内的激波形状，利用已知激波求解流场的二维逆向特征线法对各微吻切平面的流场进行求解，随后将各离散微吻切平面内获得的压缩型线进行组合，获得能够生成指定复杂三维曲面激波的压缩型面。研究结果表明，基于非轴对称吻切技术的三维激波逆向求解方法能够很好地求解出生成指定三维激波曲面的乘波构型，激波形状吻合度较高；运用该逆向设计方法求解获得的截面流场信息与数值模拟结果相似度较高，其流场分布表现出相同的规律，但在增压比的预测上存在一定的误差，主要集中在远离对称面三维效应明显处，最大误差约为8.4%，可见该逆向乘波设计方法的精度基本满足要求。另外，针对特定的二次锥面激波方程，设计截面内激波曲线越高，乘波体的升阻比越低。

关键词：非轴对称吻切技术；特征线；乘波体；升阻比；逆向设计

中图分类号：V211.3文献标识码：ADOI：10.19452/j.issn1007-5453.2020.11.005

基金项目：国家自然科学基金（51276151，91441128）；基础科研项目（B1420133058）；中央高校基本科研业务费（20720140540）；江西省教育厅科技项目（GJJ190523）

近空间飞行器的设计和开发是目前国际航空航天领域的研究热点，同时也是各国竞相争夺空间技术的焦点之一[1-5]。尽管国内外大量学者就此类问题进行了详细的研究，但是仍然存在一些亟待解决的问题。其中飞行器的升阻力特性因其直接决定了飞行器的气动性能，并且通过传统的设计方法很难构造出具有高升阻力特性的飞行器机体，因此该问题长期以来受到国内外学者的广泛关注。区别于亚声速流动，超声速流动最显著的特征在于流场中存在激波，流场中的激波一方面给飞行器带来了额外的激波阻力，但另一方面，气流经过激波后能够实现减速增压的特点，从而为飞行器部件中需要高压的区域提供合适的气流，如飞行器机体下表面或吸气式推进系统的内流通道等。

利用超声速流动具有激波的优势，Nonweiler[6]首次提出运用二维斜激波波后流场设计飞行器机体，该机体前缘产生的斜激波能够完全附着于机体下表面，从而抑制下表面气流于机体边缘处上洗至上表面导致升力缺失，较好地解决了飞行器升阻力的问题。该设计方法被命名为楔导乘波理论，生成体即为楔导乘波体。运用类似的理论，乘波理论的激波形状由最初的斜激波[7]发展至圆锥激波、椭圆锥激波等三维回转激波生成体，基于该方法生成的乘波体即为锥导乘波体[8]，锥导乘波理论在很大程度上拓宽了乘波理论的应用范围，但由于激波形状较固定，使得该理论仍然存在着较强的约束。为解决以上问题，Sobieczky[9]等提出了吻切锥导乘波理论，该理论将激波的平面形状离散成激波微元段，每段微元均可根据当地曲率半径获得相应的圆锥激波，从而拓宽了乘波理论在激波形状方面的限制，使得任意曲率中心连续过渡的激波曲线均能作为生成乘波体所需的激波曲面。

经过长期的发展，乘波理论已经逐渐趋于成熟，近年来国内外学者逐渐将重心转移至乘波体气动性能的优化[10-12]、乘波體与推进系统的一体化[13-17]等问题。作者认为，虽然各类乘波理论经过长期的发展已经基本趋于完善，但在实际工程运用中仍然存在着许多约束，其中最重要也是最难解决的问题仍然是传统乘波设计理论，特别是吻切乘波理论在激波形状设计中的局限性。传统吻切乘波理论在设计过程中首先给定的是设计截面内的激波形状，通过截面激波形状可反向推导出该激波曲面的三维结构，因此，激波曲面的三维形状是无法预先设计的，而是根据激波的截面形状与激波生成体唯一确定。此外，在实际工程应用中，受到激波曲线的形状约束，设计前缘型线时通常会遇到较强的约束，这主要是因为吻切乘波体设计理论中，前缘型线必须位于激波曲线的曲率中心与曲线之间，若前缘型线超越激波曲线曲率中心则无法构造乘波体下表面。因此，发展一套基于任意复杂三维激波曲面的气动反设计方法从而进一步拓宽乘波理论的适用范围，是亟待解决的关键问题。

气动反设计的本质是预先给定激波反求激波生成体的过程，基于二维有旋特征线法，国内外学者开展了大量的二维气动反设计研究，包括给定二维激波形状反向求解压缩型面[18]，给定压力分布反求压缩型面[19]以及给定马赫数分布求解压缩型面[20]等。以上研究较全面地解决了二维环境下的反设计问题，但是对于三维复杂激波的求解目前公开文献较少，一些学者尝试使用三维特征线法对其进行逆向求解[21-22]，但由于特征线法的复杂性以及三维求解过程中特征线易相交等问题，未能很好地解决三维激波曲面逆向求解的问题。

基于以上分析，在传统吻切乘波体理论基础上，本文进一步发展了基于非轴对称吻切技术的乘波体设计方法。不同于常规吻切乘波理论，该方法将三维激波曲面在横向上进行离散的同时，在流向上也离散为若干微小平面，本文将其命名为微吻切平面。利用改良后的二维逆向特征线法对各微吻切平面的流场进行求解，进而组合各微吻切面内的型线，获得能够生成指定复杂三维曲面激波的压缩型面。本文针对不同的三维激波曲面，构建了三种乘波构型，对比其激波形状与流场信息，用于验证本文所述逆向求解技术的可行性和正确性。为了便于对比，在本研究中，乘波体前缘型线的二维形状均指定为直线。针对激波方程为二次锥面的情况，本文探讨了设计截面内激波曲线高度对乘波体升阻特性的影响。

1非轴对称吻切技术

1.1吻切流基本原理

吻切乘波理论由Sobieczky[9]等于1990年率先提出。该理论认为：一般三维超声速的运动方程都可以在二阶精度范围内用一个轴对称流的运动方程来逼近。这个轴对称的轴线位于通过该点流线的吻切平面（osculating plane， OP）内。于是，当地的三维流动就能够由局部的二维流动（轴对称流动也是二维的）来描述。在吻切面中，非轴对称的激波后流动处理为锥形流。在出口平面内，沿激波曲线使用一系列吻切面来定义，每个吻切面内激波角为常数，以保证展向的连续性，而每个面中的锥形流顶点则由激波角和当地的曲率半径确定。因此，此技术可将完整的激波形状分解成一系列半径不同的圆锥形激波“切片”，而激波后的流场则为一系列锥形流场的耦合流场。由此可见，锥导乘波理论只是吻切乘波理论的一个特例：锥形激波相当于吻切方法中半径固定为一恒定值的情况。图1给出了吻切乘波体的设计原理图，图中ICC为进气道捕获曲线，FCT代表前缘捕获型线。

1.2非轴对称吻切技术

直接对三维曲面激波进行逆向求解显然是极为复杂的，因此需要进行降维处理，将三维的求解过程简化为二维过程。参照吻切流理论，本文将三维曲面激波在横向上进行离散，进一步地在流向上也离散为若干微吻切平面，在各微吻切平面内对激波曲线进行二维的逆向求解。本文将其命名为非轴对称吻切技术，并认为一般的三维超声速流动可由这些微吻切平面内二维流动的耦合来近似。该技术的核心问题是如何离散三维曲面激波得到一系列微吻切平面，本文将采用步进法根据前一点所得数据逐步对流场进行离散和逆向求解。基于非轴对称吻切技术的三维激波逆向求解过程可分为4步（见图2）。

（1）设计激波曲面的形状，提取该激波曲面的方程

激波作为超声速流场具有的最显著的特征在很大程度上影响着飞行器及推进系统的气动性能，对于外流部件，保证激波封闭飞行器下表面抑制上洗流能够实现飞行器的高升阻力比特性。而以同样的方式将推进系统的进气道进口封闭能够有效提高进气道流量捕获能力，从而进一步增加发动机提供的推力。因此激波的形式与分布方式的研究是实现超声速飞行的关键。分析乘波理论的设计过程可以发现，无论是二元楔导乘波理论、锥导乘波理论，还是具有较高任意性的吻切锥导乘波理论，在设计之初均为预先指定激波的分布形式。然而三者在设计过程中均存在着局限性，即仅在某一截面（设计截面）内构造激波曲面的二维几何形状，随后假设流动沿流向为平面流动（楔导乘波）或标准圆锥流动（锥导、吻切乘波）进行求解，对于三类乘波理论，其激波的二维几何形状分别为直线形、圆形和具有二次连续特征的曲线形。因此，乘波理论对激波三维形状的考虑显然是不充分的，从另一个角度分析，乘波理论实际上放弃了激波曲面三维特征这一自由度。

对于本文的研究，激波曲面的设计不再局限于某一平面，而是对全流场所具有的激波形状进行设计。理论上任意满足设计要求的三维曲面方程均能作为预输入激波曲面。需要说明的是，激波曲面沿流向和横向的曲率方向，也就是曲率的正负决定着气流的压缩方式（轴对称外压缩或轴对称内压缩），针对以上两种压缩方式求解过程显然是不同的，由于求解过程的复杂性，本文暂不将曲率方向发生变化的激波曲面作为研究对象。因此，激波的设计过程需保证激波曲面各点的曲率方向相同。而曲面方程的确定对于求解特征点的曲率半径显然是必需的。

（2）设计物面前缘型线的三维构型

为保证所生成物面边缘不产生溢流以提供较高升阻比，在前缘型线的设计中需要保证前缘型线与激波曲面相交。因此在已获得激波曲面的前提下，前缘型线的三维构型将由其横向投影形状或流向投影形状唯一确定。若在横向截面内设计前缘型线则将其沿流向投影至三维激波曲面，若在流向截面内设计则将其沿横向投影至三维激波曲面，从而获得前缘型线的三维构型。随后将三维前缘型线离散与三维激波曲面共同作为逆向求解的输入条件。

（3）求解前缘离散点对应初始微吻切平面，运用逆向特征线法解得该微吻切平面的压缩型线

初始微吻切平面可根据来流方向与对应激波点的法矢量唯一确定。如图2（a）所示，本文将蓝色箭头所表示的气流方向矢量与An点对应法矢量（normal vector）构成的平面定义为初始微吻切平面，该微吻切平面与给定激波曲面相交能够获得如图2（a）中AnB所示的初始激波微元段。本文假设若An点与B点之间间距足够小，则能够将复杂的三维流动简化为二维轴对称流动。该假设在后续的数值模拟验证中将得到证实。获得以上条件后，运用基于激波曲线的逆向特征线法能够较容易地求解出所需的压缩型线，如图2（a）中AnD所示。

（4）求解后续微吻切平面及对应压缩型线，运用相同的方法获得全流场参数

通过步骤（3）获得了初始微吻切平面内的流场参数，为后续流场求解提供了可能。初始微吻切平面之后的微吻切平面不再将来流流动方向作为构造依据，而是由上一微吻切平面内最后一条左行马赫线（图2（b）中DB）与激波点B对应的法矢量唯一确定，考虑到后续微吻切平面的流动将主要受前部已扰动来流的影响，这是本文提出的非軸对称吻切技术对吻切乘波理论在流动方向上的主要发展。前文分析中指出吻切平面将三维流动在周向内离散成轴对称流动从而解决了三维流动的简化求解问题，本文进一步拓展，将复杂三维流动在周向和流向内同时离散从而实现对更加复杂的三维流动的简化求解。流向流场求解思路将在下一节沿流向曲率中心可变的逆向特征线法中进行详细的分析。

将以上方法同时运用于其他离散前缘点（如图2中An-1，An+1），并将离散流场进行三维组合后完成复杂三维流场的逆向设计。

2沿流向曲率中心可变的逆向特征线法

实现三维激波逆向求解的另一核心手段是逆向特征线方法（method of characteristic， MOC），运用该方法的逆置激波点过程能够在二维环境内有效求解出指定激波所需的压缩型面。传统的二维特征线法多为正向求解，以来流条件为初始参数，以压缩型线为边界条件，根据已知两点发出的异簇特征线，联立特征线方程和相容性方程求解特征线交点的参数，最终获得整个流场信息。而逆向特征线法的求解思路则与之相反，先给定来流参数和预想激波形状，再采用特征线法逆向求解压缩壁面型线的几何位置参数，将后验的参数变成可设计的参数。换言之，利用逆向二维特征线法可以根据所需的激波形状反设计得到能够产生该波系结构的压缩型面几何构型，但该逆向求解方法显然无法实现任意三维曲面激波的求解，为此本文发展了一种沿流向曲率中心可变的逆向特征线法，该方法将在下文进行详细研究。

特征线法是一种求解双曲型偏微分方程的精确步进型方法。在定常超声速流场中，由于其控制方程为双曲型偏微分方程，流场中任一点的流动具有仅取决于上游流场中有限区域的性质，因此可以使用特征线法求解该流场。所谓特征线是指沿着该曲线积分将偏微分方程简化为易于求解的全微分相容性方程。在超声速流场中由于马赫线就是特征线，因此可将控制方程简化为以下两个全微分方程组。特征线方程组：

通过联立特征线方程组和相容性方程组求解特征线交点的参数，最终可获得整个流场信息。

由于特征线理论具有较高的求解精度与较短的求解时间，长期以来被广泛运用于超声速领域的飞行器设计问题。此外，基于该理论，一系列逆向求解方法也得到了发展。其中包括给定壁面压力分布求解超声速流场、给定壁面马赫数分布求解超声速流场以及给定激波曲线求解超声速流场等问题。针对本文的研究内容，给定激波曲线逆向求解超声速流场将成为主要的实现手段。传统的给定激波反向求解轴对称超声速流场的求解原理如图3所示，可以发现，图中具有两条橙色曲线，其中平行于x轴的橙色虚线为轴对称流场的回转中心。其求解过程首先需要确定来流条件与预设计激波曲线形状（图3中红色实线所示）并获得激波起始点An，然后将该激波曲线离散为一系列微元段（如图3中AnB、BB），进而运用特征线理论的逆置激波点过程求解激波点B处的左行马赫线并与起始点An发射的流线相交得到壁面点D的初始值，随后使用校正迭代法对D点进行修正直至精度满足要求。同理可对下游点D进行求解进而求解整个超声速流场，在此不再赘述。

传统的给定激波求解流场的特征线程序能够很好地实现二维流动或准二维流动（轴对称流动）的精确求解，但该理论对于三维流动的求解显然是不适用或者不精确的，其本质的原因在于各激波离散点的曲率中心受到三维效应的影响发生了变化。图4为圆锥曲面与椭圆锥曲面沿程各激波点曲率中心分布对比，显然圆锥流动的求解可以完全依赖于二维特征线法，但该方法无法满足具有三维效应的椭圆锥流场的求解。无法求解的内在原因在于圆锥曲面与椭圆锥曲面沿流向曲率中心位置的不同。

对于圆锥曲面（包括母线为曲线的曲锥面），曲面上各点对应曲率中心均落于回转轴线上，即曲面的物理中心与曲率中心重合，因此流場在周向位置内存在相似性能够满足二维流动对轴对称流场的求解要求，然而对于椭圆锥曲面，曲面上各点的曲率中心不再落于某一轴线，而是由当地曲线曲率决定，因此不存在周向相似的特点，无法使用传统的特征线法进行求解。经过以上分析可以发现，两类流动的本质区别是曲率中心的问题，进一步而言是曲率半径的问题，因此若在特征线法的求解中引入曲率半径这一变量将有望解决特征线法解决三维流动的问题。

结合曲率半径的变化，本文进一步发展了二维特征线理论，得到如图5所示的曲率中心可变逆向特征线求解原理。求解过程首先需要获得离散激波点当地曲率半径，由于三维效应，该曲率中心将偏离回转中心（图5中X轴），且向下游偏离程度逐渐增大，故将呈现出图5中橙色实线所示的偏离趋势，该趋势在图4（b）中椭圆曲面的曲率中心分布中也有同样的体现。曲率半径的不同对应特征线求解公式中即为y值的不同，分析式（5）可以发现当y值趋于0时流动为标准的轴对称流动，而当y值趋于无穷时流动将向二维流动转变，因此当y值落于0与无穷之间时求解的流场为具有曲率半径为y的圆台的轴对称流场（即前缘起始点不为原点），该流动将介于二维流动与轴对称流动之间。

上述对y值分析仍然是针对曲率中心固定的流场求解，对于类似图4（b）所示的变曲率中心流场的求解则需要在各个微元段内多次运用特征线法的以上特点，以图5为例，求解激波微元段AnB流场时，以An点对应曲率半径作为y值带入特征线方程，而当求解BB时对应y值将采用B点的曲率半径，随后以相同的思路对下游流场进行求解获得能够满足类似图4（b）中所示的变曲率中心流场的逆向求解。上述曲率半径选择方式能够成立是因为对于超声速流动，气流参数仅受上游参数影响，因此当微元段足够小时以微元起始点的曲率半径作为微元段的曲率半径进行求解是可行的。此外，需要说明的是，图5中所示的求解过程对应于图4（b）中椭圆曲面短轴对应母线的求解，该母线沿流向曲率半径逐渐增大，曲率中心偏离物理中心程度也同时增大，因此流场将出现沿流向趋于二元流动的趋势。而对于长轴对应母线，观察图4（b）可以发现沿程各点曲率中心虽然偏离物理中心，但是其偏离方向与短轴对应母线相反，因此流场沿流向将更加趋于圆锥流动，即呈现出比“圆锥流动更圆锥”的流动特点。

因此，本文提出的基于非轴对称吻切技术的三维激波逆向求解方法的实质是指在利用非轴对称吻切技术对三维激波进行离散的每张微切面内，采用曲率中心可变的逆特征线法求解流场。

3复杂激波曲面逆向求解对比分析

本节将对上文所述的基于非轴对称吻切技术的三维激波逆向求解方法进行实例分析，并与数值模拟结果进行对比分析以验证该方法的正确性与可行性。需要特别指出的是，本文验证逆向求解方法正确性的步骤为：首先设计三维激波生成体，随后运用计算流体力学（CFD）对其进行求解，进而提取所具有的激波曲面进行对比验证。因为运用CFD求解获得的激波曲面存在一定厚度，故在提取激波曲面时将不可避免地存在误差。

构造三维激波曲面时，本文采用图6所示的一般锥面方程作为初始输入。如图6所示，过定点且与定曲线相交的所有直线构成的曲面称为锥面，其中定点称为顶点，定曲线为锥面的准线，构成锥面的每一条直线叫作母线，而准线所在的截面则称作设计截面即为乘波体的结尾截面。因此，给定顶点坐标与准线方程即可获得锥面方程。值得注意的是，椭圆锥实际上是顶点为坐标原点，准线方程为椭圆方程的特殊锥面。

根据上述锥面定义，本文选取了三种锥面方程，并将乘波体前缘型线的二维形状指定为直线，构建了三种不同的乘波构型，用以对比激波形状及流场信息。三种不同锥面方程及乘波体设计的具体参数见表1。其中模型A与模型B的激波曲面实际上为椭圆锥面，而模型C的激波曲面在设计截面中的方程是二次曲线，因此，该锥面又称为二次锥面。

在来流马赫数Ma=6，H0=27km的条件下，根据前缘捕获型线与预先给定的三维激波面方程，运用基于非轴对称吻切技术的逆向求解方法便可获得乘波体的下表面及其流场信息。

本文使用ANSYS Fluent求解基于有限体积法的欧拉方程来验证基于非轴对称吻切技术的三维激波逆向乘波设计方法的准确性。选择基于密度的求解器来进行无黏计算。对流矢量选择二阶AUSM格式，CFL数设定为0.5以确保计算的稳定性。假设来流气体为理想气体，定比热比γ为1.4。对于边界条件而言，将乘波体壁面均设置为绝热无滑移固体壁面。计算域入口选用压力远场边界条件，出口设置为压力出口边界。网格方面，使用ICEM对乘波体流场划分结构化网格，在乘波体壁面布置C型网格，在壁面和激波位置附近均进行网格加密处理，各算例网格数目均在350万左右。

图7为模型A乘波体结尾截面激波示意图。其中左半为通过数值模拟获得的乘波体结尾截面激波形状与预设计截面激波形状（图中黑色虚线所示）的对比图，右半为乘波体的俯视图。此外，图7左半上图表示基于非轴对称吻切技术，采用传统的特征线法求解得到的结果，左半下图表示特征线求解技术引入曲率中心变化获得的结果。对比可以发现，在引入曲率中心变化之前，数值模拟得到的激波形状与预设计激波形状存在较大差异。而考虑曲率中心的变化之后，能够很好地求解出生成指定椭圆锥激波曲面的乘波构型。该结论验证了沿流向曲率中心可变的逆向特征线法的正确性，并实现了非轴对称吻切技术在三维流场中的应用。

除横向截面激波形状，在流向截面內，本文所述逆向求解技术同样表现出了较高的精度。图8分别提取了模型A的4个流向截面（Y=0，Y=0.2，Y=0.4，Y=0.6）内的激波形状与设计结果进行对比，其中红色实线为预设计激波形状，黑色实线为逆向求解获得的压缩型线，蓝色实线为基于该压缩型面运用数值模拟方法获得的流向激波曲线。可以发现，蓝色激波曲线基本被预设计的红色激波曲线覆盖，仅在某些局部位置出现微小偏差，该现象说明在流向截面内逆向求解方法能够很好地复现预设计的激波形状，求解精度基本能够满足要求。

图9为模型B数值模拟结果与预设计激波形状对比。可以发现，逆向求解结果与数值模拟结果同样具有较高的相似度。模型A与模型B的激波曲面均为椭圆锥面，但不同于模型A，模型B的结尾截面激波形状所具有的曲率中心显然已经超越前缘捕获型线，即激波曲线与曲率中心位于前缘型线的同一侧。因此，运用已有的吻切乘波设计理论显然无法获得具有图9所示激波形状的乘波构型，这是因为乘波理论成立的客观条件是前缘捕获型线位于截面激波曲线与激波曲线的曲率中心之间。而非轴对称吻切技术能够实现的原因是其不仅在横向上对三维流场进行离散，在流向截面内同样对流场进行离散，此外，沿流向曲率半径变化的引入也为此类流场的求解提供了可能。

模型C的激波曲面在设计截面中的方程是二次曲线，其数值模拟结果与预设计激波形状对比结果如图10所示，与模型B类似，其结尾截面激波形状所具有的曲率中心同样已经超越了前缘捕获型线。可以看到，逆向求解结果与数值模拟结果仍然具有较高的相似度。以此为基础，本文在下文中将会详细探讨不同的二次锥面对乘波体升阻特性的影响。

图11进一步提取了模型A的Y=0，Y=0.2，Y=0.4三个流向截面内的压力分布信息，并与数值模拟结果进行对比分析。由图11可知，运用逆向求解方法获得的截面流场信息与数值模拟结果具有较高的相似度，各截面内的流场分布规律均与CFD结果表现出相同的趋势。但在增压比的预测上存在一定的误差，其中对称面内增压比基本相同，而在Y=0.4截面内的增压比误差最大为8.4%。这主要是因为基于非轴对称吻切技术的逆向求解方法实际上是将三维复杂流动简化为一系列微吻切平面内的二维流动进行求解的过程，假设所有流动均在这些微吻切平面构成的流面内进行。在Y=0.4截面附近，激波曲率变化剧烈，流场中的横向流动效果显著，导致部分气流并未沿微吻切平面流动，这在一定程度上将带入一定的误差。

由此可见，基于非轴对称吻切技术的三维激波逆向求解方法能够运用于三维流场的求解，且在激波形状求解上具有较高的精度，激波吻合度较高，而在流场信息预测上，逆向求解获得流场分布与CFD结果具有相同的规律，但在增压比预测上存在一定误差，误差主要集中在远离对称面的三维效应较为明显的部分。总体而言，本文提出的三维激波逆向求解方法能够以较小的误差完成复杂三维激波曲面的逆向求解。

4激波曲面对乘波体升阻特性影响

在上文模型C的研究基础上，本文着重探讨了不同二次锥面的激波形状对乘波体升阻特性的影响。本文选取了4种不同的二次锥面，其顶点坐标均为（0，0，0），设计截面为X=5.1；设计截面内的准线方程均为二次曲线：

鉴于B的Z坐标的绝对值实际表示设计截面内激波曲线距离压缩面的高度，故用H表示B的Z坐标的绝对值，根据H值由小到大分别将乘波构型命名为构型1～构型4。

圖12为4种乘波体在设计截面内的激波曲线及曲率中心分布图。可以发现，随着H的增大，激波曲线对应的曲率中心整体不断下移。当H=2.0000时，已经有部分曲率中心位于FCT和ICC之间，即激波曲线与曲率中心位于前缘型线的同一侧。此时现有的吻切乘波理论显然并不适用，而这正是本文提出的基于非轴对称吻切技术的三维激波逆向求解方法的优势之一，其能够求解更为复杂的激波曲面，拓宽了乘波理论的应用范围。

4种乘波体的几何模型如图13所示。乘波体的下表面由本文所述求解方法逆向求解指定激波所得，上表面则直接将前缘型线延来流方向直接拉伸至设计截面所得。4种乘波体的宽度相同，均为W=3.0000，长度L基本相同，但随着H的增大，L值也逐渐增加，其中构型4的值最大，为4.1344。乘波体高度方面，可以看到随着H值的增大，乘波体的高度也越大，整体的迎风面积也逐渐增加。

图14给出了4种乘波体数值模拟结果与预定激波形状对比图。可以看到，在设计截面上4种乘波体所产生的激波与预设激波形状具有较高的相似度。这说明本文提出的基于非轴对称吻切技术的三维激波逆向求解方法在激波预测上具有较高的精确度，这与前文的分析是一致的。此逆向求解方法可用于探讨不同高度的二次锥面激波对乘波体升阻特性的影响。

图15为4种乘波体升力系数与阻力系数随H值的变化曲线。其中，CL是升力系数，CD表示乘波体的阻力系数。左上坐标图内σ表示升阻力系数用各自H=1.8时的值无量纲化后的结果，用于对比升阻力系数随H值的变化速率。可以看到，升力系数与阻力系数均随着H值的增大而增大，阻力系数与H值成线性增长关系，增长速率基本不变，表现为阻力系数成直线分布；而升力系数则成非线性增长关系，升力系数的增长速率随着H值的增大而减缓，表现为升力系数成曲线分布。同时，从左上坐标图中可以发现，阻力系数随着H值的增长速率明显大于升力系数。

图16给出了4种乘波体升阻比和迎风面积随H值的变化曲线图。图中，L/D代表乘波体升阻比，Aw表示乘波体的迎风面积。参照图16，乘波体的升阻比随着H值的增加而减小，但并未成线性分布。这是因为4种乘波体的升力系数与阻力系数均随着H值而增长，但阻力系数的增速更快，因此导致整体升阻比随着H值而降低，且升阻比的降低速率逐渐减缓。同时可以发现，4种乘波体的迎风面积随着H值成线性增长，这与前文的分析结果是一致的（见图13）。

综上所述，针对二次锥面激波，随着H值的增加（即设计截面内激波曲线的高度增加），乘波体的迎风面积将成线性增加，导致其阻力特性同样成线性增长。同时，随着高度的增加，乘波体的升力系数也会增长，但升力系数的增长速率则逐渐减缓，且其增长速率低于阻力系数。最终导致其升阻比随着高度的增加而降低，且降低速率逐渐减缓。可见，针对生成指定二次锥面激波的乘波构型，迎风面积的变化规律在其升阻特性中占主导地位，迎风面积越大，整体的升阻比越低。

5结论

借鉴吻切乘波原理，本文提出了一种非轴对称吻切技术，并在此基础上，对传统的二维逆向特征线法进行修正，发展了一种沿流向曲率中心可变的逆向特征线法，将其与非轴对称吻切技术结合可用于逆向求解复杂的三维激波曲面。利用此逆向求解方法，本文设计了三种生成不同锥面激波的乘波体，并与CFD结果进行对比，用以验证此逆向求解方法的正确性与可行性。同时，针对二次锥面激波，本文着重探讨了设计截面内激波曲线高度对乘波体升阻特性的影响。研究结果表明：

（1）吻切乘波理论实际上是将三维流动在周向上离散成一系列等波强的二维轴对称流动，在吻切面内，所有的激波角是恒定的。而非轴对称吻切技术在此基础上，进一步对流向同样进行了离散，此时，每张微吻切面内的激波强度不再相同。因此，非轴对称吻切技术本质上是将复杂的三维流动简化为二维流动的过程，是对传统吻切乘波理论的进一步拓展。

（2）本文提出的沿流向曲率中心可变的逆向特征线法在求解的过程中引入了曲率半径这一变量，以此带入激波曲面的三维效应，可用于求解具有三维效应的基准流场。因此，本文提出的基于非轴对称吻切技术的三维激波逆向求解方法的实质是在微吻切平面内运用曲率中心可变的逆特征线法。

（3）对照模型A、模型B和模型C可以发现，本文提出的三维激波逆向求解方法在激波形状求解上具有较高的精度，激波吻合度较高，而在流场信息预测上，逆向求解获得流场分布与CFD结果具有相同的规律，但在增压比预测上存在一定误差，误差主要集中在远离对称面的三维效应较为明显的部分，其能够以较小的误差完成对复杂三维激波曲面的逆向求解。同时，由于本文提出的方法相比传统的吻切乘波理论具有更高的自由度，其计算效率不可避免地有所降低。

（4）针对二次锥面激波，设计截面内激波曲线的高度越高，乘波体的升力系数与阻力系数均增大，且阻力系数的增长速度要快于升力系数，因而乘波体的升阻比呈下降趋势。同时，可以看到，乘波体的迎风面积变化规律在升阻特性中占主导地位，迎风面积越大，整体的升阻比越低。

参考文献

[1]Anderson B H. Design of supersonic inlets by a computer program incorporating the method of characteristics[R]. NASA TN D-4960，1969.

[2]Alexander K，Alexey K. Atmospheric cruise flight challenges for hypersonic vehicles under the ajax concept[J]. Journal of Propulsion and Power，2008，24（6）：128-131.

[3]劉济民，沈伋，常斌，等.乘波体设计方法研究进展[J].航空科学技术， 2018， 29（4）：1-8. Liu Jimin， Shen Ji， Chang Bin， et al. Review on the design methodologyofwaverider[J].AeronauticalScience& Technology， 2018， 29（4）：1-8.（in Chinese）

[4]王骥飞，蔡晋生.形面渐变内收缩进气道设计方法研究[J].航空科学技术， 2017， 28（1）： 30-35. Wang Jifei， Cai Jinsheng. Research on design of shape morphing inward turning inlet[J]. Aeronautical Science & Technology， 2017， 28（1）： 30-35.（in Chinese）

[5]尤延铖，梁德旺，郭荣伟.高超声速三维内收缩式进气道/乘波前体一体化设计研究评述[J].力学进展，2009，39（5）：513-525. You Yancheng， Liang Dewang， Guo Rongwei. Overview of the integration of three-dimensional inward turning hypersonic inlet and waverider forebody [J]. Advances in Mechanics，2009，39（5）：513-525.（in Chinese）

[6]Nonweiler T R F. Aerodynamic problems of manned space vechile[J]. Journal of the Royal Aeronautical Society，1959，63（585）：34-40.

[7]Naruhisa T，Mark L. Wedge-cone waverider con-figuration for engine-airframe interaction[J]. Journal of Aircraft，1995，32：1142-1144.

[8]Naruhisa T，Mark L，Mary L. Waverider configu-ration development for the dual fuel vehicle[C]// Space Plane and Hypersonic Systems and Technology Conference，1996.

[9]Sobieczky H，Dougherty F C，Jones K D.Hypersonic waverider design from given shock waves[C]// Proceedings of the First International Hypersonic Waverider Symposium，1990.

[10]Rodi P E. Waverider vehicle optimization with volumetric constraints for sonic boom[C]// AIAA Aerospace Sciences Meeting，2018.

[11]耿永兵，刘宏，姚文秀，等.锥形流乘波体优化设计研究[J].航空学报， 2006， 27（1）：23-28. Geng Yongbing， Liu Hong， Yao Wenxiu， et al. Viscous optimized design of waverider derived from cone flow [J]. Acta Aeronautica etAstronautica Sinica， 2006， 27（1）：23-28.（in Chinese）

[12]汤波，高云峰，李俊峰，等.基于CFD计算和遗传算法的乘波体优化[J].动力学与控制学报， 2008， 6（3）：270-274. Tang Bo， Gao Yunfeng， Li Junfeng， et al. Optimization of waverider configurations by CFD calculation and genetic algorithm[J]. Journal of Dynamics and Control， 2008， 6（3）：270-274.（in Chinese）

[13]You Y，Zhu C，Guo J. Dual waverider concept for the integration of hypersonic inward-turning inlet and airframe forebody[C]//16th AIAA/DLR/DGLR International Space Planes and Hypersonic Systems and Technologies Conference，2009.

[14]Li Y，You Y，Integration methodology for waverider-derived hypersonicinletandvehicleforebody[C]//19thAIAA International Space Planes and Hypersonic Systems and Technologies Conference，2014.

[15]贺旭照，秦思，周正，等.一种乘波前体进气道的一体化设计及性能分析[J].航空动力学报， 2013，28（6）： 1270-1276. He Xuzhao， Qin Si， Zhou Zheng， et al. Integrated design and performance analysis of waverider forebody and inlet[J]. Journal ofAerospace Power， 2013，28（6）： 1270-1276.（in Chinese）

[16]范晓樯，李桦，易仕和，等.侧压式进气道与飞行器机体氣动一体化设计及试验[J].推进技术， 2004， 25（6）： 499-502. Fan Xiaoqiang， Li Hua， Yi Shihe， et al. Experiment of aerodynamic performance for hypersonicvehicle integrated with sidewall compression inlet[J]. Journal of Propulsion Technology， 2004， 25（6）： 499-502.（in Chinese）

[17]易军，肖洪，商旭升.两种高超声速一体化构型的气动性能对比分析[J].航空工程进展，2011，2（3）： 305-311. YiJun，XiaoHong，ShangXusheng.Aerodynamic performanceresearchoftwointegratedhypersonic configurations[J]. Advances in Aeronautical Science and Engineering，2011，2（3）： 305-311.（in Chinese）

[18]贺旭照，倪鸿礼.密切曲面锥乘波体——设计方法与性能分析[J].力学学报， 2011， 43（6）：1077-1082. He Xuzhao， Ni Hongli. Osculating curved cone （OCC） waverider： design methods and performance analysis[J]. Chinese Journal of Theoretical and Applied Mechanics， 2011， 43（6）：1077-1082.（in Chinese）

[19]李怡庆，韩伟强，尤延铖，等.压力分布可控的高超声速进气道/前体一体化乘波设计[J].航空学报， 2016， 37（9）：2711-2720. Li Yiqing， Han Weiqiang， You Yancheng， 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）

[20]钟启涛，张堃元，庞丽娜.出口马赫数分布可控的二维高超进气道双重反设计[J].推进技术， 2015， 36（7）：982-988. Zhong Qitao， Zhang Kunyuan， Pang Lina. Double inverse design of 2D hypersonic inlet with controllable exit Mach number distribution[J]. Journal of Propulsion Technology， 2015， 36（7）：982-988.（in Chinese）

[21]蓝庆生，赵玉新，赵一龙.降维的三维特征线方法与应用[C]//中国航天空天动力联合会议，2017. Lan Qingsheng， Zhao Yuxin， Zhao Yilong. Method and application of 3D characteristic line of dimension reduction[C]// Joint Conference ofAerospace Propulsion， 2017.（in Chinese）

[22]蓝庆生，赵玉新，赵一龙，等.三维超声速流动的压力反问题[J].空气动力学学报， 2017， 35（3）：429-435. Lan Qingsheng， Zhao Yuxin， Zhao Yilong， et al. Pressure inverse problem of three-dimensional supersonic flow[J]. Acta Aerodynamica Sinica， 2017， 35（3）：429-435.（in Chinese）

[23]Center K B，Sobieczky H，Dougherty F C. Interactive design of hypersonic waverider geometries[C]// AIAA 22nd Fluid Dynamics，Plasma Dynamics and Lasers Conference，1991.

[24]Ding F，Shen C B，Liu J，et al. Influence of surface pressure distribution of basic flow field on shape and performance of waverider[J].ActaAstronautica，2015，108：62-78.

（责任编辑王昕）

作者简介

张涛（1997-）男，硕士研究生。主要研究方向：高超声速空气动力学、计算流体力学。

Tel：15256555449E-mail：1457745879@qq.com

郑晓刚（1994-）男，博士研究生。主要研究方向：高超声速推进技术研究、内外流流体力学。

Tel：18959216039E-mail：xiaogangzheng@stu.xmu.edu.cn汤祎麒（1998-）男，硕士研究生。主要研究方向：高超声速气体动力学、弯曲激波理论、计算流体力学。

Tel：18850013602E-mail：tangyiqi2019@163.com

李怡慶（1989-）男，博士，助理教授。主要研究方向：高超声速推进技术研究、内外流流体力学、计算流体力学。

Tel：13306019011E-mail：yiqingli@nchu.edu.cn

尤延铖（1981-）男，博士，教授。主要研究方向：高超声速空气动力学、内流流体力学、高超声速进气道设计、复杂湍流数值模拟（LES/DES）和CFD计算数值方法研究等。

Tel：18060979961E-mail：yancheng.you@xmu.edu.cn

Inverse Waverider Design for 3D Shock Wave Based on Non-axisymmetric Osculating Cones Method

Zhang Tao1，Zheng Xiaogang1，Tang Yiqi1，Li Yiqing2，You Yancheng1，*

1. Xiamen University，Xiamen 361102，China

2. Nanchang Hangkong University，Nanchang 330063，China

Abstract： Based on the principal of traditional osculating cones method， a new non-axisymmetric osculating cones method was proposed. On this basis， the traditional 2D MOC was improved， which can be used to solve the 3D complicate shock wave. By specifying the appearance of three-dimensional shock wave in advance， the local micro osculating plane corresponding to the discrete shock points can be obtained according to the direction of the air flow and the local curvature direction of the shock surface. The improved two-dimensional inverse MOC is used to solve the flow field in each micro osculating plane， while the intersection line between the micro osculating plane and the shock surface is defined as the main input. Subsequently， the compression surface， designed to generate specified 3D shock wave， are formed with all the compression lines gained from each micro osculating plane. The results show that the inverse technology for solving 3D shock wave via non-axisymmetric osculating cones method shows high agreement in the shape of shock wave， and the flow field information obtained by this inverse technique is quite similar to the numerical simulation results. However， there is still a certain deviation in the prediction of the pressure ratio， while the maximum deviation is about 8.4%. In addition， for the specific quadric conical shock wave， the higher the shock curve is in the design section， the lower lift-drag of waveriders is.

Key Words： non-axisymmetric osculating cones method； method of characteristic （MOC）；waverider； lift-drag ratio； inverse design