0 引言

稀薄气体电离现象是航天器再入速度持续增大面临的新问题，直接导致通信黑障等传统难题向稀薄区域大幅延伸。传统的再入物理^[1]通过求解含化学反应的N-S方程组^[2]得到再入体绕流电子数密度分布，但是该方程组对于稀薄气体流动失效。近年来发展的稀薄气体数值模拟方法中，基于Boltzmann方程的分析和计算^[3-4]尚无法包含电离过程，只有基于粒子随机统计模拟的DSMC(Direct Simulation Monte Carlo)方法^[5-9]是最有可能解决稀薄气体电离数值预测难题的手段。

第一宇宙速度下的稀薄气体电离度极低，通常在10^-5量级，相对于N₂、O₂等来流气体组分以及反应生成的NO、N和O，电子和离子属于稀有组分。使用传统DSMC方法模拟计算网格中的一个电子，将需要10⁵个其它粒子，计算代价过高而不可接受。基于这一原因，早期稀薄气体电离过程的DSMC模拟^[10-11]主要针对10 km/s以上的低维再入流动，此时电离度在1%量级，一维驻点线和轴对称计算的稀有组分统计涨落可以通过增大仿真分子数目控制在可以接受的范围内。文献[12]建立了三维复杂外形航天器再入稀薄电离特性分析的DSMC方法和并行计算平台，但是文中结果表明，RAM-C Ⅱ的电子数密度分布存在极为显著的统计涨落。

电子等稀有组分粒子数密度的统计波动带来若干不利影响。首先，真实的电子数密度分布应该较为光滑，存在较d大波动的电子数密度分布在物理上是不真实的，这种统计波动是数值原因造成的，应该可以用数值的方式予以抑制。第二，电子和离子在网格中统计数目的波动，给网格中分子碰撞对数目的确定带来了一定程度的不确定性，而分子碰撞对数目的正确计算是含电离化学反应过程正确模拟的基础，亦即稀有组分粒子数密度的统计波动可能导致物理上不真实的化学反应计算。第三，从应用角度说，比如电子数密度对通信设计的预测分析，如果电子数密度分布本身波动较大，在应用上的指导意义也大打折扣。

现有稀薄气体DSMC模拟的稀有组分处理方法可以笼统地分为三种类型：间接方法、权重格式和权重格式的改进格式。间接方法不直接模拟稀有组分粒子，如Bird的工作^[13]并不真正计算NO⁺和电子，只是在N和O碰撞时标识其发生电离的比例，根据这一比例和N、O数密度计算电子数密度分布，然而这一处理方式仅考虑N和O的联合电离，导致其预测的电子数密度低估了2到10倍^[14]，而且这一方法无法处理离子的电荷交换过程；Shevyrin的工作^[15]不将电子作为一种粒子模拟，其数密度由网格中的离子之和得出，记录电子能量用以计算电子温度，对联合电离逆反应单独处理，这一方法显然未包括电子-原子碰撞导致的直接电离。

权重格式最早由Bird提出^[16]，用每个仿真粒子代表的真实粒子数F_Ni作为粒子权重，权重较小的粒子对应稀有组分。两个不同权重的粒子发生碰撞时，将权重较大的粒子一分为二，其中一个与稀有组分具有相同权重；让两个具有相同权重的粒子二体碰撞，然后将碰后非稀有组分粒子与分裂得到的另一个粒子融合。Boyd指出^[17]，上述格式不能保证每次碰撞动能守恒，进而提出了守恒格式，即将上述每次碰撞的动能损失记录下来，添加到下一个时间步长的常规组分碰前相对平动能中。随后，Boyd等将其守恒权重格式发展用于模拟OH紫外辐射的变权重形式^[18]。受当时的研究局限，上述权重格式均未拓展到含电离的化学反应过程。

