轴对称短舱自然层流优化及转捩敏感性分析

曹凡，张美芳，胡骁，唐智礼

南京航空航天大学 航空学院，非定常空气动力学与流动控制工业和信息化部重点实验室，南京 210016

发动机是飞机的核心部件之一，为了保护发动机，降低阻力和噪声，短舱应运而生。短舱主要由进气道、整流罩和内、外涵道喷管组成，是推进系统的重要组成部分[1]。随着涡扇发动机的发展，涵道比不断增大，短舱直径也越来越大，势必带来阻力的上升，影响飞行效率。面积越大，短舱的设计难度就越高，大涵道比发动机短舱是综合了气动设计、结构设计、材料和工艺等多种学科的复杂系统。为了减小短舱面积增大带来的摩擦阻力，通过气动设计的手段对大涵道比发动机短舱进行外形层流化设计是近年来兴起的有效手段[2]。

相同来流条件下，层流的摩阻远小于湍流摩阻，如果将层流技术应用于机翼、尾翼和短舱，能够减少整机阻力15%以上，因此，层流减阻技术在新型民机的气动设计中应用前景广阔。自然层流技术是一种层流主动控制方法，该技术通过对短舱型面的设计，使得短舱外表面保持较大的顺压梯度区域，从而抑制流向T-S(Tollmien-Schlichting)波的增长，推迟转捩的发生保持大范围层流区域[3]。目前，机翼和机身的优化减阻技术比较成熟，但对短舱的研究较少，在现有的经验和技术基础上对短舱进行优化，实现有效的减阻仍然是主要攻关问题。

近年来，国内学者对NLF短舱气动外形设计也开展了一系列研究，但相关的风洞和飞行试验却很少。西北工业大学的何小龙等[12]研究了带动力的自然层流短舱气动优化设计；孟晓轩等[13]基于双eN方法研究了CRM(Common Research Model)短舱流动转捩的影响因素，认为马赫数和攻角对转捩位置影响较大,而雷诺数和湍流度相对较小。2019年,复旦大学Wang等[14]应用差分进化算法优化NLF短舱，将层流范围提升了16.64%，获得良好减阻效果;杜玺等[15]针对设计的自然层流短舱在德国ETW(European Transonic Wind)跨声速风洞进行了对比验证，结果表明低雷诺下转捩预测结果吻合较好，而高雷诺数下转捩提前发生与CFD(Computational Fluid Dynamics)结果差别明显。目前来看，针对NLF短舱的气动外形设计已经取得较好的进展，但实际设计过程中，识别对层流影响显著的变量更加重要，因此需进行参数灵敏度分析。高雷诺数下，保持大范围的层流流动极其困难，在短舱层流化设计中，外形几何控制参数对转捩位置的影响和气动参数与层流范围的关系，尚有待研究。

以CRM短舱为研究对象，在高雷诺数跨声速状态满足几何约束条件下利用发展的进化/确定性混合算法进行气动外形优化，分别以总阻力最小化和层流最大化设计轴对称NLF短舱，验证层流短舱减阻的可能性。为了分析几何控制参数、流量系数和马赫数对层流范围的影响，基于层流最大化的短舱构型开展转捩敏感性研究，分析层流与阻力系数的关系，以期为后续自然层流短舱的设计和鲁棒优化提供依据。

1 层流转捩预测

(1)

Langtry和Menter[17]采用的关联函数为

(2)

(3)

(4)

式中：U为当地速度；θ为动量厚度；s为流线弧长；Ui为速度张量分量；xi为笛卡尔坐标张量分量；μ为层流黏性系数；μt为湍流黏性系数；S为流线弧长；k为湍动能；Tu为湍流度；λθ为压力梯度；ρ为密度；Pγ、Eγ为生成项、耗散项；σr=1。

(5)

(6)

(7)

式中：Fwake为确保Fθ t在翼型下游尾流区为0的参数；y为靠近壁面最近距离；δ为边界层厚度；cθ t为源项系数；σθ t为扩散系数，取2.0；ce2为耗散项系数;t为时间。

