JFM：湍流高超声速边界层中的热化学非平衡效应

通过直接数值模拟，对马赫数 12.48、壁面冷却条件下的高超声速空间发展湍流边界层进行了分析。在选定的工况下，大量动能向内能的转换引发了热非平衡和化学非平衡现象。空气被假定为含五种组分的反应混合物，并采用双温模型来描述振动非平衡特性。壁面冷却在一定程度上抵消了摩擦加热的影响，边界层内的温升激发了振动能量模式，同时导致氧气发生轻微的化学离解。振动非平衡主要由分子氮驱动，其弛豫速率慢于混合物中的其他分子。结果表明，湍流混合维持了热非平衡状态，导致壁面附近形成振动欠激发状态，而边界层外层则出现振动过激发状态。通过定义相互作用指标，量化了湍流与热效应之间的强耦合关系。基于振动湍流普朗特数的定义，提出了一种振动能量湍流通量的建模策略。此外，还评估了热非平衡条件下强雷诺类比的有效性。强可压缩性效应促进了平动 - 振动能量交换，但与无约束湍流构型的观测结果不同，在膨胀 / 压缩过程与振动过激发 / 欠激发状态之间未发现明显的相关性。

关键词：可压缩边界层、高超声速流动、湍流反应流

       

高速流动是众多构型研究中极具挑战性的课题，涵盖进入行星大气层的物体以及大气层内高超声速飞行等场景（Gnoffo 等人，1999）。在这类流动中，大量动能向内能的转换会导致流场温度急剧升高。高温引发的效应包括化学反应和振动弛豫现象，其特征时间尺度与流动时间尺度相当。此类条件下产生的复杂热化学状态，可能对高速飞行器设计中的关键参数产生显著影响（Bertin 与 Cummings，2006；Candler，2019）。过去，针对高超声速条件下的热化学效应，研究多集中于相对简单的构型。多数研究聚焦于驻点流动，这与钝体再入大气层的场景密切相关（Fay 与 Riddell，1958；Armenise 等人，1996；Bonelli等人，2017；Colonna、Bonelli 与 Pascazio，2019）。在再入轨迹的后期阶段，或在自由流密度更高的亚轨道飞行条件下，流动描述因层流向湍流的转捩而进一步复杂化。这种转捩可能由自由流扰动、热防护系统（TPS）的侵蚀与烧蚀，或表面缺陷引发，导致壁面摩擦阻力和热流急剧增加。因此，层流 - 湍流转捩已被确认为精确预测高超声速飞行物体周围气动热流场的核心问题。文献中已有大量研究探讨化学非平衡条件下平板边界层的线性和非线性稳定性，部分研究甚至涉及转捩的初始阶段。Malik 与 Anderson（1991）、Hudson、Chokani 与Candler（1997）以及 Perraud 等人（1999）的开创性研究指出，化学反应倾向于促进所谓的 “第二模态不稳定性”，这与强冷却壁面的作用类似。近年来，Marxen 等人（2013）以及 Marxen、Iaccarino与 Magin（2014）通过直接数值模拟（DNS），评估了在存在化学反应的情况下，大振幅或小振幅波所引发的最大流向速度扰动。Kline、Chang 与 Li（2019）以及 Knisely 与 Zhong（2019）研究了热化学非平衡（TCNE）对边界层稳定性的影响；烧蚀的影响也被纳入研究范畴（Mortensen 与 Zhong，2016；Miró Miró 与 Pinna，2021）。

另一方面，尽管多年来人们对低焓条件下（即空气表现为量热完全气体）的壁面约束湍流进行了细致研究，但热化学非平衡条件与湍流之间的相互作用却鲜少受到关注。研究人员采用不同的湍流注入技术和壁面处理方法，对可压缩湍流边界层（TBL）构型展开了研究；关于自由流马赫数相对较低（M∞≤3）的湍流边界层直接数值模拟的详细综述，读者可参考 Wenzel 等人（2018）的研究及其中的参考文献。众多研究致力于将针对不可压缩构型推导的标度律推广至可压缩构型（Morkovin，1962）。Duan、Beekman 与 Martín（2011）率先通过高保真模拟对高超声速 regime 展开研究，他们对名义自由流马赫数高达 12 的时间发展湍流边界层进行了直接数值模拟，以证明即使在高速流动中，Morkovin 假设依然具有稳健性。Franko 与 Lele（2013）对马赫数 6 下的转捩空间发展边界层进行了直接数值模拟，而 Zhang、Duan 与Choudhari（2018）则构建了冷却边界层的数据库，其自由流马赫数高达 14，流动条件与超高速风洞中的情况具有代表性。Huang 等人（2020）对零压梯度高超声速湍流边界层进行了直接数值模拟，并对不同雷诺平均纳维 - 斯托克斯（RANS）模型的性能进行了先验和后验分析。

仅有少量研究关注高焓流动，即自由流条件下会引发离解与复合反应的流动（Martin 与 Candler，2001；Duan 与 Martín，2009）。Duan 与 Martín（2011b）对低焓和高焓、具有反应性的时间发展边界层开展了对比研究，发现两者极为相似，因为在低焓条件下验证有效的标度律在高焓应用中仍然适用，或可推广至高焓场景。

近期，Passiatore 等人（2021）以及 Di Renzo 与Urzay（2021）研究了有限速率化学反应与转捩及湍流空间发展边界层的相互作用，前者重点分析了化学非平衡对拟绝热湍流边界层的影响，后者则探讨了其对冷却湍流边界层的作用。在这两项研究中，热非平衡效应均被忽略，因为所选的边界条件预计不会引发显著的振动弛豫现象。迄今为止，关于这些过程对湍流流动行为的影响，仅在无约束流动构型中进行过研究，例如各向同性湍流衰减（Neville 等人，2014；Khurshid与 Donzis，2019；Zheng 等人，2020）和混合层（Neville等人，2015；Fiévet 等人，2019）。

本研究旨在通过对同时处于热非平衡和化学非平衡条件下的湍流边界层进行研究，填补这一知识空白。为此，我们对平板边界层进行了直接数值模拟，该边界层处于 6° 尖楔以马赫数 20 飞行时的激波后条件。边界处的热力学状态使得边界层内的振动弛豫现象显著，同时伴有轻微的化学活性。

