新一代战斗机非定常流动数值研究综述

2020-07-08肖志祥崔文瑶刘健罗堃宇孙元昊

航空学报 2020年6期

肖志祥，崔文瑶，刘健，罗堃宇，孙元昊

清华大学航天航空学院，北京 100084

新一代战斗机强调4“S”或5“S”标准，其中“S”指：超机动能力(Super-maneuverability)、隐身(Stealth)、超声速巡航(Supersonic Cruise)、超视距空战(Super Beyond Visual Region)或短距起降(Short Taking-off and Landing)；除第4个“S”外，其余所有“S”均与空气动力学密切相关，第1和第2个“S”则对非定常空气动力学提出严峻挑战，如大攻角静态/动态失速、内埋弹仓、大S弯进气道等。

要获得上述非定常气动性能，风洞试验、模型飞行和真实飞行都是不错的研究手段，但是要么尺度太小(风洞试验和模型飞行)，要么成本高数据少(飞行试验)。在方案设计和初步设计阶段，采用计算流体力学(Computational Fluid Dynamics, CFD)方法开展相关数值仿真研究，具有很强的紧迫性。除了获得非定常流场外，更重要的是获得非定常压力脉动数据，为结构设计提供输入。

1 数值方法

1.1 非定常RANS-LES混合方法

要准确预测非定常特征，真正的非定常湍流预测方法必不可少，比如大涡模拟(LES)、直接数值模拟(DNS)，至少要求解非定常雷诺平均Navier-Stokes(URANS)方程组。传统的RANS无法满足非定常湍流预测精度要求，其中:① URANS 虽然所需计算资源少、效率高，但对于小尺度非定常流动结构的预测能力有限，无法准确给出压力脉动、动载荷等近场非定常湍流特征，几乎无法满足非定常流动预测的精度要求；② LES 能准确计算大尺度运动，对小尺度运动采用亚格子应力模式进行模化，具备真正的揭示复杂流动机理的能力，但在模拟高雷诺数壁湍流等时，计算量过大；③ DNS所需计算资源较LES还多得多，更难以应用于工程高雷诺数湍流的精细预测。Spalart根据解析Kolmogorov尺度的需求和充分考虑计算机技术发展后，认为要实现准确预测百万量级雷诺数的流动，LES在2045年可变成现实，而DNS则需要等到2080年[1]。

因此，工程上迫切需要发展兼顾计算精度和效率、适于高雷诺数、考虑更多物理机制、高精度、高效率的非定常湍流预测方法。目前看来，在未来10～20年里唯一能满足上述要求的仍只有RANS-LES混合方法：它利用RANS模化近壁区域流动节约计算网格、提升计算效率；利用LES解析远离物面区域的大尺度流动，兼备RANS和LES方法的优点、且利用对方的优点克服自身的不足；该类方法自20年前提出起，就得到了高度重视，并获得广泛应用。

广义而言，RANS-LES混合方法指的是所有同时采用RANS和LES方法以获得解析湍流的计算方法。经过20余年努力，RANS-LES混合方法也发展出了较多分支，主要有雷诺应力混合模型，脱体涡模拟(DES)、分区RANS-LES模型，RANS限制LES模型及表现出LES特性的新一代URANS方法，如部分平均Navier-Stokes(PANS)、尺度自适应模拟(SAS)等。

DES类方法因其构造方式简单、对复杂外形的适应能力强是目前应用最广泛的RANS-LES混合方法。DES类方法通过引入网格过滤尺度，不同的计算区域采用不同的湍流模型，分界面动态变化，在壁面附近用RANS模拟小尺度湍流，分离区域采用类似Smagorinsky的亚格子应力模型。Spalart等[2]1997年提出原始DES方法；Strelets于2001年基于剪切应力输运(SST)模式提出DES方法的一般形式[3]；为解决网格诱导分离现象，Menter等于2003年[4]和Spalart等于2006年[5]提出延迟DES，即DDES(前者由Menter提出)；2008年[6]后者团队为解决对数区不匹配问题，通过引入壁面LES模型提出改进的DDES，即IDDES方法。