σθ t的大小决定了临界动量厚度雷诺数对流动历史和当地流动参数的依赖程度，σθ t越大,转捩对流动历史越不敏感，更多地取决于湍流度和压力梯度等当地流动参数。扩散系数和源项系数的变化虽然会带来转捩位置的改变，但与气动参数和几何参数的变化相比，其影响程度相对较小[18-19]。

混合函数Fθ t在边界层外为0，在边界层内为1，起到关闭Pθ t的作用。为了耦合k-ωSST湍流方程，引入有效间歇因子γeff与湍动能方程进行交互，表达式为

γeff=max(γ,γsep)

(8)

Fθ t

(9)

式中：γsep用于捕获并模拟分离诱导转捩；Freattach为一种开关函数，表示当黏度比大到足以引起层流再附着时，禁止其修改；Rev为涡量雷诺数；Reθc为临界动量厚度雷诺数。

(10)

(11)

图1 NLF(1)-0416摩擦阻力系数曲线

表1 NLF(1)-0416翼型转捩结果对比[20-21]

DLR-F5机翼前缘后掠角20°，采用超临界对称翼型。计算工况:Ma=0.82，α=2°，Re=1.5×106，参考弦长0.15 m，湍流强度0.5%，湍流黏性比10，物面表面网格y+<1，第1层网格高度L=1.5×10-6m，网格量约为247万，如图2所示。

图2 DLR-F5计算网格

图3 DLR-F5试验结果与数值模拟结果对比

2 高雷诺数下轴对称自然层流短舱优化设计

大型民用运输机短舱的流动雷诺数一般高于2×107，为了实现高雷诺数下的层流流动，需要保持更加强烈的顺压梯度，这远比在低雷诺数下的优化更加困难。通过精确的参数化方法、高精度流场求解、可靠的转捩预测结合高效的优化算法进行高雷诺数下的自然层流短舱优化设计。

2.1 轴对称短舱几何设计

为了设计满足几何约束的短舱，同时对比优化前后的层流范围，首先对CRM标模短舱[22]进行自然通气状态下仿真。设计工况：巡航高度11 000 m，基于短舱参考弦长雷诺数3.5×107，自然来流马赫数Ma=0.85，攻角α=0°，湍流强度Tu=0.05%，湍流黏性比为10，流场计算采用有限体积方法，空间离散为二阶精度，通量计算采用 Roe 格式，时间推进为 LU-SGS (Lower Upper Symmetric Gauss Seidel) 隐式方法。

计算网格如图4所示，划分六面体结构网格共850万，其中物面流向和周向各201个网格节点，第1层网格高度4×10-6m。

图4 CRM短舱物面网格

计算结果如图5所示。CRM短舱在自然通气状态下，其层流范围很小，且流动在靠近唇口的位置便开始转捩发生流动分离。

图5 CRM短舱计算结果

经测量CRM短舱，其最大弦长6.06 m，半径1.75 m，为了满足几何约束(见表2)，对CRM短舱剖面进行CST(Class Shape Transformation)参数化[23]。

表2 短舱设计约束条件

CST方法由一个类别函数和一个形状函数构成，具有精度高和设计参数少等特点，其设计参数具有明确的物理意义，能够直接反映翼型的变化过程。其数学表达式为

(12)

式中：c为翼型弦长；Rle为翼型前缘半径；Δzte/c为后缘点纵坐标；β为后缘倾角；bi为形状函数控制参数(i=1,2,…,n-1)。采用基于Bernstein多项式的CST方法时，为了提高参数化精度,且防止过高的阶数使得参数化趋于病态，取n=6。此时，参数化短舱外剖面共需要8个设计变量(Rle、Δzte、β、b1、b2、b3、b4、b5)。

对于流量系数的约束，采用了闫海津和杜玺[24]提出的在短舱出口处添加堵锥的方式代替动力短舱实现流量系数的连续调节，初始流量系数(Mass Flow Ratio,MFR)MFR=0.7，最终形成的满足约束的轴对称短舱半模如图6所示。

图6 带堵锥的轴对称短舱

