高超声速边界层转捩模型横流效应修正与应用

朱志斌，尚庆，沈清

中国航天空气动力技术研究院，北京 100074

在高超声速条件下，边界层转捩会对飞行器摩擦阻力、热流密度和流动分离再附特征等产生显著影响。特别是转捩后的湍流边界层会引起物面气动热载荷的大幅升高，进而对高超声速飞行器热防护系统带来严峻考验。高升阻比飞行器外部流场中，普遍存在三维边界层流动，出现显著的横流效应，其转捩过程同时受到第二模扰动、驻定横流模态和行进横流模态的共同作用。目前，对高超声速边界层横流转捩过程尚未认识清楚，对横流转捩的准确预测难度较大。

椭圆锥外形在0°迎角时存在显著的横流效应，并且与真实的高超声速飞行器相似，是高超声速边界层横流转捩研究的典型外形。Kimmel等对椭圆锥外形横流转捩现象做了开创性的研究工作，对比了不同横纵比的尖前缘椭圆锥边界层横流转捩特征。HIFiRE5(Hypersonic International Flight Research Experimentation)飞行试验研究项目选择横纵比为2∶1的椭圆锥作为试验载荷，针对三维边界层转捩问题开展了大量的风洞试验和数值预测以及两次飞行试验。HIFiRE5地面风洞试验主要通过缩比模型表面热流测量，观测到了椭圆锥表面转捩形态特征。Juliano等在普渡大学马赫数6静音风洞中(Boeing/AFOSR Mach-6 Quiet Tunnel，BAM6QT)，通过温敏漆对模型表面温度分布的整体测量，系统研究了来流噪声、壁面粗糙度、迎角和雷诺数对边界层横流转捩的影响。对HIFiRE5构型进行转捩预测的数值手段主要是流动稳定性分析,但线性稳定性分析在三维外形应用中存在固有缺陷，需要采用基于准平行性假设和曲率效应的抛物化稳定性分析方法进行转捩预测。经历了首次飞行试验失败之后，在2016年开展的第2次飞行试验(HIFiRE5b)取得了总体成功，该飞行试验控制有效并且记录完整，可作为校准地面测试和转捩预测的有效算例。以上3种技术手段对高超声速边界层的转捩研究各具有独特的优势和明显的不足。近期国外采用直接数值模拟等精细数值模拟方法对HIFiRE5风洞试验模型边界层转捩现象开展了研究，对高超声速边界层横流转捩现象获得了进一步的认识，但存在计算量过大、难以工程应用等问题。

基于雷诺平均方程的转捩模型预测方法能够考虑边界层转捩机制，并具有计算稳定、代价低等优势，在高超声速复杂外形边界层转捩预测中具有广泛应用前景。目前，国内外已发展了多种转捩模型，包括- 模型、--模型以及--模型等，但这些模型主要针对二维边界层流向行波不稳定性引起的转捩现象进行构造，并没有考虑横流效应。对横流诱导转捩的预测主要基于各类横流判据，构造横流转捩模型,这些横流转捩模型在特定流动中取得了较好的效果，但在不同外形中的普适性仍有待提高。此外，横流转捩预测主要以亚声速流动为研究对象，对高超声速中的应用研究相对较少。周玲等基于传统横流雷诺数判据对--转捩模型进行了改进, 基本预测了HIFiRE5椭圆锥表面出现的双肺叶状转捩阵面。Zhang等在- 转捩模型基础上，采用新定义的横流雷诺数，通过有效间歇因子的形式实现横流转捩预测，得到的HIFiRE5实验状态转捩位置与测量结果符合较好。已有横流转捩模型所采用的转捩判据大都需要计算边界层外缘信息以及边界层内横流速度，因此求解过程复杂、计算效率较低，并且非当地量的求解使得该判据在大规模并行计算中受到限制。

本文基于- 模型，通过构造新的横流雷诺数判据，提出一种基于当地流场变量的转捩模型横流效应修正方法，并以椭圆锥外形为研究对象，开展转捩模型横流效应修正方法的分析、验证及应用。

1 转捩模型预测方法

1.1 γ-Reθ t转捩模型

(1)