最近几年，为了解决RANS-LES混合方法由RANS向LES过渡的“灰区”问题，人们也逐渐采用重叠方法、引入合成湍流方法等，取得了一定效果，但是仍未完全解决。本文作者通过修改DDES方法中的网格尺度常系数(CDES)为自适应系数(DDES-AC)，有效提升了强剪切流动到充分发展各向同性湍流的预测精度[7]；还采用考虑涡矢量方向的长度尺度，有效提升了跨声速射流近场流动和远场噪声的预测精度[8]。

1.2 自适应耗散格式

究其本质，RANS-LES混合方法(耦合自适应耗散格式)是一类解析的非定常湍流预测方法，它要求计算网格足够密，格式精度足够高，耗散足够小；然而，当流场中存在激波、壁面或远场时，网格往往不足以解析当地小尺度流动结构，也不足以克服当地数值误差形成的扰动，要求耗散不能太小，否则极易导致计算发散。

因此，必须构造适于非定常流动、与RANS-LES混合方法匹配的自适应耗散格式，否则其计算结果将严重偏离实际。本文作者研究了原始耗散、直接降低耗散和自适应耗散在串列双圆柱流动中的表现，发现原始耗散严重高估物面的压力脉动[9]；在研究高超声速凹腔诱导转捩时，发现如果自适应耗散中不考虑激波影响，计算极易发散；而通过引入激波探测器，则该困难迎刃而解[10]。

1.3 其他数值方法

除上述RANS-LES混合方法及与之匹配的自适应耗散格式外，含多种物理机制的基准湍流模式，尤其是转捩/湍流一体化模式[11]、网格数目足够且根据流动特征生成高质量计算网格、高精度时间推进方法(时间步长足够小、子迭代步数足够多)，都是必不可少的。

本文作者开展的研究均基于具有完全自主知识产权的、自研CFD软件UNITs，它采用任意多块结构化网格，多种空间离散格式、时间推进方法和RANS-LES混合方法。该软件历经多个具有详细实验数据的标模检验和验证，兼具流动机理揭示和工程实际应用的特点，已经在众多工程型号中得到应用。

本文主要采用两套新一代战斗机标模对UNITs 程序进行验证，分别命名为早期的“标模-1”[12]和近期的“标模-2”[13]，见图1和图2。需要说明的是，两套标模均为不通气模型，且已进行过风洞模型试验，具有试验的集中气动力系数CL和CD与力矩系数Cm，但是没有表面压力及脉动压力试验数据，如图1所示，数据来源包括试验(Exp)和两方程k-g模型和代数B-L模型，图中α为攻角。

无黏通量采用旋转Roe格式计算，变量的重构采用三阶MUSCL插值；扩散项采用二阶中心格式；时间推进采用隐式LU-SGS方法，非定常计算采用伪时间子迭代，可达二阶精度。

常规URANS方法可高效获得全机气动力，如孙元昊[13]和崔文瑶等[14]采用SST模式、3套网格(210万、2 230万和4 500万)获得了标模-2的气动力特征，并与风洞试验进行对比：在0°～90°攻角范围内，相对稀疏的网格(210万)获得的升力、阻力和力矩可以被接受，见图2。

图1 标模-1的气动力系数对比[12]

图2 标模-2的气动力系数对比[13]

本文涉及新一代战斗机典型非定常流动数值预测，拟从如下3方面进行探讨：大攻角静态失速及控制、动态失速、内埋弹仓动载荷及其控制。

2 战斗机大攻角静态失速特性及控制

动态失速较静态失速更难，计算量更多，故下面先介绍战斗机大攻角静态失速流动及控制。

2.1 大攻角静态失速特性

