湍流边界层与高温气体效应耦合的直接数值模拟

吴正园, 莫 凡, 高振勋,2,*, 蒋崇文,2, 李椿萱,2

(1.北京航空航天大学, 北京 100191 ;2.中国航天科技集团 航天飞行器气动热防护实验室, 北京 100191)

0 引 言

现有的高超声速飞行器，当高温气体效应显著时其表面绕流多为层流，故已有对高温气体效应的研究多针对层流流态。然而，随着飞行器设计理念的发展和设计能力的提升，未来的高超声速飞行器将会以较高马赫数在相对低的高度如25～40 km高度上飞行。此时，由于来流雷诺数仍较高，飞行器表面将转捩形成湍流边界层。同时，由于飞行马赫数较高，在激波和黏性共同作用下，大量动能转化为内能，使得激波后边界层内的总温较高，边界层内将出现高温气体效应。即随着飞行马赫数的提高，空气将逐渐出现振动能激发、气体分子离解/复合等现象[1]。高温气体效应引起的热力学/输运性质变化等物理化学过程，将对湍流边界层产生影响，进而改变表面压强和气动力/热分布。另一方面，湍流引起的温度、压强脉动也会影响高温空气的化学反应进程和空气热力学/输运性质的变化。可见，上述情况下边界层内存在湍流与高温气体效应的耦合作用，该耦合作用将对飞行器的气动力/热、结构负载等产生影响。

受限于研究手段和应用需求，现阶段对湍流边界层与高温气体效应耦合现象的研究非常少。数值模拟是开展高温气体效应湍流边界层研究的主要手段之一。对湍流的数值模拟，工程上常采用雷诺平均N-S方程(RANS)方法，然而RANS方法对存在高温气体效应的流场的适用性未知。相比之下，DNS方法由于不需要引入任何经验性模型，并且可以提供流场内详尽的统计平均和脉动信息，更适合开展湍流边界层与高温气体效应耦合现象的机理研究。

然而，目前已有的针对空气介质的高超声速湍流边界层DNS研究绝大多数都是在低焓、无高温气体效应的条件下进行的。童福林等[2]采用DNS方法对高超声速激波-湍流边界层干扰流动进行了研究。Duan等[3-4]开展了马赫数5～12条件下的DNS研究，讨论了壁温和马赫数的影响。谭杰等[5]开展了高超声速平板/空气舵热环境数值模拟研究。韩宇峰等[6]开展了高超声速横流转捩的DNS研究。Zhang等[7]开展了高超声速湍流边界层数据库构建的DNS研究。当前，仅有少部分相关学者初步了开展了高温气体效应湍流边界层的高精度数值模拟研究[8]。

Duan和Martin[9]以楔形激波后的湍流边界层为背景流动，设计了低焓和高焓算例，采用一致的高温气体效应模型开展DNS研究。低焓和高焓流动来流马赫数分别为3.4和9.4，分别反映较弱的高温气体效应和较强的高温气体效应。Duan等人对高温气体效应下的湍动能输运、梯度输运假设、可压缩效应和雷诺比拟等进行分析，发现低焓边界层中成立的关系式如范-德列斯特变换和雷诺比拟等理论在高焓条件下仍然成立。同时，Duan等人还就高温气体效应对湍流统计量和气动力/热的影响进行了研究。然而，由于其低焓算例和高焓算例的来流马赫数、雷诺数等条件并不相同，对比湍流统计量结果时无法排除来流条件的差别产生的影响，即该条件下低焓和高焓算例湍流统计量结果的差异可能由来流条件的差别引起，而不是源自高温气体效应的影响。因此，Duan等人的研究方式无法精确地分析高温气体效应对湍流边界层湍流统计量和气动力/热(表面压强、摩阻系数和热流)等的影响。

Duan等人还在上述研究的基础上对湍流与化学反应干扰(TCI)现象开展了研究[10]，发现一方面湍流使各组分源项明显增大，但另一方面TCI对湍动能、平均和脉动温度及密度和气动力/热等的影响较小。此外，Chen和Li[11]对来流马赫数为6和10的槽道湍流开展了考虑高温气体效应的时间发展DNS研究，发现高马赫数条件下Morkovin假设仍然适用，这一发现与Duan的结论一致，但未就高温气体效应对湍流统计量的影响开展研究。综上可见，目前尚未有研究就高温气体效应对湍流边界层统计性质和气动力/热的影响规律开展深入的分析，关于高温气体效应对湍流边界层的统计性质和对气动力/热的影响规律尚不清楚。

本文将开展高温气体效应下高超声速湍流边界层的DNS研究。首先开展高超声速湍流边界层DNS数值模拟技术的研究。在此基础上，设计高温气体效应湍流边界层算例，分别采用量热完全气体模型和化学反应模型开展相同来流条件下的DNS研究，分析高温气体效应对湍流统计量和气动力热等的影响规律。