本文结构如下：第 2 节介绍控制方程、数值方法、模拟概况及主要参数；第 3 节呈现数值结果，重点关注完全湍流状态，阐述动态场的行为并对比流向速度的不同转换形式，随后通过强调与湍流混合机制的强耦合，分析热非平衡效应，还评估了强雷诺类比得出的相关性以及雷诺平均纳维 - 斯托克斯（RANS）方程中新增项的经典封闭模型，最后研究了可压缩性效应与热非平衡之间的耦合关系；第 4 节得出结论。

2、高速实验深度强化学习（DRL）控制框架

2.1. 控制方程与热力学模型

本研究中所考察的流体是高温空气，建模为包含 N₂、O₂、NO、O 和 N 五种组分的混合物。当考虑处于热非平衡条件下的气体时，混合物中分子的振动能级会被部分激发，且不再与转动 - 平动能级处于平衡状态。其直接后果是，即便假设粒子布居遵循玻尔兹曼分布，使用单一静温来描述所有能量模式的做法也不再适用（Anderson，2006）。

文献中处理此类条件的经典方法被称为多温度模型，其核心是通过为每种分子引入额外的 “振动” 温度，单独考虑振动能级。为使需积分的方程数量保持在合理范围内，本研究采用 Park（1988）提出的双温模型进行计算。该模型在以往研究中被广泛应用（如 Hudson 等人，1997；Franko、MacCormack 与 Lele，2010；Bitter 与 Shepherd，2015），其假设每种分子的振动能态遵循玻尔兹曼分布，且混合物中所有双原子组分（即 N₂、O₂和 NO）共用一个振动温度 Tᵥ

因此，还需要一个额外的总振动能量eV守恒方程。因此，这类流动的特性由多组分化学反应和热松弛气体的可压缩纳维 - 斯托克斯方程控制，其表达式如下：

在上述公式中，ρ为混合物密度，t为时间坐标，xj为笛卡尔坐标系中第j方向的空间坐标，uj为同一方向上的速度矢量分量，p为压力，δij为克罗内克符号，τij为粘性应力张量，其模型表达式为：

其中，μ为混合物的动力粘度。在式（2.3）中， qjTR和qjV分别为转动 - 平动和振动对热流的贡献；unjD表示扩散速度，hn为第n种组分的比焓。在组分守恒方程（2.4）中，ρn=ρYn为第n种组分的分密度（Yn为质量分数），ω˙n为第n种组分的生成率。分密度之和等于混合物密度，即ρ=∑n=1NSρn，其中NS为总组分数。为保证总质量守恒，求解混合物密度和NS−1个组分守恒方程，第NS种组分的分密度通过ρNS=ρ−∑n=1NS−1ρn计算得出。在后续研究中，我们将分子氮设为该组分，因为其在整个计算域内含量最高。在式（2.5）中，eV为混合物的振动能，其表达式为：

其中，eVm 为分子的振动能，NM 为分子的总数。在同一方程中表示振动模式与平动模式之间的能量交换（这是由分子碰撞引起的，且与能量弛豫现象相关），而 ∑m=1NMω˙meVm 表示由分子消耗或生成导致的振动能量损失或增加。

假设每种物质均表现为热完全气体；根据道尔顿分压定律，可得到热状态方程。

   

其中，Rn和Mn分别为第n种物质的气体常数和分子量，为通用气体常数。高温空气物质的热力学性质计算需考虑平动、转动和振动模式的贡献；具体而言，内能可表示为：

这里，hf,n0 是第n种物质在参考温度（Tref=298.15 K）下的生成焓，cp,nT 和 cp,nR 分别是第n种物质的等压热容中平动和转动模式的贡献，其计算方式为：

其中，eVn 为第n种物质的振动能，其表达式为：

其中，θn为每种分子的特征振动温度（对于N2、O2和NO，分别为3393 K、2273 K 和 2739 K）。对守恒方程进行数值积分后，直接根据比内能（不含振动贡献部分）计算转动 - 平动温度T，而振动温度TV则采用牛顿 - 拉夫逊迭代法根据式（2.7）计算。

热通量通过傅里叶定律建模，即qRjT=−λTR(∂T/∂xj)和qjV=−λV(∂TV/∂xj)，其中λTR和λV分别为转动 - 平动热导率和振动热导率。纯物质的黏度和热导率分别采用 Blottner、Johnson 和 Ellis（1971）提出的曲线拟合公式以及 Eucken 关系式（Hirschfelder & Curtiss 1969）计算。相应的混合物性质通过 Wilke 混合法则（Wilke 1950）进行评估。质量扩散通过菲克定律建模。

其中，右侧第一项表示有效扩散速度，第二项是质量修正项；在处理物种扩散系数非恒定的情况时，为满足连续性方程，需考虑这一修正项（Poinsot & Veynante, 2005）。具体而言，Dn是物种n在混合物中的等效扩散系数，其计算遵循赫希菲尔德近似（Hirschfelder & Curtiss, 1969），以 Gupta 等人（1990）通过曲线拟合得到的二元扩散系数为基础。需要注意的是，式（2.12）中忽略了分子量梯度的影响，因此它代表了一种相对简化的模型（尽管该模型允许质量扩散系数可变且刘易斯数非恒定）。这五种物种通过包含五个可逆化学反应步骤的反应机理相互作用（Park, 1990）。

其中，M为第三体（所考虑的五种物质中的任意一种）。解离和复合过程由反应R1、R2和R3描述，而置换反应R4和R5则代表重排过程。第n种物质的质量生成率由质量作用定律控制。

其中，νnr和νnr′分别为第r个反应中第n种物质作为反应物和产物的化学计量系数，NR为反应的总数。此外，kf,r和kb,r分别表示第r个反应的正向和逆向反应速率，其建模采用阿伦尼乌斯定律。由于此类流动中化学非平衡与热非平衡同时存在，二者存在强耦合关系，这通过对计算反应速率时所用的温度值进行适当修正来体现。对于式（2.13）中的解离反应R1、R2和R3，采用几何平均温度。 

