王 浩，徐晶磊
(1.北京航空航天大学 能源与动力工程学院，北京 100191；
2.上海交通大学 航空航天学院，上海 200240；
3.北京航空航天大学杭州创新研究院（余杭）航空学院，杭州 310023)
高超声速转捩对高超声速飞行器的飞行状态有显著影响，主要体现在气动热/力与动力方面。一方面，转捩后湍流的壁面摩擦力和传热速率远高于层流，大约是3～5 倍[1]，使得高超声速飞行器热防护系统的设计非常具有挑战性；
转捩后的气动载荷分布也会发生显著变化，影响气动配平、带来飞行安全问题。另一方面，湍流边界层的抗分离能力更强，可以有效地克服超燃冲压发动机入口处的边界层分离，并且湍流边界层的掺混能力也更强，可以增强发动机中燃料与空气的掺混，使燃料可以充分燃烧。如果能够准确预测到转捩的位置和区域，就更有把握通过优化外形、改变壁面条件等多种手段来延迟或促进转捩的发生，以满足种种工程需求。在真实飞行器的三维边界层中，受压力梯度和后掠角的共同作用，边界层内近壁区往往会出现与边界层外势流平面内流线方向垂直的流动分量，称为横流流动[2]。当横流速度出现拐点时，将会产生流动失稳并容易发展为三维边界层转捩。因此，高超声横流转捩的研究有着迫切的实际需求，发展相应的预测或计算方法具有重要的工程意义。
高超声速横流转捩是集来流湍流-激波相互作用、可压缩湍流、横流、转捩等复杂流动因素于一身的、长期困扰科研人员的难题。转捩本身的难点在于，其起始位置对扰动的演化十分敏感，只有精确计算出扰动的演化过程才能得到准确的起始位置。而演化过程涉及到上述4 个问题的前3 个，即使这3 个问题得到准确建模，建模公式的同时使用也可能产生耦合误差。目前在工程应用层面，基于雷诺平均N-S方程（RANS）的湍流/转捩模型是预测转捩较为实际和有效的选择。2000 年Suzen 和Huang[3]构造出间歇因子输运方程，它与湍流输运方程合并成为转捩模型。2006 年Menter 等[4-5]构造出γ-Reθ-SST 四方程转捩模型，只采用当地流动变量，从而与现代CFD 程序兼容。随后该模型集成在ANSYS 商业软件，达到工业级成熟度，成为转捩模型的标杆。该模型没有考虑转捩机理，而是基于经验关系式构造转捩判据，而经验关系式是通过低速湍流边界层、压气机等标定的。2009 年王亮和符松[6]提出k-ω-γ转捩模型，对超声速平板边界层和超/高超声速圆锥边界层的预测结果与试验结果定量符合。直到2016 年，周玲在高超声速转捩模型中引入了横流模态时间尺度和横流转捩准则，对原模型预测的HIFiRE-5 横流转捩图像有所改善，是高超声速横流转捩的先驱建模工作。2016 年张毅锋等[7]对γ-Reθ-SST 转捩模型的低速经验关系式进行间歇因子的马赫数修正和湍流普朗特数修正，显著改进了原模型的性能，对多个高超声速尖锥的计算结果与试验符合较好。2017 年张毅锋等[8]进一步引入以横流雷诺数为基础的转捩准则，对HIFiRE-5 横流转捩标模，在不同雷诺数下都取得了与噪声风洞实验非常吻合的计算结果。2021 年向星皓等[9]使用eN方法对高超声速横流转捩数据进行拓展，结合横流强度与表面粗糙度构造当地化的高超声速横流转捩判据，构造的新模型对多状态下的高超声速尖锥进行横流转捩预测，取得了与实验结果符合较好的预测效果。
从上述间歇因子转捩模型的进化路线来看，为了捕捉高超声速横流转捩，在保证常规转捩的预测精度的同时，还需引入高超声速（或可压缩）效应和横流判据等，建模难度和工作量很大。基于RANS 的“转捩模型”还有一条古老的路线，那就是直接用湍流方程（k、ω等）计算转捩。2019 年徐晶磊等将KDO 模型[10-11]推广至转捩的计算[12-13]，具体做法是，以输运变量r=μt/μ重新标定模型参数，使模型整体达到流动结构自适应。而模型捕捉转捩的机制是，通过流动结构的自适应和湍流输运特性的保存，自发演化出转捩。这样，不用引入特定的转捩机理，却有可能捕捉到特定类型的转捩，具有潜在的多机理普适性。目前，该模型已成功预测了T3A 与T3B 平板边界层旁路转捩、T3A-平板边界层自然转捩、Aero-A 翼型分离泡转捩、DLR-F5 横流转捩、6 ∶ 1 椭球横流转捩、超声速平板边界层转捩以及超声速尖锥转捩。本文将该模型推广至高超声速横流转捩的预测，采用经典标模HIFiRE-5，评估其捕捉高超声速横流转捩的潜力与不足，为高超声速横流转捩模型的发展提供参考。
1.1 CKDO 湍流模型
精确的可压缩湍动能控制方程[14]为：
对于可压缩湍流，可压缩效应的建模是十分必要的。对于可压缩修正的KDO 模型，即CKDO 模型，方程右端各项依次分别模化为：
为了维持计算的稳定性，ε2的建模被分为两层建模，分别黏性子层和外层:
若Rek＜10，
若Rek＞10，
模型中的许多经验数值，只是为了精准刻画从DNS 数据得到的r～Rb与Rek～Aε分布曲线而采用的多层结构标定方法，导致的多层系数[12]。本质上，模型的经验参数只有2 个，即雷诺应力本构关系的Bradshaw 参数Rb及ε方程毁灭项系数Aε，这两个经验参数经平板边界层的Reθ=4 060 剖面标定后，就再未改变。
1.2 CKDO 转捩模型
注意到，雷诺应力本构关系中的Bradshaw 系数，在公式（9）被标定为关于Rek的函数，Rek包含壁面距离，而无论流动形态如何，壁面距离是固定不变的，不反映流动结构，这就是为什么原始 KDO 在T3B 平板算例中展现转捩行为，但转捩过早。KDO 以及CKDO 转捩模型引入湍流黏性与层流黏性之比r，它传达湍流的强度，是一个输运变量。再次使用平板边界层的Reθ=4 060 截面把Rb标定为：
CKDO 和KDO 模型的经验参数标定于充分发展平板湍流边界层的Reθ=4 060 剖面，但又能用于充分转捩区域甚至横流区域，这一点是令人疑惑的。解释如下：首先，理论上允许这样的CFD 计算系统存在，正如N-S 方程一样，并不考虑任何的转捩机理，却能捕捉转捩，CKDO-tran 计算转捩的机制与N-S 方程是一致的。其次，可以注意到，CKDO 和KDO 捕捉转捩的关键在黏性比r=μt/μ＜ 1 的涵盖层流区域和层流-湍流转捩区域（“pre-transition”区域）。只要r＜1 时，模型采用的相似性规律是普适的，那么模型就能普适预测转捩。实际情况是，当r＜ 1 时，公式（14）在Reθ=150～900 的“pre-transition”区域仍然普适成立[12]。为什么普适?r＜ 1 的区域常常对应于y+＜20 的非常近壁区，即黏性子层和部分过渡层，而该区域的流动属于壁面主导的流动，该区域的流体微团只感受到壁面的影响而基本无视外界的影响，形成一个普适的“局部平板边界层”。将采用公式（14）的CKDO 模型记作CKDO-tran，将采用公式（14）的KDO 模型记作KDO-tran。
为了展示对于超声速转捩流动，考虑可压缩效应的CKDO-tran 模型相对于KDO-tran 模型的优势，本节使用超声速平板边界层的DNS 数据作为基准验证算例。德州大学的Jiang 等[15]进行了来流马赫4.5 的平板直接数值模拟（U∞=705.093 m/s，T∞=61.111 K），壁温恒定等于层流流态下的绝热温度。文献用oblique-mode 扰动激发湍流，这意味着来流的湍流度未知，因此不同的研究者采用不同的来流湍流度，来使他们的模型获得最佳表现。对于KDO-tran 和CKDO-tran，来流湍流度分别取1%和0.73%，计算结果见图1。CKDO-tran 预测的湍流区摩阻大小与DNS 结果更接近，且捕捉到了转捩后的摩阻波动；
而KDO-tran 预测的摩阻偏低且过于平滑，原因在于未考虑可压缩效应。低估超声速流动的摩阻，这一现象在文献[16]也有体现，就是因为各转捩模型搭载的湍流方程忽略了可压缩效应的缘故。然而KDOtran 计算的转捩过渡段好于CKDO-tran，这与KDOtran 的一贯表现不一致（KDO-tran 的转捩过渡段总是很短），推测是KDO-tran 的湍流黏性过小，对过渡段造成了缓冲，歪打正着。
图1 平板表面摩阻分布Fig.1 Skin friction distribution on the flat plate
HIFiRE（Hypersonic International Flight Research Experimentation）是由美国空军研究实验室和澳大利亚防卫科研组织联合举办的高超声速飞行器测试项目，其目的是发展和验证下一代高超声速系统的关键技术。其中HIFiRE-5 外形设计为长短轴之比为2 ∶ 1 的椭圆锥，使得三维效应和横向流动更加明显，被专门用于研究高超声速三维边界层转捩。本文计算采用Juliano 等[17-20]风洞试验中的HIFiRE-5 缩比模型，总长为328 mm，横截面长轴半径为82 mm，短轴半锥角为7°，短轴所在截面头部半径为0.95 mm，模型三维示意图如图2 所示。HIFiRE-5 模型的计算网格采用从短轴中心线平面划分的二分之一模型，计算采用三重网格并行计算技术，共分8 个网格块，每个网格块的各边网格点数均控制为4n+1，经过网格无关性验证后，最终模型计算网格流向共181 个网格点，法向共113 个网格点，周向共129 个网格点，总网格量约为252 万。法向第一层网格高度约为5×10-4mm，并在激波位置附近加密。HIFiRE-5 结构网格三维示意图如图3 所示。
图2 HIFiRE-5 模型三维示意图Fig.2 3D model of HIFiRE-5
图3 HIFiRE-5 结构网格示意图Fig.3 Structured mesh of HIFiRE-5
本节主要模拟HIFiRE-5 模型在不同雷诺数、湍流度和迎角的工况下的转捩结果，并与Juliano[17-20]等的风洞实验结果进行对比分析。实验来流马赫数Ma=5.8，来流总温T0=410 K，等温壁面温度Twall=300 K。变化的单位雷诺数Re、迎角条件α、以及来流湍流度Tu，其工况在表1 中给出。
表1 HIFiRE-5 模型实验工况Table 1 Experimental setup of HIFiRE-5
边界层发生转捩，最显著的特征以及影响是壁面摩擦阻力和热流系数的急剧上升，这也是目前判断边界层转捩与否的常见方法之一。Juliano 等使用温敏漆（temperature-sensitive paint，TSP）技术进行测量分析，温敏漆显示的强度与温度变化成正比，所以可以用来判断壁面热流系数的相对大小，如果观察到热流急剧升高，即为发生边界层转捩现象。图4（a）显示了HIFiRE-5 模型在单位雷诺为10.2 × 106/m 工况的壁面温度变化ΔT分布，对应计算算例Case3。Case3正是各文献研究较多的工况，本文亦围绕这一工况开展多种效应研究。
图4 Case3 转捩图像Fig.4 Transition pattern of Case3
3.1 可压缩效应
高超声速状态下，湍流的可压缩效应十分显著，需要对湍流模型进行可压缩修正。本节主要评估可压和不可压两个版本的KDO 转捩模型，即CKDO-tran和KDO-tran 的预测能力。图4（b）和图4（c）显示分别的是CKDO-tran 和KDO-tran 预测的Stanton（St）分布，可通过St计算出实际热流。可见，CKDO-tran 预测的双肺叶转捩图像与实验符合较好，而在此湍流度下KDO-tran 也预测出双肺叶转捩图像，但转捩位置显著推迟了，造成转捩形状、大小与实验差距很大。此外，对KDO-tran 还进行了变来流湍流度计算，无论怎么调整来流湍流度，都无法取得与实验符合的转捩图像。CKDO-tran 和KDO-tran 的计算结果表明，对于用湍流模型直接计算转捩的方法，可压缩效应的引入对高超声速横流转捩的预测至关重要。
3.2 多雷诺数工况评估
Juliano 的风洞实验在保证多数来流参数的条件下，进行了多个变雷诺数工况研究。文献多是对Case3 的雷诺数工况进行研究，本文将转捩模型拓展至更多雷诺数的研究，一方面通过更多工况验证模型的适用性、另一方面为他人的对照研究提供参考。图5、图6 给出Case1、Case2 工况下，实验温敏漆显示的转捩图像和CKDO-tran 计算的St云图。联合Case1、Case2、Case3 的实验转捩图像来看，随着雷诺数的增加，转捩形成的高温双肺叶区域面积增加，并且形状变得更加瘦长。CKDO-tran 的计算结果符合这种趋势，且预测的转捩图像与实验总体符合较好。
图5 Case1 转捩图像Fig.5 Transition pattern of Case1
图6 Case2 转捩图像Fig.6 Transition pattern of Case2
3.3 迎角效应
尽管实验给出了带迎角的HIFiRE-5 高超声速转捩工况，文献中却缺乏相应的转捩模型研究，可见带迎角的高超声速转捩具有一定挑战性。本节在此进行Case4 计算评估，用以展示迎角效应、评估模型。
图7 给出Case4 工况下，实验温敏漆显示的转捩图像和CKDO-tran 计算的St云图。转捩图像来自于HIFiRE-5 模型的迎风面，此时转捩图像从原来的双肺叶形状变成顶部带尖角的钟形。CKDO-tran 预测的转捩图像与实验差距较大，主要是因为转捩位置整体提前。然而，CKDO-tran 捕捉到的中轴线上的蓝色细长泡、蓝色细长泡上下方的细长泡（低温区域）、以及随后的钟形及尖角，都与实验图像一一对应，表明CKDO-tran 预测的转捩图像与实验具有相似性。向星皓等[9]对高超声速尖锥横流转捩的计算，同样发现了随着迎角增大，计算的转捩位置靠前、转捩图像相似但有所变形的现象。可见，真实情况的迎角效应是转捩图像变化；
而湍流/转捩模型的迎角效应，是随着迎角增加，预测精度下降。
图7 Case4 转捩图像Fig.7 Transition pattern of Case4
鉴于CKDO-tran 计算转捩的机制是流动结构自适应自发捕捉，并已考虑可压缩效应，变迎角导致的迎风面压缩性变化不应导致如此大的计算误差。鉴于各向同性湍流-正激波这种最简单的湍流过激波问题，各种标准湍流模型的计算误差都很大[21]，本文推测，误差应来自来流湍流与头部脱体激波的相互作用。两者间的作用在迎角0°时导致的误差尚可接受，一旦有迎角，湍流-激波作用的各向异性造成的波后湍流度不均匀性就会增强，反映在迎风面上的转捩面的剧烈变化。因此，为了精准预测高超声速横流转捩，需加强湍动能或间歇因子过激波过程的建模。通过改进来流湍流过激波的建模，也能提高转捩模型对迎角0°工况的计算精度。
3.4 湍流度的影响
图8 给出了湍流度对转捩的影响。在前文，分析了湍流-激波作用的影响，认为湍流-激波作用建模的误差，导致飞行器头部感受到的来流湍流度失真,失真的湍流度造成失真的转捩图像。这一观点成立的前提是，来流湍流度是影响转捩的主要因素之一。事实上，转捩图像对湍流度的变化十分敏感，本节在此展示。仍然采用带迎角的Case4 算例，仅变化其中的湍流度，形成工况Case5 和Case6，计算的转捩图像分别见图8（a）和图8（b）。可见，随着湍流度的降低，CKDO-tran 计算的转捩位置逐渐后移；
转捩图像逐渐偏离实验结果，与实验图像的相似性最终消失。
图8 湍流度对转捩的影响Fig.8 Effect of turbulence intensity on transition
本文采用HIFiRE-5 经典标模，评估了CKDO-tran模型预测高超声速横流转捩的性能，并研究了可压缩因素、雷诺数因素和迎角因素对转捩的影响。结论如下:
1）CKDO-tran 的高超声速横流转捩预测精度显著高于KDO-tran 模型，表明可压缩效应的捕捉，对基于流动结构自适应的转捩模型至关重要；
2）采用CKDO-tran 计算了迎角0°、多雷诺数的HIFiRE-5 标模，丰富了雷诺数工况，为他人研究提供了参考。计算得出的横流转捩与实验符合较好，表明该模型对0°迎角的横流转捩具有较好的适用性；
3）4°迎角下，CKDO-tran 计算的横流转捩图样结构与实验相似，但转捩线过于提前，表明CKDO-tran对非0°迎角的高超声速横流转捩的预测具有一定潜力，但仍需进一步完善该模型。
4）基于间歇因子γ的转捩模型，用湍流模型提供全湍流解，用γ指示层流区域和转捩区域。然而，γ方程里显含k与μt，这正是来流k与μt影响转捩位置的原理。基于γ的转捩模型是有可能在k方程过激波计算失误的情况下，通过调整γ方程的其他参数而获得准确的预测结果的。因此，间歇因子转捩模型的转捩位置同时取决于湍流方程和γ方程。
本文完全通过湍流方程的演化捕捉转捩，排除了γ方程的因素，所得的结论完全针对于湍流方程，可为各模型的湍流方程的改进提供可靠的信息，如需加强来流湍流过激波的研究等。
猜你喜欢横流来流迎角两种典型来流条件下风力机尾迹特性的数值研究能源工程(2022年2期)2022-05-23连续变迎角试验数据自适应分段拟合滤波方法北京航空航天大学学报(2021年6期)2021-07-20横流热源塔换热性能研究煤气与热力(2021年3期)2021-06-09不同来流条件对溢洪道过流能力的影响水利规划与设计(2020年1期)2020-05-25基于横流风扇技术的直升机反扭验证北京航空航天大学学报(2017年5期)2017-11-23脊下横流对PEMFC性能影响的数值分析电源技术(2015年11期)2015-08-22弹发匹配验证试验系统来流快速启动技术研究弹箭与制导学报(2015年1期)2015-03-11失速保护系统迎角零向跳变研究科技传播(2014年4期)2014-12-02横流中多孔射流的稀释特性实验研究实验流体力学(2011年5期)2011-01-14来流边界层拟序结构对超声速混合层流场结构的影响研究实验流体力学(2011年5期)2011-01-14