钝体头部热化学非平衡效应的CFD与DSMC仿真分析

任 翔，袁军娅，蔡国飙

（北京航空航天大学，北京，100191）

0 引 言

高速流动广泛存在于战略战术打击、导弹系统防御、航天器发射与再入、地外天体探测等方面。高速流动一般指流速Ma＞5的流动。流动中往往存在薄激波层、壁面附近强的熵梯度和粘性耗散，并产生强烈的高温现象。高温条件下气体的振动能、电子能被激发，且当分子的振动能级和电子能级达到一定程度时发生分解、电离等一系列化学反应，这称为高温真实气体效应。同时高速飞行往往出现在上层大气或者低密度环境中，流动往往伴随着稀薄效应。此时稀薄效应和高温真实气体效应共同造成强烈的热化学非平衡，这为高速流动和飞行器的气动特性的准确预测提出挑战。

考虑到高速流动的跨流域特性，CFD和DSMC方法广泛用于高速流动的机理研究和工程应用，尤其针对热化学非平衡过程，CFD和DSMC方法分别从宏观和微观角度发展了多种模型。CFD方法自20世纪60年代提出Lighthill离解模型、Landau-Teller振动松弛模型，到20世纪80年代的Park双温模型。而DSMC方法，主要使用L-B（Larsen-Borgnakke）或者量子L-B模型来模拟分子转动、振动态的激发与松弛，先后发展了总碰撞能（Total Collision Energy，TCE）、振动态激发解离（Vibrationally Favored Dissociation，VFD）、广义碰撞能（Generalized Collision Energy，GCE）和Q-K（Quantum Kinetic）等模型来模拟化学反应过程，并逐步考虑振动-离解耦合效应。

迄今为止，中国少有同时采用CFD和DSMC方法模拟高速流动中的热化学非平衡效应，并进行对比分析。另一方面，考虑到高速流动的地面试验困难，综合使用CFD和DSMC方法恰好相互验证。本文以不同来流速度的高速圆柱绕流为研究对象，分别采用基于Park双温模型的CFD方法以及基于L-B内能松弛模型和Q-K化学反应模型的DSMC方法进行仿真，获得热化学非平衡效应对流场和气动热的影响，为高速气动热的预测提供参考。

1 模型和计算方法

1.1 物理模型和网格

采用来流工质为N2的二维高速圆柱绕流[1,2]，圆柱直径为0.3048 m，壁温为500 K。来流克努森数为Kn=0.01，来流速度取为Ma=5～45，其他来流参数保持恒定，具体见表1。

表1 来流条件 Tab.1 Conditions of Freestream Flow

因为圆柱后侧为亚声速流动，这将导致模拟收敛速度缓慢，而本文主要关注驻点处的压强和热流，所以仅对圆柱前半段进行仿真，采用贴体四边形网格，计算域和网格结构示意如图1所示。对于不同的来流速度工况，依据文献[1]中的 Recell,w=1准则设定靠近壁面附近的法向网格尺寸。最大网格长宽比不超过100，具体的网格数量和网格尺寸见表2。

图1 计算域及网格示意 Fig.1 Computational Domain and Grid Structure

表2 不同来流速度下的网格尺寸设置 Tab.2 Grid Size Setting for Different Freestream Flow Velocities

续表2

1.2 数值求解方法

CFD方法是在宏观上通过数值求解N-S方程，而DSMC方法是在微观上基于概率方法直接模拟分子的运动和碰撞。

开源程序OpenFOAM是一个高度可扩展的CFD软件开发包，并包含了Lagrange粒子库可实现DSMC方法。为了模拟高速流动，基于OpenFOAM框架，Casseau等[3,4]开发了适用于高速的双温CFD求解器hy2Foam；White等[5]发展了基于L-B过程的内能松弛过程和基于Q-K的化学反应过程dsmcFoamPlus求解器。本文同时使用这两个求解器进行对此圆柱绕流开展数值模拟研究。

1.2.1 CFD方法

CFD 方法采用有限速率化学反应的Park T-Tv 双温松弛模型[6]。双温模型中平动-振动的能量传输是由Landau-Teller方程决定的[7]，而它的松弛时间采用 Millikan-White-Park模型[8]。化学反应采用基于Arrhenius有限速率模型，且其速率系数取自DSMC的Q-K模型[9]。本文只考虑氮气的离解过程 （N2+N2—2N+N2、N2+N—2N+N），未考虑电离过程。

组分的粘性系数采用幂律模型，与DSMC的变径硬球分子的碰撞过程保持一致。热扩散系数与粘性系数间符合Eucken关系。混合气体的粘性系数和热扩散系数采用Wilke混合规则。质量扩散采用Fick定律，且其二元扩散系数则遵循Gupta的曲线拟合[10]。使用基于第一激发能级的简单谐振子模型和多电子激发态模型共同来确定气体的热力学参数。

不考虑壁面上的速度滑移，且认为是非催化壁。由于来流Re较小，流动属于层流，不用考虑湍流模型。

1.2.2 DSMC方法

DSMC方法使用L-B模型模拟各内能态的激发与松弛，与CFD的双温模型中假设分子的平动和转动能间是平衡相比，DSMC将平动能、转动能和振动能都独立分开。N2离解反应过程采用Q-K模型，设定VHS分子碰撞模型与CFD方法保持一致，并使用NTC方法模拟分子的运动和碰撞过程。除采用与CFD相同的网格外，DSMC算例还保证最少网格粒子数不少于7，时间步长小于最小平均碰撞时间。为了减小DSMC的统计误差，采样时间取为105倍的平均碰撞时间。