新一代战斗机强调超机动能力，要求60°攻角可控。在如此大攻角下，全机流动充满了各种旋涡(机翼涡、边条涡、鸭翼涡、机身涡)，旋涡与下游部件(尤其是立尾)的干扰非常严重。当攻角逐渐增加时，旋涡破裂顺序一般为机翼涡、边条涡、鸭翼涡和机身涡等；破裂位置则随攻角增加而靠前。在特定角度下，破裂的旋涡如果刚好位于斜立尾(隐身考虑)前方，此时低频、高能量的非定常气流将周期地冲刷立尾，造成不可忽视的动载荷，形成立尾抖振，久而久之，立尾会出现不同程度的裂纹，影响飞行安全和寿命。

以标模-2为例，图3[13]给出了标模-2在不同攻角下边条涡破裂点情况，在边条涡所在的纵截面处画出无量纲流向速度U的云图，去掉U>0部分，留下的云图即为旋涡破裂区域。可看出，攻角为20°与24°时，边条涡并没有发生破裂；当攻角为28°时，首次出现了流向速度U<0的区域；随着攻角增大，U<0的区域范围增加。

在2005年，本文作者团队[12]就采用两方程k-g模式数值分析了标模-1的大攻角流动(12°、24°和40°)，发现40°时破裂的边条旋涡涡团正好位于斜立尾正前方，见图4。由于当年的计算网格数目有限(全模1百万)，预测精度有待改善。

然而，URANS方法无法获取准确的频谱特征，也无法获得小尺度的旋涡结构。要获得上述特征，必须采用真正的非定常湍流预测方法。Morton等[15]对比了SST、SA和DES-SA预测得到的F/A-18C战斗机在攻角30°时的流场，DES-SA能够捕捉到更精细的流场结构(图5)，且通过对比数值与飞行功率谱密度，发现计算的误差可接受。Forsythe和Woodson[16]用DES方法模拟了F/A-18E战斗机在大攻角时的流动，证明DES方法在预测表面压力分布时也比URANS方法也更具优势。Jeans等[17]使用DDES方法研究了翼身组合体在攻角23°时前缘涡破裂现象，与实验上观测到的涡破裂行为一致。Rizzi和Luckring[18]使用RANS-LES混合方法研究了F-16XL战斗机低速大攻角流动，证明了RANS-LES混合方法可更加深入地揭示流动机理。

图3 标模-2的不同攻角旋涡破裂位置[13]

标模-2是近年来多个课题组研究较多的标准算例，外形、流动参数保持一致。其中孙元昊[13]用IDDES预测了最大俯仰力矩系数所在的32°攻角时的流场，如图6(a)所示，清晰地揭示了旋涡结构以及它们与下游部件的相互作用；Zhang等[19]在非结构网格基础上，用IDDES方法预测攻角为36°时的流动结构、旋涡破裂以及气动力，如图6(b)所示；孟德虹等[20]使用IDDES方法研究了攻角40°时旋涡对斜立尾的非定常作用(图6(c)[20])；Xu等[21]使用DDES探究了攻角40°(图6(c)[20])、50°(图6(d)[21])及60°(图6(e)[21])时的大范围流动分离现象，并分析了其升力系数的功率谱密度(PSD)。尽管计算方法不完全相同，计算网格也不一样，计算平台迥异，但是先进的DES类方法均能获得丰富的流动结构，展示旋涡破裂特征；上述方法也可获得破裂的旋涡与倾斜的立尾之间的非定常干扰，获得其频率和幅值特征，为结构分析和设计提供输入参数。

图4 标模-1在3个攻角下的机翼涡和边条涡[12]

图5 F/A-18C战斗机SST、SA、DES-SA预测得到的旋涡结构及功率谱密度对比(α=30°, Re=13×106)[15]

对于具备双斜立尾的战斗机，实验上已经证明机翼涡或边条涡的破裂会导致垂尾抖振[22]，并且会导致垂尾的结构疲劳[19]。当垂尾浸没在涡破裂后的尾迹中时，垂尾表面的压力脉动会急剧增加[12]。因此，垂尾抖振的问题必须尽可能地控制。

2.2 大攻角流动控制

