力学学报  2018 , 50 (4): 722-733 https://doi.org/10.6052/0459-1879-18-056

Orginal Article

高温气体热化学反应的DSMC微观模型分析

杨超12, 孙泉华12*

ANALYSIS OF DSMC REACTION MODELS FOR HIGH TEMPERATURE GAS SIMULATION 1)

Yang Chao12, Sun Quanhua12*

中图分类号:  O354.7, O362

文献标识码:  A

通讯作者:  *孙泉华, 研究员, 主要研究方向: 稀薄气体与非平衡流动. E-mail:qsun@imech.ac.cn*孙泉华, 研究员, 主要研究方向: 稀薄气体与非平衡流动. E-mail:qsun@imech.ac.cn

版权声明:  2018 力学学报期刊社 力学学报期刊社 所有

基金资助:  中国科学院战略性先导专项(XDAI7030100)和国家自然科学基金(11372325, 91116013) 资助项目.

展开

摘要

热化学耦合的非平衡现象一直是高温气体热化学问题研究的难点, 制约了诸如爆轰波胞格结构、低温点火速率等现象的分析. 本文以高温氮气离解和氢氧燃烧中的链式置换反应为例, 从微观反应概率、振动态指定的反应速率、热力学非平衡态的宏观反应速率、碰撞后的能量再分配等角度, 分析了直接蒙特卡罗模拟中的典型化学反应模型(TCE, VFD, QK模型)的微观动力学性质. 研究发现, 无论是高活化能的高温离解反应还是低活化能的链式置换反应, 实际参与反应的分子的振动能概率分布都偏离了平衡态的Boltzmann分布, 包含较强振动能额外影响的VFD模型可以很好地模拟高温离解反应, 而TCE (VFD的一个特例)和QK模型对活化能较低的链式置换反应的预测效果相对更好. 此外, 化学反应碰撞后的能量再分配应遵循微观细致平衡原理, 细微的偏差都可能造成平动能和振动能难以达到最终的平衡状态. 直接蒙特卡罗模拟的应用评估结果表明, 化学反应的振动倾向对热化学耦合过程产生了明显的影响, 特别是由于高振动能分子更多地参与了化学反应, 气体平均振动能的下降将影响后续化学反应的进行.

关键词: 热化学非平衡 ; 直接蒙特卡罗方法 ; 反应模型 ; 微观分析

Abstract

The non-equilibrium phenomenon of thermochemical coupling has been a difficult problem in high temperature aerothermal dynamics, and hinders to analyze phenomena such as cell structure of detonation wave and ignition speed of low temperature combustion. In this paper, typical chemical reaction models (TCE, VFD, QK models) employed in the direct simulation Monte Carlo (DSMC) simulation are analyzed using two examples (namely, N2 dissociation at high temperature, and chain displacement reaction in H2‒ O2 combustion) from microscopic reaction probability, vibrational state specific reaction rates, total reaction rate under thermal nonequilibrium condition, and post-collision redistribution of internal energy. It is found that the probability distribution of vibrational energy of reacted molecules deviates from the equilibrium Boltzmann distribution for both the high temperature dissociation reaction having high activation energy and the chain displacement reaction having low activation energy. The VFD model with strong vibrational favored contribution can predict well the high temperature dissociation reaction, whereas the TCE model (a special case of VFD model) and QK model are better for the chain displacement reaction. Besides, the post-collision redistribution of internal energy should follow the principle of detailed balance, as small deviations may cause inequality between the translational and vibrational energy under final equilibrium state. The DSMC simulation results also show that the vibrational favor of chemical reactions has an obvious effect on the thermochemical coupling process. Particularly, because molecules having high vibrational energy are more easily to have chemical reactions, the decrease of the average vibrational energy of the gas will affect the subsequent chemical reactions.

Keywords: thermochemical nonequilibrium ; DSMC method ; reaction model ; microscopic analysis

0

PDF (7207KB) 元数据 多维度评价 相关文章 收藏文章

本文引用格式 导出 EndNote Ris Bibtex

杨超, 孙泉华. 高温气体热化学反应的DSMC微观模型分析[J]. 力学学报, 2018, 50(4): 722-733 https://doi.org/10.6052/0459-1879-18-056

Yang Chao, Sun Quanhua. ANALYSIS OF DSMC REACTION MODELS FOR HIGH TEMPERATURE GAS SIMULATION 1)[J]. Acta Mechanica Sinica, 2018, 50(4): 722-733 https://doi.org/10.6052/0459-1879-18-056

气体分子在高温状态下往往发生复杂的物理化学变化, 既有分子内能状态的激发和松弛等热力学过程, 也会出现离解、复合和置换等化学反应过程. 这些热化学过程通过分子间的碰撞来实现, 在分子碰撞过程中可能发生热化学过程的耦合, 在宏观上表现为热力学非平衡和化学非平衡以及两者之间的相互耦合. 热化学非平衡耦合现象在高超声速流动中比较常见, 如再入飞行器头部激波层内的空气在超过5 000 K的环境中会经历强烈的热化学变化, 可显著影响波后的流场参数[1-3] (如温度、压力等). 近年来的研究还发现在燃烧问题的常见工况范围内(如1 000 ~4 000 K之间)热化学过程的耦合可能对火焰面[4]、 爆轰波5-6]等结构有不同程度的影响. 但由于热化学过程的微观机理复杂, 文献中有关理论及数值模拟结果还难以与实验观测结果完全符合, 高温热化学问题依然是目前研究的难点和热点.

热化学非平衡过程的宏观处理一般是在数值计算时引入Landau-Teller (L-T)方程求解内能松弛过程, 并通过双温[7]或三温[8]反应模型修正热非平衡态下的化学反应速率. L-T方程通过特征松弛时间来描述转动或振动温度的松弛速度, 无法给出内能在松弛过程中的概率分布或者相对平衡态分布的偏离程度. Park双温反应模型利用平动温度和振动温度的加权平均值 (T1-ξTvξ)来计算热力学非平衡态的反应速率, 其中的修正系数 ξ取0.3 ~0.5时定性上能较好地预测高温空气的化学反应过程. 但近年来Voelkel 等[9]的研究发现, 对于氢氧燃烧相关的一些链式反应, Park 模型中的修正系数要很小时才能使双温反应模型的反应速率与计算化学的准经典轨迹理论(QCT)结果较为接近, 其他含振动非平衡修正的反应模型如CVCV 模型也与计算化学结果存在明显偏差.