其中，τm 是第 m 种分子的弛豫时间，其计算采用 Millikan 和 White（1963）提出的表达式。具体而言，第 m 种分子相对于第 n 种物质的弛豫时间可表示为：

其中，p为压力，patm=101325 Pa，amn 和 bmn 为 Park（1993）中给出的系数。由于该表达式在温度高于 5000 K 时往往会低估实验数据，因此 Park（1989）提出了一种高温修正方法。

2.2.数值方法

通过高阶中心有限差分格式对纳维 - 斯托克斯方程进行数值积分（Sciacovelli 等人，2021）。对对流通量采用十阶中心差分进行离散，并辅以高阶自适应人工耗散项。该人工耗散项由两部分组成：一部分是基于守恒变量十阶导数的九阶精度耗散项（用于抑制网格尺度振荡），另一部分是低阶激波捕捉项。采用一种高选择性传感器 —— 基于杜克罗（Ducros）对詹姆森（Jameson）压力传感器的扩展（Ducros 等人，1999），在流动不连续区域的紧邻范围内启动激波捕捉，适用于除振动能量方程之外的所有方程。对于振动能量方程，为确保对虚假振荡的适当抑制，优先采用基于振动温度二阶导数的激波传感器。粘性通量采用标准四阶差分格式。时间积分采用三阶总变差减小（TVD）龙格 - 库塔格式（Gottlieb与 Shu，1998）。关于当前数值技术的更多细节，以及其在包括化学反应高超音速边界层在内的多种高可压缩流动问题中的完整评估，可参见 Sciacovelli 等人（2021）的研究。      

2.3.计算设置

研究的构型为空间发展的零压力梯度平板边界层，如图1 所示。给定的边界条件为：来流马赫数Me=12.48、来流温度Te=594.3 K、来流压力pe=4656 Pa，这些条件模拟了 6° 尖楔在马赫数M=20、海拔约 36km 飞行时产生的激波下游状态。边界层边缘的滞止焓为，该值与 Duan 和 Martín（2011b）以及 Di Renzo 和 Urzay（2021）所研究的高焓案例相当。在此来流条件下，空气被假设处于热化学平衡状态（XN2=0.79，XO2=0.21，其中Xn为第n种物质的摩尔分数）。值得注意的是，Kline 等人（2019）在稳定性研究中也考虑了类似场景。计算域为激波层内的矩形区域，在图 1 中以红色突出显示。域的范围为Lx×Ly×Lz=3000δin×120δin×30πδin，其中x、y、z分别表示流向、壁法向和展向，δin为入口截面处边界层的位移厚度，其计算公式为δ=∫0δ(1−ρu/(ρeue))dy（其中δ为达到来流速度 99% 处的边界层厚度）。计算网格为Nx×Ny×Nz=9660×480×512，总网格点数约为 24 亿。流向和展向的网格间距均匀，而壁法向采用 0.7% 的恒定网格拉伸率。计算从平板前缘距离x0处的层流区域开始。其轮廓为

入口处给定的守恒变量是通过求解二维化学非平衡且振动平衡边界层的局部自相似理论生成的（Sciacovelli 等人，2021）；基于入口位移厚度的入口雷诺数为Reδin=6054。采用热平衡假设简化了数值设置，且预计不会改变后续分析中关注的湍流区的定性行为。在入口边界下游设置了一个海绵层，范围直至(x−x0)/δin=20，以防止边界层发生突变。

顶部和右侧边界采用基于特征线的边界条件，展向设置为周期性条件。假设壁面为非催化壁面（即∂Yn/∂y=0）且为等温壁面，温度T=TV=1800 K。第一个条件意味着壁面不参与化学过程。这是对实际飞行条件下情况的理想化处理：飞行器热防护系统（TPS）的材料可能具有催化性，从而促进近壁区域混合物中原子的复合。有限速率催化作用的研究超出了本文的讨论范围，但在未来的研究中，探究壁面催化作用对热应力的影响可能具有重要意义（Bonelli、Pascazio 和Colonna，2021）。另一方面，壁面热平衡的假设主要是由于缺乏关于振动温度最适用条件的认识。在先前对空心圆柱 flares 周围层流的研究中（Kianvashrad 和 Knight，2017，2019），振动能量既采用了绝热边界条件，也采用了等温边界条件，在传热方面得到了相似的结果。此外，过去考虑热非平衡的边界层稳定性研究均假设壁面处为热平衡状态（Bitter 和 Shepherd，2015；Kline 等人，2019；Knisely和 Zhong，2019）。遗憾的是，目前尚无关于湍流流动的相关信息。尽管如此，振动传热对总壁面热通量的贡献相对于转动 - 平动传热而言较小（如第 3 节后续详述），这为将狄利克雷边界条件作为一阶近似的选择提供了后验合理性。另一方面，所选的壁面温度值代表了实际飞行数据（Park，2004）。

通过一种 suction-blowing 策略诱导转捩，该策略与 Passiatore 等人（2021）采用的策略相似。具体而言，以下壁面法向速度在靠近入口的壁面条带上施加了以下壁面法向速度：

其中，g(x)=(x−xforc)/(δin2σ)，且2σ=0.85(2π/ω)，xforc是调制强迫函数的类高斯分布中心。强迫位置处的雷诺数为Reδforc=104，对应xforc=7.4×10−2 m。在式（2.19）中，无量纲脉动ω=ω~δin/c∞（c∞为来流声速）对应的有量纲频率为f~=ω~/2π=75 kHz，β=0.4/δin为展向波数。最后，A为强迫振幅，已设定为 5%。吸吹注入条带中存在强振幅的定常模式，这模拟了粗糙点的影响，用于在规定计算域内加速流动向湍流的破裂过程。

2.4.数据收集与分析

在下文呈现的结果中，待初始瞬态过程消除后，通过在时间和展向均匀方向上进行平均来计算流动统计量。本文将呈现并讨论多种流动量的一阶矩和二阶矩。对于给定变量f，我们用fˉ=f−f表示标准的时间和展向平均，其中f为相应的脉动值；而f~=f−f表示密度加权的法夫雷平均，其中f为法夫雷脉动值，且f~=ρf/ρˉ。采样时间间隔是恒定的，恰好对应于每个强迫谐波周期 300 个样本；具体而言，Δtstats+=Δtstats(uτ2ρw/μw)=0.16，其中uτ=τw/ρw是基于计算域末端平均壁面剪应力τw的摩擦速度。数据收集时间超过三个周转时间，对应于Tstats(uτ2/(μw/ρw))≈6700，共计约 40000 个时间快照。除流动统计量外，还分别以强迫函数基频的 20 倍和 40 倍频率提取瞬时平面和特定网格线，样本数量分别为 3200 个和 6400 个。