由图2可知，战斗机在32°攻角附近，其俯仰力矩系数为正(抬头)且达到最大，如何有效地降低大攻角下的俯仰力矩从而避免机头持续上仰[23]，或如何在大攻角下迅速低头，改出失速，最直接的方法是矢量喷管，如F-22[24]，当矢量喷管的倾角为15°时，可额外产生0.4的俯仰力矩系数。然而，目前多数战斗机暂不具备矢量喷管；可行的方法是采用措施控制压力分布，即增大重心前压力，同时减小重心后压力。与之对应，使机头涡或边条涡提前破裂，或机翼流动分离延迟，就可实现上述控制需求。通常地，控制措施可分为被动式或主动式。

图6 标模-2在攻角为32°、36°(非结构网格)、40°、50°和60°时的Q等值面云图

通过在前体安装吹气装置，破坏F-16的机头涡和边条涡，可有效地降低俯仰力矩[25]，这是典型的主动控制措施；但是主动控制措施通常要增加额外的系统和装置，复杂度高，应用较少。

目前常用被动控制措施包括机身外形修型[26]、前体边条[27]、前缘延伸[28]等；其中，机身上安装的减速板虽然主要用于减速，但同时也是一种减小俯仰力矩的被动控制措施。Anderson等[29]在低速风洞中对F/A-22 V-9模型进行了测试，得到了垂尾抖振的频率特征；当马赫数为0.2、攻角为32°时，其Strouhal频率与结构响应的频率恰好一致，极有可能引起垂尾的强烈抖振，翼刀则可有效地降低这种抖振响应。类似地，Sheta[30]预测了F/A-18战斗机前缘延伸处安装流向翼刀时，能够明显地减弱攻角30°时斜立尾的抖振(见图7)。齐孟卜[31]和Dong[32]等已通过风洞试验，证实了减速板对俯仰力矩的控制效果，但他们只得到了积分形式的力或力矩，并未提供流场细节，很难对俯仰力矩减小的物理机理进行详细分析。

崔文瑶等[14]在标模-2的前机身安装减速板，对最大俯仰力矩系数(α=32°)进行控制，发现当减速板打开角度为60°时，可降低整机俯仰力矩的60.22%；减速板的存在有效地增加了边条附近的压力，使得抬头力矩得到了有效的控制。而且，减速板不仅可以降低抬头力矩，还极大地降低了垂尾表面的压力脉动，声压级最大可降低11.8 dB(图8)。

3 战斗机大攻角动态失速

本节主要介绍机动动作中的大攻角动态失速特性，较静态失速更贴近真实飞行，因而更具挑战性。一般来说，飞机从小攻角上仰至失速过程，仅需要几秒，其气动特性依次由附着流、前缘分离及旋涡、旋涡破裂和大范围分离等主导。与静态失速最大的区别在于其迟滞效应，且受到初始攻角、俯仰频率及振幅的严重影响。新一代战斗机兼备三代机大攻角气动特性，其机翼平面形状类似于双三角翼；因此，研究双三角翼的动态失速，也可以获得新一代战斗机的动态失速特征。

图7 F/A-18边条流向翼刀及其立尾压力脉动控制效果[30]

真实飞行器的动态失速其机动动作复杂，运动规律多变。采用风洞或数值研究完全自由飞行的大攻角动态失速难度很大，为了降低研究难度，更多的是采用限制自由度的动态失速，如大攻角强迫俯仰振动或强迫滚转等。

Rizzi等[33]通过求解层流Navier-Stokes方程模拟了后掠角70°的三角翼在振荡状态下的非定常流动，同时考虑了时间积分参数效应、空间离散方法等的影响。Visbal等[34-35]探究了基于层流假设的低雷诺数三角翼在上仰过程中涡破裂位置的响应，并发现用临界点理论判断涡破裂位置的方法同样适用于振荡状态下的三角翼。Ekaterinaris和Schiff[36]用基于Baldwin-Lomax涡黏假设湍流模型的URANS方法探究了后掠角为76°/40°的双三角翼在大幅度振荡过程中的流动现象，捕捉到了涡破裂点的延迟。

