考虑化学非平衡效应的高温湍流边界层统计特性分析

2022-08-23刘朋欣傅亚陆袁先旭

空气动力学学报 2022年4期

刘朋欣，李辰，孙东，傅亚陆，袁先旭，*

(1. 中国空气动力研究与发展中心空气动力学国家重点实验室，绵阳 621000；2. 中国空气动力研究与发展中心计算空气动力研究所，绵阳 621000)

0 引言

Morkovin假设[1]是否成立是可压缩湍流边界层研究的一个重要部分。该假设指出，对于中等马赫数的自由来流（Ma≤5），可压缩湍流与不可压缩湍流的区别可以通过引入流体特性的平均变化来消除。这个理论是van Direst变换的基础，可压缩边界层内的流向速度经过变换后，其分布能与不可压缩边界层较好地符合在一起。较多的试验和数值模拟结果[2-3]都证实了van Direst变换后边界层速度型分布的一致性，且在高超声速湍流边界层中也基本成立[4-7]。另外一个在不可压缩湍流边界层中常用的分析是强雷诺比拟关系（strong Reynolds analogy，SRA），它表征了温度脉动与速度脉动之间的关系。SRA关系式也已经在较为宽泛的马赫数范围和不同的壁面温度下得到了评估和验证[7-11]，并发展了一些修正的SRA关系式，比如GSRA[9]、RSRA[10]、HSRA[11]。

高超声速飞行器在中低空以极高马赫数飞行时，头部会产生较强的激波并压缩来流空气，使其温度急剧升高，导致化学非平衡效应的出现。基于量热完全气体假设的一些湍流边界层经典假设和关系式在高温化学非平衡湍流边界层中是否仍然成立，有待于进一步的验证。Duan等[12]研究了焓值对高超声速湍流边界层的影响，在高焓时考虑空气化学反应模型，而低焓时则使用量热完全气体模型。他们发现，低焓流动中大部分的尺度关系，比如边界层平均速度剖面、SRA关系式，通过消除量热完全气体假设后仍然能够成立，并提出了更通用的关系式GHSRA；高焓效应带来的可压缩效应较小，对湍流结构的影响不明显。Kim等[13]进一步考虑了热非平衡效应的影响，对比了O2或N2单组分的完全气体模型和化学反应模型，研究表明，热化学非平衡效应下，基于Morkovin假设的相关尺度关系仍然能够成立。

但是，目前关于热化学非平衡湍流边界层特性的研究还比较有限。且Duan等[12]在计算时，使用高焓化学非平衡气体模型与低焓完全气体模型对比，没有在同一工况下评估两种气体模型之间的差异；而Kim等[13]虽然在同一工况下对比了两种气体模型，但采用的是O2或N2单组分工质，没有考虑包含多种组分的复杂反应。有必要选择更多的来流工况，更加广泛地研究高温化学非平衡湍流边界层特性。我们之前[14]研究了高温化学非平衡湍流边界层的标度律，并发现经典标度律理论[15]仍然成立。本文进一步研究高温化学非平衡效应对基于Morkovin假设的湍流边界层特性关系式的影响。

1 计算模型与方法

采用含化学反应的三维Naiver-Stokes方程进行直接数值模拟，具体求解过程及算例验证可参考文献[14，16]。无黏通量的离散采用七阶WENO-Z格式[17]，黏性通量的离散采用四阶中心差分格式，时间推进采用显式的三阶TVD（Total Variation Diminishing）型Runge-Kutta法。空气化学反应模型采用5组分（N2，O2，N，O，NO）6基元反应模型[18]，只考虑空气组分的离解。

计算状态示意图见图1。在30 km高度、以马赫数20进行飞行的楔形体，会在头部产生一道斜激波。选择斜激波后的气体状态作为湍流边界层外缘流动状态：马赫数Me=4.5，温度Te=3400K，密度ρe=0.10248kg/m3。此时温度已经足够激发空气发生化学反应，产生非平衡效应。

图1 计算模型状态示意图Fig. 1 Schematic of computational condition

在考虑高温湍流边界层中的化学非平衡效应时，来流组分的质量分数设置为f(N2) = 0.7325，f(O2) =0.2675，记为Multi算例。作为对比，本文还使用了同一来流状态的完全气体模型进行计算，以评估完全气体模型的适用性及化学非平衡效应的影响，并记为Perf算例。

初场采用相同来流状态下SST-RANS[19]（Shear Stress Transfer-Reynolds Averaged Naiver-Stokes）计算结果的剖面生成。湍流入口脉动的生成采用数字滤波合成湍流方法[20]；出口为特征边界条件；上边界设置为初始值，并保持不变；壁面设置为黏性等温，壁温Tw等于Te；展向为周期边界。

边界层厚度和雷诺数见表1。此时的边界层名义厚度 δ约为5 mm，基于动量厚度的雷诺数Reθ约为2600。分析所采用流向位置如图2中黑色线所示，且如无特殊说明，下文中分析默认基于此流向截面处的流场。并在下文的表述中约定：q¯代表变量q的雷诺平均，对应的脉动量为（或）表示其Favre平均，对应的脉动量为。

