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]进一步考虑了热非平衡效应的影响，对比了O₂或N₂单组分的完全气体模型和化学反应模型，研究表明，热化学非平衡效应下，基于Morkovin假设的相关尺度关系仍然能够成立。

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

1 计算模型与方法

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

计算状态示意图见图1。在30 km高度、以马赫数20进行飞行的楔形体，会在头部产生一道斜激波。选择斜激波后的气体状态作为湍流边界层外缘流动状态：马赫数 $ {M_e} = 4.5 $ ，温度 $ {T_e} = 3\;400{\text{ K}} $ ，密度 $\;{\rho _e} = 0.102\;48{\text{ kg/}}{{\text{m}}^{\text{3}}}$ 。此时温度已经足够激发空气发生化学反应，产生非平衡效应。

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

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

边界层厚度和雷诺数见表1。此时的边界层名义厚度 $\delta $ 约为5 mm，基于动量厚度的雷诺数 $R{e_\theta }$ 约为2600。分析所采用流向位置如图2中黑色线所示，且如无特殊说明，下文中分析默认基于此流向截面处的流场。并在下文的表述中约定： $ \bar q $ 代表变量 $ q $ 的雷诺平均，对应的脉动量为 $ q' = q - \bar q $ ； $ \tilde q = {{\overline {\rho q} } \mathord{\left/ {\vphantom {{\overline {\rho q} } {\overline \rho }}} \right. } {\overline \rho }} $ （或 $ < q > $ ）表示其Favre平均，对应的脉动量为 $ q'' = q - \tilde q $ 。

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

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

2 结果与讨论 2.1 平均量分析

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

2.2 Walz关系式

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

$ \frac{{\tilde T}}{{{T_\delta }}} = \frac{{{T_w}}}{{{T_\delta }}} + \frac{{{T_{aw}} - {T_w}}}{{{T_\delta }}}\left( {\frac{{\tilde u}}{{{u_\delta }}}} \right) + \frac{{{T_\delta } - {T_{aw}}}}{{{T_\delta }}}{\left( {\frac{{\tilde u}}{{{u_\delta }}}} \right)^2} $

(1)

其中， $ {T_{aw}} $ 为恢复温度， ${T_{aw}} = {T_\delta }\left[ {1 + r\dfrac{{\left( {\gamma - 1} \right)}}{2}M_\delta ^2} \right]$ ，恢复因子r = 0.9。

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

为了消除热力学状态、化学反应对温度-速度关系的明显依赖，可以引入了一个无量纲参数—恢复焓 $ h_r^* $ ^[12]，并建立恢复焓与速度之间的关系。 $ h_r^* $ 定义为： $h_r^* = \dfrac{{{{\tilde h}_r} - {h_w}}}{{{h_{aw}} - {h_w}}}$ ，其中 ${\tilde h_r} = \tilde h + r\dfrac{{{{\tilde u}^2}}}{2}$ ，恢复因子仍设置为r = 0.9。在壁面处， $ {\tilde h_r} = {h_w} $ ，则 $ h_r^* = 0 $ ；在边界层边缘 $ y = \delta $ ，有 ${\tilde h_r} = {h_\delta } + r\dfrac{{\tilde u_\delta ^2}}{2} = {h_{aw}}$ ，则 $ h_r^* = 1 $ 。图5给出了 $ h_r^* $ 与速度之间的关系。可以看到 $ h_r^* $ 与 $ {{\tilde u} \mathord{\left/ {\vphantom {{\tilde u} {{u_\delta }}}} \right. } {{u_\delta }}} $ 的值接近，并且与文献[12]符合较好。这说明存在如下关系式， $h_r^* = {{( {{{\tilde h}_r} - {h_w}} )} \mathord{/ {\vphantom {{( {{{\tilde h}_r} - {h_w}} )} {( {{h_{aw}} - {h_w}} )}}} } {( {{h_{aw}} - {h_w}} )}} = f( {{{\tilde u} \mathord{\left/ {\vphantom {{\tilde u} {{u_\delta }}}} \right. } {{u_\delta }}}} )$ ，或者

$ \frac{{\tilde h}}{{{h_\delta }}} = \frac{{{h_w}}}{{{h_\delta }}} + \frac{{{h_{aw}} - {h_w}}}{{{h_\delta }}}f\left( {\frac{{\tilde u}}{{{u_\delta }}}} \right) - \frac{r}{2} \frac{{u_\delta ^2}}{{{h_\delta }}}{\left( {\frac{{\tilde u}}{{{u_\delta }}}} \right)^2} $

(2)

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