1 数学物理模型和数值方法

首先介绍本文数值模拟所采用的数学物理模型和数值方法。

1.1 化学非平衡流动的瞬时N-S方程

在笛卡儿坐标系下，含化学反应流动的瞬时Navier-Stokes(N-S)方程可写成如下守恒形式：

(1)

(2)

关于组分质量分数的边界条件，选择了完全非催化壁面边界条件，即不考虑壁面对化学反应的催化作用，其具体数学形式为:

(3)

其中η表示壁面法向方向，Y为组分质量分数。

本文开展DNS数值模拟时，无黏通量的离散采用7阶WENO格式，黏性通量离散基于4阶中心差分格式，时间推进采用双时间步迭代方法。

1.2 入口湍流边界条件

采用DNS方法模拟充分发展湍流时，为了节约计算量，常引入上游入口边界条件，即给定当前数值模拟入口处的湍流脉动条件。目前常用的湍流入口脉动生成方法包括白噪声随机扰动方法、转捩诱导方法、湍流库方法和“回收/调节”方法[14-18]等。

本文DNS入口湍流脉动生成采用湍流库方法的思路，即将已有DNS研究的湍流脉动数据库通过合理的标度之后，应用到当前数值模拟的入口边界条件中诱发边界层失稳，并使之经过一段流向距离的调整后形成充分发展的湍流边界层。具体的标度方法参考了Martin[19]的研究。Martin[19]通过对Spalart 不可压DNS[20]的数据采用无量纲化和强雷诺比拟关系式进行转换后，得到了可压流DNS所需的湍流入口脉动数据。本文的湍流数据库取自课题组已开展的超声速(M∞=2.92)平板湍流边界层DNS研究。对超声速的湍流数据库根据Duan的方法进行转换，获得本文DNS的入口湍流脉动数据。具体转换公式如下：

(4)

(5)

本文假设入口边界为二维平面(y-z平面)，其中y为壁面法向方向，z为展向。速度脉动V′(u′、v′和w′)根据式(4)进行转换，其中uτ=(τw/ρ)1/2为壁面摩擦速度。z+=z/yτ，在边界层内层，y*=y/yτ(yτ为边界层内层特征长度尺度)；在边界层外层，y*=y/δ(δ为边界层厚度)。温度和密度脉动根据式(5)获得，平均温度、密度和马赫数和uτ等结果通过RANS得到。同样，入口处的平均流场亦用RANS计算得到。则入口边界流场基本量瞬时值可写为：

(6)

式(5)中的系数b用于保证瞬时温度T为正值,式(6)中q表示u、v、w、ρ、T。

图1给出了该方法的示意。

图1 入口湍流边界生成方法示意

2 数值方法验证

本文的数值模拟采用课题组开发的ACANS有限差分CFD软件，该软件的准确性已经得到大量数值模拟算例的验证[21-22]。本节基于ACANS软件开展了马赫数7.3低焓平板湍流边界层的DNS研究，通过与实验结果的定量对比，验证本文所采用的数值方法用于开展高超声速湍流边界层DNS研究的可靠性。

2.1 计算设置

以来流条件作为初始计算流场，入口处的边界条件由平均值〈q〉和湍流脉动值q′两部分组成。〈q〉取相同条件下RANS收敛流场中边界层厚度δ=7.1×10-3m处的结果。q′基于湍流库方法得到。非定常计算时全局时间步长为5×10-8s，推进10Tref物理时间后开始进行结果的统计(Tref流动的特征时间Lx/u∞，其中Lx为特征长度，取0.1 m，u∞为自由来流速度)。选择统计15Tref物理时间后的结果进行分析，此时统计的流场平均值和脉动均方根已基本稳定，不再随统计时间的增加而变化。

2.2 结果分析

图2展示了DNS得到的截面瞬时温度云图。在近入口区域，未观察到明显的流场失稳现象。而在下游区域，温度分布则呈现出了较明显的湍流脉动特征，说明本文的DNS能够模拟出流场在入口边界湍流脉动诱发下失稳形成湍流边界层的过程，并且可获取到较精细的湍流边界层特征结构。

图2 瞬时流场温度云图

图3和图4分别给出了DNS得到的平均密度和平均速度结果与实验测量结果的对比。实验[23]测量了x=0.065 m处沿边界层法向的平均密度、平均流向速度u和u的脉动均方根urms的分布。可见DNS的平均结果与实验符合较好。

图3 平均密度型结果对比

图4 平均速度型结果对比

脉动速度结果的对比如图5所示。实验测量的urms用u*进行无量纲化。u*定义为u*=(τw/ρ)1/2。可以看到，DNS的urms峰值结果略低于实验的，整体而言，DNS的脉动速度结果与实验结果较为一致。

