高超声速CO2稀薄气体振动激发与DSMC方法实现

2023-07-03王学德

弹道学报 2023年2期

谢东,王学德

(南京理工大学能源与动力工程学院,江苏南京 210094)

作为地球轨道外的第一近邻,火星的地质条件等与地球较为接近[1],是人类地外行星探测活动中的重要对象。探测器进入火星、减速过程穿越火星稀薄大气的历时很短,但由于探测器的速度极高,需要经受严酷的气动加热过程,因此该过程是探测任务中最危险也是最重要的环节[2]。

火星的大气成分与地球显著不同,主要由95.7%的CO2、2.7%的N2和其它稀有气体组成[3],其地面平均大气密度、大气压力约为地球的1%。美国早期的火星探测器飞行任务模拟主要采用以CO2大气环境修正的地面实验方法,但是其费用高昂、试验难度大。随着CFD技术的快速发展,美国从探路者号探测器(MPF)之后开展的气动评估,几乎全部基于计算流体力学和数值模拟[4]。目前国内主要采用求解连续Navier-Stokes方程和模拟稀薄气体流动的直接模拟蒙塔卡洛方法(DSMC)进行火星大气高温真实气体效应的研究,其中化学反应模型大多采用5组分17反应[5]或者8组分54反应[6]。

探测器进入火星大气层需要经历以CO2气体为主导的高超声速化学反应流动过程。由于CO2具有多个振动模态,热力学模型比地球大气更复杂,且CO2相较于其它双原子分子更容易离解,并在高温下发生置换、电离等复杂的化学反应,这使得火星大气环境下的热化学非平衡效应比地球再入过程更为严重。

因此本文旨在探究CO2分子在高温下的振动激发模式、化学反应模型以及在DSMC框架下的实现方式。构建了针对CO2气体环境下5组分17反应的热化学非平衡流动模型,编制了计算程序,并进行了验证,证明了该物理化学模型与计算方法的可行性。

1 CO2的振动激发

1.1 CO2分子的振动模型

CO2分子为直线型三原子分子,常温下具有3个平动自由度和2个转动自由度。CO2分子具有4种振动模态,其激发程度由温度决定。如表1所示,其中第1、2模态为退化的弯曲模态,特征温度均为960 K;第3模态为对称拉伸模态,特征温度为1 920 K;第4模态为反对称拉伸模态,特征温度为3 380 K;CO2的离解特征温度为66 000 K。低温下CO2的大部分振动能都储存在弯曲模态;在高温下,振动模式会被部分或完全激发,所以CO2分子的实际振动模式是这几种模态的组合。

表1 CO2分子的振动模态Table 1 Vibration modes of CO2molecules

1.2 CO2分子的振动能级与振动能

经典模型下,BORGNAKE和LARSEN[7]提出的L-B方法可以通过给每个分子分配连续振动能级,进而求得有效振动温度。然而此方法涉及基于碰撞总能的隐式温度,需要迭代求解。为此在DSMC方法中一般根据L-B模型给振动模态m分配离散的振动能级im:

im=⎣-ln(fR)T/Θm」

(1)

式中:⎣」表示取整,fR为(0,1)的随机数,T为分子热力学温度,Θm为对应振动模态的振动特征温度,对应的振动能为

εm=imkΘm

(2)

式中:k为玻尔兹曼常数。

由于实际的热化学反应流场温度可达到10 000 K,分子的振动激发十分激烈,分子的原子间距增大,因此本文考虑了振动能的非谐性效应,从MORSE势推导得到非谐振子模型[8]下各振动能级对应的特征温度:

(3)

式中:Θd为离解特征温度,η为非谐性因子,取值范围为[0,1]。当η=0时,上式对应简单谐振子模型。本文取η=0.90,则各模态的最高振动能级为

(4)

式中:Edis为离解能阈值。

图1为当η=0.90,离解温度为Θd/4时,CO2分子各模态的振动能级与相应能级的振动特征温度。由于CO2的光谱图并未显示出很高的振动能级,本文假设各振动模态相互独立,结合内能按自由度均分原则[9],将各模态振动能上限均设为Edis/4,也即各模态的离解特征温度为Θd/4。这一方法既保证了振动能级不会太高而超出实际,也考虑到了非谐性效应,各能级的能量间隔随着特征温度升高而变小。