分析将主要集中在五个选定的流向位置：一个位于层流区域（用于对比），一个位于转捩区域，最后三个位于域内的湍流部分。表 1 列出了选定位置处基于距前缘距离的雷诺数Rex=ρeuex/μe以及一些边界层特性。在完全湍流区域，基于当地动量厚度的雷诺数Reθ=ρeueθ/μe（其中θ=∫0δ(ρu/(ρeue))(1−ρu/(ρeue))dy）达到接近 6000 的值，对应的Reθinc=Reθμe/μw约为 3000。基于位移厚度的雷诺数Reδ=ρeueδ/μe和摩擦雷诺数Reτ=ρwuτδ/μw分别达到约16×104和 1100。值得注意的是，网格间距确保了各处均具有类直接数值模拟（DNS）的空间分辨率。此处，符号 “・+” 表示采用粘性长度尺度lv=μw/(ρwuτ)进行归一化。除非另有说明，统计量的壁法向演化均以内部半局部单位y+=ρˉuτy/μˉ展示，其中uτ=τw/ρˉ。

3、结果

3.1.整体流动特性      

首先讨论壁面处选定物理量的流向演化。图 2（a）和图 2（b）分别给出了摩擦系数Cf以及热流系数Cq和CVq的分布，其定义为：

总净热通量由这两部分贡献之和给出。从(x−x0)/δin≈1100处开始的上升段导致摩擦系数Cf曲线出现急剧超调，壁面剪应力峰值几乎是其相应层流值的五倍。流动从(x−x0)/δin≈1800处开始进入湍流状态，摩擦系数Cf平稳下降直至计算域末端。结果表明，壁面热通量的转动 - 平动贡献CTRq远大于振动贡献CVq（为匹配图 2b 的范围，振动贡献已乘以 10），这主要是因为振动热导率λV的值约比λTR小一个数量级。尽管CTRq紧随摩擦系数曲线变化，但CVq却呈现出不同的演化趋势：在强迫条带之后，它迅速下降，并在流动破裂为湍流之前基本保持恒定；此后，湍流区域的CVq值几乎是层流状态下的两倍，并在整个计算域末端保持恒定。热通量的这种不同趋势与边界层内振动温度梯度的演化有关，这将在 3.3 节中详细讨论。

为了分离平均场和脉动场对表面摩擦的贡献，我们采用了 Renard 和 Deck（2016）提出的分解方法（后来 Li 等人在 2019 年将其推广到可压缩流动）。值得注意的是，该分解方法的基本假设使其即使在存在热非平衡和化学非平衡效应的情况下也能直接应用。因此，表面摩擦系数可以重新表示为：

其中，Cf,1、Cf,2和Cf,3分别代表平均场分子耗散、湍流耗散以及与边界层空间增长相关的效应。这三项的总和及其各自的贡献如图 3 所示。在类层流区域和完全湍流区域，计算结果与直接计算得到的Cf吻合极佳，仅在转捩区域存在微小偏差。Cf,3的贡献在各处均可忽略不计，但在流动向湍流破裂的位置除外 —— 此处由于边界层的突然增厚，Cf,3呈现负值。Cf,3的减小被雷诺应力相关项的大幅增加所抵消，而平均场贡献项Cf,1的增长则略滞后于其他两项。在湍流区域，Cf,1和Cf,2几乎重合，这与 Passiatore 等人（2021 年）的观测结果不同，他们的研究中前者占主导地位。这种差异可归因于当前构型中达到的更大摩擦雷诺数，从而导致湍流贡献增强（Fan、Li 和 Pirozzoli，2019 年）。      

3.2.平均流动分析

针对平均流向速度剖面，已测试了多种尺度变换方法。首先是范德瑞斯特（Van Driest，1956）和特雷特尔与拉尔森（Trettel & Larsson，2016）提出的变换方法。      

分别在图 4（a）和图 4（b）中展示了三个湍流位置的这两种尺度变换结果。前者的尺度变换已被证明在绝热边界层中效果较好，即使在

高速流动中（Passiatore 等人，2021），前者的尺度变换已被证明效果较好。但另一方面，当考虑强冷却构型时，无论是边界层（Zhang 等人，2018；Huang 等人，2020）、管流（Ghosh、Foysi和 Friedrich，2010）还是槽道流（Modesti 和 Pirozzoli，2016；Sciacovelli、Cinnella 和 Gloerfelt，2017），这种尺度变换都会变得不准确。图 4（a）证实了这一趋势 —— 线性区和对数区均相对于解析定律出现偏移。值得注意的是，图 4 中对数区由(1/κ)logy++C和(1/κ)logy∗+C描述，其中κ=0.41，C=5.2。

如预期所示，Trettel 和 Larsson 的半局部尺度变换改善了近壁区的预测结果，因为它明确考虑了整个内层的应力平衡条件。然而，在对数区观察到较大的离散性，其雷诺数依赖性与 van Driest 尺度变换类似。尽管这种尺度变换在内流中效果尚可，但在外部构型中表现欠佳，这很可能是由于与尾迹区的相互作用 —— 图 4（b）中显示尾迹区被过度拉伸。

最近，Griffin、Fu 和 Moin（2021）提出了一种新的基于总应力的变换方法（以下称为 Griffin-Fu-Moin 尺度变换，uGFM），该方法基于恒应力层假设，表达式为：