2.2 自然层流短舱优化设计方法

通常更大范围的层流要求短舱表面维持较大范围的顺压梯度，这会导致短舱后半段的压力快速恢复产生强激波，增大压阻和波阻从而对总阻力的减小产生负面效应。本文着眼于自然层流短舱的优化设计，期望在短舱的表面实现更大范围的层流，但对于具体的可优化程度没有经验参考，伴随摩擦阻力减小对激波阻力与总阻力的影响也未知。因此，为了更好地分析自然层流短舱的局部敏感性，分别以转捩位置最大化和总阻力最小化为优化目标进行优化，分析不同参数影响下的层流范围，验证层流减阻的可能性。希望获得一款能够保持较大层流范围的短舱进行敏感性分析，有利于研究层流范围变化规律。

2.2.1 进化/确定性混合优化算法

解决气动优化问题可以用3种不同的方法：① 无导数搜索算法(Derivative-Free Search Algorithns,DFSAs)；② 基于梯度的算法(Gradient-Based Algorithms,GBAs)；③ 进化算法(Evolntionan Algorithms,EAs)。每一种算法都有各自的优缺点，可以结合不同的算法[25]，把其优点同时发挥出来，比如:EAs的鲁棒性和全局搜索,以及GBAs的快速收敛能力。通常使用基于梯度的搜索方法，计算目标函数的导数以确定最佳的下降方向，梯度法已与遗传算法(Genetic Algorithm,GA)结合用来设计翼型和机翼[26-27]。

以GA为主算法，用拟牛顿方法(BFGS算法进行)辅助加速收敛，发展了一种新的多级混合优化算法。文献[28]对该混合算法的原理进行了详细介绍。进化和确定性算法在这种混合优化方法中交替进行，并成功应用于高超声速飞行器气动外形优化，提升了优化效率。GA/BFGS算法流程如下：

Step1初始化种群，并计算适应值。定义初始种群数量Np进化代数Ke，局部搜索选择频率为f。

Step2执行GA算法，产生较优种群。产生一个随机整数Ns(Ns≤Np)，并从当前种群中选中并标记Ns个个体。

Step3对标记的Ns个个体执行确定性优化算法BFGS，BFGS将这Ns个个体每个迭代Kb步，产生Ns个新个体,并将其加入种群替代原有个体。

Step4最优选择(即从GA种群中选择前Ns个最优个体)，新的搜索空间简化为收敛的较小空间。随着搜索空间的缩小，种群规模和搜索频率也可以减少。

Step5继续对简化的搜索空间和局部搜索频率进行优化，直至收敛。

多级方法在这个交替的过程中耦合。为了保持基因的多样性，可以在进化过程开始时使用大种群规模，以便选择更好的基因并传给下一代。该算法交替过程如图7所示。

图7 GA/BFGS混合优化算法示意图

为验证算法的有效性，对(GA/BFGS算法)与纯梯度(BFGS)算法[29]、差分进化算法(Differential Evdution,DE)[30]、GA算法[31]在Ackley多峰函数上进行了对比验证，如图8所示。种群个体均为100、代数5 000、交叉概率0.9、变异概率0.1，差分算法中缩放因子取0.5。

Ackley多峰函数为

(13)

Ackley函数取值范围为x∈[-32, 32]，维度d=30，具有很多局部极小值，虽然梯度算法进化效率很高，但容易陷入局部最优解，而GA、DE这种进化算法能够利用交叉、变异算子跳出局部最优解，进化效率相对较慢。本文所使用的混合算法在进化前期效率不如梯度算法，如图8所示，但避免了陷入局部最优解，很好地平衡了进化效率和寻优能力。

图8 Ackley 函数收敛曲线

2.2.2 轴对称短舱自然层流设计

自然层流短舱设计过程中，目标是在激波前保持顺压梯度，抑制T-S波转捩，从而获得大范围层流区域。采用GA/BFGS混合优化算法结合六阶CST(形状类别变化函数)参数化方法，以总阻力最小化和转捩位置最大化为目标分别进行单点优化。

