双锥/双楔流动中的高温气体效应仿真模拟

袁军娅，任 翔，蔡国飙，王伟宗，贺碧蛟

(北京航空航天大学宇航学院，北京 100191)

引 言

在大气层内或跨大气层高超声速飞行的飞行器，与周围空气剧烈摩擦并对前方空气强烈压缩，由于黏性耗散效应和激波强烈压缩，气体损失的巨大动能转变为内能，造成流场温度急剧升高，形成高温高压的气动热环境，引起空气分子一系列复杂物理化学反应现象如离解、电离和内部能态激发等，从而造成高温真实气体效应[1]。

在高温真实气体描述方面，计算流体力学方法和直接模拟Monte Carlo(direct simulation Monte Carlo，DSMC)方法都得到适当发展并在科学研究和工程实践中得到了广泛应用[2-3]。CFD和DSMC方法对高温真实气体效应的模型化包括热力学参数、输运参数、热力学非平衡过程、化学非平衡过程的描述。目前广泛使用的各类模型主要是从常温向高温扩展、平衡态向非平衡态扩展，在常温和平衡态下相互符合的不同模型在高温非平衡态下会产生明显的差异，这些差异会影响流场和气动参数的模拟预测。

目前热化学非平衡流CFD模拟，自20世纪60年代提出Lighthill离解模型、Landau-Teller振动松弛模型[4]，到80年代的Park双温反应模型[5]，这些物理模型由于某些局限性或者缺陷而不断推旧出新。Park双温模型主要关注振动态的松弛过程，并假设转动态与平动态时刻达到平衡，从而使用两个温度来描述流场的热运动状态。同时CFD中大气组分的化学反应速率包含Dunn等[6]，Gupta等[7]和Park[8]机理数据。

DSMC方法主要使用L-B(Larsen-Borgnakke)[9]和量子L-B模型分别模拟分子转动和振动态的激发与松弛，并先后发展了总碰撞能量(total collision energy，TCE)模型[10]、振动控制离解(vibrationally favored dissociation，VFD)模型[11]、广义碰撞能量(generalized collision energy，GCE)模型[12]和量子动力学(quantum-kinetic，Q-K)反应模型[13-15]来描述化学反应过程。其中Q-K模型不同于其他依赖于宏观物理量的化学反应模型，是仅基于两个碰撞粒子的基本属性的分子级化学模型。Q-K模型对于分子振动能采用简单谐振子模型，并综合考虑分子的离解与振动能的激发。Liechty[16]对Q-K模型进一步发展并囊括了电子能级和电离反应。

双锥和双楔模型是验证激波-激波和激波-边界层干扰、流动分离与涡旋等现象的主要模型，目前学者对该模型进行了广泛的理论分析、风洞测量和数值模拟研究[17-25]。1999年，Olejniczak等[18]采用带化学反应的N2流进行双楔实验，并采用量热气体模型进行数值仿真，模拟的分离区比实验值小。2003年，Nompelis等[26]对高超声速双锥实验的数值模拟将模型的传热速率高估了约20%。2007年，Hash等[27]研究了化学反应速率对CUBRC LENS设备中高焓空气实验的影响。2012年，Holden等[22]介绍了在LENS超高速地面试验装置中为获得基本测量数据而进行的实验方案，并以此来评价和发展用于预测非平衡流动化学、边界层过渡和激波/湍流边界层相互作用的模型。2017年，Knight等[28]采用Mach数为7.1的双楔形模型，在2.1 MJ/kg和8 MJ/kg 的空气和氮气中，对高超声速激波层流边界层相互作用的CFD预测能力进行了评估。

本文通过可描述内能松弛和化学反应的CFD和DSMC对3个具有试验数据的双锥双楔流动进行仿真分析，研究高温真实气体效应对流动预测的影响。

1 物理模型和数值方法

1.1 物理模型

使用CFD和DSMC方法分别对焓值为2.23，8.17，26.1 MJ/kg这3组双锥或双楔流动进行仿真[18，21，28]，来流工质均为N2，3组几何模型见图1，来流参数见表1。仿真过程中考虑稀薄效应和高温真实气体效应的影响。

表1 来流参数和壁面条件Table 1 Incoming flow parameters and wall conditions

(a) Double cone

1.2 数值方法

1.2.1 CFD方法

采用Park双温模型来模拟高超声速流动中的热化学非平衡效应，控制方程除包含连续性方程、动量方程、总能量方程外，还添加了振动能守恒方程

式中，U为守恒量，通量分为无黏通量Fi，inv和有黏通量Fi，vis,W为源项，下标i指代三维空间，对于Park双温模型