本文除了采用上述的可模拟热化学非平衡过程的DSMC方法外，还采用忽略离解反应的DSMC方法进行对比。结果中使用的简写记号的含义见表3，其中Tt、Tr、Tv分别表示在DSMC中独立考虑的气体分子平动、转动和振动的温度。而CFD采用Park T-Tv双温松弛假设，其认为平动态和转动态、振动态和电子态分别是平衡的，即采用2个温度Ttr和Tve表征。

表3 数值方法设置 Tab.3 Numerical Method Settings

1.3 计算方法验证

为了验证本文的CFD和DSMC程序的有效性，使用了美国Holden等[11]做的“CUBRC”第7次试车实验数据进行方法验证。“CUBRC”实验测量稀薄超声速来流对双圆锥体的气动力和气动热值，是一个典型的稀薄气体动力学实验。

图2为CFD和DSMC验证计算中双圆锥体表面压强、热流分布与实验对比。双锥体表面压强和热流密度计算值与实验值相符，最大值和最小值的位置一致，准确再现了双锥流动中的分离涡和激波-边界层干扰。

图2 双锥的气动参数分布 Fig.2 Distribution of Aerodynamic Parameters on the Double Cone

2 计算结果分析

图3为来流速度Ma=25且考虑化学反应的DSMC流场结果。结果表明，高超声速来流在圆柱头部形成弓形激波，经过激波后速度快速下降，而压强和温度快速上升，其中压强在壁面的驻点处达到最大值。由于壁面设定为恒温壁面，平动、转动和振动温度是先增加后减小，并考虑到内能松弛效应，平动、转动和振动温度存在明显差别。

图3 Ma=25时的DSMC模拟的流场 Fig.3 Flow Field from DSMC Simulation at Ma=25

图4为DSMC模拟的驻点线上平动、转动、振动温度和氮原子质量分数分布，温度使用来流静温无量纲化。结果表明，平动和转动温度间的非平衡效应主要出现在激波位置处，而振动温度在整个波后位置均表现出与前二者的非平衡。考虑氮气离解反应时各模式温度的极大值明显降低。

图4 Ma=25时的DSMC模拟的驻点线上的温度分布 和氮原子质量分数分布 Fig.4 Temperature and N% Distributions on the Standing Line for the DSMC Simulation at Ma =25

图5为各模式温度和N质量分数的最大值随Ma的变化，为了统一不同来流速度下的温度，这里的温度使用了来流总温无量纲化。图5a表明，DSMC与CFD在考虑氮气离解时，获得的氮原子质量分数与来流速度的关系基本相同，并且当Ma＜15时，氮气离解是可以忽略的，此时图5b中DSMC有无离解反应的平动温度rT以及转动温度tT是重合的。由于在DSMC中考虑了转动能的松弛，图5b中的转动温度tT的分布曲线总是稍低于平动温度rT，而平动温度rT会出现大于考虑平动-转动平衡的理想总温值的现象。当Ma＞15时，随着马赫数的增加，氮气离解越来越多，考虑化学反应的平动温度rT和转动温度tT与无化学反应得到的温度差距越来越大，DSMC在考虑氮气离解时所获得的平动温度rT和转动温度tT要比CFD的下降的快。图5c中DSMC在考虑氮气离解时所获得的振动温度vT也比CFD低，随着来流速度的增加，振动温度vT呈现先减小后增大再减小的变化规律，这是由于振动温度vT的最大值是受驻点线上流动特征时间、振动态的弛豫特征时间以及弛豫过程中离解反应对振动能消耗的综合影响。

图5 各模式温度和N质量分数的最大值随Ma的变化 Fig.5 Maxima of Temperture and N% with Freestream Velocity

续图5

图6为DSMC和CFD方法获得的驻点处压强系数Cp和热流系数CH。

图6 驻点处压强和热流的预测 Fig.6 Prediction of Pressure and Heat Flux at the Stationary Point

分别使用下式进行无量纲化：

式中wp为壁面压强；U∞为来流速度；wq为壁面热流。

由图6可知，高速条件下，驻点压强系数一直维持在1.80～1.85之间，呈现出马赫数无关效应，也不随化学反应而变化。对于热流系数，当Ma＜15时，两种方法预测的结果基本一致，随马赫数的增加而增大，并且化学反应效应不明显；在15≤Ma＜40时，化学效应显现，随着来流速度的增加，考虑化学反应的热流系数先减小后增大，且CFD和DSMC方法的规律基本一致，而忽略化学反应会造成热流预测偏高。

对比CFD和DSMC预测的表面参数，两者的结果基本是一致的，并表现出随来流速度相同的变化规律，但是两者对流场的预测出现一定偏差，需要进一步细致研究两者在物理模拟方法的差异。

3 结 论

本文采用可描述热化学非平衡过程的CFD和DSMC方法对不同来流速度下的高速圆柱绕流进行了数值模拟，获得了驻点线上的各模式温度、驻点气动参数随来流速度的变化关系，并讨论了热化学非平衡效应对气动热的影响，得到了如下结论：

a）CFD和DSMC方法分别在宏观和微观上发展了模拟热化学非平衡的模型，在相同的热物性设定基础上可以获得基本一致的气动参数分布；

b）在稀薄来流条件下，分子内能的激发与离解过程存在松弛效应，并且随着来流速度增大，热化学非平衡效应对气动热的准确预测影响越显著；

c）在充分考虑热化学非平衡过程时，热流系数随来流速度的增加呈现先增后减再增的变化规律，而忽略化学反应会造成热流预测偏高。