面向超声速民机层流机翼设计的转捩预测方法

聂晗，宋文萍，韩忠华,*，陈坚强，段茂昌，万兵兵

1. 西北工业大学 航空学院 环保型超音速客机研究中心/气动与多学科优化设计研究所，西安 710072

2. 西北工业大学 翼型、叶栅空气动力学国家级重点实验室，西安 710072

3. 中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000

对超声速民机而言，减阻是继低声爆技术基本攻克之后下一个亟需突破的核心关键问题[1]，是决定其能否持续商业运营的“拦路虎”，而层流减阻技术则是超声速民机气动减阻中的最重要途径之一[2]。层流减阻技术是一种通过精细化外形设计、流动控制等手段推迟飞行器表面流动转捩的发生，从而维持大范围稳定层流的气动减阻技术。由于层流边界层的摩擦阻力远小于湍流边界层，采用层流减阻技术可显著降低飞行器表面摩擦阻力，从而大幅改善全机气动性能。例如，对于典型超声速民机构型，在巡航状态下黏性阻力占总阻力40%左右，通过外形设计在机翼上表面实现40%范围的自然层流，可使全机升阻比提升约5%[3]。由于巨大的减阻潜力，层流技术被国际航空运输协会[4]在2050年飞机技术路线图中评选入下一代飞机最有发展前景的技术。

转捩过程的物理机理复杂、影响因素众多，是流体力学中尚未完全解决的前沿问题之一。为了满足工程设计应用的需求，需要结合现有经验和理论研究，发展鲁棒、准确的转捩预测方法。现有转捩预测方法可分为以下几类：① 基于稳定性理论的方法。稳定性理论是最早用来研究转捩现象的理论，其具备物理基础，并经过了长期发展，在工程问题中得到广泛的应用。这类方法包括：基于LST[5-9]/Bi-global[10]/PSE[11-12]的eN方法、DNS方法[13-14]等。② 转捩经验关系。转捩经验关系是通过理论研究、试验分析等归纳得到一系列经验性的转捩判断准则，包括：包络线拟合法[15]、数据库方法[16]、代理模型方法[17]、AHD准则[18]、Gleyzes准则[19]和C1准则[20]等。③ 转捩模式。转捩模式是以输运方程等形式将转捩关系式和CFD求解进行高效耦合的方法，包括：γ-Reθ模型[21-22]、k-ω-γ模型[23-24]、C-γ-Reθ模型[25]、基于AHD/C1准则的扰动放大因子输运模型等[26]。

针对三维后掠层流机翼设计，现有转捩预测方法中，目前在工程上应用最广泛的是基于线性稳定性理论的eN方法[27]。德国宇航院(DLR)[28-31]、法国宇航院(ONERA)[32]、荷兰宇航院(NLR)[33]等均发表了比较成熟的耦合eN方法和RANS求解器的转捩自动判断方法，并在三维机翼、翼身组合体等构型上进行了应用。国内以周恒[34]为代表的天津大学力学系早在20世纪60年代就开展了边界层稳定性和转捩的基础研究。唐登斌[35]开展了后掠机翼可压缩边界层稳定性分析研究。张坤[36-37]、朱震[38]和左岁寒[39]等分别基于线性稳定性理论和线性抛物化稳定性方程开展了三维机翼/翼身组合体边界层转捩预测研究。徐国亮和符松[40]研究了可压缩流动中横流扰动的失稳及其控制。黄章峰等[41-42]研究了后掠机翼边界层的横流不稳定特性以及转捩预测，并总结了稳定性理论在复杂外形流、非平行流、局部突变流、强三维流等工程问题中的应用。徐家宽等[43-44]通过线性稳定性理论得到转捩关系式，并发展了相应的转捩模式方法，用于预测亚/跨/超声速边界层中流向和横流不稳定性诱导的转捩。杨体浩等[45]通过eN方法对自然层流机翼翼套飞行试验中的机翼转捩位置进行了预测。韩忠华等[46-51]基于耦合双eN方法的RANS求解器，开展了自然层流后掠机翼优化设计的工作。