(2)

式中：间歇因子输运方程的生成项的构造形式为

=()05(1-)

(3)

间歇因子输运方程的耗散项的构造形式为

=(-1)

(4)

转捩动量厚度雷诺数输运方程中生成项 的表达式为

(5)

式(4)和(5)中:为涡量值;为湍流、层流黏性系数比值的函数; 为边界层指示因子;模型常数和 取值分别为50.0和0.03。

为实现对分离诱导转捩现象的模拟预测，以涡量雷诺数与临界动量厚度雷诺数的比值作为增加湍动能生成的表征参数，通过分离间歇因子反映流动分离时出现的湍动能快速增长现象，其构造形式如下：

=

(6)

最终的有效间歇因子为

=max(,)

(7)

本文采用-SST湍流模型与- 转捩模型进行耦合,利用有效间歇因子对原SST湍流模式中湍动能输运方程的生成项和耗散性进行修正：

=

(8)

=min[max(,01),1]

(9)

湍动能比耗散率输运方程与原SST湍流模型中一致。

1.2 横流效应修正

在边界层转捩横流效应分析中，常采用流向涡强度来表示横流的强度。流向涡强度将无量纲的速度矢量与涡矢量点乘以获得流线方向的涡量大小，又称为Helicity参数，记为，其表达式为

=|·|

(10)

式中：

=

从的表达形式可以看到，该参数为正值，对于二维边界层，的值自然退化为0。但为有量纲量，不便于模化横流效应的局部相对特性。考虑到流向涡强度量纲为，与应变率绝对值的量纲一致，本文借鉴涡量雷诺数的构造形式，提出一种新的横流雷诺数构造形式，具体表达式为

(11)

式中：、为当地流场密度及动力学黏性系数;为距壁面距离。从式(11)中可以看到，新构造的横流雷诺数完全基于当地流场变量求解，并且计算过程简单，便于融合于非结构求解框架及大规模并行计算。

借鉴分离诱导转捩修正构造形式，进一步构造横流间歇因子，为便于对比与的作用强度，模型常数3.235仍保留，具体表达形式为

(12)

式中：和为常数。与分离诱导转捩一致，取为常数2。横流效应修正模型参数反映横流效应的作用强度，值越大，横流效应将表现得越为显著，可通过风洞试验数据或精细模拟结果进行标定。

在椭圆锥高超声速绕流流场中，同时存在中心线处的流向分离涡、模型中部横流涡以及二者间的相互干扰，主要体现三维涡结构的效应，而主要反映横流效应。三维剪切流动对边界层转捩的效果可认为由二者中较强的效应来体现。因此，定义间歇因子输运方程的生成项中的有效间歇因子为

=max(,,)

(13)

1.3 数值求解方法

采用基于结构化网格的有限体积法对RANS流动控制方程、湍流模型以及转捩模型输运方程进行耦合求解。采用时间相关法在空间离散后迭代推进得到方程组的收敛解。在空间离散中，流动变量和湍流变量在网格单元界面处的对流通量分别采用Roe格式和Lax-Friedrichs通量分裂格式进行构造，黏性通量计算采用传统的二阶中心差分格式，时间推进采用LU-SGS(Lower-Upper Symmetric-Gauss-Seidel)隐式方法。本文发展的转捩模型预测方法已在高超声速典型流动中得到了验证和应用。

2 横流效应修正分析

2.1 计算模型及工况

以在普渡大学静音风洞开展的椭圆锥风洞试验为研究对象，该试验模型为38.1%比例的HIFiRE5外形，模型长度328 mm，头部鼻尖短轴方向半径为0.95 mm。本文主要针对单位雷诺数10.2×10/m噪声来流工况开展数值模拟研究，风洞试验工况如表1所示。风洞试验噪声条件下来流湍流度达到约3%，本文转捩模型计算中，来流湍流度设为3%，壁面取等温无滑移边界条件，壁温设为300 K。

表1 基准风洞试验工况Table 1 Conditions of benchmark wind tunnel experiment

2.2 网格影响