U=(ρ，ρs，ρu1，ρu2，ρu3，Eve，E)
W=(0，ωs，0，0，0，ωv，0)

式中，ρ为密度;u为速度;E和Eve分别为总能量和振动-电子能;ωs为组分s的质量变化，取决于化学反应过程;ωv为振动能的变化量，由振动-平动松弛过程和化学反应共同决定

ωv=QV-T+QC-V

式中，振动-平动项QV-T采用Landau-Teller公式计算[4]

振动-平动松弛时间τV-T取为Millikan-White拟合与Park修正项的和值[29-30]。

化学反应带来的振动-电子能量变化为

其中，A,β,Ta分别为指前因子、温度指数、活化温度，Tc，f为正向反应的控制温度，由Ttr和Tve共同决定。

这里只考虑N2的离解反应(N2+ N2= 2N + N2，N2+ N = 2N + N)且其反应速率取自DSMC的Q-K模型[31]，忽略电离反应。使用基于第一激发能级的简单谐振子模型和多电子激发态模型共同来确定气体的热力学参数。

1.2.2 DSMC方法

DSMC方法是由Bird提出的基于概率的直接分子模拟，并将分子运动与碰撞过程解耦处理[32]。本文使用量子L-B模型模拟分子振动态的激发与松弛过程。DSMC分别模拟转动能和振动能的松弛，而CFD的双温模型中假设分子的平动和转动能间是平衡的。

采用Q-K模型模拟N2的离解反应[14]，即对于一对碰撞分子，由碰撞能量确定可达到的最大振动能级

如果此能级高于分解能级，则发生离解。

1.2.3 模型设置

本文分别在CFD和DSMC方法中考虑高温条件下振动能和化学反应过程的模拟，设置了如表2中所示的几种模式，并且为了在结果中使用方便，约定了简写标记。受限于存储容量和计算能力，仅对低焓流动使用DSMC方法。

表2 数值方法设置Table 2 Numerical method settings

1.3 网格无关性验证

CFD和DSMC均采用贴体四边形网格，依据文献[33]中的Recell，w=1准则分别设定双锥/双楔靠近壁面附近的法向网格尺寸。CFD方法具有成熟的离散方法和数值格式，而DSMC方法对网格的依赖严格，这里专门对DSMC方法进行网格无关性验证。

分别采用网格数量为7.8×104，3.12×105，1.248×106且贴壁网格尺寸为0.5，0.2，0.2 mm的3套网格对双锥工况进行无关性验证。图2展示了获得的壁面热流分布。结果表明，几套网格主要在分离区和激波-边界层干扰位置处存在一定差异，并且后两套网格的结果相近，因此在后续仿真中采用第2套网格。另外两个双楔流动也进行了网格无关性验证，这里不再赘述。

图2 网格无关性Fig.2 Grid independence

2 结果与讨论

2.1 低焓双锥流动

美国 Holden 等的CUBRC第7次试车实验数据作为稀薄流仿真软件的验证算例。实验测量稀薄超声速来流对双圆锥体的气动力和气动热值，是国际上一个典型的稀薄气体动力学实验。

图3展示了使用考虑振动能松弛及氮气离解反应的CFD方法对双锥流场的模拟结果。从第一锥的头尖部产生一道附体激波，因第二锥而形成非附体激波，且形成三波点;透射激波入射到第二锥的锥面上，使流动分离，且产生一个很大的当地压升，这个压升使上游流动分离并制造出一个超声速的欠膨胀射流，射流在靠近第二锥锥面处流向下游，分离区尺寸取决于激波入射点位置及激波强度。在双锥体的拐角上形成的分离带和与第一个锥角上的斜激波相互作用产生了分离激波。

(a) Ma

注意，此风洞实验将经过Laval喷管加速的 N2作为来流气体，由于膨胀后的N2特别稀薄所以振动态松弛明显，采用喉部的振动温度作为此来流的振动温度，此时气体Tvib>>Ttr，并且在双锥附近激波内的流动，依然表现出强烈的振动态弛豫效应。另一方面，虽然考虑了N2离解反应，但是激波层内生成的N原子密度极低，化学反应过程可以忽略。

图4展示了在CFD和DSMC方法中采用多种高温气体模型获得的双锥壁面压强和热流分布。CFD和DSMC的结果都表明，考虑振动态激发和松弛过程条件所获得的壁面热流升高而壁面压强无变化，因为压强主要由来流气体动量决定;而化学反应十分微弱，所以化学反应过程是否考虑对热流和压强均无影响。需要强调，如果考虑振动态激发，但是认为它总是与平动-振动保持平衡，会造成分离涡区域减小而二次反射波不明显。忽略振动态激发会造成第一锥面上的热流预估值偏小，这是由于双锥流动实验中，来流的振动能远高于平动能。CFD和DSMC在不同设置下获得的壁面压强和热流的分布，具有很好的一致性，且都与实验值符合。这说明，CFD和DSMC方法分别从宏观和微观发展了对热化学非平衡现象的描述方法，均能很好地描述此复杂流动，另一方面，由于此流动比较稀薄，所以DSMC方法是适用的。后续的中焓双楔-1和高焓双楔-2流动十分稠密，而DSMC方法要求网格尺寸要严格小于局地的分子平均自由程，受限于计算能力，后续模拟只使用了CFD方法。