对于马赫数1.3～2.0超声速民机的层流机翼设计，基于NTS-NCF双N因子积分策略的eN方法还存在如下问题[31]：双N因子方法主要考虑边界层二维Tollmien-Schlichting(TS)波和横流(Crossflow, CF)驻波诱导的转捩，而超声速机翼边界层还需要考虑三维TS斜波和CF行波不稳定性。高雷诺数、大后掠等特点，使得超声速民机层流机翼和传统低速、跨声速层流机翼在边界层稳定性上存在很大区别。一方面，随着马赫数的增加，二维TS波的扰动放大率大幅降低，三维TS斜波成为更不稳定模态[7-8]，其中最不稳定三维TS斜波波角约为45°～60°[52]；另一方面，随着马赫数和后掠角的增加，CF不稳定性也有所增强，且除了频率为0 Hz的横流驻波不稳定性，一般还要考虑横流行波不稳定性[53]。牛海波等[54]在马赫数为6的三角翼静音风洞试验中发现了机翼表面存在横流行波主导的转捩。万兵兵等[55]对飞行试验测得的低频信号进行分析，首次证实了在真实飞行条件下高超声速飞行器边界层存在横流行波。此外，日本宇宙航空研究开发机构(JAXA)[53]和美国国家航空航天局(NASA)[56]等将eN方法应用于超声速大后掠层流机翼设计，并通过风洞试验、飞行试验等开展了转捩判据标定，也提出了在转捩预测中要考虑TS斜波和CF行波。

针对超声速民机机翼层流减阻设计中考虑TS斜波和CF行波的转捩自动判断需求，本文发展了耦合RANS求解器的eN转捩预测方法，主要工作包括3方面：① 改进了扰动放大因子积分策略，考虑了超声速机翼边界层TS斜波和CF行波诱导的转捩；② 开展了NASA超声速试验机翼的边界层稳定性分析，验证了本文方法的正确性；③ 开展了马赫数2.0、雷诺数1.39×107、后掠角60°的无限翼展层流机翼初步设计，探索了超声速大后掠机翼边界层横流不稳定性的抑制方法，验证了本文方法对于超声速民机层流机翼设计的适用性。

1 耦合RANS求解器的eN转捩预测方法

针对超声速机翼边界层转捩预测需要考虑TS斜波和CF行波的问题，本文发展了耦合RANS求解器的eN转捩预测方法。

1.1 eN转捩预测方法

eN方法进行转捩预测的思想是，在边界层中引入一个微小扰动，并通过线性稳定性分析，得到扰动的空间演化情况；当扰动幅值增长到初始幅值的eN倍时，认为转捩发生。其中，线性稳定性理论假设的三维小扰动形式为

(1)

对于三维扰动应用线性稳定性理论，为了求解特征关系式，需要补充扰动展向放大率的计算或选取方式。Mack提出[57]，可以假设扰动沿群速度方向(简化为无黏流线方向)增长，引入流线方向和弦向的夹角θ0，得到展向扰动放大率为βi=αitanθ0，展向距离为dz=dxtanθ0。对于无限翼展后掠机翼，该条件可进一步简化为βi=0，即假设扰动在展向的增长可以忽略[58]。Arnal等[59]通过数值试验发现βi是否为0不会给扰动放大因子N的计算带来较大差异。因此，为了简化稳定性方程的求解，本文对于任意的机翼构型均假设βi=0成立。

考虑到超声速机翼边界层中可能导致转捩的TS斜波和CF行波，本文发展了固定波角和固定频率方法来寻找不稳定TS和CF波扰动，进而计算中性曲线，再通过固定展向波数/频率方法或包络线方法积分求解扰动放大因子N。下面将详细介绍TS波和CF波扰动的区分、中性曲线计算、扰动放大因子积分和转捩判据Ntr等。

1.1.1 TS和CF波扰动的区分和中性曲线计算

中性曲线是扰动放大率为0的位置，即不稳定扰动幅值增长的起点。为了计算中性曲线，首先应区分超声速机翼边界层中可能存在的TS和CF波扰动，包括二维TS波(2D TS wave)、TS斜波(Oblique TS wave, OTS)、CF驻波(Stationary CF wave, SCF)和CF行波(Traveling CF wave, TCF)等，在流线坐标系下给出其频率/展向波数分布示意图(见图1)。在低速和跨声速层流机翼设计中主要考虑的是二维TS波和CF驻波，可以看出，两者的区分方式是：二维TS波沿无黏流线方向，波角为0°，而频率f不为0 Hz；CF驻波的频率为0 Hz，而波角不为0°。因此对于二维TS波，可通过固定波角φw=0°计算中性曲线，并得到不稳定TS波扰动的频率范围；对于CF驻波，则通过固定频率f=0 Hz计算中性曲线，从而得到不稳定CF波扰动的展向波数范围。