权重格式的另一种形式是对所有组分使用相同的F_N，但是对每个组分设置不同的权重W_i；常规组分的权重为1，稀有组分的权重小于1.这一方法在稀薄气体电离DSMC模拟中的应用可以追溯到Bartel等人的工作^[19]，但他们并未考察不同权重粒子之间的碰撞细节。Boyd^[14]在Bird^[16]的基础上，人为增大电离化学反应概率R倍，同时复制电离反应产物S次，即赋予带电粒子权重为1/(R×S)；这一权重用于确定正确的碰撞对分子数目，两个不同权重的粒子碰撞过程仍采用Bird的分裂思想。

权重格式的改进包括Boyd等的重叠方法^[20]，而樊菁等的TSS算法^[21]可以看成是权重格式的最新发展。重叠方法是一种两步流场计算方法，第一步(基础模拟)使用较大权重，不考虑稀有组分；第二步(重叠模拟)使用较小权重，常规组分的产生受到抑制。在重叠模拟中，常规-稀有组分碰撞中的常规组分粒子由根据网格中常规组分属性抽样产生的虚拟粒子代替，常规-常规组分碰撞生成稀有组分粒子数目由网格中的反应物粒子数密度等参数计算得到。重叠方法主要用于辐射模拟，后来又发展处理表面过程^[22]，但尚未见到处理稀薄气体电离过程的报导。TSS方法与常规DSMC的主要区别是稀有组分的产生不依赖碰撞过程，而根据网格中反应物的粒子数目计算，稀有组分的内能状态根据当地宏观参数抽样得到。

本文工作基于稀有组分权重格式处理稀薄气体含电离化学反应过程。与前述方法不同的是，本文在前期工作^[12]的基础上，对不同权重粒子之间的碰撞依概率改变内能状态，对不同权重粒子之间的化学反应依概率删除或保留反应物粒子，同时对生成物中的稀有组分粒子进行复制。该方法在对成熟DSMC模拟框架和程序作极小化修正的前提下，可显著提高电子等稀有组分数密度等值线光滑性和计算精度。

1 稀薄气体电离反应的DSMC方法 1.1 化学反应模型

再入过程的纯大气组分二体碰撞化学反应，通常可以写成如下形式

$ {\rm{A + B}} \to {\rm{C + D + E}} $

(1)

涉及电离的化学反应包括五种：双原子分子的离解反应、置换反应、联合电离反应及其逆反应、直接碰撞电离反应、离子交换反应。对于离解反应，E代表原子；对于直接电离反应，E代表电子；而对于其它反应，E为空。

化学反应速率通常表示为如下Arrhenius形式的方程^[23]

$ k\left( T \right) = \mathit{\Lambda }{T^\eta }{\rm{exp}}({\rm{ - }}{E_a}/\kappa T) $

(2)

其中Λ和η为常数，E_a为反应活化能，κ为Boltzmann常数，T为温度。当碰撞总能量E_c＞E_a时，上述常数Λ和η通过下式与化学反应概率P相关联^[24]

$ \begin{array}{l} P = \frac{{{\pi ^{1/2}}\varepsilon \mathit{\Lambda }T_{{\rm{ref}}}^\eta }}{{2{\sigma _{T, {\rm{ref}}}}\left( {\kappa {T_{{\rm{ref}}}}} \right)}}\frac{{\mathit{\Gamma }\left( {\bar \xi + 5/2 - {\omega _{AB}}} \right)}}{{\mathit{\Gamma }\left( {\bar \xi + \eta + 3/2} \right)}}{\left( {\frac{{{m_r}}}{{2\kappa {T_{{\rm{ref}}}}}}} \right)^{1/2}} \cdot \\ \;\;\;\;\;\;\;\;\;\frac{{{{({E_c} - {E_a})}^{\eta + \bar \xi + 1/2}}}}{{{E_c^{\eta + \bar \xi + {\omega _{AB}}}}}} \end{array} $

(3)