图1 振动能级与对应的振动温度Fig.1 Vibration energy level and corresponding vibration temperature

1.3 CO2分子的振动弛豫

分子间的碰撞分为弹性碰撞和非弹性碰撞,非弹性碰撞会导致总能量在分子的平动模态和内部模态之间重新分配。在碰撞过程中一般会引入振动松弛因子Zv和转动松弛因子ZR来调整弹性碰撞和非弹性碰撞的比例,实现控制内能变化的速率,确保各参数不会失真的目的,各能量模式松弛概率满足:

(5)

式中:v为碰撞频率,τv为振动松弛时间。

CO2分子的转动松弛因子可取常数ZR=6.5,其它双原子分子为5[8]。而振动松弛因子需要根据流场的局部环境参数决定。MILLKAN-WHITE[10]在高温激波管内通过干涉仪观测了分子振动松弛过程,拟合了在10 000 K温度范围内,振动松弛时间τv与压强p和温度T的关系式:

(6)

式中:μr为分子碰撞对的折合质量。

为保证上式在高温下的松弛时间切合实际,PARK提出了高温修正因子τs[11]。而在DSMC方法中,振动松弛因子不允许小于转动松弛因子,因此本文设置了当Zv

此外由于CO2分子具有多个振动模态,为了简化多振动弛豫时间的处理,本文依据CAMAC[12]的实验结果,将CO2分子的4种振动模态以相同的速率松弛:

(7)

图2为CO2分子3种模态的弛豫时间和压力的乘积与温度的关系图。由于弯曲模态的振动温度最低,因此低温下大部分的振动能量都以这种模态存在,这导致低温下弯曲模态的弛豫速率和CAMAC的弛豫速率相差很小;而在高温下,3种模态的弛豫时间将会趋于一致,它们与CAMAC的弛豫速率相差也不太大。

图2 CO2各模态振动弛豫时间Fig.2 Relaxation time of various modes of CO2 vibration

2 碰撞后内能分配方法

本文根据L-B模型,将各内能模式逐一与平动能进行重新分配,这一过程中通过调整弹性碰撞与非弹性碰撞的松弛因子,来控制弛豫过程使得结果与实验相符合。下面以两个CO2分子碰撞后内能的分配为例,描述具体实现方法。

设ET为两个分子的碰前相对平动能;EVL和EVM分别是分子L和M的碰前总振动能,EVLm和EVMm为振动模态m(m=1,2,3,4)的振动能;ERL和ERM为碰前转动能;Mm为1.2节中各模态的最大振动能级。碰撞后依次考虑分子L和M的能量分配:

先考虑分子L(振动模态m=1,2,3,4):

①首先判断当Zv>fR时,才会发生能量分配,否则跳到步骤⑨;

②从m=1开始,设定初始碰撞总能量为EC=ET+EVLm;

③按式(3)确定振动能为EC时的振动能级,向下取整为imax;

④由im=fR(imax+0.999…),求得此振动模态的新能级和振动能E′>VLm;

⑤若im>Mm,返回步骤④;

⑥,根据分子对黏性-温度指数ωLM,计算概率P=(1-E′>VLm/EC)(1.5-ωLM)若P

⑦将剩余平动能ET返回第一步分配给其余模态的振动能,ET=EC-E′>VLm;

⑧所有振动模态能量分配结束,按照上述步骤逻辑进行T-R能量分配;

3 化学反应理论及实现

本文以纯CO2气体作为初始组分,高温下,气体分子将会发生离解反应、置换反应和复合反应等,主要涉及了CO2、CO、O2、C、O五种组分,共计17种基元反应,化学反应式如下,各反应的反应常数参见文献[13]。

CO2+GCO+O+G
CO+GC+O+G
O2+GO+O+G
CO2+OO2+CO
CO+OO2+C

(8)

式中:G代表第三体分子或者催化粒子。

3.1 离解反应

对于简单谐振子模型,一般认为当分子的碰后振动能级大于离解能所对应的振动能级时,离解发生,即:

imax>Θd/Θv

(9)