其中，μ+表示平均粘度相对于其壁面值的归一化处理。这种尺度变换已被证明能够成功整合槽道流、管流和边界层构型的结果，即便在大马赫数情况下也是如此。图 4（c）展示了在Reθ=6200处提取的速度剖面中，uTL与uGFM随y的变化对比。在内层，粘性应力占主导地位，且大致等同于总剪切应力，因此两种尺度变换的结果均能很好地重合。相反，在对数区，湍流剪切应力成为主导因素，而基于总应力的尺度变换能最准确地考虑这一影响。这使得uGFM与通用对数剖面的重合度优于uTL，后者对数区的斜率被大幅高估。不过，与经典不可压缩流动的数值相比，uGFM变换仍预测出对数区具有更大的斜率和截距，这与 Lee、Martin 和Williams（2021）近期的研究结果一致。最后，我们验证了基于恒应力层假设的变换（即边界层内τ/τw≈1，从而得到式（3.4））与基于实际总剪切应力（从 DNS 数据中提取）的变换之间不存在任何可分辨的差异，这证实

恒应力层假设的有效性。化学活性和热弛豫过程都不会显著改变这些变换的有效性，因为正如 Passiatore 等人（2021）以及 Di Renzo 和 Urzay（2021）先前所观察到的，热化学过程与湍流活动之间存在一定程度的解耦。

图 5 展示了在湍流位置Reθ=6200处，以半局部摩擦速度uτ归一化的雷诺应力剖面；并将结果与从 Zhang 等人（2018）（M14Tw018和 M8Tw048 案例）以及 Xu 等人（2021）（M8T1 案例）中提取的低焓构型结果进行了对比。尽管来流马赫数、摩擦雷诺数、绝对壁温以及壁面冷却速率的值存在显著差异，但湍流强度剖面具有可比性，且未表现出任何可归因于热化学非平衡（TCNE）效应的明显影响。与先前的观察结果（Duan 等人，2011；Lagha 等人，2011；Zhang 等人，2018）一致，较大的流向分量和较小的横流分量可能表明存在强烈的可压缩性效应，这将在 3.6 节中详细讨论。

在壁面约束可压缩湍流中，另一个经常研究的尺度变换是零压力梯度边界层中平均速度与平均温度剖面之间的关系。由 Walz（1969）推导的修正克罗科关系式为：

其中，Taw/Te=1+r((γ−1)/2)Me2，恢复因子r设定为 0.9；该关系式如图6（a）所示。在先前针对高马赫数、壁面冷却、高焓边界层的研究中（Duan & Martín，2011b），发现此关系式与从 DNS 数据中提取的精确T~/Te剖面存在偏差；实际上，在0.2<u~/ue<0.7的范围内显示出显著差异。这并不奇怪，因为关系式（3.5）是在量热完全气体假设下推导得出的。为了消除对热模型和化学模型的显式依赖，Duan 和 Martín（2011b）提出了一个类似的基于焓的方程，其表达式为：

其中，haw=he+21rue2。函数f(u~/ue)必须与来流条件、壁面温度以及表面催化作用（若存在）无关。对于量热完全气体，      

(3.6)式可简化为 (3.5) 式，只是经典公式中的项f(u~/ue)在此时等于u~/ue。通过对 DNS 数据进行曲线拟合，Duan 和 Martín（2011b）得出了关于f(u~/ue)的如下关系式：      

这种变换的结果如图 6（b）所示，与经典的基于温度的公式相比，其重合度有所提升。

3.3.热非平衡效应

通过将振动弛豫时间与流场的特征时间尺度进行对比，分析热非平衡的存在情况。为此，考虑以下振动达姆科勒数：

其中Da=O(1)表示非平衡效应起重要作用的状态。对于给定的流向位置，每个达姆科勒数的分子是固定的，且对三种分子而言均相同；因此，这三个双原子物种中每个参数的不同趋势完全由物种的分子弛豫时间tm决定，而tm取决于气体的局部状态。这些无量纲数的演变如图 7 所示。（3.8a–c）中的第一个比值将粘性底层（VS）内流体运动的特征时间与分子m的振动弛豫时间联系起来。图 7（a，d，g）展示了混合物中三种分子的这一量的演变。结果表明，对于N2和O2，各处的DavVS均远大于 1，而仅对于NO，其数量级为 1，但NO的浓度几乎可以忽略不计（见 3.5 节）。这两个时间尺度的解耦表明，在边界层内部（较冷）区域存在一个厚度为lv的层，其中分子处于振动冻结状态。在层流位置观察到DavVS远大于 1，这是由于粘性应力占主导地位。

值得关注的是DavLE，它表示大涡（LE）周转时间与弛豫时间的比值，通过位移厚度δ和uτ构建。对于分子氮，这两个时间尺度的数量级相当（图 7b），在湍流区域流动演变过程中，当y≈10时，DavLE的峰值范围为 0.4 至 0.7。另一方面，图 7（e）和 7（h）显示，对于O2和NO，DavLE远小于 1，这意味着这两种分子相对于大尺度湍流运动的特征时间更快地达到振动平衡。对于DavRT也有类似的结论，它将振动弛豫时间与局部流动停留时间（RT）进行比较。随着流体沿平板流动，O2和NO有足够的时间达到热平衡，而N2在所有位置都明显处于非平衡状态。这一结果是分子氮特有的高振动温度（以及较慢的弛豫速率）导致的。值得注意的是，由于湍流混合的贡献，在三个湍流位置，DavLE和DavRT均表现出高得多的值。

归一化的平均平动温度和振动温度的壁面法向剖面分别如图 8（a）和 8（b）所示。平动温度在壁面法向位置y≈10处达到约[2.3−2.5]Tw的峰值，对应于湍流强度的峰值。在整个边界层中，振动温度TV比T~小 20% 至 30%，这证明了在所研究的热力学条件下存在振动弛豫层，并且其峰值位于缓冲层和对数下层之间。在对应于Cf峰值的Reθ=2960处，可观察到最大的差异。此外，壁面梯度∂TV/∂y对峰值的大小和位置不太敏感，这解释了先前在图 2（b）中观察到的CVq近乎恒定的趋势。

归一化温度差(T~−TV)/Tw分别以y和y/δ为变量展示于图 9（a）和图 9（b）中。该量在转捩区域增大，在摩擦系数Cf达到峰值处达到约 2000K 的最大值，随后在湍流区域趋于减小。

朝着湍流区域移动时，壁面法向最大值从边界层的外部（y/δ≈0.6）移动到缓冲层，并稳定在y≈8处。热力学非平衡状态一直保持到计算域末端，这一点已通过对振动达姆科勒数的观察得到证实。