上式中σ_T为总截面；ε为对称因子，当A和B相同时，ε取值为2，否则为1；T_ref为参考温度；Γ为伽马函数，ξ为平均内自由度，ω_AB为黏性指数平均值；m_r为折合质量，取值为m_r=m_Am_B/(m_A+m_B)。

1.2 弱等离子体效应简化处理

电子与离子及常规组分的巨大质量差异是稀薄气体电离过程DSMC模拟的本质困难，电子真实质量的使用和运动轨迹的捕捉要求极短的时间步长，进而导致不可接受的计算代价。受Boyd工作^[14]启发，本文采用增大电子质量三个数量级的方式对上述等离子体效应简化处理。电子的真实质量为9.11×10^-31 kg，增大三个数量级至9.11×10^-28 kg。以氧原子O为例，其质量为2.656×10^-26 kg，调整后的电子质量约为氧原子质量的3.4%。在温度相同的情况下，调整后电子的随机热运动速度均值比氧原子大不到一个数量级。同时，为保证化学反应过程质量守恒，增大电子质量须同步减小离子质量。以氧离子O⁺为例，真实质量为2.656×10^-26 kg，调整后的质量为2.5649×10^-26 kg，质量变化约为3.4%，亦不会对其热运动速度造成明显影响。

电子质量增大导致的一个结果是涉及电子碰撞的折合质量m_r增大，由式(3)知其将导致涉及电子化学反应概率的增大。为保持真实反应概率，可调整涉及电子的化学反应速率系数Λ，简单而直接的做法是将Λ减小$ \sqrt {1000} $倍。对于空气中的11组元(N₂、O₂、NO、N、O、N₂⁺、O₂⁺、NO⁺、N⁺、O⁺和电子e)，采用Gupta模型^[25]，涉及电子的相关化学反应常数如表 1所示，其中Q_rea为反应热。

2 含电离化学反应的稀有组分权重因子方法 2.1 权重因子对粒子碰撞的影响

本文采用权重格式的第二种形式，即对每个组分设置不同的权重W_i，N₂、O₂、NO、N、O等常规组分的权重因子设置为1，离子和电子等稀有组分的权重因子小于1。

如Boyd^[16]指出，粒子权重在两个方面影响DSMC的碰撞计算：碰撞对数目和二体碰撞机制。基于NTC(Non-Time-Counter)方法^[26]，文献[12]提出将电子作为一个组P，而将除电子以外的其它组分作为另一个组Q；不考虑P组内电子的相互碰撞，于是碰撞可以分为两类，即Q组内粒子之间的相互碰撞以及P组与Q组粒子之间的相互碰撞。同时，统计P组的电子数目N_P和Q组的粒子数目N_Q，于是网格i中Q组内的粒子碰撞对数为

$ \frac{1}{2}{N_{\rm{Q}}}{{\bar N}_{\rm{Q}}}{F_N}{({\sigma _T}{c_r})_{{\rm{QQ, max}}}}\Delta t/{V_c} $

(4)

P组与Q组粒子之间的分子碰撞对数为

$ \frac{1}{2}{N_{\rm{P}}}\bar N_Q{F_N}{({\sigma _T}{c_r})_{{\rm{PQ, max}}}}\Delta t/{V_c} $

(5)

权重因子的引入使得式(4)和式(5)中的N_P、N_Q和N_Q将变大。传统的做法是用一个数组CS(2, i, j)来记录碰撞网格i中的分子数目，这一数组在排序模块中更新，其做法是遍历仿真粒子并记录在网格i中组分j的仿真粒子数目。需要注意的是，由于网格i中不同组分的仿真粒子具有不同的权重，式(4)中的N_Q须包含组分权重，即