图2 Multi算例瞬时三维流场和展向中间截面参数分布Fig. 2 Instantaneous 3D flow field and parameter contours in the spanwise middle plane for the Multi case

表1 边界层厚度、雷诺数和网格分辨率参数[14]Table 1 Parameters for the boundary layer thickness, Reynold number and grid resolution[14]

计算域设置为流向约20δ，法向约4δ，展向约2δ。湍流充分发展段的网格分辨率见表1。在之前的研究[14]中，进行了同等工况下的网格无关性验证，最终采用的网格为901×161×301（流向×法向×展向），可以满足直接数值模拟的要求。

图2（a）给出了瞬时流场的Q准则等值面（Qcr=0.05），并用O组分质量分数着色。可见经过一段距离后，湍流边界层得到充分发展。O的生成主要集中在壁面附近温度较高的区域。展向中间截面的瞬时温度和O组分质量分数分布见图2（b）和图2（c），边界层的温度在局部已经可以达到5000 K以上，且壁面为3400 K的冷壁，较大的温度梯度和高温来流条件会使得壁面产生较大热流，使飞行器面临恶劣的热环境；在湍流充分发展区，已经有大量的O2分子发生离解反应并生成O原子，O原子的质量分数最大已经超过0.12，这表明此时发生了较强的化学非平衡效应。图2（d）给出了展向中间截面的瞬时无量纲密度梯度大小，也可以清楚地看到，本文分析所采用的流向位置处的湍流边界层已经充分发展。

2 结果与讨论

2.1 平均量分析

课题组前期工作中[14]分析了相同计算设置下高温非平衡湍流边界层的平均速度剖面、温度剖面和雷诺应力强度。这里只给出了平均温度剖面，见图3，以方便后续分析讨论。研究发现平均速度剖面仍然存在明显的对数区，对数区内斜率基本保持不变（为1/k，k为0.41），截距有所升高C= 6.2。相较于完全气体模型，考虑空气化学非平衡效应会使得边界层内温度显著降低，这是由于空气发生了较强的吸热离解反应，使得组分内能转换为化学能。采用摩擦速度无量纲化的雷诺应力分布与文献中超声速湍流边界层[21-22]、高焓湍流边界层[12]结果一致。这说明采用归一化方法后，高温化学非平衡湍流边界层中的平均速度剖面和雷诺应力仍然与完全气体具有相似的分布。

图3 平均温度的法向剖面[14]Fig. 3 Wall-normal profiles of the averaged temperature[14]

2.2 Walz关系式

对于零压力梯度的平板边界层，Walz提出的修正Crocco关系式[23]（称为Walz关系式）来建立平均温度和平均速度间的关联，具体表达式为：

其中，Taw为恢复温度恢复因子r= 0.9。

图4给出了DNS得到的平均速度-平均温度的关联曲线与Walz关系式之间的对比，可以看到两者存在很大的差异，尤其是对于空气化学反应模型得到的结果差异更大，这说明高温化学非平衡效应对平均温度-速度的关联性有相当大的影响。这可能是由于高温气体效应的吸热反应对平均温度分布有明显的降低作用，由此导致边界层内温度分布发生显著改变。

图4 Walz关系式Fig. 4 Walz relationship

为了消除热力学状态、化学反应对温度-速度关系的明显依赖，可以引入了一个无量纲参数—恢复焓并建立恢复焓与速度之间的关系。定义为：其中恢复因子仍设置为r= 0.9。在壁面处，在边界层边缘y=δ，有则。图5给出了与速度之间的关系。可以看到与的值接近，并且与文献[12]符合较好。这说明存在如下关系式，或者

图5 无量纲恢复焓与速度之间的关系Fig. 5 Relation between the non-dimensional recovery enthalpy and the velocity

经过对DNS数据的拟合可得：

将式（3）代入式（2），并给出修改后的平均静焓与速度之间的关系（图6）。可见采用式（2）得到的曲线与DNS数据吻合较好。这说明在边界层内采用恢复焓，并使用拟合得到的关系式去修正平均静焓与速度的关系式是有效的，可以消除自由来流的马赫数、壁温、化学反应等因素的影响，且对不同工况具有一定的普适性。

图6 修正后的∼ 关系式Fig. 6 Modified ∼relation

2.3 强雷诺比拟分析

可压缩边界层内各参数脉动量之间的关系也很重要，这可以通过SRA来表征。MorKovin[1]在1962年提出了5个强雷诺比拟关系式，其中4个为：

对应于上述三种改进，c=1.0、c=1.34、c=Prt。其中Ma为当地马赫数，为总温。前期的研究中显示[2,7,24]，HSRA对低焓情况下的绝热或非绝热湍流边界层都有较好的适用性。以上三种修正是基于量热完全气体假设，另外还有一种修正关系式GHSRA[12]，其消除了完全气体假设并适用于考虑化学非平衡的工况，具体为：