式中:imax是碰撞能EC(平动能+振动能)对应的最高振动能级。

由于本文采用的非谐振子模型,所以imax需要根据式(3)进行截断处理,并且由于CO2分子有4种振动模态,需要按照第2节中的逻辑,当总振动能εv>kΘd时,离解才发生。

3.2 复合反应

复合反应可以看作离解反应的逆反应,反应过程是一个三体碰撞过程。复合反应发生的概率,依据平衡碰撞理论[14],可简化为

(10)

式中:σR为化学反应截面,σT为碰撞截面,nM为第三体分子的数密度,ϖ和γ为不同反应的常数值[15]。利用网格累计概率法判断,当网格单元中的累加概率大于1时反应发生。

3.3 置换反应

置换反应是典型的双分子反应,其可根据反应能量的变化方向,分为吸热置换反应和放热置换反应两类。本文对吸热的置换反应采用TCE模型[16],最终简化得到反应发生的概率表达式为

(11)

式中:Ea为反应的活化能,Pf为吸热的置换反应发生概率。

对于放热的置换反应,由于其本身的活化能较低,反应过程又向周围散发热量,反应发生的概率也相对较大。本文简化得到的化学反应抽样几率公式[17]为

(12)

4 计算方法验证

4.1 计算条件

本文选取火星探路者号(MPF)为研究对象,其外形是一个顶角为140.38°的球形钝头大底,后接一个46.63°的锥形,探测器最大直径为2.65 m。其参数和几何外形如图3所示。来流参数见表2,参考文献[18],选取该探测器在65 km处来流为模拟状态。

图3 MPF探测器外形Fig.3 Shape of MPF detector

表2 来流参数Table 2 flow Parameters of incoming flow

来流气体为纯CO2气体,壁面条件为完全漫反射,不考虑催化。计算网格采用具有更好的贴体性的非结构网格[19]。分子模型采用变径硬球(Variable Hard Sphere,VHS)模型。化学反应模型采用TCE模型,转动松弛数设为5,振动松弛数依据1.3小节设置。

4.2 结果分析

计算结果如图4所示。图4(a)、(b)为流场速度分布云图和压力云图,MPF探测器迎风面出现了一道强烈的弓形激波面,气体的压缩效应明显。图4(c)～(f)分别为流场总温、平动温度、转动温度和振动温度的分布云图。激波后,气体被强烈压缩,分子数密度增加导致分子间的碰撞频率加快,激波附近出现局部高温区,且平动、转动和振动温度差异很大,流场中出现了显著的热力学非平衡现象。

图4 探测器流场分布云图Fig.4 Cloud map of detector flow field distribution

图5为沿驻点线附近的4种温度分布图。可以看出,由于远端来流中分子密度低、温度低,流场处于平衡态;平动温度的峰值远高于转动、振动温度,转动温度峰值高于振动温度,且各温度的峰值逐步靠近壁面,这是由于平动松弛速率高于内能松弛速率,转动弛豫速率高于振动弛豫速率导致的;随着流场靠近壁面,分子密度变大,能量交换速率加快,3种能量模式的温度逐步趋于一致。这4种温度的变化趋势及大小关系初步验证了本文将CO2离解能均分到4种振动模态,从而处理分子振动激发、化学反应及能量计算方法的正确性。

图5 驻点线温度Fig.5 Temperature of stagnation point line

图6和图7分别给出了MPF探测器物面压力和迎风面热流密度沿探测器上半壁面(从前驻点沿上物面至后驻点)的变化图,横坐标S表示壁面沿前驻点的距离。由图中可以看出,本文的计算结果与文献[18]中给出的结果基本一致,在前驻点迎风位置和探测器肩部拐角处的压力、热流密度与参考文献几乎一致,且在探测器尾部出现了0.1 Pa左右的近乎“真空”区域,符合实际情况。

图6 物面压力Fig.6 Surface pressure

图7 迎风面热流密度Fig.7 Heat flux density on the windward side