图5 脉动速度结果对比

以上结果验证了本文采用的数值方法包括数值格式和湍流入口生成方法等的准确性。在此基础上，将开展高温气体效应下的湍流边界层的DNS研究。

3 高温气体效应湍流边界层DNS

本节构造了包含高温空气化学反应的高超声速平板湍流边界层算例，分别选择量热完全气体模型和空气化学反应模型开展DNS研究。以量热完全气体模型DNS结果为对照，分析高温气体效应与湍流的耦合作用。

3.1 算例设置和基本结果

设计的高温气体效应下的湍流边界层自由来流条件为M∞=7、P∞=50 000 Pa、T∞=3000 K。考虑到本文的研究目的是研究高温气体效应与湍流边界层相互作用机理，因此选择了简单的平板外形。自由来流条件的设计则假设一个20°的斜楔以马赫数20在30 km高度飞行，选择与此条件下斜楔波后的流场相近的流动参数作为DNS的自由来流参数。在本文选取的来流条件下，流场中的空气组分以N2、O2、NO、N、O为主，因此各组分间的化学反应过程采用Gupta的5组分6反应模型进行模拟。

采用结构网格，流向和展向网格均匀分布，间距分别为Δx=6.7×10-5m、Δz=6.7×10-5m。壁面法向第一层网格的间距为Δy1=2.5×10-5m。对应的无量纲距离分别为Δx+=6.2、Δz+=6.0、Δy1+=0.45。网格总量为4000万。采用等温壁面条件，壁面温度Tw=2400 K。全局时间步长为2.5×10-8s。非定向计算推进15Tref物理时间后开始进行结果的统计，选择统计30Tref物理时间后的结果进行分析。

展向截面流场瞬时温度云图如图6所示。可以看出此时边界层内的温度最高在6000 K左右。O组分的瞬时质量分数分布如图7所示，可以看出当前的来流条件下，已有大量的O2组分分解成O原子，O原子浓度最高已到0.17左右，说明本文所设计的来流条件已足以产生较明显的高温气体效应。

图6 化学反应模型DNS 瞬时流场温度云图

图7 DNS研究O原子组分的瞬时质量分数云图

图8展示了下游x=0.06 m处沿边界层法向的各组分质量分数分布(如无特别说明，下文展示的沿y方向各物理分布结果均取自该流向位置处)，以说明边界层内空气组分间化学反应的进行程度。该站位处边界层内层区域大部分O2分子已分解，N2分子只分解少量。

图8 x=0.06 m处沿y方向平均组分浓度分布

3.2 统计平均量和脉动量结果分析

图9展示了沿y方向平均温度结果的对比，图中采用单温度模型(本文的温度均指平动-转动温度)，图中Calorically表示量热完全气体模型，Chemical表示化学反应模型。可以看出，在边界层内，对化学反应模型，空气组分的分子振动能激发使部分平动-转动能量转换为分子振动能。同时，以分解反应为主的化学反应的吸热效应亦吸取平动-转动能量。故高温气体效应使流场平均温度显著降低。在边界层外，由于温度相对较低，高温气体效应不明显，高温气体效应对平均温度影响较小。

图9 平均温度结果对比

平均密度结果的对比如图10所示，在边界层内，采用化学反应模型计算的平均密度值较采用量热完全气体模型计算的结果明显增大，即高温气体效应使平均密度增加。

图10 平均密度结果对比

图11和图12分别展示了量热完全气体和化学反应气体模型的温度、密度脉动均方根值分布结果的对比。本文采用脉动量的相对值(脉动均方根值与当地平均值的比值)，来对比分析高温气体效应对温度和密度脉动的影响。显然，在整个边界层内，采用化学反应模型计算得到的温度脉动均方根明显小于量热完全气体模型的对应结果，表明高温气体效应对温度脉动有较强的抑制作用。对密度脉动而言，在边界层的外层区域，高温气体效应对密度脉动的影响作用同对温度的相同，即使密度脉动减弱。然而，在边界层靠近壁面区域，化学反应模型计算的脉动密度值却要高于量热完全气体模型的，说明高温气体效应对密度脉动的影响趋势在边界层内并不一致。

图11 脉动均方根温度结果对比

图13展示了流向、壁面法向和展向速度脉动结果。在近壁面区域，化学反应模型计算的urms值较高，且峰值较大。然而过峰值后，urms值随y增大的衰减趋势更明显，故在边界层外层区域，化学反应模型的urms值小于量热完全气体模型的结果。高温气体效应对壁面法向速度脉动的影响趋势与对流动速度脉动的一致，并且化学反应模型的vrms值与量热完全气体模型的差别更加显著。对展向速度脉动而言，在整个边界层范围内，化学反应模型的wrms值均小于量热完全气体模型的。可见，高温气体效应使边界层内的展向速度脉动减弱。