从粘性底层（我们记得壁面处T=TV=1800 K）到对数区的第一部分，均观察到T~−TV为正值。然而，随着流体沿平板流动，出现了一个以T~−TV<0为特征的区域，这意味着流体在振动上过度激发。由于边界层增厚，符号发生变化的高度从边界层边缘平稳移动到y/δ≈0.45处。这在图 10 中清晰可见，其中我们展示了纵向切面温度差的快照（顶部）以及由(T−TV)着色的速度梯度张量第二不变量的等值线。

Fiévet 等人（2019）在对振动非平衡与湍流混合之间的耦合进行分析时发现，大多数非平衡状态出现在激发不足的条件下。少数存在过度激发非平衡状态的区域是由转动 - 平动能量与机械能之间的联系所决定的。静态温度通过状态方程与其他变量耦合；因此，湍流强度的降低（朝着边界层边缘移动）会导致平动温度的降低。与此同时，振动能量几乎保持恒定，因为其弛豫速度较慢，这会在外部层中导致过度激发，直到两种温度都达到来流平衡状态。

图 11 展示了边界层内平动和振动温度的法夫雷脉动的归一化均方根（r.m.s.）。在图 11（a）中，在T~max位置（y≈10）的两侧各出现一个峰值，这在壁面冷却的湍流边界层（TBLs）中是常见现象（Duan、Beekman和 Martín，2010；Zhang 等人，2018；Di Renzo 和 Urzay，2021）。在湍流区域，外部均方根峰值占主导地位，这主要是由于T~向外边界方向降低，而内区和外区的绝对温度脉动具有可比性。相比之下，振动温度脉动没有表现出局部最小值（图 11b），而是单调增加到边界层边缘，在那里突然下降。这与层流区域通常观察到的剖面不同（见图中用于对比的Reθ=1600处的剖面），这可归因于热非平衡：振动弛豫时间比平动弛豫时间慢，并且与大尺度涡旋的周转时间相当

由于振动 - 平动能量交换引起的振动温度（TV）变化，其速率与由壁面法向湍流运动引起的变化速率相当，这是因为振动弛豫时间比平动弛豫时间慢，且与湍流时间尺度相当。因此，后者在边界层内重新分配振动温度方面起着主要作用，从而抹平了局部最小值。

图 12 给出了脉动温度的预乘能量谱，该谱以归一化展向波数λz为函数，进一步证实了这一趋势。对于两种温度而言，在λz≈1500和y≈2500处都能看到一个大尺度峰值，而仅在转动 - 平动温度谱中观察到一个内层（稍弱的）峰值，这与温度脉动的均方根是一致的。此外，发育良好的谱表明流动已达到完全湍流状态，且强制扰动策略的影响已完全消失。

边界层内热非平衡状态的持续存在与向湍流状态的转捩密切相关。在相同构型下未进行边界层触发的初步数值实验表明，对于层流状态，在平板末端，转动 - 平动温度与振动温度之间的差异几乎可以忽略不计（约为 50K）。正如 Urzay 和 Di Renzo（2021年）通过先验考虑所预期的那样，湍流情况下观察到的热非平衡与扫掠和喷射湍流运动密切相关，从而导致混合增强。在图13 中，我们报告平均值。

在Reθ=6200条件下，针对多个壁面法向位置，图中展示了基于扫掠事件（特征为u>0且v<0）和喷射事件（特征为u<0且v>0）的两种温度平均值。在粘性底层内，扫掠会从上方区域卷入高温流体，因此在近壁区域，其平均温度高于整体平均值，而喷射事件则呈现相反特征。由于转动 - 平动温度T~存在更显著的峰值，扫掠与喷射的条件平均温度差异在转动 - 平动模式中最为明显。

在温度峰值位置以外的区域（y≈10），观察到相反的趋势。此时，扫掠将外层的冷空气带入温度较高的内层区域，平动模式会迅速达到平衡，而振动模式则需要更长的弛豫时间，因此振动温度TV始终低于转动 - 平动温度T~。另一方面，喷射运动将高温且热非平衡的混合物带到边界层外层，由于该区域环境温度较低，振动弛豫时间迅速增加，因此在这一区域TV>T~

扫掠与喷射事件之间的最大温差出现在对数区的起始部分（直至y≈100），其中振动模式的温差高达 800K，转动 - 平动模式的温差则达到 1000K。图 14 所示的两种温度的条件概率密度函数（p.d.f.s）进一步证实了这一现象。三个壁面法向位置

考虑了三个壁面法向位置：记录到最大温差的位置（图 a 和 b，y≈10）、对数区（图 c 和 d，y≈100）以及边界层边缘（图e 和 f，y≈3000）。在靠近边界层边缘处，两种事件的概率密度函数均呈现出明显的右尾特征，且振动温度TV的值略大于平动温度T的值。向边界层内部移动，在y≈100处，两种概率密度函数均趋于对称。最显著的差异出现在缓冲层，平动温度T的分布呈左偏态，而振动温度TV的分布呈对称态，其分布形状与图 11 中所示的平动温度脉动平方（T2）存在局部最小值以及振动温度脉动平方（TV2）存在小平台相关。

3.4.热模型与雷诺类比

对于热非平衡流动，通过法夫雷平均振动能量方程（2.5）会出现一个新的未封闭项，即振动能量的湍流输运项ρujeV。该术语与总能量方程中出现的湍流输运项类似。后者通常通过引入湍流普朗特数Prt来建模，其定义为：

振动能量的湍流通量可通过引入振动湍流普朗特数来建模

此处，振动湍流普朗特数（如图 15b 所示）在这一相对较大的区域内取值接近 1，该区域中湍流输运占主导地位；与Prt类似，它在边界层边缘和近壁处会出现发散。

通过在图 16 中分别给出流向速度脉动u、法向速度脉动v与平动温度T、振动温度的相关性，我们可以获得更深入的认识。正如在低焓和高焓边界层中所观察到的（Zhang 等人，2018；Di Renzo 和 Urzay，2021），在对数区和外层区域，u与T呈负相关（图 a），v与T呈正相关（图 b）。这是由于：（i）扫掠事件将温度较低的流体卷入缓冲层；（ii）喷射事件将温度较高的流体卷向边界层边缘.在缓冲层下方，出于相同原因，这两种相关性均会改变符号，这与先前在图 13 中讨论的趋势一致。