短舱计算工况：巡航高度11 000 m，基于短舱最大长度雷诺数Re=2.5×107，流量系数MFR=0.7，湍流度Tu=0.05%，黏性比为10。划分结构网格如图9所示，第1层网格高度5×10-3m，网格前缘用C型网格保证正交性，并进行加密处理，轴向分布200个网格点，网格数量约11万。

图9 轴对称短舱剖面网格

多层混合交替初始阶段种群规模Np=40，局部搜索选择频率为f=2，进化代数100代，保证优化过程中的基因多样性，随着搜索空间的缩小，种群规模也可以减少，这里减少到20个，局部搜索频率减少到1，继续进化10代，确保最后几代个体能够收敛到全局最优解。

计算结果如表3所示，Baseline为短舱设计工况，minCD是以总阻力系数最小为优化目标的工况，maxXtr是以转捩位置最大为优化目标的工况。ΔCD为短舱内外剖面的总阻力系数相对于基准外形的变化量，由于此次设计不改变内型线，其值变化可近似代表外剖面的阻力变化；Xtr为短舱转捩位置；Cf-up为短舱外剖面的摩擦阻力系数，参考面积为383.7 m2；Dmax为短舱外罩最大直径；Rle为短舱前缘半径。

表3 轴对称短舱优化前后计算结果对比

对比分析发现在总阻力系数最小时，转捩位置推迟了3.6%，摩擦阻力也略有降低，但短舱外罩最大直径大大降低，短舱剖面面积减小约28.5%，前缘半径变小，如图10所示。当层流面积最大化时，短舱外剖面的摩擦阻力减小24%，表明延迟转捩，增大层流范围极大地降低摩擦阻力，但也使得总阻力明显增加。这表明合理的平衡层流面积最大化和总阻力最小化2个目标，才能实现真正的短舱减阻目的。

图10 优化前后短舱外形对比

归一化后Cp、Cf曲线如图11、图12所示，原始外形由于唇口吸力峰的存在过早转捩，层流范围最大化的外形在前半段维持了良好的顺压梯度，可达到总长的27%左右，但是在短舱尾部出现由激波导致的强逆压梯度，压力快速恢复，出现相对强烈的激波，这是导致阻力略有增加的原因。阻力最小化外形后半段外剖面迅速向内剖面靠近，压缩了短舱面积。

图11 优化前后压力系数曲线对比

图12 优化前后摩阻系数曲线对比

3 轴对称自然层流短舱局部敏感性分析

基于单点优化设计，获得了2种不同的短舱构型，研究表明如果以总阻力为优化目标最后的优化结果的确会使总阻力减小，但是层流转捩位置推迟较小，层流范围增大不明显。单点优化的设计结果在偏离设计点时，往往不能保证飞行性能。同时，为了提取层流短舱的设计规律，分析了几何参数和飞行参数的敏感性。为了更好地分析单变量影响下轴对称自然层流短舱的转捩变化规律，选取层流范围最大的短舱构型为初始构型，较大范围的层流更容易对比分析不同参数带来的影响。

3.1 局部单变量几何敏感性

由分析可知，短舱剖面采用六阶CST参数化方法进行参数化，考虑几何约束去掉后缘厚度Δzte后还有7个参数，分别为前缘半径Rle、后缘倾角β、短舱剖面形状函数的控制参数(b1、b2、b3、b4、b5),控制范围如图13所示。在保持飞行状态不变的情况下，对CST参数化设计变量初始值的(-15%, 15%)区间如表4所示，取6组数据进行单一变量的几何敏感性分析，并给出Cp以及归一化后转捩位置,如表5、图14所示。

表4 设计变量取值范围

通过对比表5的数据，按对转捩位置影响大小进行排序，转捩位置随b1、b2、β的变化最为明显，其次为b5、Rle，影响最小的参数为b3、b4。其中b1对层流范围的影响最为显著，当其降低10%时，层流范围降至10%以下。由图13可知,b1主要控制短舱截面的头部区域，在短舱的设计加工中需要重点考虑。b2、β变化下，层流范围在10%、20%内波动较大，是次要考虑的设计点。而b3、b4的变化对转捩位置的影响不敏感，层流范围都可以保持在20%以上。