在马赫数1.3～2.0的超声速机翼边界层中，还需要考虑三维TS斜波和CF行波。从图1可以看出，在流线坐标系下，TS斜波和CF行波的展向波数和频率均不为0，仅从展向波数和频率上无法作出显式区分。一种区分TS斜波和CF行波的方式是：观察扰动放大率随波角的变化(图2)，CF行波只在波角为正或负的单一方向上增长(αi<0)，且最不稳定CF行波波角相对较大，一般大于70°；TS斜波在波角为正和负的2个方向上均有所增长，且最不稳定TS斜波波角相对较小，一般为45°～60°。该方法适用于由TS斜波或CF行波主导的典型边界层流动，而当两种扰动均存在时，可进一步结合波角等信息来判断放大因子包络线上的不稳定扰动类型。上述方法是一种隐式区分最不稳定扰动波类型的方法，其对于一般问题的普适性还需后续验证。

图1 超声速边界层TS和CF波扰动频率/波数分布示意图

图2 超声速边界层中TS斜波和CF行波扰动放大率随波角的变化

在此基础上，对于超声速机翼，可分别选取固定波角φw和固定频率f两种方法来确定不稳定扰动。首先，计算不同波角上的扰动放大率分布情况，并判断扰动类型；然后，为了确定不稳定扰动的频率范围，采用固定波角φw方法来计算扰动中性曲线，在其上寻找不同频率的TS斜波和CF波扰动；最后，为了找出各频率下不稳定扰动的波数范围，采用固定频率f的方法来计算扰动中性曲线，在中性曲线上找出各频率下不稳定TS斜波和CF波扰动的展向波数，从而确定不稳定扰动的频率/波数组合。

1.1.2 扰动放大因子N的计算

对于可能导致转捩的不同波数/频率的扰动，从中性曲线出发，可采用两种方法积分计算扰动放大因子。一种是固定展向波数/频率方法，该方法最早由Mack提出[58]，是通过对扰动波引入无旋假设，并采用无限翼展机翼条件简化得到的。该方法的好处是可以根据扰动波角、频率等区分不同类型的扰动，如TS斜波、CF行波、CF驻波等，但由于要针对许多不同频率/波数组合的扰动进行分析，计算量相对较大。类似的方法还有固定波长/频率、固定波角/频率方法等[7]。另一种方法是包络线方法，假定扰动总是以各方向上扰动放大率的最大值进行累加，即在求解放大因子时通过改变展向波数或波角得到扰动放大率的最大值。该方法计算量相对较小，但不能判断是哪种扰动导致转捩发生。采用固定展向波数/频率和包络线方法计算扰动放大因子的公式分别为

(2)

(3)

在马赫数1.3～2.0超声速民机的机翼边界层中，由于马赫数的增加，二维TS波受到显著抑制，应主要考虑三维TS斜波；而边界层转捩由CF驻波还是由CF行波占主导，则取决于来流湍流度、壁面粗糙度等条件[40]。例如，在来流湍流度大或机翼前缘粗糙时，转捩一般由CF行波主导；在来流湍流度小、机翼表面光滑的情况下，转捩一般由CF驻波主导。这一结论在超声速流动下还需进一步确认。静音风洞和飞行试验表明，在低湍流度、高马赫数环境下，飞行器表面也存在横流行波主导转捩现象[54-55]。而在超声速层流机翼设计中，考虑到CF行波的扰动放大率一般高于CF驻波，按照CF行波进行设计更加保险[60]。

1.1.3 转捩临界N因子的选取

对于计算得到的不同频率/波数扰动的放大因子，结合相应的转捩判据，就可以预测得到机翼边界层转捩位置。eN方法是一种半经验方法，通常需要借助风洞试验和飞行试验等手段对转捩判据进行标定。

针对不同类型的TS和CF波扰动，给出相应的转捩判据。根据与DLR-F4机翼风洞试验的结果对比[38]，二维TS波和CF驻波的临界扰动放大因子可分别选用Ntr_2DTS=10.5,Ntr_SCF=7.5。根据NASA的静音风洞试验结果[61]和JAXA的飞行试验结果[3]标定得到，在采用eN-包络线方法时，可选用Ntr_ENV=14[61]或Ntr_ENV=12.5[3]作为超声速机翼转捩预测的临界N因子；根据NEXST-1超声速民机构型的飞行试验转捩位置与稳定性分析结果的对比可得[53]，在采用固定展向波数/频率方法时，可选用Ntr_OTS=9.5和Ntr_TCF=9.8分别作为TS斜波和CF行波临界N因子。应注意的是，超声速边界层转捩受到来流扰动条件、来流雷诺数、壁面粗糙度等许多因素的影响，而ONERA和JAXA[60]等仅对来流雷诺数、来流扰动的影响进行了定性分析。为了保证eN方法转捩判据的泛化能力，还需要针对各种因素的影响开展深入研究。