$ f\left( {\frac{{\tilde u}}{{{u_\delta }}}} \right) = 0.18454{\left( {\frac{{\tilde u}}{{{u_\delta }}}} \right)^2} + 0.83602\left( {\frac{{\tilde u}}{{{u_\delta }}}} \right) $

(3)

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

2.3 强雷诺比拟分析

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

$ \frac{\dfrac{T_{{\text{rms}}}^{''}}{\tilde T}}{{\left( {\gamma - 1} \right)Ma^2\left(\dfrac{ {u_{{\text{rms}}}^{''}}}{\tilde u} \right)}} \approx 1 $

(4)

$ {R_{{u^{''}}{T^{''}}}} \approx - 1 $

(5)

$ {R_{{u^{''}}{v^{''}}}} \approx - {R_{{v^{''}}{T^{''}}}}\Bigg( {1 - \frac{{\overline {v^{''}T_t^{''}} }}{{\overline {v^{''}T^{''}} }}} \Bigg) $

(6)

$ {Pr_t} = \frac{{\overline {\rho {u^{''}}{w^{''}}} \left( \dfrac{{\partial \tilde T}}{{\partial y}} \right)}}{{\overline {\rho {w^{''}}{T^{''}}} \left( \dfrac{{\partial \tilde u}}{{\partial y}} \right)}} \approx 1 $

(7)

图7（a）中黑色曲线给出了采用式（4）得到的温度脉动 $ T_{{\text{rms}}}^{''} $ 与流向速度脉动 $ u_{{\text{rms}}}^{''} $ 之间的关系。由于此关系式没有考虑壁面热流的影响，DNS的结果在0.1～0.5之间，与SRA预测值相差较大。为了考虑壁面温度和热流的影响，一些学者提出了修正的SRA关系式，比如GSRA^[9]、RSRA^[10]、HSRA^[11]，形式如下：

$ \frac{{\dfrac{T_{{\text{rms}}}^{''}}{\tilde T}}}{{\left( {\gamma - 1} \right)M{a^2}\left( {\dfrac{u_{{\text{rms}}}^{''}}{\tilde u}} \right)}} \approx \frac{1}{{c\left[ {1 - \left( {{{\dfrac{\partial {{\tilde T}_t}} { { {\partial \tilde T}}}}}} \right)} \right]}} $

(8)

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

$ T_{\text{rms}}^{''} = - \frac{1}{{P{r_t}}}\frac{{\partial \tilde T}}{{\partial \tilde u}}u_{{\text{rms}}}^{''} $

(9)

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

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

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