图13 形状控制参数控制范围

表5 不同几何参数变化下的转捩位置

通过分析层流短舱的初始构型压力分布如图14 中所示，在短舱截面前半部分保持较长范围的顺压梯度，紧接着出现由激波导致的强逆压梯度，压力迅速恢复。跨声速状态下，流动从翼型前缘开始，穿过一系列膨胀波加速到超声速状态，而后保持超声速流动直到后缘附近。由于压力需要恢复到后驻点处的压力值，同时由超声速状态变为亚声速状态，因此在层流短舱截面后缘附近产生一道激波，通过逆压梯度的大小可以评估激波的强度。

图14 自然层流短舱几何参数变化下的Cp曲线

随着b1的增大，曲线吸力峰值减小，外形表面的压力变化更为平缓，转捩位置变化较小，但b1减小时吸力峰变强使得提前发生转捩，层流范围变小；b2对前缘吸力峰的影响也较大。同时由于转捩位置变化较大，导致后缘处激波位置和强度变化也较大；b3、b4、b5、β在吸力峰后的压力变化趋于平缓，对压力分布的影响较小；随着前缘半径Rle增大前缘吸力峰变强，对激波强度的影响较小。

3.2 流量系数敏感性

MFR是指进入进气道的实际空气质量流量与同马赫数下进入进气道的理论最大质量流量之比[24]。本文优化后的轴对称层流短舱MFR=0.699，在(0.67，0.73)范围内通过移动堵锥实现流量系数的调节，转捩位置变化如图15所示。

图15 转捩位置随流量系数变化关系

流量系数增加的过程中，转捩位置的大小波动明显，MFR=0.678时，转捩位置前移到0.4 m以下。可以发现，巡航状态下进气流量对层流影响很大，此构型下尽量保持MFR>0.69才能获得较大范围层流。

3.3 马赫数敏感性

短舱能否在不同马赫数下产生较为强烈的顺压梯度以保持稳定的层流范围以适应飞行状况的变化是非常关键的。此次马赫数的单变量影响分析围绕巡航点依次在(0.8, 0.9)内取了10个工况，并在亚声速区间取4个工况,分析该自然层流短舱的变化规律，结果如图16～图18所示。

图16 转捩位置随马赫数变化关系

图17 Cp随马赫数变化曲线

图18 Cf随马赫数变化曲线

压力和阻力系数曲线都对马赫数变化体现出了明显的敏感性。需要特别指出的是当Ma=0.87 时，转捩位置出现强烈的波动，经重复计算排除了计算意外。除此之外，马赫数越大顺压梯度越强，吸力峰小幅降低，激波增强且后移。尤其是在飞行速度相差较大的情况，转捩位置在亚音速区间比跨声速区间更加敏感。跨声速状态下，随着马赫数增大，短舱整体上能够保持层流范围大于20%，层流范围波动较小。

4 结 论

在几何约束条件下,针对短舱剖面开展轴对称自然层流短舱优化设计研究，并利用CST参数化后带有物理意义的几何参数，分析了几何扰动、流量系数和马赫数的转捩敏感性。主要得到以下结论：

2) 基于层流最大化的轴对称短舱构型，开展了局部几何控制参数敏感性研究。外形控制参数b1、b2对压力分布和层流范围影响最大，而前缘半径Rle、后缘倾角β对压力分布影响很小，是层流范围的次要影响因素。b3、b4、b5的敏感性相对较低，其主要影响短舱外表面湍流段的压力分布。为保持大范围层流，要重点关注短舱剖面头部外形的设计。形成了基于CST控制参数的短舱剖面外形设计准则。

3) 转捩位置会随着流量系数和马赫数的增大呈波动状态，敏感性较高。在设计点(MFR=0.7，Ma=0.85)附近，层流范围能够保持在20%以上，降低了摩擦阻力，除了Ma=0.87以外，设计点附近均能保持较大范围的层流，为后续鲁棒性优化设计提供指导。