1.2 耦合流场求解的转捩自动判断流程

本文方法将eN方法和RANS求解器相结合，其总体思路是通过从Navier-Stokes方程的层流边界层解中提取边界层信息，找出各位置的不稳定扰动，计算相应的扰动放大因子，并结合一定的转捩判据预测转捩位置，最后将预测的转捩位置返回流场求解器，实现闭合循环。流程框架如图3所示。下面具体介绍该方法中的各个步骤。

图3 面向超声速民机层流机翼设计的转捩预测方法框架

步骤1给定预估的转捩位置，用于激活湍流模型中的生成项，从而保证预估转捩位置的上游为层流边界层，以便开展稳定性分析。

步骤2在给定转捩位置下，进行流场求解，直至收敛。

步骤3提取机翼层流边界层信息，用于稳定性分析。一种提取方式是首先从流场解中得到机翼表面压力分布、边界层外边界速度等流动参数，然后求解可压缩边界层方程，从而获取光滑的边界层速度、温度剖面及其一阶、二阶导数。另一种方式是直接从流场解中提取层流边界层信息。根据NLR的研究结果[33]，采用后一种方式的话，至少需要在边界层内法向布置约60～70个网格点才能保证准确提取边界层信息。

步骤4采用固定波角φw和固定频率f方法，求解线性稳定性方程，寻找不稳定TS和CF波扰动，计算得到不同频率/波数扰动的中性曲线，作为扰动幅值增长积分计算的起点。

步骤5从中性曲线上选取不稳定扰动，采用固定展向波数/频率方法或包络线方法，求解线性稳定性方程，得到不同有量纲波数/频率的扰动放大率，积分得到相应的扰动放大因子N。

步骤6根据计算得到的不同波数/频率扰动的放大因子N，并结合相应的转捩判据，预测新的转捩位置。

步骤7将预测得到的转捩位置返回流场求解过程，并重复进行步骤2～步骤7，直至预测的转捩位置收敛。

本文的转捩自动判断方法框架，和DLR[28]、NLR[33]等发表方法的主要区别在于，对其中eN方法的扰动放大因子积分策略进行了改进，考虑了超声速机翼边界层TS斜波和CF行波诱导转捩的预测。

2 超声速机翼边界层转捩预测验证算例

为了验证所发展的转捩预测方法的正确性，本文选取了NASA的Owens等[62]研究的超声速后掠机翼模型，在不同来流雷诺数下进行了流场求解和边界层稳定性分析，并将计算得到的CF波扰动N因子和文献中的计算结果进行了比较。

2.1 NASA试验机翼模型和计算状态

所选用于验证转捩预测方法的机翼模型源自NASA的Swept-Wing Laminar Flow项目[62]，目的是研究超声速机翼边界层横流不稳定性及控制。NASA近年内依次开展了CFD分析、风洞试验和飞行试验研究。风洞试验在NASA Langley研究中心的20 in(1 in=2.54 cm)超声速风洞中进行，飞行试验在NASA Armstrong飞行研究中心的F-15B载机平台上进行。

该机翼的平面形状和剖面翼型如图4[62]所示，机翼前缘后掠角为65°，翼尖切角35°。根弦长croot=566.8 mm，展长b=266.7 mm，垂直前缘的弦线长度cn=239.5 mm,xn为垂直前缘距离。机翼剖面翼型为双弧形翼型，沿垂直机翼前缘方向布置。本文所用翼型由文献[62]结果拟合得到。

图4 NASA试验模型[62]的平面形状和剖面翼型

NASA开展风洞试验和飞行试验的来流马赫数为Ma=2.0，来流单位雷诺数为Re=4.6×106～13.2×106m-1，攻角为α=0°。风洞试验中观测到机翼表面的转捩前沿呈光滑状，因此认为是横流行波诱导的转捩；飞行试验中观测到的转捩前沿则是锯齿状，因此认为是横流驻波诱导的转捩。

2.2 NASA试验机翼流场求解

为了获得用于稳定性分析的基本流，在Re=6.2×106m-1和Re=11.5×106m-1两个试验状态下，对试验机翼的绕流流场和边界层流动进行分析。流场求解采用课题组自主开发的结构多块网格RANS方程求解程序——“PMNS3D”[63]。计算总网格量390万，机翼表面沿流向布置400个网格单元、展向布置60个网格单元、法向布置100个网格单元。第1层网格高度5×10-6m，法向增长率1.1。上下表面初始转捩固定在90%弦长位置，进行流场求解，然后根据从流场中提取的机翼表面压力分布、边界层外边界速度等信息，求解三维可压缩边界层方程，获取边界层速度剖面、温度剖面及其一阶、二阶导数。