$ {{N'}_{\rm{Q}}} = \sum\limits_j {CS\left( {2, i, j} \right) \times {W_j}} $

(6)

W_j为组分j粒子对应的权重。式(5)中的N_P亦须包含电子的权重，即

$ {{N'}_{\rm{P}}} = CS(2, i, e) \times {W_{\rm{e}}} $

(7)

W_e为电子e的组分权重。用式(6, 7)中的N′_Q、N′_P及对应N′_Q代替式(4、5)中的N_Q、N_P和N_Q即可得到真实的分子碰撞对数目。

粒子权重对DSMC碰撞计算的第二个方面是二体碰撞机制。当两个粒子权重相同时，碰后的粒子速度为

$ {{\mathit{\boldsymbol{U'}}}_1} = {\mathit{\boldsymbol{U}}_m} + \frac{{{m_2}}}{{{m_1} + {m_2}}}\mathit{\boldsymbol{g}} $

(8a)

$ {{\mathit{\boldsymbol{U'}}}_2} = {\mathit{\boldsymbol{U}}_m}{\rm{ - }}\frac{{{m_1}}}{{{m_1} + {m_2}}}\mathit{\boldsymbol{g}} $

(8b)

U_m为碰撞粒子对的质心速度，g为碰撞粒子对的相对速度。当两个粒子权重不同时，碰撞后较低权重粒子的速度由式(8a)或式(8b)确定，较高权重粒子则仅有一部分速度发生变化(由式(8a)或式(8b)确定)；较高权重粒子速度发生变化的概率为

$ P = \left\{ \begin{array}{l} {W_1}/{W_2}, {\rm{若}}{W_1} < {W_2}, \\ {W_2}/{W_1}, {\rm{其他。}} \end{array} \right. $

(9)

粒子的碰撞后的内能状态类似确定。

上述方法并不保证每一步的动量和能量守恒，但是在足够大量的碰撞时可以保持动量和能量近似平均守恒。

2.2 权重因子对化学反应的影响

对于含电离的11组元空气模型，涉及稀有组分的化学过程有四种类型：常规组分与常规组分反应生成稀有组分，如联合电离反应N+O→NO⁺+e；常规组分与稀有组分反应生成稀有组分，如直接电离反应N+e→N⁺+e+e；稀有组分与稀有组分生成常规组分，如联合电离逆反应NO⁺+e→N+O；常规组分和稀有组分反应生成常规组分和稀有组分，如电荷交换反应N₂+N⁺→N+N₂⁺。为处理方便，假定所有稀有组分取相同权重因子W_r。

对于联合电离反应N+O→NO⁺+e，将生成物NO⁺和电子e复制1/W_r份。对于直接电离反应N+e→N⁺+e+e，本质上的反应物只有原子N，生成物为N⁺和电子e，其处理方式是对生成的N⁺和电子e复制1/W_r份，参与反应的电子e只改变其速度。对于联合电离逆反应NO⁺+e→N+O，依概率W_r保留生成物N和O。对于电荷交换反应N₂+N⁺→N+N₂⁺，其处理方式是反应物中的常规组分N₂依概率W_r消除，生成物中的常规组分N依概率W_r保留。

采用上述权重因子法的根本目的是提高稀有组分的仿真分子数目，进而抑制其统计涨落，得到稀有组分较为光滑的分布等值线。依据上述处理方式，稀有组分的仿真分子数目扩大到1/W_r倍。

3 典型再入的稀有组分权重因子效果验证 3.1 RAM-C Ⅱ飞行试验稀薄电离过程计算验证

RAM-C Ⅱ飞行试验是目前可查到的仅有飞行试验电子数密度测量公开案例，常作为检验电离计算模型和程序的经典算例。文献[12]考察了不使用稀有组分权重因子方法的电离过程模拟，结果显示尽管在测点方向的最大电子数密度与试验测量值吻合较好，球锥体周围流场的电子数密度分布存在很大统计波动。采用与文献[12]完全相同的计算条件，取稀有组分权重因子W_r为0.01。图 1给出了是否采用稀有组分权重因子方法得到的宏观流场参数比较情况。很明显，使用稀有组分权重因子方法与否，得到的流场宏观参数几乎不存在差异，意味着电离导致的稀有组分对宏观流场参数几乎不存在影响，也从一个方面印证了第一宇宙速度下常规气动特性计算通常不考虑电离过程在工程意义上的合理性。

然而，若考察电子数密度分布等再入物理领域极为关心的核心问题，权重格式使用与否，对于能否得到光滑的流场电子数密度分布，有着极为重要的影响。图 2给出了使用稀有组分权重因子方法与否得到的电子数密度分布云图。不使用权重因子方法时，在电子密度较高的头部附近，电子数密度等值线较为光滑，但是电子数密度较低的身部和尾部区域，巨大的统计涨落直接淹没了电子数密度等值线。使用权重因子方法时，上述情形得到了极大改善，电子数密度较高和较低的区域均能得到光滑的等值线。

图 3给出了使用稀有组分权重因子方法与否得到的沿测点方向最大电子数密度比较情况。该图表明，在不使用稀有组分权重方法时，尽管总体上DSMC计算值与试验测量值较好吻合，在中部两个测点仍存在一些偏差；使用稀有组分权重方法以后，上述偏离能得以较好矫正，计算值与试验测量值明显更为靠近。甚至可以认为，不使用稀有组分权重因子时的偏差是因为统计涨落造成的，事实上图 2上半部分也确实存在极其明显的统计涨落，这一涨落在使用权重因子方法时可以得到有效抑制，进而得到与试验测量值更为靠近的计算结果。

3.2 彗星探测器Stardust再入稀薄段电离过程计算验证

彗星探测器Stardust^[27]实现了迄今人造飞行器的最高速度, 它以12.8 km/s的速度再入大气层，必将引起强烈的化学反应和电离过程。采用与文献[12]完全相同的计算条件，取稀有组分权重因子W_r为0.1。图 4给出是否采用稀有组分权重因子方法得到的宏观流场参数比较情况，与RAM-C Ⅱ算例的情况一样，稀有组分权重因子的使用不会给宏观流场参数分布带来影响。

图 5给出了核心关注的Stardust在80 km高度的电子数密度分布情况。该图表明，对于12.8 km/s量级的再入速度，不使用权重因子方法时，DSMC方法计算得到的电子数密度分布等值线在流场大部分区域已经较为光滑，但在尾流的低密度区仍存在明显的统计波动。权重因子方法的使用，使得上述波动被进一步抑制，得到的等值线光滑性更好。事实上，相对于文献[12]，使用权重因子方法得到的电子数密度与文献[27]相比，数据吻合性更好，如图 6所示。如引言所述，电子等稀有组分数密度的统计波动可能导致不真实的化学反应计算，其原因是统计波动给网格中碰撞分子对计算带来一定程度的不确定性，这也是图 5中壁面附近电子数密度显示出一定差异的可能原因。

4 结论

在前期工作^[12]的基础上，提出了一种处理含电离化学反应DSMC模拟的稀有组分权重因子方法。该方法对DSMC二体碰撞的影响包括碰撞分子对数的确定和碰撞机制两个方面，前者的修正需考虑权重因子，后者则根据权重因子之比确定的概率进行状态确定。对于空气11组元的含电离化学反应，归类为四种情况分别处理，基本思想是根据权重因子对生成的稀有组分粒子进行复制、对反应物中的常规组分按概率保留或删除。数值测试表明，在诸如RAM-C Ⅱ飞行条件的极弱电离情况下，权重因子方法的使用使得电子数密度等值线光滑性大幅改善；在Stardust飞行条件的强电离情况下，权重因子方法的使用也可以改善低电子数密度区域等值线的光滑性，得到的电子数密度分布与参考文献值吻合较好。同时，权重因子方法的引入不会给宏观流场参数计算带来影响。