图8为沿驻点线上5种组分的摩尔分数变化图。CO2的摩尔分数沿驻点线逐渐减小,物面驻点位置的含量仅为0.1%,离解较为彻底。CO分子和O原子随着CO2的离解,摩尔分数随之升高,趋势较一致。但随着分子来到激波高温区,CO发生离解,使得O原子摩尔分数继续升高,CO分子的摩尔分数随之降低,C原子得以产生。当分子靠近壁面时,由于冷壁效应(壁面的温度远低于流场温度),CO分子的复合概率增加,C原子和O原子在壁面附近生成CO,致使CO的摩尔分数重新回升,C原子和O原子的比例随之下降。整个过程中,由于高温,很难发生O2的复合反应,所以探测器前端O2含量几乎为0。本文中各组分比例的变化规律与理论相符,且O原子的CO分子的最大摩尔比例分别为0.48和0.41,与文献[18]给出的结果相近,验证了本文计算方法的正确性。

图8 驻点线各组分摩尔百分数Fig.8 Molar percentage of each component on the stationary point line

5 不同参数对化学反应的影响

5.1 密度影响

为了进一步模拟火星探测器进入过程中面临的真实气动环境的影响,选取文献[18]中85～55 km海拔高度下的来流计算条件,探究密度变化对流场特性的影响,来流参数见表3。

表3 不同高度计算条件Table 3 Calculation conditions for different heights

图9为进入火星大气过程中,不同海拔高度下沿探测器前驻点线的流场总温变化趋势图。各个高度下的总温变化趋势一致,高温区域都集中在激波与物面之间的过渡区域。随着海拔下降,从75～55 km高度处,总温峰值温度逐渐下降,而85 km处总温分布却高于75 km处,这是因为总温是平动温度、转动温度和振动温度与各温度自由度乘积之和的平均值。

图9 不同高度下流场沿驻点线总温分布Fig.9 The total temperature distribution of the flow field along the stationary line at different heights

图10为不同高度驻点线上CO2分子的摩尔分数,其中55 km与65 km海拔下,CO2分子在激波内部几乎离解完全,在到达驻点位置时,含量在1%以下,并且由于55 km海拔处的分子压缩效应更强,所以该驻点线上CO2分子开始离解的位置比65 km海拔处更靠近物面。海拔75 km和85 km处激波后方的CO2分子已经无法完全离解,75 km海拔驻点处CO2分子的摩尔分数为15%左右,这些未离解的处于高振动激发态的CO2分子随流场运动到探测器肩部以后,将导致探测器尾部的流场总温升高。

图10 不同高度下驻点线CO2摩尔分数Fig.10 CO2 molar fraction at different altitudes of the stationary line

图11 不同高度下物面特性参数Fig.11 Surface characteristics parameters at different heights

5.2 速度影响

火星大气中的音速仅为240 m/s,选取65 km海拔处的火星大气条件,来流马赫数取25、27、29、31,计算模型仍为MPF探测器,其余参数见表3。

图12为不同马赫数下驻点线上CO2分子的摩尔比例。随着马赫数的增大,流场中的能量增加,提供了更多的能量供CO2分子发生离解,所以当来流数密度相同时,马赫数越大,流场中发生的化学反应越激烈,CO2分子的离解率越高。所以在激波区域相同的横坐标位置,马赫数越高的流场,CO2分子的摩尔比例越低。图13为不同马赫数下物面压力与迎风面热流密度变化曲线图。

图12 不同马赫数下驻点线上CO2的摩尔比例Fig.12 Molar ratio of CO2 molecules on the stationary line at different Mach numbers

图13 不同马赫数下物面特性参数Fig.13 Different Mach number object surface characteristics parameters

随着马赫数增加,物面压力和热流密度均增大,这与流场总温的变化趋势一致。但是在肩部往后的背风区域,4种流场中的物面压力及热流密度迅速减小且趋于一致。从驻点处的热流峰值来看,随着马赫数的增加,物面热流密度的增幅几乎保持一致。图13(a)中,4种流场的物面压力系数几乎没有差别,这符合牛顿流体在超音速下的流动特性;从图13(d)物面热流密度系数变化图来看,除了头部驻点附近的热流密度系数有些许偏差外,4种马赫数流场的迎风面热流密度系数的变化趋势以及数值大小均保持一致,这说明马赫数对物面的热流密度分布情况以及分子与物面的传能效率影响较小,仅仅对物面热流密度影响较大。