3.3 湍流马赫数及雷诺应力结果分析

湍流马赫数Mt和声速a定义为:

(7)

计算的湍流马赫数结果的对比如图14所示。可以看到，在边界层内Mt最大值接近0.8。一般认为Mt>0.3时可压缩效应已不可忽略。可以认为当前算例的可压缩效应已经比较明显。

图14 湍流马赫数结果对比

雷诺应力<-ρu′v′>结果的对比如图15所示。可以看出，在边界层内层y+<500的区域内，采用化学反应模型计算得到的<-ρu′v′>显著大于采用量热完全气体模型获得的结果，并且峰值结果更大。然而，过峰值点后，化学反应模型的<-ρu′v′>值随y+增加降低得更为剧烈，导致在边界层外层区域内，高温气体效应反而使雷诺应力<-ρu′v′>减小。可发现，高温气体效应对雷诺应力的影响趋势与对密度脉动和速度脉动的影响规律一致。

图15 平均雷诺应力<-ρu′v′>结果对比

3.4 气动力热结果分析

本节将在3.3节对雷诺应力分析的基础上，讨论高温气体效应对表面压强、摩阻和热流的影响规律。壁面压强结果的对比如图16所示。整体上，高温气体效应使壁面压强明显增加，增幅约为10%。

图16 平均壁面压强结果对比

壁面压强脉动均方根结果的对比如图17所示。化学反应模型计算的壁面压强脉动略高于量热完全气体模型的。高温气体效应对壁面压强脉动有一定的增强作用。

图17 壁面压强脉动结果对比

图18给出了壁面处速度梯度、黏性系数结果的对比。化学反应模型计算的无量纲速度梯度∂u/∂y高于量热完全气体模型，其差别约为10%左右。本算例采用等温壁边界条件，对量热完全气体模型，μw采用Sutherland公式计算，只与壁面温度有关，因此其壁面黏性系数沿流向不变。化学反应模型计算的壁面处黏性系数μw值略高，与量热完全气体模型结果的差别在1%作用。

(a)平均速度梯度∂u/∂y

图19给出了平均摩阻系数结果的对比，化学反应模型的摩阻系数结果明显高于量热完全气体。由上述分析可知，高温气体效应使壁面处速度梯度增大，对壁面处黏性系数影响很小，故最终使摩阻增大。

图19 平均摩阻结果对比

图20展示了壁面处平均温度梯度和热传导系数结果的对比。高温气体效应使平均温度降低。尽管高温气体效应增强了湍流引起的动量和能量交换，但吸热效应对温度的影响更为明显，因此高温气体效应最终使壁面处平均温度梯度降低，如图20(a)所示。平均壁面热流结果见图21，化学反应模型计算的热流值显著低于量热完全气体模型的结果。此外，采用化学反应模型计算的壁面热传导系数略低于量热完全气体模型。高温气体效应对壁面热流的影响主要体现在对壁面温度梯度的影响上。

(a)平均温度梯度∂T/∂y

图21 平均物面热流结果对比

4 结 论

本文开展了高超声速高温空气化学反应湍流边界层的直接数值模拟(DNS)研究，基于模拟结果分析了边界层内高温气体效应与湍流的耦合作用机理。主要结论如下：

1)本文所采用的数值方法以及基于湍流数据库的入口脉动量生成方法能够实现高超声速湍流边界层的DNS高精度模拟。DNS结果与马赫数7.3的平板湍流边界层实验数据吻合良好，验证了DNS数值技术的可靠性。

2)基于所设计的高温气体效应湍流边界层DNS数值研究发现，高温气体效应使边界层内平均温度显著降低，平均密度显著升高。对于湍流脉动量，高温气体效应对温度脉动有较强的抑制作用，而对流向及壁面法向速度脉动和密度脉动的影响趋势在近壁面区域和边界层外层是相反的。在近壁面区域，流向及壁面法向速度脉动和密度脉动都被增强，而边界层外层，高温气体效应则会抑制流向和壁面法向速度脉动和密度脉动。但展向速度脉动在整个边界层内均被降低。

3)对于边界层内的雷诺应力，在y+>500的区域，高温气体效应反而使雷诺应力值减小; 而在y+<500的近壁区域内，高温气体效应明显增大了雷诺应力值，进而增强了湍流引起的动量交换，从而使壁面处速度梯度增加，最终使摩擦阻力增加。

4)高温气体效应使壁面平均压强增大，并且增强了湍流引起的壁面压强脉动。高温气体效应使壁面处温度梯度大幅降低，同时使壁面热传导系数减小，因此使平均壁面热流降低。