现代中远途民用客机和高空高速长航时无人机巡航状态的阻力主要由摩擦阻力和升致阻力构成，分别占总阻力的55%和35%左右^[1]。飞机的总阻力每减小1%，直接使用成本可以降低0.2%或者增加1.6 t的有效载荷^[2]。根据Bregust公式分析得出，在飞机航程不变的情况下，150座民机每降低一个阻力点(ΔC_d=0.0001)就可以多承载3或4名乘客^[3]。目前，在跨声速状态下，飞机表面层流化是能够大幅度减小摩擦阻力的有效方法。研究表明：飞机机翼表面层流区域从10%扩大到90%弦长时，航程可以增加近50%；或者减少一半以上的起飞重量^[3]。超临界翼型的发展，使得飞机在跨声速状态下能够实现较高的升阻比，所以跨声速大展弦比飞机都采用中等后掠角的超临界机翼。然而在飞行雷诺数约10⁷条件下，后掠翼极易诱发前缘附着线转捩及边界层展向流的不稳定性，从而使机翼过早地发生转捩进入到湍流状态，大大增大了飞机表面的摩擦阻力^[4]。大展弦比超临界机翼的设计是高空高速长航时无人机的需要，拥有成熟的跨声速自然层流翼型/机翼设计技术，将会极大地延长未来高空长航时无人机的留空时间，使其具备持久的情报收集和战场监视能力；而对于民航客机而言，层流化技术不仅可以减小CO₂的排放量，亦可提高运营效率^[5]。

自然层流(Natural laminar flow, NLF)翼型设计的关键就是设计翼型使其表面保持大范围的顺压梯度。然而在跨声速状态下，层流翼型表面流向大范围顺压梯度的存在，使得其后缘因压力恢复产生较强的激波，在减小摩擦阻力的同时又增加了激波阻力，反而削弱了层流化的效果，因此有必要对翼型/机翼后缘处的激波进行控制。在过去的十几年中，多种不同的方法被应用于跨声速状态下的流动控制。为了控制翼型/机翼表面的强激波，提出了许多流动控制方法，例如吹吸气技术、涡流发生器、震荡片等技术。近年来，安装于跨声速机翼激波位置处控制激波强度的鼓包技术得到大量的研究和发展，因未极大地增加结构重量及结构复杂性，展现出其减小总阻力的潜力。二维激波控制鼓包概念最初是由Ashill等人于1992提出的^[6]。为了提高二维鼓包技术的鲁棒性，Wong等人^[7-8]提出了三维激波控制鼓包的概念。Birkemeye r等人^[9-10]研究鼓包耦合边界层吹吸气技术对激波控制的效率。Ogawa等人^[11]应用楔形鼓包、圆角楔形鼓包、等轮廓鼓包等不同构型的鼓包来研究不同外形对激波控制的效果。后缘装置技术(Trailing edge device, TED)是通过增加后缘后加载减小上翼面压力恢复，进而减小流动穿过激波的总压损失来控制机翼表面激波强度的。后缘装置是安装于翼型/机翼后缘处的机械控制装置，其尺寸小于传统的控制面，通常只有4%或者更小的弦长比例^[12]。后缘装置包含多种控制外形，其中Gurney后缘、分裂式后缘、离散式后缘是较常见的几种装置^[13]。Gurney后缘是由Liebeck^[14]在1978年第一次提出来的。离散式后缘是基于1985年Boyd^[15]的专利提出并应用于流动控制之中。Henne等人^[16]通过实验手段和数值模拟在跨声速状态进行了带有分离式后缘装置的跨声速翼型的稳定性设计。Li等人^[17]在超临界翼型上开展了Gurney后缘装置的流动控制研究工作。