以Re=11.5×106m-1状态为例，给出试验机翼表面压力云图、边界层外边界势流流线如图5 所示。可见，势流流线首先沿机翼前缘的方向，然后迅速恢复到接近自由来流方向，且和机翼表面压力梯度方向存在一定的夹角。沿垂直机翼前缘的剖面(见图5)，分别给出机翼剖面压力分布、当地后掠角分布，如图6所示，其中c为机翼弦长。此处，当地后掠角φsw定义为边界层外边界势流速度方向和当地压力梯度方向的夹角。从图6 可以看出，试验机翼表面存在较强的顺压梯度，而顺压梯度会增强横流不稳定性；机翼表面当地后掠角从机翼前缘的90°迅速减小到65°左右,接着缓慢减小，逐渐接近横流最不稳定的45°当地后掠角。给出试验机翼不同流向位置的边界层横流速度形(如图7所示，横流速度w为边界层内弦向和展向速度投影到垂直外边界势流方向的速度分量，ue为边界层外边界流向速度，x1、x2和x3为不同的流向位置)，可见机翼上表面横流速度方向均为正(从翼根到翼梢)，且其最大值沿向下游的方向而增大；不同流向位置的边界层流向速度形如图8所示(图中u为流向速度)，可以看出边界层沿向下游的方向厚度呈增加的趋势。下面将对NASA试验机翼的边界层扰动特性进行具体分析。

图5 NASA试验机翼表面压力系数Cp分布和边界层外边界势流流线(Ma=2.0,Re=11.5×106 m-1,α=0°)

图6 NASA试验机翼剖面压力系数和当地后掠角分布(Ma=2.0, Re=11.5×106 m-1, α=0°)

图7 NASA试验机翼不同流向位置的边界层横流速度形(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%)

图8 NASA试验机翼不同流向位置的边界层流向速度形(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%)

2.3 NASA试验机翼边界层稳定性分析

作为对转捩预测方法的验证，分别在Re=6.2×106m-1和Re=11.5×106m-1两个试验状态下对NASA试验机翼进行边界层稳定性分析。首先，给出不同频率下NASA机翼边界层有量纲扰动放大率随着波角的变化情况(见图9)，可见扰动仅在波角为正的方向增长。结合机翼表面压力分布、横流速度型分布等信息分析可得，NASA试验机翼边界层转捩由CF波占主导。然后，比较了不同频率、不同雷诺数下的CF波扰动放大因子，并和文献的计算结果进行了对比。

图9 NASA试验机翼边界层内不同频率和波角下的扰动放大率(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%, xn=23.5 mm)

2.3.1 不同频率/波数下的CF波扰动放大因子

以Re=11.5×106m-1、风洞试验状态为例，试验测得的转捩位置在垂直前缘距离xn=23.5 mm 附近。在0～50 kHz的频率范围内，通过固定频率方法找出中性曲线，并采用固定展向波数/频率方法计算扰动放大率。图10为不同频率下的CF波扰动放大因子包络线(图中Nenvelope为各流向位置下，固定展向波数的扰动放大因子曲线族的最大值)，可见随着频率首先从0 Hz增加到33 kHz，各流向位置的CF波扰动放大因子均有所增长，最不稳定横流行波频率约为33 kHz；而随着频率进一步增加到50 kHz，在xn=35 mm附近下游的扰动放大因子开始下降，而在xn=35 mm 附近上游的扰动放大因子略有增长。

图10 NASA试验机翼边界层内不同频率的CF波N因子包络线计算结果(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%)

图11分别为频率0 Hz的CF驻波和频率33 kHz 的CF行波的扰动放大因子，可见最不稳定的CF行波放大因子高于CF驻波，且最不稳定CF行波的展向波数约在2 424～2 727 m-1之间，低于最不稳定CF驻波的3 540～4 215 m-1之间。

图11 NASA试验机翼边界层内不同频率/波数的CF波N因子计算结果(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%)

2.3.2 不同雷诺数下的CF波扰动放大因子

对Re=6.2×106m-1和Re=11.5×106m-1两个状态下的CF驻波与CF行波扰动放大因子进行比较，如图12所示。可以看出，保持频率为0 Hz，最不稳定CF驻波的扰动放大因子NSCF的最大值均随着雷诺数的增加而增加；给出不同频率的CF波N因子包络线，得到最不稳定CF行波的放大因子包络线Nenvelope的最大值也随着雷诺数的增加而增加。