为分析计算网格对转捩模型预测结果的影响，划分了3套结构化多块网格。初始网格以半模模型为对象，流向、展向和法向网格点数分别为326×151×131，网格单元总量6 175 000，壁面法向第1层网格尺度为0.001 mm。改进网格以1/4模型为对象，流向计算域延伸至350 mm，并在初始网格流场计算结果基础上，对网格点分布进行调整，提高计算网格与激波的匹配性，改进后的网格流向和展向网格点数分别为361和91，法向网格点数及最小网格尺度与初始网格相同。加密网格是将网格流向和展向网格点数增大为1.5倍，网格单元总量达到9 213 750，法向网格分布不变，划分的计算网格如图1所示。

计算网格第1层网格尺度根据网格雷诺数准则及以往转捩湍流模拟经验设定。图2显示了改进后的中等网格下采用- 转捩模型计算得到的无量纲法向网格尺度分布云图，从图中可以看到，模型大部分区域<0.2，已满足边界层转捩湍流模拟<1的要求，表明第1层法向网格尺度设置是合理有效的。

图1 计算网格示意图Fig.1 Schematic of computational grid

图2 y+分布云图Fig.2 Contour of y+

采用未引入横流效应修正转捩模型在3套网格下计算得到的热流分布如图3所示，对比可发现，不同网格得到热流分布特征总体相似，在模型中后部形成三角形的两个高热流区域，二者被中心线低热流条带隔开。随网格加密，中后部区域高热流条带结构愈加清晰，反映出对流场旋涡结构的捕捉更加精细。

图3 不同网格热流分布Fig.3 Heat flux distribution on different grids

典型位置处的热流曲线对比如图4所示，从中可以看到，不同网格计算得到的热流曲线分布趋势一致。沿中心线(=0 mm)热流曲线，改进网格和加密网格结果更为接近；流向截面位置(=300 mm)处，中心区域外侧计算热流曲线近乎重合, 网格尺度差异的影响主要体现在对称中心线附近。这是由于椭圆锥中心线附近形成的汇聚流向旋涡结构丰富、流动空间及时间尺度差异大，网格越密，流动结构刻画越细致。转捩模型对于该处分离流动的失稳转捩判定是基于当地旋涡流动剪切特性的流动信息，这就导致转捩模型对该处流动状态及热流分布的预测出现网格差异，难以达到严格意义的网格无关性。

综合以上结果可以得到，不同网格获得的热流分布形态一致，在中心线区域外侧由改进后的中等网格和加密网格计算得到的热流曲线相互吻合较好，可认为对横流影响区域热流预测达到了网格收敛性。转捩模型预测方法核心是基于时均流场对转捩形态特征进行模化，改进网格已能够反映流场结构及热流分布特征，并且计算量相对较小，在后续计算分析中均采用此网格。

图4 不同网格得到的热流曲线对比Fig.4 Comparison of heat flux curves for different grids

2.3 横流流动特征

高超声速椭圆锥边界层会出现典型的三维流动特征。图5为椭圆锥绕流层流流场，以流向截面无量纲压强等值线和壁面摩阻系数云图及壁面极限流线展示。从图中可以看到，在高超声速来流条件下，椭圆锥模型沿周向激波强度不同，长轴方向波后压强显著大于中心线处，形成的压力梯度使得流动由再附线向中心线汇聚。

图5 椭圆锥绕流层流流场Fig.5 Laminar flowfield around elliptical cone

图6展示了=200 mm截面处的层流流场。从图中可发现，模型侧缘流动向中心线汇聚形成蘑菇状的速度、温度分布形态；在中心线区域外，边界层内部和外部展向速度分量方向相反；横流涡强度在靠近壁面的边界层内数值较大，在模型侧缘肩部附近横流作用最为显著。图7对比了边界层外缘和底层流线分布(上半部分为边界层外缘，下部分为边界层底部)。从中可以看到，在边界层外部，流线基本平行于轴方向；在边界层底部，紧邻中心线的一窄条区域，流线仍平行轴，而在模型中心线和侧前缘间的中部区域，流线方向与轴方向呈一夹角。边界层底层和外缘速度方向的差异进一步反映了边界层流场中存在显著的横流效应。