合理的鼓包外形和安装位置可以弱化激波、减小激波阻力，然而应用于层流翼型设计中后缘处激波的控制时，由于鼓包安装于激波位置处，激波处的流动发生改变使得流动提前发生转捩，从而减小翼型表面层流区域。另外鼓包的流动控制鲁棒性较差，当飞行状态偏离设计点时，鼓包控制激波的效果将减弱，严重时会产生两道激波增加总阻力^{[8, 12]}。扩展到三维层流机翼优化设计中后缘激波控制时，由于一系列离散鼓包的存在，需要在鼓包位置处的展向和弦向布置足够数量的网格来准确模拟其对流场的影响，将大大增加优化设计中的流场分析时间，从而延长三维层流机翼设计周期。因此，本文开展后缘装置对激波阻力控制的研究，探讨后缘装置不同高度、长度下对激波控制的规律，最后应用于层流翼型优化设计中后缘处强激波的控制。

文中使用基于线性稳定性理论的eN方法作为转捩判断的工具，采用后缘装置对层流翼型后缘处激波进行控制，基于此，应用多岛并行进化算法耦合均衡Pareto策略在推迟翼型表面转捩位置降低摩擦阻力的同时优化后缘装置(高度、长度)对翼型进行多目标优化设计。

1 流动计算和层流优化方法 1.1 流动计算和转捩预测

本文求解基于有限体积法的控制方程对流场进行数值模拟，见式(1)。

$ \frac{\partial }{{\partial t}}{\smallint _\mathit{\Omega} }\mathit{\boldsymbol{U}}{\rm{d}}\mathit{\Omega} + \oint_{\partial \mathit{\Omega} } {({\mathit{\boldsymbol{f}}_{\rm{c}}} - {\mathit{\boldsymbol{f}}_{\rm{v}}}){\rm{d}}s} = 0 $

(1)

式中U, f_c, f_v分别代表守恒相矢量、对流通量矢量和黏性项通量矢量。

采用有限体积格式对流场网格进行离散，S-A湍流模型用来模拟边界层内湍流流动。Roe通量差分分裂格式作为空间离散格式用来计算无黏通量相，黏性通量的离散采用二阶中心差分格式，时间离散采用四步Runge-Kutta推进格式，远场边界条件采用黎曼不变量。由于流动高速可压，将气体视为理想气体，黏性系数由Sutherla nd公式计算得到。计算时通过适当减小松弛因子来提高数值计算的稳定性，隐式残差光顺和当地时间步长方法用来加速流动收敛。网格采用C形网格，上下翼面均布置150个网格点，尾迹区和法向均布置100层网格。

转捩位置的预测一直是流体力学的一个难题，根据目前的理论，边界层的转捩方式可以划分为以下3种：(1)“模态转捩”。通常出现在背景扰动的湍流度比较低的情况下。其过程通常为，边界层对于外部初始小扰动而诱发形成T-S波过程，随后是经历了较长的扰动线性增长，最后进入非线性演化阶段，当流场中高频广谱扰动形成后，流动最终达到湍流。(2)“旁路转捩”。Klebanoff等^[18]通过边界层转捩实验研究观察到当初始扰动足够大，扰动的线性阶段将被直接跳过，边界层扰动出现突变式的增长，这种转捩被称为旁路转捩。(3)Schmid-Henningson^[19]发现边界层转捩中存在另一类重要的转捩现象即“斜波转捩”，运用直接数值模拟方法首次研究了槽道流动中一对三维扰动模态诱发转捩的时间发展过程。

线性稳定性理论(eN方法)是基于小扰动理论建立起来的转捩判断方法，通过数值模拟流场中的微小扰动随时间在流场中的演变，观察这个扰动随着时间的推移是放大还是衰减来判断流动是否发生转捩。如果扰动衰减，就认为流动保持层流状态；如果扰动振幅随着流动不断增加，则流动就会发生转捩并发展成湍流^[20-21]。线性稳定性方程是从Navier-Stokes方程推导得出的，将方程中独立的物理量分解为平均量和扰动量之和${q_i} = {\bar q_i} + {\tilde q_i} $，并进行线化处理后得到可压缩线性扰动方程组^[21-22]。定义u，v，p，T和ρ为流场中瞬时速度、瞬时压力、瞬时温度和瞬时密度项等，包含流动平均相和波动相，其表达式为