然而，尽管近壁处的流向速度脉动与平动温度脉动完全相关，但其与振动温度脉动的相关系数（图 c）较弱在 0.7 到 0.85 之间）。这一现象再次归因于振动模式的弛豫时间更长与平动温度T相比，振动温度需要更多时间来适应周围温度。

结合湍流普朗特数的演变以及温度 - 速度脉动相关性可知，强雷诺类比（SRA）假设

正如文献中已观察到的那样，式（3.11）中的关系式与强雷诺类比（SRA）的估计结果（未展示）存在较大偏差，因为它没有考虑壁面的非绝热特性。已有多种修正形式被提出，其中最为成功的一种是 Huang、Coleman 和Bradshaw（1995）提出的修正形式。

通常称为 HSRA（修正强雷诺类比），其中\(T_t\)为滞止温度。通过摒弃量热完全气体的假设（Duan & Martín，2011b），可以得到广义修正强雷诺类比（GHSRA）：

式（3.13）左右两边比值的壁面法向演变如图 17（a）所示，在边界层的大部分厚度范围内（至少在速度脉动和温度脉动均不为零的区域），观察到了令人满意的一致性。通过采用与处理\(Pr_t\)相同的方法，广义修正强雷诺类比（GHSRA）也可扩展到振动温度，即：      

图 17（b）所示结果证明，广义修正强雷诺类比（GHSRA）在热非平衡条件下仍然有效，并且对振动能量同样适用。除了湍流输运，还应重点关注振动能量方程源项所引发的非线性湍流 - 热非平衡相互作用的建模问题，这些源项与平动 - 振动能量交换以及化学离解 / 复合过程相关。虽然热非平衡状态可通过的大小进行合理量化，但对源项的分析能够分离出产生这种非平衡状态的物理过程。

首先，此类非平衡状态的产生与这些物理过程密切相关。在该特定构型中，化学活性极弱，因此几乎可以忽略不计。

图 18（a）展示了沿壁面法向的分布，其分布趋势与图 9（a）和图 9（b）中的分布趋势高度一致，尽管在外层区域的局部放大图中可见的负值比内层区域的负值微弱得多（低两个数量级）。这种显著差异是由于以下综合效应所致：绝对温度值较低，且（在较小程度上）的差值较小，这分别导致式（2.15）中的分母增大和分子减小。

在当前情况下，尽管物理意义相似，但对后者进行分析主要是出于建模角度的考虑；也就是说，只有对源项进行恰当建模，才能准确描述非平衡状态。在这方面，借鉴 Duan 和 Martín（2011a）用于评估湍流 / 化学相互作用的方法，并为了量化湍流与热弛豫之间的耦合关系，我们构建了一个指标，其定义为：

偏离零值的程度反映了与湍流脉动相关的振动能量变化量。图 18（b）显示，利用平均量计算会低估其在整个边界层内的值，在内层中，偏差超过15%。这些数据进一步证实了本节所讨论的热非平衡与湍流之间的紧密耦合关系，并强调必须重视合适湍流闭合模型的开发，尤其是在雷诺平均 Navier-Stokes（RANS）建模中。

3.5.化学活性

本节研究高焓边界层中由摩擦加热引发的化学活性。图中展示了存在物种的法夫雷平均质量分数分布

混合物中各物种的法夫雷平均质量分数分布如图 19 所示。在本研究的流动条件下，分子氧是解离程度最高的物种，会生成原子氧。混合物中原子氧的最大质量分数约为 1%。边界层内的最高温度过低，即使在温度峰值处，也无法使氮发生显著解离，因此原子氮的量与其他物种相比可忽略不计。尽管如此，在当前条件下，通过正向置换反应 R4 和逆向置换反应 R5 会生成少量一氧化氮（约 0.7%）。

化学活性较低的部分原因是强烈的壁面冷却作用，这将最高温度限制在略高于氧解离阈值的水平（T>2000K）。热非平衡状态进一步限制了化学反应。高振动激发态分子的解离速度更快，因为它们在碰撞过程中所需交换的能量更少。相反，在激发不足的状态下，解离发生的可能性较低。这一机制在反应速率系数的计算中得到了数值上的体现（计算时采用两种温度的几何平均值，即\(T^{0.7}T_V^{0.3}\)）；结果表明，当存在强烈的非平衡效应时，反应速率系数会更低。另一方面，吸热化学反应往往会消耗转动 - 平动能量，不过这种效应会被振动模式的弛豫部分抵消 —— 振动弛豫会提供能量，从而在一定程度上抑制温度的降低。

对于所有物种，基于大涡周转时间的化学达姆科勒数最大仅为10^{-2}量级，这表明化学反应与湍流运动之间的相互作用较弱，呈现出近乎冻结的化学行为。关于绝热湍流化学非平衡边界层的详细讨论，可参见 Passiatore 等人（2021 年）的研究。尽管条件不同，但在本研究中可以得出类似的定性观察结果

3.6. 可压缩性效应

本分析通过讨论可压缩性效应的作用及其与热非平衡条件的相互作用来完成。图 20 展示了湍流马赫数和脉动马赫数的演变情况。尽管外部马赫数较高，但湍流马赫数在各处均小于 1，在湍动能产生峰值位置的最大值约为 0.8–0.9。不过，这样的值足以形成涡旋激波，粘性底层内声速线（图 20a 中三条竖线所示）的位置也证明了这一点。相应地，开始，脉动马赫数的值大于 1。外层的脉动马赫数峰值出现在（其值与 Zhang 等人（2018 年）在量热完全气体高速边界层中得到的值一致），这与边界层边缘沿湍流 / 非湍流界面存在的热力学性质大幅脉动有关。