图12 NASA试验机翼边界层内不同雷诺数的CF驻波和CF行波N因子计算结果(Ma=2.0, α=0°, z/b=50%)

2.3.3 计算结果和文献结果的比较

将Re=6.2×106m-1和Re=11.5×106m-1两个状态下，计算得到的CF驻波与CF行波扰动放大因子，与NASA文献[62]给出的计算结果进行比较，如图13所示，图中xtr_TCF和xtr_SCF分别为不同试验中测量的横流行波和横流驻波诱导转捩的位置，Ntr_CAL和Ntr_REF则是相应的计算与文献给出的N值。可见，在不同来流雷诺数下，本文计算的CF行波的放大因子和文献结果吻合良好，CF驻波的放大因子也基本吻合，验证了本文发展的转捩预测方法的正确性。其中一些微小差异可能是由采用的边界层信息提取方法或扰动放大因子积分、稳定性方程求解方法的差异所导致，对转捩预测的影响不大。此外，还将飞行试验状态Ma=1.98、Re=10.9×106m-1、α=0°下预测的转捩线和试验测量的转捩线进行了比较(见图14)，飞行试验中转捩由横流驻波诱导发生。可以看到，当选取横流驻波临界转捩放大因子Ntr_SCF=7.5时，预测的转捩位置和试验转捩位置吻合良好，验证了本文方法对于超声速机翼转捩预测的适用性。

图13 NASA试验机翼边界层内CF驻波和CF行波N因子计算结果与文献[62]结果的比较(Ma=2.0, α=0°, z/b=50%)

图14 预测的NASA试验机翼转捩位置和飞行试验测量转捩位置[62]的比较(Ma=1.98, Re=10.9×106 m-1, α=0°)

3 超声速层流机翼设计算例

为了进一步验证本文发展的转捩预测方法对于超声速民机层流机翼设计的适用性，将其用于后掠角φsw=60°、来流马赫数Ma=2.0、弦长雷诺数Rec=1.39×107的无限翼展层流机翼初步设计。首先，基于可压缩三维边界层相似解，探讨了有利于抑制横流不稳定的理想压力分布；然后，通过气动反设计方法获得了满足理想压力分布特性的机翼外形；最后，基于固定展向波数/频率方法，给出了设计机翼的边界层TS波和CF波扰动稳定性分析结果与预测的转捩位置。

3.1 理想压力分布的探讨

在高雷诺数、大后掠超声速民机的层流机翼设计中，为了避免扰动在机翼前缘快速增长并诱发转捩，关键在于抑制边界层横流不稳定性。对此，通过对不同当地后掠角φsw、压力梯度因子βh和边界层外边界马赫数Mae的三维可压缩边界层相似解(即Falkner-Skan-Cooke解)进行比较分析，探讨了适合于抑制横流不稳定性的设计目标。由于横流不稳定性源自边界层横流速度剖面上的拐点，是一种无黏不稳定性，因此在分析不同的边界层剖面时，认为横流速度最大值wmax越大，即横流速度形越饱满，则横流不稳定性越强。

针对上面提到的3种因素，分别开展了3组Falkner-Skan-Cooke方程解分析，其流动参数分别为：①Mae=2.0,βh=0.3,φsw=0～90°；②Mae=2.0,βh=-0.06～0.5,φsw=60°；③Mae=0.5～4.0,βh=0.3,φsw=60°。第1组解如图15所示，其中η为无量纲法向坐标。可以看出当地后掠角从0°增加到45°时，横流速度最大值wmax随着后掠角的增大而增大；而随着当地后掠角的进一步增大到90°，wmax反而逐渐减小；说明当地后掠角越接近45°，横流不稳定性越强。第2组解如图16所示，可见，当压力梯度为0时，不存在横流速度分量，横流不稳定性最弱；而不管压力梯度因子>0还是<0，即边界层为顺压或逆压情况下，横流不稳定性都随着压力梯度的绝对值增大而增强。第3组解如图17所示，可以看出，wmax随马赫数增加而增加，横流不稳定性也随之增强。

图15 不同当地后掠角下Falkner-Skan-Cooke解的横流速度形(Mae=2.0，βh=0.3，φsw=0°～90°)

图16 不同压力梯度因子下Falkner-Skan-Cooke解的横流速度形(Mae=2.0，βh=-0.06～0.5，φsw=60°)

图17 不同马赫数下Falkner-Skan-Cooke解的横流速度形(Mae=0.5～4.0，βh=0.3，φsw=60°)