高温气体的热化学反应在微观上一般由反应物分子间的高能碰撞引起, 涉及分子的内能状态和碰撞方位等因素, 理论上可以通过计算化学方法得到任意反应在态与态层面的完整信息, 但是具体计算从建模到应用都十分繁琐, 目前阶段还难以用于实际流动问题的模拟. 文献中报道的微观层面的热化学反应研究, 大多是在直接蒙特卡罗模拟(DSMC)方法框架下对大量模拟粒子进行碰撞模拟, 通过统计分析可以得到反应系统(无论系统的热力学状态是否达到平衡)的宏观反应速率、热力学状态以及流动发展等信息. DSMC 方法模拟化学反应流的关键是建立合理的微观模型来处理分子运动和碰撞等过程, 主要的模型包括分子碰撞模型、化学反应模型等, 其中化学反应模型是模拟热化学问题的关键.

DSMC 方法是一种模拟气体粒子微观运动的统计模拟方法[10], 在稀薄气体领域[11]得到了广泛使用, 并在热化学非平衡问题以及流动的微观机理等研究中有重要应用. DSMC 方法通过模拟一定数量的粒子来代表真实气体分子, 假设分子运动和碰撞在极短时间内可以解耦, 并对模拟粒子的微观信息进行统计获得宏观信息. 模拟粒子通过碰撞模型(如VHS和VSS模型)来处理分子间的碰撞, 碰撞频率由分子运动论决定, 碰撞过程中是否发生热化学反应由热化学模型确定, 其中内能交换概率由转动能和振动能的松弛碰撞数 ZrZv确定, 而化学反应概率通过反应概率 Preac(或反应截面 σR)定量描述.

DSMC 模拟中的化学反应模型一般基于化学反应碰撞理论建立, 其中具有代表性的模型有TCE模型[10]、VFD模型[12]等. Bird早期提出的TCE模型, 假设反应概率 Preac是总碰撞能 Ec的函数, 函数中的待定参数可通过宏观反应速率的Arrhenius 公式确定. TCE 模型简单易用, 但在模拟振动与离解耦合问题(coupled vibration dissociation, CVD)时低估了振动能对离解反应的影响. 为此, Haas 和Boyd 提出了具有振动倾向(vibrational favored)的VFD 模型; Boyd 等[13]还将VFD 模型进一步发展成为同时考虑转动能和振动能影响的GCE 模型. 此外Bondar等[14]将VFD 模型推广到了存在多个振动模态的多原子分子. 其他考虑振动能对反应概率影响的模型还有KSS 模型[15]和Bias 模型[16]等. 近年来反应模型又有了新的发展[17-22], 如Bird 基于离散的振动能级提出了原理上不需要宏观反应速率标定系数的QK模型[17].

这些化学反应微观模型都是唯象模型, 其反应动力学性质并不一致, 模型的对比及验证一直是相关研究的重点内容[23-26]. 如Boyd[23]发现VFD 模型模拟得到的氮气离解过程的DSMC 结果能够与激波管实验结果较好地吻合; Kim和Boyd的研究[24] 揭示了TCE 模型明显低估了N 2分子在高振动能级的离解反应概率(以QCT结果为参考). Wysong 等[25]在研究飞行器再入稀薄气体条件下的生成NO的置换反应时, 发现QK模型由于忽略转动能对反应概率的影响导致预测的NO浓度过高, 而TCE模型预测结果与QCT数据较为吻合. 可见不同化学反应微观模型对于具体问题的有效性差异明显. 文献研究大多针对一个具体算例的模拟结果开展比较, 尚没有一个相对完整的化学反应微观模型评估研究的梳理和总结.

实际上, 以往针对化学反应微观模型动力学性质的研究较多地关注了热力学状态对反应速率的影响, 而忽视了化学反应对热力学状态影响的研究. 以分子离解过程为例, 由于高振动能的分子更容易发生离解反应, 剩余未参与反应分子的平均振动能会降低, 并使得系统出现振动温度低于平动温度的准稳态(QSS)[16-27]. 除了反应物, 生成物分子的内能分布也会对系统热力学状态产生影响, 而这种影响对多步反应系统如燃烧问题的微观模拟十分重要. 尽管Bird[28]指出碰撞后的能量再分配需要保证化学平衡时系统能够达到热力学平衡, 但如何实现碰撞过程的微观细致平衡尚未有明确报道. 此外, 对于化学反应微观模型, 以往研究主要关注空气组分的热化学反应, 对于燃烧机理的微观模拟研究[29-31]仍缺乏相应的动力学分析研究.

本文在文献研究工作基础上, 选取几种代表性的化学反应微观模型, 分析反应模型的微观动力学性质, 讨论能量再分配的技术细节, 并通过氮气离解的热化学耦合问题和氢氧预混自燃的DSMC 模拟对化学反应微观模型进行了应用评估, 以提升高温气体热化学流动的DSMC 反应模型的认识.

1 化学反应DSMC微观模型介绍

目前在DSMC模拟中常用的化学反应微观模型主要有TCE模型、VFD模型和QK模型, 其中VFD模型是对TCE模型的修正, 而QK模型相对较新但有待发展.

1.1 TCE模型

Bird在20世纪70年代提出的TCE (total collision energy)模型假设分子碰撞发生化学反应的反应概率(即反应截面/碰撞截面)是总碰撞能的函数. 这里的总碰撞能包括发生碰撞的分子对的相对平动能、分子转动能和振动能. TCE模型假设仅当总碰撞能 Ec超过某一阈能 εr时才可能发生化学反应, 其反应概率与总碰撞能满足如下的函数关系

Preac=C(Ec-εr)mEcn

式中 C,m,n为待定系数.

根据反应概率式(1), 在( εr, + ∞)区间上对自变量 Ec积分可以得到平衡态反应速率 kTCE

kTCE=vεrPreac(Ec)fEckBTdEckBT

式中 ν为分子碰撞率, f(Ec/kBT)为总碰撞能的概率分布函数, kB 为玻尔兹曼常数. 通过与实验数据拟合的反应速率Arrhenius公式 k=aTbe-Ea/(kBT)对比, 即可确定 C,m,n, 具体结果如下(采用VHS 模型确定碰撞率 ν)