图7（b）给出了基于空气化学反应模型算例并采用四种不同SRA修正关系式的分布，可以看到考虑了化学非平衡过程的GHSRA的值基本上在1附近，表现最好。靠近壁面附近的突变和符号的改变是由于突变附近的温度梯度等于0引起的，见平均温度分布图（3）。图7（a）还对比了采用原始SRA公式（4）与GHSRA公式（9）时，两种不同气体模型的温度脉动-流向速度脉动关系式的分布，可以看到两种气体模型下GHSRA的值都基本上在1附近，这也说明此式在消除量热完全气体假设时的有效性。

图7 温度脉动与流向速度脉动的强雷诺比拟关系Fig. 7 Strong Reynolds analogy relation between the temperature fluctuation and the streamwise velocity fluctuation

图8 给出了流向速度脉动与温度脉动之间的关系式Ru′′T′′、法向速度脉动与温度脉动之间的关系式Rv′′T′′、流向速度脉动与法向速度脉动之间的关系式Ru′′v′′的分布。可以看到在边界层外区，u′′和T′′并不是完美的反相关性，Ru′′T′′大部分为−0.6～−0.7之间，这与文献[2，12，24-26]计算结果一致。式（6）表明Ru′′v′′与Rv′′T′′呈现相反的相关性，图8中两者的分布曲线证实了这一关系式的正确性。

图8 不同参数脉动量之间的雷诺比拟关系Fig. 8 Strong Reynolds analogy relation between the fluctuations of different parameters

式（7）给出了湍流Prandtl数的定义，它表征了平均运动中动量交换与热交换之比。强雷诺比拟中认为该值近似于1。图9中给出了两种气体模型得到的Prt的分布，可以看到Prt约在0.8附近波动。这与Pirozzoli等[3]、Duan等[8]和Li等[27]计算较高来流马赫数的工况结果一致。另外，图9中还给出了表征湍流质量扩散的特征参数Prρ的分布，其定义为：

图9 Prandtl数的法向分布Fig. 9 Wall-normal distributions of the Prandtl number

Pirozzoli等[3]认为Prt≈Prρ，即在Morkovin假设下，湍流边界层中动量交换、热量交换和质量交换过程处于近似相当的状态。图9的结果表明，在边界层中的大部分区域，采用完全气体模型时，Prρ的大小与Prt基本相当。但当考虑高温化学非平衡效应而采用空气化学反应模型时，Prρ与Prt之间的差异稍大；且越靠近边界层外层，两者的差异越明显，说明三种交换过程的相对大小发生了变化。

以上分析表明，当前计算的工况下的高温化学非平衡湍流边界层中，强雷诺比拟关系式及其修正基本上仍然能够成立。

2.4 湍动能生成耗散机制

湍动能输运方程常被用来分析湍动能生成、耗散等机制。湍动能的定义为：

其输运方程为：

图10给出了湍动能输运方程各项沿壁面法向的分布，并采用壁面参数进行归一化，其中。从图10（a）可以看到，高温化学非平衡湍流边界层内湍动能的生成项P、输运项T、扩散项D和黏性耗散项ε占主导地位，其他各项的值都比较小。图10（b）对比了两种气体模型的结果。在采用壁面尺度归一化的法向坐标y+下，考虑化学非平衡效应时得到的生成项P、输运项T、扩散项D的峰值更靠近于壁面。这可能是由于两种气体模型计算得到的壁面附近参数不同所导致的。

图10 壁面尺度下湍动能输运方程各项的法向分布Fig. 10 Wall-normal distributions of the turbulent kinetic energy budget terms scaled in wall units

图11 半当地尺度下湍动能输运方程各项的法向分布Fig. 11 Wall-normal distributions of the turbulent kinetic energy budget terms scaled in semi-local units

2.5 可压缩效应

本文计算的工况马赫数较高，且考虑了化学非平衡效应，因此有必要分析可压缩效应的影响。湍流马赫数Mt常用来表征可压缩效应的大小[26-27]，其定义为：

且通常认为Mt＞ 0.3时就要考虑可压缩效应对湍流特性的影响。图12给出了Mt在边界层内的分布。可以看到，在近壁区域，两种气体模型得到的Mt均大于0.3，且空气反应模型的值略大。这可能是由于化学反应使得Multi算例中气体的比热降低、声速降低导致。

图12 湍流马赫数的法向分布Fig. 12 Wall-normal distributions of the turbulent Mach number

进一步可以通过胀量项来研究可压缩效应的影响，包括压力-胀量项 Πd和膨胀-耗散项 εd，这两项都是由于可压缩性导致速度散度不为零引起的。图13给出了压力做功项 Π中三个分项 Πp、 Πt和 Πd沿壁面法向的分布，并采用生成项P进行归一化。可以看到压力相关项 Π主要来自于压力输运 Πp的作用，其余两项所占比例很小。在边界层中，压力-膨胀项的大小不足湍动能生成项的3%。膨胀-耗散项与螺旋耗散项之比 εd/εs，可以用来表征黏性耗散中压缩性部分与不可压缩部分的比例。图14给出了两者之比的分布，可以看到在边界层中 0