刘健等[37]使用基于SST模式的URANS方法，结合刚性动网格技术模拟了振荡过程中标模-2的流动现象，详细分析了俯仰运动对非定常流动(特别是涡结构)的影响，其Ma为0.2，平衡攻角为40°，振幅为20°，图9给出了其升力系数和俯仰力矩系数在不同减缩频率k下的迟滞曲线。

动态失速中时间延迟是重要特征。Srinivas等[38]评估了破裂点的时间延迟，ΔτU∞/c，对于后掠角大于70°的三角翼来说，其范围为1～2。Gursul[39]提出了一种基于波传播特性的新解释，试图统一描述时间延迟机制，这也适用于稳定状态下的类似现象。螺旋模态的不稳定行为也引起了人们的关注，但相关的研究却不如时间延迟效应丰富。Coton等[40]发现在俯仰运动过程中螺旋模态不稳定的中心激励频率随流向位置的变化与静止状态时不同，且受振荡频率的影响。工程上最关心的是流动结构的时间延迟与动态俯仰稳定性之间的关系。目前一致认为动态稳定性取决于大攻角下的减缩频率和振幅[41]。然而，它们的定量关系及物理机制尚不清楚。

受到计算资源的限制，关于动态大攻角失速的研究大都基于层流假设或使用URANS方法，只能获得气动力等积分量的迟滞效应。而大攻角下旋涡破裂及其沿时间的历程效应，是动态失速的主要研究内容，只能采用更先进的数值计算方法才能获得。

刘健等[42]采用刚性动网格技术和DDES方法，模拟了80°/65°双三角翼静态36°和以36°为平衡攻角的正弦俯仰运动(振幅6°)过程中的非定常流动，见图10,t为无量纲时间。重点分析了涡破裂点位置、螺旋模态不稳定性、压力脉动和动态俯仰稳定性。其中涡破裂点位置的定义为：(cr-x)/c，c为弦长。与实验结果一致，涡破裂点的运动几乎接近一个简谐运动，并与上仰运动的频率锁定，同时还伴随着相位延迟。此外，还建立了破裂点的一阶和二阶微分数学模型，这些模型能预测一定范围内减缩频率下旋涡破裂位置的迟滞效应。

图8 标模-2 Q等值面及压力系数、均方根压力系数[13-14]

图9 标模-2不同减缩频率对升力和俯仰力矩系数迟滞曲线的影响[37]

Xu等[43]采用DDES方法研究了标模-2在0°～80° 攻角范围内大振幅俯仰运动时的气动特性，主要研究振荡减缩频率f*的影响，分别为0.4 Hz 和0.6 Hz。随着振荡频率的增加，非对称涡结构的发展将受到抑制；与此同时，还能够观察到飞机上表面的左右两侧分别出现气泡式和螺旋式涡破裂，见图11。但他们未能针对此流场进行深入分析旋涡破裂相关的物理机制及频谱特征。

然而上述大攻角下的正弦振荡运动仅仅是动态失速的一种典型状态，其他诸如滚转、偏航、侧滑或其耦合机动动作，还需要将空气动力学、飞行力学、飞行控制、矢量发动机等结合起来。只有充分认识和理解动态失速中的旋涡发展、演化、迟滞以及与飞机部件的相互作用，才能真正发挥新一代战斗机的气动潜力。

图10 双三角翼上的Q等值面及旋涡破裂点的动态响应[42]

图11 标模-2做俯仰运动过程中左右两侧非对称涡破裂现象[43]

4 内埋弹仓动载荷及其控制

对于新一代战斗机而言，隐身无异于生存力。为了有效降低雷达散射截面积，所有武器均需内埋。但是，当投放内埋武器时，弹仓开启后必定导致非常剧烈的动载荷(同时破坏了隐身)。通常地，内埋弹仓可简化为空腔流动，包含非常复杂的物理机理，如边界层分离、剪切层失稳、压力振荡以及对空腔后壁的冲击和噪声等。旋涡与腔体的相互作用会导致腔体内部设备或结构损坏，前缘剪切层会对存储武器投放产生不可预期的不利影响。为了内埋武器的安全存储、投放及弹仓结构的安全，非常有必要揭示弹仓非定常流动机理，并开展有效控制。