文献[6-7]风洞试验采用温敏漆测试技术获得了模型表面整体热流分布；此外，前期研究工作采用大涡模拟方法得到了时均热流分布，参考的热流分布如图8所示。这些结果反映了三维高超声速边界层转捩形态特征，可以为横流效应修正的参数标定确立基准。

图6 流向截面流场(x=200 mm)Fig.6 Flowfield on flow direction section (x=200 mm)

图7 边界层外缘及底部流线对比Fig.7 Comparison of streamlines on outer edge and bottom of boundary layer

图8 热流分布参考结果Fig.8 Reference results of heat flux distribution

2.4 横流效应表征参数

利用转捩模型流场计算数据，可得到椭圆锥绕流流场中的横流雷诺数流向截面及空间等值面(=350)分布(见图9)。从中可以看到，新构造的横流雷诺数只在边界层范围内有值，在边界层以外的空间流场中值为0，并且在椭圆锥中心线和侧缘间的中部区域横流雷诺数数值较大，能够反映高超声速三维边界层流动的横流效应。

图9 横流雷诺数分布Fig.9 Distribution of crossflow Reynolds number

与分离诱导转捩修正方法相似，将横流转捩雷诺数与临界动量厚度雷诺数的比值作为判定转捩模型横流效应修正的表征参数。不同比值等值面在空间的分布如图10所示，从图中可以看到，该比值等值面分布于模型中心线和侧缘间的中部区域，随增大，等值面覆盖范围逐渐减小；=1.75时的等值面分布形态与风洞试验获得的转捩形态相似。

图10 不同Recf/Reθc比值等值面Fig.10 Isosurface of different ratios of Recf to Reθc

2.5 横流修正参数标定

从两方面来对横流效应修正参数进行分析和标定：一是利用未加横流修正的转捩预测流场数据，直接计算横流间歇因子，得到横流效应影响作用区域；二是采用不同取值的横流效应修正模型参数开展转捩模型计算，获得模型表面热流分布。结合计算结果与试验数据及精细数值模拟结果对比，最终确定转捩模型横流修正参数取值。

图11给出了基于转捩预测流场数据采用不同横流参数取值得到的横流转捩间歇因子等值面分布。其中认为间歇因子取2时发挥作用，并因 在边界层内部近似为常数0.995，因此将横流转捩间歇因子等值面数值取为1.95。从横流转捩间歇因子等值面分布范围可以看到，取值较小时(=1)，横流间歇因子未达到2，无法体现横流效应；而取值较大时(=9)，等值面覆盖全部侧缘范围，横流效应过为显著；在4～5范围内取值时，等值面位置与风洞试验数据及精细数值模拟结果转捩起始位置接近，其中=45，等值面形态与参考转捩起始位置最为相似。

图11 不同Acf取值横流间歇因子等值面Fig.11 Isosurface of crossflow intermittent factor with different Acf

图12给出了不同取值下横流效应修正转捩模型计算得到的热流分布。从图中可以看到，较小和较大的数值(=1和=9)预测结果与试验数据有明显偏差，较小值不能反映横流转捩影响，而较大值高估了湍流范围；在4～5范围内取值时，热流分布与试验数据及精细模拟结果接近，其中=45能够较清晰地反映出肺叶状热流分布形态，同时侧缘附近未显著升高。因此，综合上述分析结果将横流效应模型参数的取值确定为4.5。

图12 不同Acf取值热流分布云图Fig.12 Heat flux distribution with different Acf

图13为有/无横流效应修正转捩模型计算得到的湍流区域对比(以湍流与层流黏性湍流比值等值面表示，比值取10)。从图中可发现，在模型中后部，引入横流效应修正计算得到的湍流区域更为饱满。结合有/无横流效应修正热流计算结果(图12(c)和图3(b))与参考数据的对比可发现，横流效应修正反映了在中心线及侧缘间的横流转捩现象，提高了转捩模型对高超声速三维边界层转捩的预测能力。

图13 湍流区域对比Fig.13 Comparison of turbulent regions

3 横流效应修正验证

3.1 风洞试验

针对椭圆锥模型风洞试验中不同来流雷诺数及迎角不为0°工况，考虑的风洞试验条件如表2所示。