计算域纵向切片上的归一化密度梯度幅值如图 21 所示，从中可清晰看到这些结构。热力学变量的这种大幅脉动并未伴随流场的强烈活动，这一点可由域末端部分的半局部标度脉动速度散度的值得到证明（见图 22）。较强的膨胀事件出现在近壁处，其幅值随离壁距离的增加而单调减小，图 23（a）展示了这一情况，同时还给出了 Xu 等人（2021 年）的结果以供对比。图 23（b）显示了在壁面法向位置处计算得到的的两个概率密度函数。两种分布都略微偏向负值，这表明强压缩的发生略多于强膨胀，尽管在内层区域这种事件更为强烈。然而联合概率密度函数（未展示）并未表现出任何偏好相关性，这意味着即使边界层内可能存在激波，其数量和 / 或强度也不足以对热力学和动力学量产生显著影响。

此前已在混合层（Fiévet 等人，2019 年）和均匀各向同性湍流（Zheng 等人，2020 年）构型中研究了可压缩性与热非平衡效应之间的相互作用，观察到

振动非平衡状态在某种程度上会因强烈的膨胀运动（无论是扩张还是压缩）而加剧。热力学量的突发性变化（如存在强密度梯度的情况），会因两种温度的弛豫时间不同而增大二者之间的差值。在此，高 |∇ρ| 区域主要存在于边界层边缘附近，该区域的平均温度值

相对较低；因此，对于边界层构型，可压缩性与热非平衡的相互作用预计不会那么显著。为此，我们在图 24 中展示了最后一个湍流位置处，与图 23（b）所示相同壁面法向位置上的概率密度函数。这些概率密度函数以强膨胀运动为条件。图 24（a）和图 24（b）中的概率密度函数在正值（负值）处达到峰值，这与图 18 中所描述的局部平均激发不足（过度激发）状态一致。正如 3.3 节中所解释的，由于激发不足状态在强度和数量上占主导地位，所有分布都呈强烈的右偏态。

在这两个位置，强压缩与略高于强膨胀的值相关；（b）图中分布之间的差异更大，这是由外层区域观察到的源项脉动小得多所致。总体而言，差异被证明是微小的，这验证了我们之前的陈述。      

4、结论

本文阐述了高超声速高焓湍流边界层（TBL）的数值模拟结果。所选的边界层边缘条件使得流动处于平动 - 振动非平衡（TCNE）状态。研究探讨了高温效应与湍流参量之间的相互作用，并对振动弛豫动力学进行了具体分析。

位于入口相似解正后方的强制扰动带引发了向湍流的突然转捩，其特征是转捩区壁面剪应力急剧增加，而在完全湍流区域则呈渐近衰减趋势。根据 Renard 和 Deck（2016）的表面摩擦分解结果，平均流场项和湍流项对摩擦系数（\(C_f\)）的贡献相当，而与边界层增长相关的项仅在向湍流破裂的区域不可忽略。

关于壁面热流，观察到其转动 - 平动贡献部分遵循与\(C_f\)相同的趋势，而振动贡献部分在整个湍流区域几乎保持恒定，数值比转动 - 平动贡献小一个数量级。范德赖斯特变换以及 Trettel-Larsson 变换在将流向速度剖面折叠到通用对数律方面效果较差。相反，Trettel-Larsson 变换和 Griffin 等人（2021）提出的基于总应力的新变换与线性律完全匹配，后者还改善了在对数区的折叠效果。

此外，Walz（1969）推导出的平均速度与平均温度剖面之间的关系式不适用于量热非完全气体；尽管存在热化学非平衡效应，Duan 和 Martín（2011b）提出的基于焓的类似方程仍显示出改进效果。

为了量化热非平衡状态，基于不同的特征流动时间尺度定义了三个振动达姆科勒数。结果表明，除了近壁区域（此处流动处于振动冻结状态）外，强烈的热非平衡条件一直持续到边界层边缘。分子氮是导致热非平衡的原因，因为其特征温度足够高，会导致弛豫速率较慢。正如边界层内扫掠和喷射事件与之的紧密耦合所观察到的那样，振动激发在很大程度上由湍流混合维持。扫掠会将较冷的空气从外部区域输送到平板壁面，而喷射则会将较热的空气带到远离内层的地方。平动模式能快速适应这种较大的温度脉动，与之相反，振动模式的重新调整则较为缓慢。在宏观层面上，直接的结果是在对数层的前半部分，振动温度滞后于平动温度；从该位置往后，流动向振动过激发状态转变，并一直保持到边界层边缘。通过一个衡量由湍流脉动引起的振动能量变化的参数进行量化后发现，湍流与热弛豫之间的相互作用在整个缓冲层内都十分显著，这证实了湍流输运发挥着主要作用。按照常规的雷诺平均 Navier-Stokes（RANS）建模策略，可以定义一个振动湍流普朗特数来封闭振动能量的湍流通量，其在边界层的大部分区域取值接近 1。与经典的\(Pr_t\)类似，在内层区域和边界层边缘会偏离 1。其余的强雷诺类比（SRA）关系式与估算结果呈现出可接受的一致性，且发现广义修正强雷诺类比（GHSRA）即使扩展到振动温度时也能较好地发挥作用。

二阶统计量显示，雷诺应力呈现出经典的自相似趋势，在自由流马赫数相当的情况下，其值与低焓条件下观察到的值相近。等温壁面条件导致温度脉动剖面呈双峰形状，其中外部峰值更为显著。相反，由于缓冲层中存在明显的非平衡条件，振动温度脉动则呈现出单调递增的剖面。

由于设定的壁面温度不足以引发显著的化学活性，因此壁面处的气体离解几乎可以忽略不计。观察到化学反应与湍流运动之间的相互作用较弱，所有物种的达姆科勒数最大仅为量级。通过半局部标度速度散度脉动对可压缩性效应进行的统计分析表明，强烈的膨胀运动主要出现在近壁处，且强压缩发生的概率略高于强膨胀。为了将可压缩性效应与热非平衡状态联系起来，计算了以局部散度值为条件的平动 - 振动能量交换源项的概率密度函数（p.d.f.）。正如先前关于热非平衡混合层和各向同性湍流构型的研究所观察到的那样，强压缩与略高的值相关。然而，对于边界层构型，这种相关性要弱得多，因为最强的密度梯度和温度梯度出现在边界层的外部区域，而该区域由于温度较低，流动呈现出振动冻结行为。