(a) Pressure

2.2 中焓双楔-1流动

图5展示了CFD方法中是否考虑振动态激发与松弛以及N2的离解过程(CFD-noVib，CFD-Vibeq，CFD-Vibneq和CFD-Vibneq+React)对中焓双楔-1温度场预测的影响。考虑激发振动态可以明显降低流场内温度的最大值，减小激波层厚度、分离区尺寸、三叉点位置等流场特征，并且CFD-Vibneq获得的平动温度介于CFD-noVib和CFD-Vibeq的温度场之间，而所获得的振动温度明显低于平动温度，在两个楔面附近的振动温度都不高于3 000 K，而平动温度主要介于3 000～6 000 K之间。

图6展示了考虑高温条件下振动态激发松弛过程以及N2的离解过程对双楔-1壁面气动参数分布的影响。可以看出CFD-noVib获得的分离涡尺寸最大，CFD-Vibneq的结果是介于CFD-Vibeq和CFD-noVib之间，而CFD-Vibneq+React与CFD-Vibneq相近。仿真获得的热流仅在第一楔面开始部分与试验结果一致，之后的分离区和第二楔面上都存在明显差别，本文的仿真结果与文献[28]中Komives等的仿真结果基本相同，这种偏差很可能是由三维效应造成的。

(a) Temperature of CFD-noVib

2.3 高焓双楔-2流动

图7展示了CFD方法中是否考虑振动态激发与松弛以及N2的离解过程(CFD-noVib，CFD-Vibeq，CFD-Vibneq和CFD-Vibneq+React)对高焓双楔-2温度场预测的影响，和前面算例一样未使用湍流模型。当忽略了N2振动态激发时，流场温度最大值将近20 000 K，考虑振动态激发的流场温度最大值约为17 200 K，进一步考虑N2离解，则流场温度最大值降低至14 400 K。第一楔面附近存在强烈的振动松弛效应而化学效应不明显，而第二楔面附近振动态与平动态间基本保持平衡且N2的离解强烈。并且考虑振动激发和离解反应所获得的斜激波层更薄、三叉点前移且分离区减小。

图6 双楔-1壁面热流分布Fig.6 Wall heat flow distribution of the double wedge-1

(a) Temperature of CFD-noVib

图8展示了考虑高温条件下振动态激发松弛过程以及N2的离解过程对双楔-2壁面气动参数分布的影响。可以看出CFD-noVib获得的分离涡尺寸最大，而CFD-Vibeq，CFD-Vibneq和CFD-Vibneq+React的分离涡尺寸较小且相近。

图8 双楔-2壁面热流分布Fig.8 Wall heat flow distribution of the double wedge-2

从图6和8可以看到，中焓双楔-1和高焓双楔-2流动的CFD计算结果与实验均存在较明显的偏差，分析可能存在以下原因:(1)双楔流动具有显著的三维效应，而本文直接进行的二维模拟，这不可避免地引入误差;(2)本文的CFD模拟均采用层流假设，高焓双楔-2流动的来流Re=2.9 ×105，会在第二楔面上发生转捩过程，后续需要采用能够预测转捩点的层流-湍流模型进一步仿真研究;(3)考虑高焓来流的稳定性以及试验测量精度，实验结果也具有很大的不确定性。

3 结论

采用可描述热化学非平衡过程的CFD和DSMC方法对三类双锥/双楔流动进行了数值模拟，并讨论了高温真实气体效应对流场结构和气动参数的影响，得到了如下结论:

(1)随着来流气体焓能的增加，高温真实气体效应越显著，数值方法对热化学非平衡过程的准确模拟会影响到流场和气动分布的预测。

(2)采用量热完全气体模型往往会高估流场温度、激波厚度，并造成拐角附近的分离涡过早分离而激波三叉点后移。

(3)不同工况的内能松弛和化学反应程度不同，比如低焓双锥流动中振动能表现出强烈的松弛现象而化学反应可忽略，高焓双楔-2流动中振动态激发且存在明显的化学反应，但是振动态松弛现象不明显。在高超声速模拟中应注意正确应用物理现象的模化模型。