层流机翼反设计中，通过调整机翼外形来获得所需要的表面压力分布。根据上述因素影响分析可知，压力梯度绝对值越小，即机翼表面压力分布越平坦，横流不稳定性越弱。由此，可以提出一种“让流动在机翼前缘迅速加速(约0.5%弦长内)，继而维持微弱压力梯度”的理想压力分布，以此来抑制横流不稳定性的增长。

3.2 两轮次压力分布反设计

针对来流马赫数Ma=2.0的无限翼展后掠机翼，为了获得在机翼上表面满足理想压力分布的气动外形，基于课题组自主开发的代理优化程序——“SurroOpt”[64]，开展了机翼气动反设计。设计机翼弦长为1 m，后掠角60°。设计状态为：Ma=2.0,Rec=1.39×107,α=2°。反设计机翼几何外形和上表面目标压力分布如图18所示。根据收敛性分析结果，选取计算网格量为520×100×16，沿翼面流向布置400个网格单元，后缘割缝处布置60个网格单元，法向布置100个网格单元，展向布置16个网格单元，并在两侧采用循环边界条件。为了实现图18中的目标压力分布，需要让流动在机翼前缘很小一段范围(约0.5%弦长)内迅速加速，即反设计对机翼前缘的几何精度具有很高的要求。而直接采用类/形函数变换(CST)参数化方法对机翼前缘的拟合精度有限。因此，在NACA0003翼型的基准机翼基础上，本文采用了两轮次设计的方法：第1轮采用直接CST参数化方法，获得基本满足目标压力分布的中间外形；第2轮在该外形基础上，采用自由变形(FFD)参数化方法获得较好满足目标压力分布的设计外形。

图18 反设计机翼平面形状和上表面目标压力分布(Ma=2.0, Rec=1.39×107, α=2°)

3.2.1 基于CST参数化方法的第1轮反设计

在代理优化框架[65-67]下，以当前压力分布和目标压力分布的差量为设计目标，以机翼上表面外形为设计对象，开展了基于直接CST方法的气动反设计。仅对机翼上表面外形进行参数化，机翼下表面维持NACA0003翼型不变。参数化方法采用8阶直接CST，并将机翼前缘的CST类函数也作为1个设计变量，得到一共10个设计变量。设计空间通过对基准的NACA0003翼型法向坐标上下扰动50%得到。优化数学模型为

(4)

式中：n为流向位置标号；n0为流向位置总数；Cp,n和Cp_target,n分别为设计机翼表面压力系数和目标压力系数；si为设计变量；si_min和si_max分别为设计变量的上下界，即设计空间。

优化收敛曲线如图19所示。首先通过拉丁超立方抽样方法，在设计空间内选取50个初始样本点，建立代理模型；然后通过改善期望(EI)和最小代理模型预测值(MSP)两种加点准则，在固定设计空间内，不断寻找新的样本点，更新代理模型；最后根据样本分布情况，自适应调整设计空间，继续加点和更新代理模型，直至收敛。得到的第1轮设计机翼表面压力分布如图20所示，可以看到机翼上表面基本满足设计目标，但压力分布仍存在微小抖动。

图19 第1轮基于CST参数化方法的机翼反设计收敛曲线和样本分布

图20 第1轮设计机翼表面压力分布和基准、目标压力分布的比较(Ma=2.0, Rec=1.39×107, α=2°)

3.2.2 基于FFD参数化方法的第2轮反设计

为了抹平第1轮设计机翼表面压力分布的微小扰动，在其基础上又开展了基于FFD方法的气动反设计。在翼型表面压力分布波动的位置布置FFD控制点，通过改变FFD控制点坐标获得新的机翼上表面翼型，一共20个设计变量。设计空间由FFD控制点坐标上下各扰动10%得到。机翼下表面同样维持NACA0003翼型不变。

优化收敛过程如图21所示，首先通过拉丁超立方抽样方法，在设计空间内选取50个初始样本点，建立代理模型；然后通过EI+MSP两种加点准则，在固定设计空间内，不断寻找新的样本点，并更新代理模型，直至收敛。得到的第2轮设计机翼表面压力分布如图22所示，可以看到设计机翼能够良好满足给定的上表面目标压力分布。

图21 第2轮基于FFD参数化方法的机翼反设计收敛曲线和样本分布

图22 第2轮设计机翼表面压力分布和基准、目标压力分布的比较(Ma=2.0, Rec=1.39×107, α=2°)

3.3 设计机翼边界层稳定性分析和转捩预测

在设计状态下，分别给出配置NACA0003翼型的基准机翼和第2轮设计机翼的上表面边界层内不同波数/频率扰动的空间演化情况以及预测的转捩位置。转捩判据采用从JAXA文献[53]中得到的Ntr_OTS=9.5和Ntr_TCF=9.8。