其中 a,b,Ea为Arrhenius公式的参数, ω,Tref,σref为VHS模型的参数, ξi为分子内能的总自由度, mr为碰撞对的约化质量, ε为对称因子. 式(1)中的 εr为活化能 Ea.

TCE模型采用了分子内能为连续分布的假设, 而DSMC程序一般将振动能作量子离散化处理(如简谐振子SHO、非简谐振子AHO等), 因此TCE模型在实际程序中复现的平衡态反应速率与理论结果有一定偏差. Gimelshein等[32]为此提出了一种针对TCE模型的系数修正办法以减小内能离散化引起的误差. 实际上, TCE模型因振动能级量子化的误差影响一般在较低温度时才明显, 因此本文没有对TCE模型采用Gimelshein的修正.

TCE模型可以方便地将每个基元反应的Arrhenius反应速率转化为DSMC模拟所需的反应概率, 因此得到了广泛应用. 但TCE模型无差别的对待平动能、转动能和振动能, 导致其低估了振动能对化学反应的影响.

1.2 VFD模型

在TCE模型的基础上, Haas和Boyd提出了VFD(vibrationally favored dissociation)模型. VFD模型在反应概率中加入了振动能 Ev的单独影响. 譬如对于采用SHO振动模型的反应概率, 其形式如下

Preac=β(Ec-Ea)b+ξi/2+1/2Ecξi/2+3/2-ωEvEcϕ

式中 ϕ是可调整的系数, Ev为振动能, β为合并的常数项

其中 ξv为振动能自由度. 当 ϕ=0时, VFD模型退化为TCE模型.

与TCE模型相比, ϕ>0的VFD模型降低了低振动能级分子的反应概率并提高了高振动能级分子的反应概率, 能更加真实地反映离解反应的实际情况. 图1比较了TCE模型与VFD模型计算得到的 N2离解反应的反应概率, 其中VFD和TCE模型中的系数通过文献[16]的Arrhenius 公式确定. 从图中可以看出, 在总碰撞能相同条件下, 当 ϕ=2时, 高振动能级(20)的分子具有比低振动能级(10)的分子大得多的反应概率.

图1   VFD和TCE计算N2离解反应的反应概率对比.

Fig.1   Comparison of Preac for VFD and TCE model on N2 dissociation

1.3 QK模型

QK(quantum kinetic)模型是Bird近年来提出的一种化学反应模型, 是从量子化的振动能级出发, 直接构造离解、复合和置换等基元反应的反应概率. 下面对离解反应和吸热置换反应进行介绍,相应的逆反应模型只需要借助正反应模型和平衡常数即可确定.

对于离解反应AB+M †A+B+M, QK模型假设只有当相对平动能 Etr和离解分子的振动能 Ev之和超过离解分子的离解能时反应才发生且一定发生, 即

Preac=1,ifEc=Etr+Ev>Ediss

实际上分子碰撞能否发生化学反应不仅与碰撞能量有关还与碰撞方位角等有关, 所以为使QK模型预测的平衡态反应速率与实验结果吻合, 常要引入折减因子(reduction factor), 即令 Preac=α<1.

对于吸热的置换反应AB+C †AC+B, QK模型假设仅当碰撞能 Ec=Etr+Ev,AB大于反应活化能 Ea时发生反应, 其反应概率为

Preac=(1-Ea/Ec)3/2-ω/i=0imax(1-ikBθv/Ec)3/2-ω

式中 θv为AB分子的特征振动温度, imaxEa

imax=Etr+EvkBθvEa=reactionheat+Eb

Eb是为了和实验结果吻合更好而添加的可人为调整的系数, 一般为0.

QK模型虽然原理简单, 但需要引入修正系数提高模型预测的准确度, 仍是一种唯象模型.

2 化学反应微观模型的动力学分析

化学反应微观模型是建立在碰撞理论基础上的唯象模型, 具有一定假设甚至带有待定参数. 在应用微观模型时, 需要通过平衡状态的反应速率Arrhenius公式或实验结果来确定待定参数. 在参数确定后, 微观模型是否在非平衡条件或者分子碰撞层次具有模拟能力需要开展相应的动力学分析.

根据文献研究[23-26], 对于离解反应的模拟, VFD等具有振动倾向的反应模型比TCE 模型的DSMC模拟结果更接近实验数据; 而对于活化能较低的置换反应, 采用TCE或QK模型的模拟结果则与实验等其他来源的结果更为吻合. 因此本文选择典型的离解和置换反应来开展动力学分析和比较(复合反应一般不存在明显的振动倾向), 并侧重梳理动力学性质在不同层次间的联系. 这里的动力学性质主要是指振动态指定的反应速率和热力学非平衡态的宏观反应速率等.