$ {{Pr} _{\rho} } = \frac{{\overline {{{u^{''}}{v^{''}}}} \left( {{\dfrac{\partial \bar {\rho} } {\partial y}}} \right)}}{{\overline {{\rho^{ '}}{v^{''}}} \left( {{\dfrac{\partial {\tilde u}} {\partial y}}} \right)}} $

(10)

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

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

2.4 湍动能生成耗散机制

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

$ \tilde k = \frac{1}{2}\frac{{\overline {\rho u_k^{''}u_k^{''}} }}{{\bar \rho }} $

(11)

其输运方程为：

$ \frac{\partial }{{\partial t}}\left( {\bar \rho \tilde k} \right) = - C + P + T + \varPi + D + M + \varepsilon $

(12)

其中 $C = \dfrac{\partial }{{\partial {x_j}}}\left( {{{\tilde u}_j}\bar \rho \tilde k} \right)$ 代表对流项； $ P = - \overline {\rho u_k^{''}u_j^{''}} \dfrac{{\partial {{\tilde u}_k}}}{{\partial {x_j}}} $ 代表生成项； $T = - \dfrac{\partial }{{\partial {x_j}}}\left(\; \overline {\dfrac{1}{2} {\rho u_k^{''} u_k^{''} u_j^{''} }} \; \right)$ 代表扩散项； $ \varPi = - \overline {u_k^{''}\dfrac{{\partial p}}{{\partial {x_k}}}} $ 代表压力做功项； $D = \dfrac{\partial }{{\partial {x_j}}}\left(\; {\overline {u_k^{''}\sigma _{kj}^{'}} } \; \right)$ 代表黏性扩散项； $ M = \overline {u_k^{''}} \dfrac{{\partial \overline {{\sigma _{kj}}} }}{{\partial {x_j}}} $ 表示脉动速度所导致的平均流黏性应力的扩散项； $ \varepsilon = - \overline {\sigma _{kj}^{'}\dfrac{{\partial u_k^{''}}}{{\partial {x_j}}}} $ 是黏性耗散项。 $ \varPi $ 又可以分解为三项之和 $ \varPi = {\varPi _p} + {\varPi _t} + {\varPi _d} $ 。 $ {\varPi _p} = - \overline {u_k^{''}} \dfrac{{\partial \bar p}}{{\partial {x_k}}} $ 代表平均运动中压力所做的功； $ {\varPi _t} = - \overline {\dfrac{{\partial \left( {{p^{'}}u_k^{''}} \right)}}{{\partial {x_k}}}} $ 代表压力输运项； $ {\varPi _d} = - \overline {{p^{'}}\dfrac{{\partial u_k^{''}}}{{\partial {x_k}}}} $ 代表压力-膨胀项。 $ \varepsilon $ 也可以分解为三项之和 $ \varepsilon = {\varepsilon _s} + {\varepsilon _d} + {\varepsilon _i} $ ，其中 $ {\varepsilon _s} = - \overline {\mu \omega _i^{''}\omega _i^{''}} $ 为螺旋耗散项， $ {\varepsilon _d} = - \dfrac{4}{3}\overline {\mu \dfrac{{\partial u_i^{''}}}{{\partial {x_i}}}\dfrac{{\partial u_k^{''}}}{{\partial {x_k}}}} $ 为膨胀-耗散项， $ {\varepsilon _i} = \varepsilon - {\varepsilon _s} - {\varepsilon _d} $ 为与脉动密度的相关项。

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

半当地尺度 $ y_\tau ^* $ ^[11]考虑了当地密度和黏性系数的变化，有助于消除不同工况的影响。用 $ y_\tau ^* = \dfrac{{\bar \mu (y)}}{{\bar \rho (y)u_\tau ^*}} $ 归一化法向距离，并将湍动能输运方程各项通过 $ \;{{\bar \rho u{{_\tau ^*}^3}} \mathord{\left/ {\vphantom {{\bar \rho u{{_\tau ^*}^3}} {y_\tau ^*}}} \right. } {y_\tau ^*}} $ 进行归一化。图11显示，当采用半当地尺度进行归一化后，两种模型的计算结果基本重合，且与文献[12]结果符合较好。这说明采用半当地尺度进行归一化可以一定程度上消除不同工况的影响。

2.5 可压缩效应

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

$ {M_t} = \frac{{{{\left( {\overline {u_k^{''}u_k^{''}} } \right)}^{\frac{1}{2}}}}}{{\bar a}} $

(13)

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

进一步可以通过胀量项来研究可压缩效应的影响，包括压力-胀量项 $ {\varPi _d} $ 和膨胀-耗散项 $ {\varepsilon _d} $ ，这两项都是由于可压缩性导致速度散度不为零引起的。图13给出了压力做功项 $ \varPi $ 中三个分项 $ {\varPi _p} $ 、 $ {\varPi _t} $ 和 $ {\varPi _d} $ 沿壁面法向的分布，并采用生成项P进行归一化。可以看到压力相关项 $ \varPi $ 主要来自于压力输运 $ {\varPi _p} $ 的作用，其余两项所占比例很小。在边界层中，压力-膨胀项的大小不足湍动能生成项的3%。膨胀-耗散项与螺旋耗散项之比 $ {{{\varepsilon _d}} \mathord{\left/ {\vphantom {{{\varepsilon _d}} {{\varepsilon _s}}}} \right. } {{\varepsilon _s}}} $ ，可以用来表征黏性耗散中压缩性部分与不可压缩部分的比例。图14给出了两者之比的分布，可以看到在边界层中 $ 0 < y < 0.8\delta $ 的范围内，膨胀-耗散项仅为螺旋耗散项的1%左右，这也说明此种工况下化学非平衡效应及高马赫数效应带来的可压缩性较弱。

3 结　论

本文对空间发展的高超声速高温湍流边界层进行了直接数值模拟，结果表明：

1）Walz关系式给出的平均温度-速度关系与DNS的结果相差较大；采用无量纲恢复焓与速度之间的关系式能够消除马赫数、化学反应等因素的影响，与DNS结果符合较好。

2）GHSRA可以较好地描述温度脉动与速度脉动之间的关系，强雷诺比拟关系式及其修正在高温化学非平衡湍流边界层中基本上成立。

3）湍动能的生成项、输运项、扩散项和黏性耗散项在湍动能输运过程中占主导地位，且使用半当地尺度归一化各输运项和壁面法向高度时，不同工况下湍动能输运方程各项的分布能够较好地符合。压力-胀量项和膨胀-耗散项的值较小，化学非平衡效应及高马赫数效应引起的可压缩效应有限。