表2 风洞试验工况Table 2 Conditions of wind tunnel experiments

从不同雷诺数工况热流计算结果(图14)可以看到，在较低雷诺数下(8.1×10/m)，椭圆锥表面边界层转捩起始位置后移，而在较高雷诺数下(11.8×10/m)，转捩位置前移。转捩模型预测结果反映了与风洞试验数据(图15)一致的雷诺数影响规律。此外，从4°迎角下转捩模型计算得到的热流分布(图16)与试验数据(图17)对比可以看到，模型迎风面中心线处出现了明显的高热流条带，表明中心线处的汇聚流动率先发生转捩，转捩模型计算得到的椭圆锥迎风面高热流区域的形态和范围与试验温敏漆数据较为接近。

图14 不同雷诺数工况热流计算结果Fig.14 Computed heat flux distribution under condition of different Reynolds numbers

图15 不同雷诺数工况热流试验结果[6-7]Fig.15 Experimental heat flux distribution under condition of different Reynolds numbers[6-7]

图16 4°迎角热流计算结果Fig.16 Computed heat flux distribution at 4° angle of attack

图17 4°迎角热流试验结果[6-7]Fig.17 Experimental heat flux distribution at 4° angle of attack[6-7]

三维边界层转捩机制复杂，对外部影响因素较为敏感。转捩模型预测得到了椭圆锥边界层总体转捩形态特征，但与试验数据还有一些差异。其原因可能在于数值模拟无法完全再现风洞试验的全部细节信息，包括来流扰动形式、壁面粗糙度、壁面温度非均匀分布等，此外风洞试验磷光图像也存在模型传热以及光线传播受边界层流动状态干扰等因素的影响。转捩模型及横流修正是对转捩现象的模化，仍有进一步改进和发展空间。

3.2 飞行试验

针对HIFiRE5b天上飞行状态，采用发展的转捩模型横流效应修正方法开展转捩预测，考虑515.10 s时刻单位雷诺数为10.0×10/m状态，此时飞行迎角为0.58°，侧滑角为-0.49°，具体飞行工况如表3所示。

计算网格拓扑形式与HIFiRE5风洞试验模型网格一致，针对全模外形，对计算网格进行整体放大，截取流向计算域至900 mm，壁面第1层网格尺度仍设为10m。采用转捩模型计算得到的迎风面热流分布如图18所示，与飞行试验数据(图19)对比可发现，转捩模型计算得到的热流分布形态及数值与飞行试验数据较为接近，预测的转捩起始位置与由实验数据勾勒的转捩阵面吻合较好，在飞行器表面出现了3个叶状的转捩形态，只是在侧滑迎风侧缘区域流动状态存在差异；此外，对称面热流计算结果与飞行试验数据的对比(图20)进一步表明，本文转捩数值模拟得到的热流数据较为准确，有效反映了HIFiRE5b表面出现的转捩现象。

表3 HIFiRE5b飞行试验工况Table 3 Conditions of HIFiRE5b flight experiments

图18 飞行工况热流计算结果Fig.18 Computed heat flow under flight conditions

图19 飞行试验转捩阵面示意图[13]Fig.19 Transition front diagram of flight experiment[13]

图20 热流计算结果与飞行试验数据对比Fig.20 Comparison of computed heat flux with flight experiment data

4 结 论

在- 转捩模型基础上，通过构造横流雷诺数，提出了一种基于当地流场变量的转捩模型横流效应修正方法。以HIFiRE5椭圆锥绕流为对象，在不同风洞试验工况及天上飞行工况进行了验证与应用，得到以下结论：

1) 构造的横流雷诺数能够有效表征高超声速三维边界层横流效应，在此基础上建立的横流间歇因子能够识别出模型中心线汇聚流向涡以及中部横流涡，可作为横流促发转捩的判据。

2) 提出的横流效应修正方法具有一定物理意义，并且完全基于当地流场变量求解，计算过程简单，便于大规模并行计算。

3) 转捩模型横流效应修正方法在不同风洞试验工况和飞行工况得到了验证和应用，获得了与参考结果接近的转捩形态，并有效预测了模型表面的热流分布特征。