$ \begin{array}{l} u = \bar u + u\prime , {\rm{ }}v = \bar v + v\prime , \\ p = \bar p + p\prime , {\rm{ }}T = \bar T + T\prime , {\rm{ }}\rho = \bar \rho + \rho \prime \\ k = \bar k + {\rm{ }}k\prime , {\rm{ }}\mu = \bar \mu + \mu \prime , {\rm{ }}\lambda = \bar \lambda + \lambda \prime \end{array} $

将瞬时变量代入到NS方程中，由于波动相u′，v′，p′，T′，μ′，k′和λ′均为小量，因此其平方项和乘积可以忽略。同时平均速度项、平均压力项和平均温度项也满足NS方程，简化后得到小扰动方程组。基于平行流近似理论：平均速度($ {\bar u}$)、平均密度(${\bar \rho } $)和平均温度(${\bar T} $)仅为y的函数。即

$ \bar u = \bar u\left( y \right), {\rm{ }}\bar \rho = \bar \rho {\rm{ }}\left( y \right), {\rm{ }}\bar T = \bar T\left( y \right), {\rm{ }}\bar v = 0 $

方程进一步得到简化。定义

$ R = \frac{{{u_{\rm{e}}}L}}{{{\nu _{\rm{e}}}}}, {\rm{ }}{M_{\rm{e}}} = \frac{{{u_{\rm{e}}}}}{{\sqrt {\gamma \bar R{T_{\rm{e}}}} }}, {\rm{ }}Pr = \frac{{{c_{\rm{p}}}\mu }}{k} $

代入方程后得到小扰动方程的量纲一形式。应用稳定性理论，即变量分析方法对以上方程组进行简化。设小扰动波(方程组的解)为一系列正弦波，则

$ \begin{array}{l} q\prime \left( {x, y, t} \right) = \hat q\left( y \right){e^{{\rm{i}}\left( {\alpha x - \omega t} \right)}}\\ \hat u = \frac{{\partial q\prime }}{{{\rm{d}}y}}, {\rm{ }}\hat v = - \frac{{\partial q\prime }}{{{\rm{d}}x}} \end{array} $

代入到小扰动方程组，整理后得到稳定性方程。

求解RANS方程得到物面压力分布，进而得到边界层外边界压力分布(边界层内$ \partial p/\partial n = 0$)。基于外边界压力分布求解可压缩边界层方程得到边界层内速度型、动量厚度等参数。然后把边界层的解作为稳定性方程的输入，通过求解稳定性方程得到空间扰动放大率。进一步通过eN方法求解得到不同流场位置处的扰动放大因子值，当其值达到给定的阈值时判定为流动转捩点。重复上述过程至转捩点收敛，见图 1。

为了验证本文转捩预测方法的准确性，选用NLF416翼型和NACA0012翼型作为算例进行验证。表 1是NLF416翼型在来流马赫数Ma_∞=0.69、迎角α=1.21°和雷诺数Re=1.2×10⁷状态下，通过eN方法计算的转捩位置与实验结果对比^[23]，表中x_upper和x_lower分别为翼型上下表面转捩位置。表 2为NACA0012翼型在来流马赫数Ma_∞=0.4、迎角α=2°和雷诺数Re=8×10⁶的状态下，eN方法计算的转捩位置和文献[24]中的结果对比。通过分析表中数据得出本方法计算得到的翼型上下表面转捩位置均接近于实验值/文献中参考值，误差在可接受的范围之内，可以作为翼型初步设计中转捩判断的依据。

1.2 自然层流翼型设计方法

层流化的目的是为了减阻，但是减阻设计却未必能实现层流化。为此目的，分别以总阻力最小化和转捩位置最大化为优化目标进行翼型优化设计。选用RAE2822翼型为初始翼型，在来流马赫数Ma_∞=0.729，迎角α=2.31°，雷诺数分别为Re=1.5×10⁷和Re=1.28×10⁷状态下，基于单目标进化算法，分别以总阻力最小化和翼面转捩位置最大化为优化目标进行优化设计，具体过程参考文献[20]。