选择的具体反应例子分别为 N2在20 000 K时的离解反应(反应#1)和氢氧燃烧在2 000 K时的一个链式反应(反应#2)

#1:N2+N2N+N+N2#2:H+O2O+H2

反应#1在氮气离解过程初期等 N2浓度较高的环境中起重要作用, 反应#2则是氢氧燃烧的三步链式反应的第一个反应. 20 000 K和2 000 K分别对应于飞行再入和燃烧的典型环境温度, 处于实际工程问题比较关注的温度范围. 两个反应的反应速率常数分别选择文献[16, 3333]中的数值, 据此确定微观模型中的参数, 特别是QK模型中的折减因子为0.47 (反应#1)和0.244 (反应#2).

当分子转动能和平动能处于平衡(多数问题均可满足)时, 宏观反应速率 k等价于振动态指定的反应速率 kiv在相应振动能级分布下的加权平均, 即

k(T)=iv[kiv(T)Pv(iv)]

式中 Pv为振动能级 iv的概率. 而 kiv的计算方法为

kiv=vmax{Ea-Ev,0}PreacfEtr+rotkBTdEtr+rotkBT

其中 f为平动能和转动能之和 Etr+rot的概率分布.

由于 Preac与碰撞能量(TCE模型)甚至与振动能(VFD和QK模型)有关, 在平动转动能总和相同时 Preac有赖于振动能的大小, 此外式(11)中的积分下限也与振动能有关, 所以 kiv与振动能大小密切有关. 图2分别对比了不同反应模型针对两个反应下的 kiv. 由图可知, 与TCE, QK模型相比, VFD模型提高了高振动能级分子的反应速率而降低了低振动能级分子的反应速率. 值得注意的是, 并不是振动能级越高的分子, 它的反应速率也越大. 如图2(b)所示, 置换反应#2的 kiviv( O2的振动能级)超过4以后就不再有明显提升, 甚至略微下降. 注意到 O2分子的 iv大于4时振动能已经大于活化能(反应#2的活化能远低于反应#1, Ea,#1约940 kJ, Ea,#2约70 kJ), 此时式(11)中的积分下限与振动能无关, 但 Preac随振动能的增大而略有减小. 以QK模型为例, 式(7)分母中大多数项随着振动能增大的速度稍大于分子中的对应项, 并且 imax随振动能的增大而增大, 总的效果是使 Preac随振动能的增大而有所减小.

图2   给定温度下反应模型对具体反应计算得到的kiviv的关系.

Fig.2   Relation of kiv and iv for given chemical reaction at specific temperature predicted with reaction model

根据图2, 不管反应#1还是反应#2, VFD模型中额外考虑了振动能的单独影响, 极大地改变了原来TCE 模型中高低振动能分子的反应速率, 当然改变幅度与振动能单独影响因子ϕ有关(ϕ=0对应于TCE 模型). QK模型要求相对平动能Etr和离解分子的振动能Ev之和超过反应的活化能时发生反应, 体现了振动能的影响, 此外在计算Preac中的imax与振动能显式相关, 说明QK 模型在一定程度上体现了振动能的优先影响. 但从图2的两个算例结果看, QK 模型预测结果与TCE模型结果比较接近, 振动能的额外影响存在但不显著.

当反应物中不同能态下的反应速率确定后, 反应物消耗的分子数就确定了, 但是碰撞能量大的分子对比较容易发生反应, 反应物中剩下分子的热力学是不平衡的. 为此, 根据 kiv进一步求出参与反应分子的振动能级分布 Pv(iv|reacted)(文中简称为反应振动能分布), 可以了解化学反应对反应物热力学状态的影响, 作为热化学耦合效应的参考.

图3给出了不同化学模型预测的两个反应在初始热力学平衡态下的反应振动能分布. 图中的Boltzmann分布是离散振动能级 iv在热力学平衡态下的分布

Pv,B(iv)=e-εv(iv)/kBTive-εv(iv)/kBT

式中 εv( iv)为能级 iv的振动能. 对于振动特征温度为 θv的SHO模型, εv(iv)=ivkBθv. 从图可以发现, 所有模型预测的结果都偏离了平衡态的Boltzmann分布. 对于离解反应#1, VFD模型预测的结果与计算化学QCT 结果[27]定性上最为吻合, 即振动能的反应分布随能级 iv先增大后减小, 且高振动能级分子对离解反应的贡献很大. 对于置换反应#2, 振动能级小于8的 O2分子占据了所有发生反应的 O2分子的压倒性绝大部分, 这一结果说明低振动能级的 kiv是影响热化学非平衡问题微观模拟结果准确性的关键. 虽然置换反应#2的反应振动能分布尚未在文献中发现可供对比的计算化学结果, 但在下文中通过双温反应速率的对比分析可以得知TCE和QK模型对反应#2 的描述较为合理, 因此图3(b)中这两个模型预测的反应振动能分布(单调递减分布)应更有参考价值. 从两个反应例子可以看出, 发生化学反应分子的振动能分布与具体化学反应密切相关. 对于高温离解反应#2, 发生化学反应的分子大多集中在高振动能级区, 而氢氧燃烧置换反应#2中的发生反应的分子几乎全部集中在低振动能级, 两者都明显偏离了平衡态的Boltzmann分布.

图3   不同化学模型预测的热力学平衡态下的反应振动能分布.

Fig.3   Pv(iv|reacted) predicted using different chemical model for reaction in a thermal equilibrium state

文献中对双温反应速率进行了较多地比较. 从式(10)也可以求出化学反应微观模型在给定平动温度 T和振动温度 Tv下对应的双温反应速率 k(T,Tv). 图4(a)和图4(b)分别为不同反应模型预测的反应#1 和#2 的双温反应速率. 对于离解反应#1, VFD模型( ϕ=2)的预测结果与经典的Park模型( ξ=0.5[7]) 以及QCT 数据[27,34] 吻合得很好. 而对于置换反应#2, TCE和QK 模型预测的双温反应速率与通过QCT 结果修正的Park模型( ξ=0.139[9])较为吻合, 其中TCE模型结果相对更好一些.

图4   不同化学反应模型预测的双温反应速率对比.

Fig.4   Comparison of k(T,Tv) predicted using different chemical model

3化学反应碰撞后的分子内能再分配

当化学反应速率相当或超过它的反应物或生成物的内能松弛速率时, 化学反应对系统热力学状态的影响将不可忽略. 在微观模拟中, 反应物消耗对热力学状态的影响是通过kivPv(iv|reacted) 等动力学性质体现的, 而生成物对热力学状态的影响则是通过化学反应后的能量再分配结合后续反应来体现. 因此化学反应后的能量再分配同样是热化学问题微观模拟的重要问题. 虽然非弹性碰撞的能量再分配已有经典的L-B方法[10], 但文献研究[12,28]指出对于VFD 和QK 模型直接使用L-B方法处理反应后的能量再分配无法保证模拟系统在化学平衡时能够趋于理想的热力学平衡状态, 即平动、转动、振动温度的相等. 为了满足化学平衡时的热力学平衡态, Boyd[12]和Bird[28]指出在化学平衡时通过正反应消耗和逆反应生成的同一组分的内能分布应完全一致. 这一结论实际上是微观细致平衡原理的必然要求, 但文献中难以找到根据这一原理研究化学反应能量再分配的具体工作, 因此本文开展能量再分配的细节工作.

对于任意一组化学反应(假设正反应为吸热反应)

A+B[reverse]forwardC+D

在化学平衡时需要满足3个必要条件(下文均以下标f和r代指正、逆反应):

(1)组分的质量守恒, 即宏观尺度对化学平衡的要求

kfnAnB=krnCnD

(2)总碰撞能为 Ec,f的正反应碰撞对的反应速率应该等于总碰撞能为 Ec,r的逆反应碰撞对的反应速率

vA,BnAnBPreac,f(Ec,f)f(Ec,f)vC,DnCnDPreac,r(Ec,r)f(Ec,r)=1

式中 Ec,f-Ec,r=Eheat,Eheat为反应热. 也就是总碰撞能为 Ec,f的A与B反应碰撞对构成的子集F和总碰撞能为 Ec,r的C与D反应碰撞对构成的子集R在反应前后存在对应关系.

(3)在条件(2)的基础上, 子集R中的组分C与D发生逆反应后变为组分A与B, 这些A, B的内能分布应该和子集F中A, B的内能分布相同才能最终确保满足细致平衡.

上述3个条件中, 条件(1)在达到宏观化学平衡时自动满足, 而条件(2)和条件(3)作为满足微观细致平衡的必要条件则需要进一步研究. 其中条件(2)能否满足是由正、逆反应的反应模型的 Preac直接决定, 条件(3) 则需要针对特定反应模型分别进行设计.

首先分析条件(2).由式(14)可知, 当正反应的 Preac由具体的反应模型确定后, 逆反应的 Preac实际也已确定. 因此DSMC模拟中常见的对正、逆反应分别使用具体反应模型的做法在理论上无法保证条件(2)的满足. 下文以对正、逆反应都采用TCE模型的做法为例进行分析.

将式(13)代入式(14), 再利用如下关系将反应速率 k进行替换

k=vEaf(Ec)Preac(Ec)dEc

就可以得到式(14)的等价形式

Preac,f(Ec,f)f(Ec,f)Ea,ff(Ec,f)Preac,f(Ec,f)dEc,f=Preac,r(Ec,r)f(Ec,r)Ea,rf(Ec,r)Preac,r(Ec,r)dEc,r

再将 f(Ec/kT)和TCE模型的 Preac的表达式代入式(16)后可得最终结果. 实际上, 式(16)的左右两边(分别对应于正、逆反应)可分别表示为

fforwardEx,fkBT=Ex,fkBTbf+ξi,f2+12e-Ex,fkBT/Ex,f>0Ex,fkBTbf+ξi,f2+12e-Ex,fkBTdEx,fkBT

freverseEx,rkBT=Ex,rkBTbr+ξi,r2+12e-Ex,rkBT/Ex,r>0Ex,rkBTbr+ξi,r2+12e-Ex,rkBTdEx,rkBT

式中 Ex,f=Ec,f-Ea,f,Ex,r=Ec,r-Ea,r.

上述推导中使用了隐含的假设 Ea,f-Ea,r=Eheat, 又由于 Ec,f-Ec,r=Eheat, 所以 Ex,f=Ex,r. 要使式(17)和式(18)相等, 只有对应的幂次 b+ξi/2+1/2也相等. 对于实际问题的正、逆反应, 这一条件一般是不满足的, 说明正、逆反应都采用TCE模型时条件(2)并不能严格满足, 但具体差异与 b,ξi的取值有关. 以上一节的反应#2 为例, 正、逆反应的 fforwardfreverseT=2000K时的结果如图5所示, 两者有差异但差别不大.

图5   fforwardfreverse的对比.

Fig.5   Comparison of fforward and freverse

针对条件(3), 不同反应模型需要分别设计对应的能量再分配方法. 对于TCE模型, 由于Preac 仅是Ec的函数, 因此在总碰撞能为Ec的碰撞对构成的集合内, Preac不会进一步影响参与反应的分子的内能概率分布, 因此TCE模型可以沿用非弹性碰撞的L-B方法. 但如果反应物中存在多个振动模态待分配能量, 每个振动模态分别和相对平动能进行能量再分配的L-B方法(serial distribution) 将会失效. 需要将所有待分配的振动模态看做一个统一的具有离散能级的振动模态来处理, 这在燃烧等涉及多原子分子(如H2O)的微观模拟中会遇到.

对于VFD模型, 由于 PreacEcEv都相关, 因此总碰撞能为 Ec的反应碰撞对的振动能分布会因为 Preac而进一步偏离平衡态, 而与之对应的逆反应必须使用对应的能量再分配方法才能保证逆反应的生成物具有与正反应的反应物相同的振动能分布(即条件(3)). 这里以正反应为离解反应为例, 给出逆反应(复合反应)的生成物分子的能量再分配方法.

考虑总碰撞能为 Ec的正反应碰撞对构成的集合, 其振动能级 iv的概率分布可以写为

P(iv|Ec)=C(Ec-Ev)5-2ω+ξrot2-1Evϕ

为了满足条件(3), 逆反应的生成物的振动能分布也应满足式(19). 在DSMC模拟中, 可以通过取舍法确定 iv所需使用的概率 Pa-r

Pa-r=1-EvEcζleft2-1EvEcϕϕζleft-1+ϕϕζleft-1ζleft-1+ϕζleft-1

式中 ζleft=5-2ω+ξrot.

为考察化学反应碰撞后的分子内能再分配的平衡情况, 采用VFD模型模拟了 T=10000K的等温热浴中的 N2的离解. 图6为离解复合达到平衡后离解反应消耗的(红色实心)和复合反应生成的(蓝色实心) N2分子的振动能级分布, 可以看出二者基本能够吻合, 化学平衡时振动温度也能够和平动温度基本一致. 这里的差别主要来源于条件(2)的误差. 但如果采用L-B方法, 复合反应生成的 N2的振动能级分布则会呈现出很大的偏差.

图6   复合反应生成的N2的振动能级分布

Fig.6   Distribution of iv of N2 through recombination reaction

4 应用评估

本节将反应模型应用于DSMC模拟中, 通过氮气离解的热化学耦合问题以及氢氧预混自燃问题来解释热化学效应对实际问题的影响.

具体问题的热化学非平衡程度在空间尺度上常采用努森数 Kn来表征. 但对于均匀各向同性问题, 热化学非平衡程度由时间尺度上的达姆科勒数 Da表征. 对于氮气离解问题中的振动能影响, Da可定义为振动松弛和离解反应的特征时间之比. 在30 000 K和20 000 K时 Da分别为0.86和0.25, 说明高温离解过程的热化学非平衡耦合会比较明显. 对于氢氧自燃问题的热化学非平衡程度, Da的定义为振动松弛特征时间和点火延迟时间之比. 以初始温度为2 000 K的自燃过程为例, 由氢气和氧气分别计算的 Da为0.96 和0.17, 说明自燃过程中也存在较强的热化学耦合.

4.1 氮气离解的热化学耦合问题

在高温时(>15 000 K)由于氮气的离解和振动松弛的特征时间比较接近, 会导致再入激波附近的非平衡区域存在所谓的振动与离解耦合现象(CVD). 这一现象对飞行再入中的高温流场有重要影响, 但由于存在复杂的微观物理和苛刻的高温环境, 实际流动问题的数值模拟结果和实验结果都存在一定误差, 二者在定量上的准确对比目前还有困难[26]. 本节的主要目的是通过空间均匀的热化学耦合问题对反应模型的动力学性质开展应用评估.

DSMC模拟的初始条件为数密度为 1025m-3的纯净氮气, 其平动、转动、振动温度分别为30 000 K, 20 000 K, 5 000 K, 初始模拟粒子总数为 1.0×106. DSMC模拟所需的转动能和振动能的松弛碰撞数 ZrZv由特征松弛时间 τrelax转化得到, 转换中都引入了相应的DSMC 微观修正[35]. 此外振动能的特征松弛时间 τv还考虑了高温修正[7]. 为了突出离解前期的振动与离解耦合过程, 只考虑了单向的离解反应, 忽略了对前期影响很小的复合反应. 另外由于QK模型和TCE模型的动力学性质比较接近, DSMC模拟只选择了TCE 和VFD模型进行评估比较.

计算结果如图7所示, 其中图7(a)给出了反应过程的 N2浓度变化, 图7(b)为平动温度和振动温度变化, 而转动温度很快和平动温度达到平衡并未在图中给出. 在反应前期, 振动温度明显低于平动温度, VFD模型预测的离解反应速率明显低于TCE模型结果, 反应吸热相对少, 因此系统整体温度高. 比较结果显示在同一时刻VFD和TCE模型预测的温度相差可达1000 K以上, 表明化学模型细节差异对于实际应用的重要性.

图7   N2离解过程中的浓度和温度变化 ||||(a) N2浓度随时间的变化;(b) 温度随时间的变化

Fig.7   Evolution of N2 concentration and temperature during N2 dissociation |||| (a)Evolution of N2 concentration;(b) Temperature history

通过微观模拟还可以分析化学反应对系统热力学状态的影响. 图8为采用VFD模型的算例在t=8.3955 ns 时N2分子的振动能概率分布. 与相同振动温度(12 698 K)的平衡态Boltzmann分布相比, 模拟结果在低振动能级与平衡态基本一致, 但高振动能级的概率明显降低. 原因是离解反应在短时间内消耗了超出平衡的高振动能级的N2分子, 而热松弛过程(非弹性碰撞)尚未补充恢复. 这说明即使采用双温反应速率的宏观化学反应模型, 在计算热化学耦合问题时仍会带来一定的误差, 而微观粒子模拟或态指定的反应速率求解方法则能包含更多反应细节、模拟结果更为可信.

图8   VFD模型计算的N2离解过程中振动能级的瞬时概率分布

Fig.8   Vibrational energy distribution of N2 during dissociation with VFD model

4.2 氢氧预混自燃问题

在氮气离解问题中的热化学耦合机理相对简单, 一方面离解产物不涉及内能分配, 另一方面逆反应以及化学平衡的影响也较弱. 但对于多步反应或反应机理较为复杂的热化学问题, 如自燃、气相爆轰等燃烧现象, 上述影响将变得重要. 其中预混气体的自燃作为燃烧中的基本现象, 不涉及流动过程, 是研究燃烧中的热化学现象以及检验微观模拟方法的理想问题. 因此选取氢氧预混气体自燃问题, 进一步分析热化学耦合对化学反应模拟的影响.

DSMC模拟的计算条件如下. 初始条件为 T=2000K, n=1025m-3, H2O2的摩尔数之比 xH2: xO2=2:1, 粒子总数为 106, 采用Maas 和Warnartz 的详细反应机理[36]和TCE模型. 气体分子的振动松弛时间 τv一般由Millikan-White 公式确定, 但对于 H2, 此公式与实验结果偏差较大, 因此针对相关的 H2- H2H2- O2的特征松弛时间进行了修正. 其中 τv,H2-O2由实验值[37]确定, 而 τv,H2-O2由于缺少实验结果, 选择了已有实验结果中分子量最接近的 τv,H2-Ar的数值.

图9给出了DSMC模拟的自燃过程的组分和温度变化情况. 总体上, 各种组分的变化和反应动力学结果类似. 但在点火过程中, H2O2的振动温度出现了先降低再升高并最终与平动温度趋于平衡的过程. 出现这一现象的原因是高振动能级的 H2O2分子更容易被链式反应消耗, 导致剩余未反应的分子的平均振动温度降低. 具体来说, 在燃烧过程中随着自由基数量(H/O/OH)的增加, 链式反应速率迅速增大; 当反应速率接近并超过 H2O2的热松弛速率时, 化学反应对热力学状态的这种影响就会显现出来; 而振动温度的降低, 会导致诱导阶段后期相关链式反应速率下降(低于平衡态的Arrhenius 式反应速率), 并最终影响到点火延迟时间的计算结果.

图9   氢氧预混自燃过程||||(a) 组分浓度;(b) 温度.

Fig.9   Auto-ignition of premixed H2-O2 mixture|||||(a) Species concentration;(b) Temperature

值得注意的是, 目前结果中组分的振动温度和平动温度在化学平衡时(4 μs)并不能完全趋于一致, 约有150 ~ 200 K (约5%)的差别, 主要原因是化学反应碰撞后的内能再分配时的细致平衡没有严格达到, 这一点在第3节中已经解释. 实际上, 采用Bird 2013年专著[28]附带的DSMC程序计算得到的振动平动能差别可达500 K以上. 显然, 分子碰撞后的能量再分配技术仍需要发展和提高.

5 结 论

本文针对高温化学反应中的热力学非平衡和化学非平衡耦合问题, 以高温氮气离解和燃烧中的链式置换反应为例, 从振动态指定的反应速率、热力学非平衡态的宏观反应速率、碰撞后的能量再分配等角度, 分析了DSMC模拟中的典型化学反应模型(TCE, VFD, QK模型)的微观动力学性质.

研究发现, 无论是低活化能的置换反应还是高温离解反应在化学反应概率的计算上都需要考虑碰撞能以及振动能的额外单独影响, 而且反应分子的振动能分布都偏离了平衡态的Boltzmann分布. 只是高温离解反应中发生化学反应的分子大多集中在高振动能级区, 而氢氧燃烧置换反应中发生反应的分子几乎全部集中在低振动能级. 由于VFD模型包含了振动能的额外影响, 在振动能单独影响因子 ϕ=2时可以较好地模拟高温离解反应, 而对活化能较低的置换反应, TCE (即 ϕ=0时的VFD)和QK模型的预测结果相对更好.

研究还发现, 化学反应碰撞后的能量再分配可对实际模拟产生较大的误差影响. 主要原因是能量再分配影响了系统的热力学平衡, 进而影响后续化学反应. 能量再分配应遵循微观细致平衡原理, 细微的偏差都可能造成平动能和振动能难以达到最终的平衡状态.

通过氮气离解和氢氧预混自燃的DSMC模拟, 可以认为化学反应普遍存在热化学非平衡的耦合作用, 而多温度宏观模型难以刻画微观非平衡细节, 应大力发展微观模拟技术, 揭示非平衡效应明显的化学反应问题的机理和特性.

The authors have declared that no competing interests exist.


参考文献

[1] 陈松, 孙泉华.

高超声速飞行流场中的最大氧离解度分析

. 力学学报, 2014, 46(1): 20-27

(Chen Song, Sun Quanhua.

Analysis of maximum dissociation degree of oxygen during hypersonic flight

.Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(1): 20-27 (in Chinese))

[2] 彭傲平, 李志辉, 吴俊林.

含振动能激发Boltzmann模型方程气体动理论统一算法验证与分析

. 物理学报, 2017, 66(20): 204703

(Peng Aoping, Li Zhihui, Wu Junlin, et al.

Validation and analysis of gas-kinetic unified algorithm for solving Boltzmann model equation with vibrational energy excitation

.Acta Physica Sinica, 2017, 66(20): 204703 (in Chinese))

[3] 张子健, 刘云峰, 姜宗林,

振动激发对高超声速气动力/热影响

. 力学学报, 2017, 49(3): 616-626

(Zhang Zijian, Liu Yunfeng, Jiang Zonglin.

Effect of vibration excitation on hypersonic aerodynamic and aerothermodynamic

.Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(3): 616-626 (in Chinese))

[4] Fiévet R, Voelkel S, Koo H, et al.

Effect of thermal nonequilibrium on ignition in scramjet combustors

.Proceedings of the Combustion Institute, 2017, 36(2): 2901-2910

[本文引用: 1]     

[5] Shi L, Shen H, Zhang P, et al.

Assessment of vibrational non-equilibrium effect on detonation cell size

.Combustion Science and Technology, 2017, 189(5): 841-85

[本文引用: 1]     

[6] 方宜申, 胡宗民, 滕宏辉.

圆球诱发斜爆轰波的数值研究

. 力学学报, 2017, 49(2): 268-273

(Fang Yishen, Hu Zongmin, Teng Honghui, et al.

Numerical study of the oblique detonation initiation induced by spheres

.Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(2): 268-273 (in Chinese))

[7] Park C.

Assessment of two-temperature kinetic model for ionizing air

.Journal of Thermophysics and Heat Transfer, 1989, 3(3): 233-244

[本文引用: 2]     

[8] Park C. The limits of two-temperature model. AIAA Paper, 2010

-

911, 2010

[本文引用: 1]     

[9] Voelkel S, Raman V, Varghese PL.

Effect of thermal nonequilibrium on reactions in hydrogen combustion

.Shock Waves, 2016, 26(5): 539-549

[本文引用: 1]     

[10] Bird GA.Molecular Gas Dynamics and the Direct Simulation Monte Carlo of Gas Flows. Oxford: Clarendon Press, 1994

[本文引用: 4]     

[11] 樊菁.

稀薄气体动力学: 进展与应用

. 力学进展, 2013, 43(2): 185-201

[本文引用: 2]     

(Fan Jing.

Rarefied gas dynamics: Advances and applications

.Advances In Mechanics, 2013, 43(2): 185-201 (in Chinese))

[本文引用: 2]     

[12] Haas BL, Boyd ID.

Models for direct Monte Carlo simulation of coupled vibration-dissociation

.Physics of Fluids A: Fluid Dynamics, 1993, 5(2): 478-489

[本文引用: 4]     

[13] Boyd ID, Bose D, Candler GV.

Monte Carlo modeling of nitric oxide formation based on quasi-classical trajectory calculations

.Physics of Fluids, 1997, 9(4): 1162-1170

[本文引用: 2]     

[14] Bondar Y, Gimelshein N, Gimelshein S, et al.

On the accuracy of DSMC modeling of rarefied flows with real gas effects

.AIP Conference Proceedings, 2005, 762(1): 607-613

[本文引用: 2]     

[15] Bondar YA, Ivanov MS.

DSMC dissociation model based on two-temperature chemical rate constant

. AIAA Paper, 2007-614, 2007

[本文引用: 2]     

[16] Wysong IJ, Gimelshein SF.

Comparison of DSMC reaction models with QCT reaction rates for nitrogen

.AIP Conference Proceedings, 2016, 1786(1): 050021

[本文引用: 5]     

[17] Bird GA.

The QK model for gas-phase chemical reaction rates

.Physics of Fluids, 2011, 23(10): 106101

[本文引用: 3]     

[18] Baikov BS, Bayalina DK, Kustova EV, et al.

Inverse Laplace transform as a tool for calculation of state-specific cross sections of inelastic collisions

.AIP Conference Proceedings, 2016, 1786(1): 090005

[19] Luo H, Kulakhmetov M, Alexeenko A.

Ab initio state-specific N2+O dissociation and exchange modeling for molecular simulations

.The Journal of Chemical Physics, 2017, 146(7): 074303

[20] Sebastião IB, Kulakhmetov M, Alexeenko A.

DSMC study of oxygen shockwaves based on high-fidelity vibrational relaxation and dissociation models

.Physics of Fluids, 2017, 29(1): 017102

[21] Ramin Z, Kamali-Moghadam R, Mani M.

A new approach for chemical reaction simulation of rarefied gas flow by DSMC method

.Computers & Fluids, 2016(140): 111-121

[22] Sebastiao IB, Luo H, Kulakhmetov M, et al.

DSMC implementation of compact state-specific N&lt;inline-formula&gt;&lt;mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml229-0459-1879-50-4-722"&gt;&lt;mml:msub&gt;&lt;mml:mrow&gt;&lt;mml:mi mathvariant="normal"&gt; &lt;/mml:mi&gt;&lt;/mml:mrow&gt;&lt;mml:mrow&gt;&lt;mml:mn&gt;2&lt;/mml:mn&gt;&lt;/mml:mrow&gt;&lt;/mml:msub&gt;&lt;/mml:math&gt;&lt;/inline-formula&gt;+O dissociation and exchange models//55th

AIAA Aerospace Sciences Meeting, 2017

[23] Boyd ID.

Analysis of vibration-dissociation-recombination processes behind strong shock waves of nitrogen

.Physics of Fluids A : Fluid Dynamics, 1992: 4(1): 178-185

[本文引用: 4]     

[24] Kim JG, Boyd ID.

Monte Carlo simulation of nitrogen dissociation based on state-resolved cross sections

.Physics of Fluids, 2014, 26(1): 012006

[本文引用: 2]     

[25] Wysong I, Gimelshein S, Gimelshein N, et al.

Reaction cross sections for two direct simulation Monte Carlo models: Accuracy and sensitivity analysis

.Physics of Fluids, 2012, 24(4): 042002

[本文引用: 2]     

[26] Wysong I, Gimelshein S, Bondar Y, et al.

Comparison of direct simulation Monte Carlo chemistry and vibrational models applied to oxygen shock measurements

.Physics of Fluids, 2014, 26(4): 043101

[本文引用: 2]     

[27] Valentini P, Schwartzentruber TE, Bender JD, et al.

Direct molecular simulation of nitrogen dissociation based on an ab initio potential energy surface

.Physics of Fluids, 2015, 27(8): 086102

[本文引用: 2]     

[28] Bird GA.

The DSMC Method

.Create Space Independent Publishing Platform, 2013

[本文引用: 5]     

[29] Bird GA.

Chemical reactions in DSMC

.AIP Conference Proceedings, 2011, 1333(1): 1195-1202

[本文引用: 2]     

[30] Bondar YA, Maruta K, Ivanov MS.

Hydrogen-oxygen detonation study by the DSMC method

.AIP Conference Proceedings, 2011, 1333(1): 1209-1214

[本文引用: 1]     

[31] Yang C, Sun QH.

Investigation of spontaneous combustion of hydrogen-oxygen mixture using DSMC simulation

.AIP Conference Proceedings, 2014, 1628(1): 1261-1267

[本文引用: 1]     

[32] Gimelshein SF, Gimelshein NE, Levin DA, et al.

On the use of chemical reaction rates with discrete internal energies in the direct simulation Monte Carlo method

.Physics of Fluids, 2004, 16(7): 2442-2451

[本文引用: 2]     

[33] Saxena P, Williams FA.

Testing a small detailed chemical-kinetic mechanism for the combustion of hydrogen and carbon monoxide

.Combustion and Flame, 2006, 145(1): 316-323

[本文引用: 1]     

[34] Bender JD, Valentini P, Nompelis I, et al.

An improved potential energy surface and multi-temperature quasiclassical trajectory calculations of &lt;inline-formula&gt;&lt;mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml230-0459-1879-50-4-722"&gt;&lt;mml:msub&gt;&lt;mml:mrow&gt;&lt;mml:mi mathvariant="normal"&gt;N&lt;/mml:mi&gt;&lt;/mml:mrow&gt;&lt;mml:mrow&gt;&lt;mml:mn&gt;2&lt;/mml:mn&gt;&lt;/mml:mrow&gt;&lt;/mml:msub&gt;&lt;mml:mo&gt;+&lt;/mml:mo&gt;&lt;mml:msub&gt;&lt;mml:mrow&gt;&lt;mml:mi&gt;N&lt;/mml:mi&gt;&lt;/mml:mrow&gt;&lt;mml:mrow&gt;&lt;mml:mn&gt;2&lt;/mml:mn&gt;&lt;/mml:mrow&gt;&lt;/mml:msub&gt;&lt;/mml:math&gt;&lt;/inline-formula&gt; dissociation reactions

.The Journal of Chemical Physics, 2015, 143(5): 054304

[35] Gimelshein NE, Gimelshein SF, Levin DA.

Vibrational relaxation rates in the direct simulation Monte Carlo method

.Physics of Fluids, 2002, 14(12): 4452-4455

[本文引用: 1]     

[36] Maas U, Warnatz J.

Ignition processes in hydrogen oxygen mixtures

.Combustion and Flame, 1988, 74(1): 53-69

[本文引用: 1]     

[37] Dove JE, Teitelbaum H.

The vibrational relaxation of H&lt;inline-formula&gt;&lt;mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml231-0459-1879-50-4-722"&gt;&lt;mml:msub&gt;&lt;mml:mrow&gt;&lt;mml:mi mathvariant="normal"&gt; &lt;/mml:mi&gt;&lt;/mml:mrow&gt;&lt;mml:mrow&gt;&lt;mml:mn&gt;2&lt;/mml:mn&gt;&lt;/mml:mrow&gt;&lt;/mml:msub&gt;&lt;/mml:math&gt;&lt;/inline-formula&gt;. I. Experimental measurements of the rate of relaxation by H&lt;inline-formula&gt;&lt;mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" id="Mml232-0459-1879-50-4-722"&gt;&lt;mml:msub&gt;&lt;mml:mrow&gt;&lt;mml:mi mathvariant="normal"&gt; &lt;/mml:mi&gt;&lt;/mml:mrow&gt;&lt;mml:mrow&gt;&lt;mml:mn&gt;2&lt;/mml:mn&gt;&lt;/mml:mrow&gt;&lt;/mml:msub&gt;&lt;/mml:math&gt;&lt;/inline-formula&gt;, He, Ne, Ar, and Kr

.Chemical Physics, 1974, 6(3): 431-444

[本文引用: 1]     

/