Lawson和Barakos[44]于2011年综述了2010年之前超过60个试验和计算的空腔流动研究。其后，还出现了许多空腔非定常流动的研究，如Liggett[45]、Temmerman[46]、Wang[47]、Arunajatesan[48]、Babu[49]、Sheta[50]、Hassan[51]等。许多研究涉及复杂外形，如真实武器舱的储存和投放问题(Lawson[52-53]、Kannepalli[54]、Khanal[55]、Chaplin[56]、Kim[57]、Barone[58]等)；自2010年后，大多数方法均为RANS-LES混合方法。罗堃宇等[59]采用特殊的网格拓扑结构，将尽可能多的网格集中于M-219凹腔区域，并采用IDDES方法耦合自适应耗散格式，可获得与风洞试验误差很小的计算结果(约1 dB)，见图12[60]。正是由于内埋弹仓具有大声压级特征，许多研究者探索了空腔流动的控制方法，也主要为被动和主动两类。无论是主动还是被动控制措施，均分为前缘装置和后缘装置两类。

图12 M-219干净空腔声压级对比[60]

主动措施主要是主动吹气(Zhang[61]、George[62]等)和等离子体(Yugulis[63]、De Jong[64]等)。主动控制措施的可变参数多，应用相对较广，但是会增加额外系统，目前尚处于研究探索阶段，距离实际使用尚存在较大距离。

与主动控制措施相比，被动措施虽不具备多种条件下的通用性，但是其实施方便，简单易行，得到了更多应用。前缘装置，通常如锯齿状[65]、块状[66]、横杆[67]和锯齿状扰流器[68]；其原理通常是为了抬高剪切层，使其远离空腔后缘，或增加其不稳定性并削弱大规模的冲击。后缘修型，如传统的后缘斜坡控制[69]，虽不常见，但却可以有效地抑制压力脉动。Saddington等[70]的试验研究表明，前缘控制技术比后缘控制技术能更有效地抑制空腔的噪声振幅。

使用先进的CFD方法(如LES或混合RANS-LES)研究流动控制主要包括前缘支杆[71]、射流[47]和吹气[61]等方面；但是对于真实空腔或武器舱外形控制装置的研究很少(如Panickar等[72]对F-35武器舱模型的研究)。

罗堃宇[60]分别研究了前缘锯齿及前缘横杆对空腔流动的影响，如图13所示。此外，还探究了不同尺寸前缘锯齿对不规则弹舱内部流动的影响。发现前缘锯齿控制装置使剪切层抬高，不再进入弹舱内部并于底部再附，而是从锯齿顶部水平向后延伸，将高速流体隔离在弹舱外部，弹舱流动形态变为典型的开式空腔流动，且弹舱内部形成一个较大的回流区，最终使得腔内声压级减小，见图14[60]。锯齿不仅降低了压力脉动，还降低了空腔阻力。此外，罗堃宇等[73]通过预测超声速凹腔，采用SPDMD方法，获得凹腔、前缘锯齿、前缘横杆、倾斜后缘的频谱特征，增强了对流动及控制机理的认识。当然，此类方法只能研究简单空腔外形，要用到复杂外形，还有很多工作要做。图15给出加装锯齿控制措施的空弹仓、待发流动图，由图可见，锯齿装置脱出的旋涡，在弹仓上空发展不充分，直接冲击到待发弹身，必定增强弹身压力脉动，削弱弹仓后缘的压力脉动[60]。

值得注意的是，内埋弹仓不会是凹腔那么简单的外形，一定会包括非常复杂的外形结构，比如进气道、内埋弹仓舱门、控制措施，存储的导弹、导弹发射架，甚至内部加强结构等。采用多块结构化网格数值分析，难度越来越大，因此必须发展非结构网格以处理复杂外形，同时提升非结构网格的格式精度，提升其解析非定常湍流流动的能力。