对比表 3, 4中RAE2822翼型、总阻力最小为优化目标得到的翼型(MinC_d)和以层流区域最大化为优化目标得到的翼型(NLF)的气动数据可以得出：以总阻力最小为优化目标得到的翼型的阻力明显减小，但是转捩位置没有推迟，维持在60%弦长(上翼面10. 9%C和下翼面48.7%C)左右。所以以总阻力最小化为优化目标不一定能实现翼型层流化。另一方面，由于摩擦阻力很难通过CFD精确计算，进一步解释了以总阻力最小为优化目标不是层流翼型设计的有效方法。相比而言，以转捩位置最大化为优化目标得到翼型的层流区域明显增加，上翼面由原来的21.02%弦长增加到优化后翼型的52.58%弦长。通过分析以上两个算例得出，应以转捩位置作为目标进行层流翼型优化设计。

分析得出在跨声速自然层流翼型的设计过程中，激波阻力会相应的增加。从表 4中层流化得到的翼型和初始翼型的气动数据可以发现，优化后翼型激波强度较初始翼型大大增加，压致阻力系数从原来的0.008 095增加到优化后翼型的0.009 403，增加了16.16%。图 2为4组翼型压力曲线和转捩位置对比图(Ite B，层流翼型优化设计中，第二代最优个体；Ite C, 第三代最优个体)，可以看出随着转捩位置的推移，翼型后缘恢复压差增加，产生更强的激波阻力，所以在设计自然层流翼型的同时也要对后缘激波强度进行控制。

2 激波控制技术 2.1 激波控制鼓包技术

激波控制鼓包技术研究的主要目的就是通过弱化激波来减小波阻力，即改变激波根部的流动结构，使强激波转变成弱的激波^[20]，见图 3。然而鼓包技术的流动控制鲁棒性较差，当飞行状态偏离设计点时，激波控制的效果将减弱，严重时会产生两道激波增加总阻力，见图 4。另外当扩展到三维层流机翼优化设计后缘激波控制时，由于一系列离散鼓包的存在，需要在鼓包位置处的展向和弦向布置足够数量的网格来准确模拟其对流场的影响，将大大增加优化设计中流场分析的时间。从表 5中看出由于网格量的增加使得单个个体的流场计算时间大大增加，进而延长三维层流机翼设计周期。因此，本文将采用另外一种方法控制翼型后缘处的激波强度。

2.2 后缘装置控制激波技术

后缘装置是一系列用来控制翼型/机翼后缘流动的机械装置。在跨声速状态下，流动从翼型前缘开始，穿过一系列膨胀波之后加速到超声速状态，而后保持超声速流动直到后缘附近。由于压力需要恢复到后驻点处的压力值，同时流动由超声速状态变成亚声速状态，因此在后缘附近产生一道激波。而安装后缘装置后使得压力不需要恢复到后驻点处的压力值，逆压梯度小于初始翼型的逆压梯度，激波强度得到抑制，波阻力小于初始翼型的波阻力。

本节开展后缘装置控制翼型后缘激波强度的数值研究。研究后缘装置不同高度和长度对后缘处激波控制效果的规律。长度分别选取2.3%C和2.5%C(C为翼型弦长)，高度从0.2%C逐渐增加到0.8%C，见图 5。计算状态为自由来流马赫数Ma_∞=0.76、雷诺数Re=7.6×10⁶。图 6为初始翼型加装不同尺寸后缘装置后压致阻力变化趋势。从图中可以发现，在给定的长度下，后缘装置的减阻效果随着后缘高度的增加开始增加，在大约0.54%弦长高度时达到最优的减阻效果，压致阻力由初始的0.011 76减小到0.011 4 5(2.3%C)和0.001 47(2.5%C)左右，见表 6。而后随着高度的增加减阻效果开始减小。其中

$ - \Delta {C_{{\rm{Dpressure}}}} = \frac{{{C_{{\rm{Dpressure(Original)}}}} - {C_{{\rm{Dpressure}}}}}}{{{C_{{\rm{Dpressure(Original)}}}}}} $

图 7为OS01和OS02两组加装后缘装置翼型和初始翼型的压力系数、转捩位置对比。数值计算结果表明加装后缘装置后，翼型后缘后加载增加，上翼面后缘处逆压梯度减小，从而弱化因压力恢复而产生的强激波，翼型的气动性能得到改善。在升力系数保持不变的情况下，压致阻力系数分别减小了2.58%和2.46%，见表 6，因此后缘装置具有良好的控制激波效果。另外从图 7中发现，后缘装置的存在不但没有减小翼面层流区域，相反增加了翼面上的层流区域。因此，后缘装置可以有效地应用于层流翼型设计中后缘处的激波控制。