对于基准机翼，从计算得到的不同频率、波角的有量纲扰动放大率分布(图23)可以看出，虽然扰动在波角为正和负两个方向都存在增长，但波角为正方向的扰动放大率明显高于波角为负的方向，且正波角方向的最不稳定扰动波角大于70°，符合CF行波特征，而负波角方向的最不稳定扰动波角在40°～50°，符合TS斜波特征，说明在波角为正方向最不稳定CF行波不稳定性要强于TS斜波模态；因此认为，在波角为正的方向计算的不同波数下最不稳定扰动放大率围成的包络线(图24)属于CF行波。根据前面给出的横流行波转捩临界因子Ntr_TCF=9.8，得到转捩位置在弦长x/c=0.27，即机翼上表面层流范围为27%。

图23 基准机翼边界层不同频率和波角下的扰动放大率(Ma=2.0, Rec=1.39×107, α=2°, x/c=0.2)

图24 基准机翼边界层内不同频率下的CF行波扰动N因子包络线计算结果(Ma=2.0, Rec=1.39×107, α=2°)

对于设计机翼，计算得到不同频率下设计机翼边界层内的有量纲扰动放大率随着波角的变化情况(见图25)，可见各频率下不稳定扰动在正负波角方向上均存在增长，可得设计机翼边界层转捩由TS斜波主导。然后，采用固定展向波数/频率方法，计算得到设计机翼边界层内不同频率的TS斜波和CF行波扰动N因子包络线(见图26)，可见CF波扰动在机翼前缘附近得到了很好的抑制，而TS斜波扰动N因子从5%弦长位置附近开始增长，在边界层终点达到最大值，但仍未超过转捩临界放大因子Ntr_OTS=9.5。最后，给出计算得到的不同频率下TS和CF波扰动放大因子如图27所示,可见，波角为负的最不稳定TS波扰动放大因子均大于波角为正的TS波；随着频率从9 kHz增加到12 kHz，在x/c=0～0.75范围内边界层最不稳定扰动放大因子Nmax增大，在x/c>0.75范围内边界层最不稳定扰动放大因子Nmax有所减小。

图25 设计机翼边界层不同频率和波角下的扰动放大率(Ma=2.0, Rec=1.39×107, α=2°, x/c=0.2)

图26 设计机翼边界层内不同频率下的扰动N因子包络线计算结果(Ma=2.0, Rec=1.39×107, α=2°)

图27 设计机翼边界层内不同频率/展向波数下TS和CF波扰动N因子计算结果(Ma=2.0, Rec=1.39×107, α=2°)

综上可得，设计机翼相比基准机翼能够更好地抑制CF行波不稳定性，在设计状态Ma=2.0、Rec=1.39×107、α=2°下机翼上表面可以接近维持全层流，具有良好的层流特性。

4 结 论

针对超声速民机机翼层流减阻设计对转捩自动判断的需求，本文发展了能考虑TS斜波和CF行波的耦合RANS求解器的eN转捩预测方法。主要结论如下：

1) 改进了eN方法的扰动放大因子积分策略，考虑了超声速边界层中可能诱发转捩的TS斜波、CF驻波、CF行波等不稳定性。采用固定波角和固定频率方法来寻找TS和CF波扰动中性曲线，通过固定展向波数/频率方法或包络线方法积分计算扰动放大因子，最后结合相应的转捩判据预测转捩位置。

2) 为验证所发展的转捩预测方法的正确性，选取NASA的Owens等发布的超声速后掠机翼模型，在不同来流雷诺数下进行了流场求解和边界层稳定性分析，结果表明计算的CF驻波和CF行波扰动放大因子N和文献结果基本吻合，验证了本文方法的正确性。

3) 将发展的转捩预测方法用于来流马赫数Ma=2.0、雷诺数Rec=1.39×107和后掠角φsw=60°的无限翼展层流机翼设计，提出了一种使机翼表面流动在前缘附近迅速加速，然后维持弱压力梯度的理想压力分布。经评估，设计机翼边界层内CF波在机翼前缘附近得到了很好的抑制，且TS斜波和CF行波扰动放大因子均未达到转捩临界放大因子Ntr，表明设计机翼上表面可以接近维持全层流，验证了该方法对于超声速民机层流机翼设计的适用性。

致 谢

感谢西北工业大学许建华、宋科、乔建领、张力文、丁玉临、许朕铭等的意见和帮助，感谢中国空气动力研究与发展中心涂国华、李晓虎、杨强、陈曦、董思卫、向星皓等的大力支持和帮助。