3 结合后缘装置的自然层流翼型多目标优化设计 3.1 多目标优化问题及设计变量定义

在设计带激波控制的跨声速自然层流翼型的过程中，层流优化的结果是转捩位置的推迟，但是激波阻力也随之增加，所以总阻力不一定有效地减小。因此在层流翼型设计过程中，需要同时对激波阻力进行控制。层流翼型设计过程中需要解决的问题有：

(1) 翼型表面实现大面积的层流区域，即推迟流动转捩的发生。

(2) 减小后缘由于压力恢复而产生的激波。

本文选用RAE2822翼型为初始翼型，以层流区域最大化和激波阻力最小化为优化目标，应用基于Pareto策略的多目标进化算法，在马赫数Ma_∞=0.729，迎角α=2.31°，雷诺数Re=1.28×10⁷的状态下进行优化设计，考虑升力约束后的优化问题如下

$ \left\{ \begin{array}{l} \mathop {{\rm{max}}}\limits_{\left( {{\rm{ }}X, Y} \right)} {J_1} = ({x_{{\rm{upper}}}} + {x_{{\rm{lower}}}})(1 + \beta ({C_{\rm{L}}} - {C_{{\rm{L0}}}})/{C_{{\rm{L0}}}})\\ \mathop {{\rm{min}}}\limits_{\left( {X, Y} \right)} {J_2} = {C_{{\rm{Dpressure}}}}(1 - \beta ({C_{\rm{L}}} - {C_{{\rm{L0}}}})/C{{\rm{}}_{{\rm{L0}}}})\\ \;\;\;{\rm{subject}}\;\;\;\;{\rm{ to}}_{\beta = 1\;\;\;\;\;{C_{\rm{L}}} < {C_{{{\rm{L}}_{\rm{0}}}}}}^{\beta = 0\;\;\;\;{C_{\rm{L}}} \ge {C_{{{\rm{L}}_{\rm{0}}}}}}\\ X = ({x_{C1}}, {x_{C2}}, {x_{{\rm{length}}}}, {x_{{\rm{height}}}})\\ Y = ({y_{{\rm{1, up}}}}, \ldots , {y_{{\rm{7, up}}}};{y_{{\rm{1, low}}}}, \ldots , {y_{{\rm{7, low}}}}) \end{array} \right. $

在优化过程中，翼型的形状由初始翼型叠加Bezier曲线来控制，即Y=Y_Ini+ΔY。其中Y_Ini为初始翼型坐标值，ΔY为Bezier函数控制的翼型纵坐标变化量。翼型控制变量和后缘装置控制变量取值空间见表 7，8。

3.2 多目标优化方法

多目标遗传优化方法是模拟达尔文的生物进化论中适者生存和近代生物遗传学机理的计算模型，属于全局搜索最优解的一种方法^[25]。通常多目标遗传算法有几种常见的处理方式：对不同的目标给定权因子后求和，再对求和的值进行优化；非支配多目标遗传算法^[26]；采用Nash均衡的遗传算法等^[27]。采用加权目标函数求解多目标问题的方法具有诸多缺点，例如，丢失信息或者如何定义目标函数加权值的优先级。1985年，Schaffer^[28]第一次提出通过矢量进化遗传算法(VEGA)求解多目标问题的基因算法。1989年，Goldberg^[29]提出了另一种基于Pareto排序和共享的解，使得多目标问题的解散布在整个Pareto阵面上。而后合作Pareto策略耦合进化算法则成为多目标优化问题求解的标准方法。

为了加速收敛，本文采用多岛并行非支配排序遗传算法(NSGA-Ⅱ)进行自然层流翼型多目标优化设计^[30]。由于Master和Slave之间的信息传递只是染色体和相应适应值的传递，因此并行多目标进化算法的流程亦是明确的。每一代种群数为150，平均分布到5个独立的子岛上，子岛之间每4代进行信息的传递。子岛每一个个体的适应值通过并行CFD求解器求解。其中，交叉概率和突变概率分别为0.8和0.1，选择算子基于锦标赛算法。图 8为多岛并行多目标进化算法流程图。

3.3 优化结果及分析

图 9为多目标进化过程历程图，经过80代进化得到该多目标问题的均衡解。图 10为该多目标问题的Pareto阵面，其上一系列的点给出了翼型气动性能的变化趋势。阵面上随着翼型转捩位置的推迟，翼型的激波阻力也随之增强；层流区域最大的翼型，其激波强度也最强。阵面上的翼型表面流动及激波强度相比初始翼型均有所改善，即翼面上层流区域增加的同时，激波强度有效地得到控制。在跨声速状态下，由于超临界翼型上表面中段比较平坦，激波不易稳定存在于此区域，转捩有可能发生在失稳点附近或者激波位置处，所以Pareto阵面出现不连续的情况。

表 9为初始翼型和阵面上PM A (Pareto Member A)，PM B和PM C 3组翼型的气动数据对比，图 11为初始翼型和选取的3组翼型的压力分布和转捩位置对比。通过对以上翼型的气动性能进行分析得出，Pareto阵面上翼型的转捩位置及激波强度明显优于初始翼型，阵面上左下缘对应翼型的激波强度较初始翼型明显减小，如PM A翼型，其压致阻力较初始翼型减小约20.46%，而翼面上层流区域保持不变；右上缘对应翼型(PM C)层流区域较初始翼型明显增大，上翼面由初始翼型的21.02%增加到优化后翼型的55.61%弦长，然而其压致阻力较初始翼型并没有增加，反而减小了11.35%。因此数值计算结果表明，优化后得到的一些列翼型的气动性能明显得到改善。

另一方面，通过图 12初始翼型、加装鼓包翼型和阵面上PM B, C两组翼型不同马赫数下的阻力曲线图发现，利用鼓包来控制翼型后缘的激波阻力，虽然在设计点可以得到理想的减阻效果，但是当飞行偏离设计状态时，其总阻力明显增加，将极大地影响飞机的飞行性能。然而应用后缘装置来控制后缘处的激波强度时，无论在设计点还是偏离设计点时，翼型均具有良好的升阻特性和鲁棒性。

4 结论

本文针对层流翼型的设计和激波减阻问题开展优化设计，采用多岛并行多目标进化算法来解决翼型层流化和激波强度之间的矛盾问题。

(1) 利用RANS方程耦合边界层方程和可压缩线性稳定性方程作为转捩位置判断的准则, 并通过NLF416翼型和NACA0012翼型验证本方法的准确性，证明此方法能够比较准确地预测流动的转捩，可以作为翼型表面流动转捩判断的方法。

(2) 通过研究后缘装置不同长度和高度下的激波控制效果表明，在给定长度下随着高度的增加激波阻力得到抑制，在0.005 4%弦长左右时控制效果最好。分析最优外形后缘装置流动控制结果发现，后缘装置不会使得翼型表面流动转捩点前移，相反会推迟转捩的发生。同时进行后缘装置应用于三维层流机翼设计中的可行性分析，证明应用后缘装置可以极大地减少网格数量，缩短流动计算时间，具有扩展到跨声速三维机翼层流设计中后缘处激波控制的潜力。

(3) 通过层流区域最大化和总阻力最小化的算例验证自然层流翼型设计的思路是通过推迟翼面转捩位置实现层流区域的最大化，而不是总阻力最小化来实现的。

(4) 分析后把层流翼型设计的单目标优化问题转变成一个多目标优化问题，并采用多岛并行多目标进化算法耦合均衡Pareto理论求解该多目标问题。优化结果表明均衡理论耦合进化算法可以快速得到该问题的Pareto阵面。数值计算结果显示阵面上翼型的气动性能较初始翼型均明显改善，即翼型表面的层流区域较初始翼型增大、激波强度得到有效控制。另外在设计点和偏离设计点后缘装置均具有稳定的激波控制效果，安装后缘装置翼型的阻力发散马赫数也得到提高。