固体火箭尾焰等离子体特性影响因素数值仿真

孙 行， 聂万胜*， 蔡红华， 陈 朋， 黄卫东， 石天一

(1.航天工程大学宇航科学与技术系，北京 101416；2.航天员科研训练中心，北京 100094)

推进剂在火箭发动机中燃烧，经拉瓦尔喷管喷出，形成高温尾焰，尾焰组分发生电离，形成稠密不均匀的弱电离等离子体[1-2]，对火箭进行测控时发现，电磁波穿过火箭尾焰时会受到严重干扰[3-8]。研究火箭尾焰等离子体的各项特征参数(等离子体浓度、等离子体频率)对于研究火箭尾焰对电磁波的干扰作用意义重大。

近年来，中外对火箭尾焰等离子体特征参数的研究不断深入。Fromentin-Denoziere等[7]描述了尾焰中碱金属元素的电离过程；Smoot等[8]描述了尾焰中自由电子的分布规律，预测了雷达信号在尾焰中的衰减；Troyes等[9]采用数值计算的方法对固体火箭推进剂是否含铝两种情况的尾焰流场进行了研究；Chen[10]对等离子体特征参数的计算方法进行了描述。清华大学的杨任刚等[11]在采用RNGk-ε湍流模型描述气相的同时采用离散相模型描述固相，研究了尾焰气固两相射流流场；石雁祥等[5-6]在考虑火箭尾焰为尘埃等离子体的基础上，对火箭尾焰等离子体特征参数进行了大量的理论研究；郭祥天等[12]建立了固体火箭发动机尾焰的气固两相流计算模型，研究了固相颗粒Al2O3对发动机尾焰流场的影响。研究发现，相较于液体火箭，固体火箭尾焰对电磁波的干扰作用更加明显。这些研究多以火箭尾焰流场参数的经验值为基础，计算等离子体特征参数，本文以某型号固体火箭为研究对象，对其不同工况下的尾焰流场进行仿真，获取准确的尾焰流场参数；在此基础上，计算尾焰等离子体特征参数分布，分析其分布规律，并讨论不同工况对尾焰等离子体特征参数的影响。

1 物理模型和计算方法

1.1 控制方程

基于Navier-Stokes方程计算火箭尾焰流场，添加组分输运方程，得到如下四个控制方程[13]。

1.1.1 质量守恒方程

(1)

式(1)中：ρ为流体的密度；v为流体的速度矢量。

1.1.2 动量守恒方程

(2)

1.1.3 能量守恒方程

(3)

式(3)中：cp为定压比热容；T为流体的温度；k为流体的传热系数；ST为流体的内热源及由于黏性作用而使流体机械能转化为热能的部分。

1.1.4 组分质量守恒方程

(4)

式(4)中：cs、ρcs及Ds分别为组分s的体积浓度、质量浓度及扩散系数；Ss为系统内单位时间单位体积通过化学反应产生的组分s的质量。

1.2 化学反应

采用文献[14]中给出的12组分17步化学反应模型计算燃烧过程，以CHEMKIN嵌入计算流体力学(CFD)的流场计算中。参与化学反应的12组分为CO、CO2、H、H2、H2O、OH、O、O2、HCl、Cl、Cl2、N2,反应速率由式(5)计算：

(5)

式(5)中：κ为化学反应速率,mol/(cm·s·k)；A为指前因子，对应于二级化学反应，cm3/(mol·s)；T为温度，K；n为温度指数1；E为活化能，cal/mol；具体反应参数如表1所示。

表1 多步反应Table 1 Multi-step chemical reaction

采用1～4步分支反应和5～8步重组/离解反应计算H2/O2系统，9～11步反应计算CO/CO2系统，12～17步反应描述抑制剂HCl的作用[14]。

1.3 湍流模型与燃烧模型

尾焰流场流动属于湍流形态，采用realizable，k-ε双方程湍流模型计算可压缩、雷诺平均N-S方程[15-16]。由于采用多步化学反应机理，燃烧模型采用涡耗散概念模型。

1.4 等离子体浓度

根据理想气体状态方程，可得气体粒子数密度为

(6)

式(6)中:ne为电离粒子的数密度,cm-3，nm为未电离粒子的数密度，cm-3;P为压强，T为温度，均采用国际单位制。

文献[10]中根据萨哈-朗缪尔方程导出处于热平衡状态的气体电离量为

(7)

式(7)中：Ui为气体的电离能(最外层电子离开原子所需的能量，单位为J)，K=1.38×10-23J/K为玻尔兹曼常数。

(8)

1.5 等离子体频率

文献[10]给出了等离子体频率的计算公式：

(9)

式(9)中：ωpe为电子的等离子体角频率；ne为电子的等离子体密度；e为电子所带电量；ε0为真空介电常数；me为电子的质量。在实际计算应用中，常采用更简洁的表达方式[10]：

(10)

式(10)中：fp为等离子体角频率对应的频率，Hz。

1.6 热力学计算

研究的固体火箭发动机使用复合推进剂，发动机燃烧室压强为9.15 MPa，燃烧室平均质量流量333 kg/s，固相燃烧过程采用热力计算简化，计算结果如表2所示[13]。

表2 热力学计算结果Table 2 Thermodynamic calculation results

1.7 网格构建及边界条件

采用CFD对固体火箭尾焰进行仿真，以获取火箭尾焰流场参数，计算区域从喷管喉部开始，轴向长度40 m，径向长度15 m，为减小仿真计算量，采用1/2对称模型。喉部采取质量流动入口边界条件，喷管壁面采取固壁边界条件，对称面采取对称边界条件，远场采取压力远场边界条件，网格构建完成后进行无关性验证，共采用15×104、44×104和62×104三套网格计算相同工况下的尾焰流场，15×104网格对应的计算结果有一定误差，44×104和62×104两套网格对应的计算结果相近，因此采用44×104网格方案，如图1所示。

图1 计算网格Fig.1 Computing grid

计算火箭发动机在四种情况下的工作情况，四种工况如表3所示。

表3 发动机工况Table 3 Engine conditions

考虑Al2O3作为离散相参与换热，粒径取50 μm，并假设Al2O3颗粒为球体，温度相同，且颗粒之间无相互作用。

2 尾焰流场计算结果

2.1 Al2O3含量对尾焰流场的影响

图2、图3是发动机处于地面，Al2O3含量10%、20%两种情况下，尾焰的温度云图与压强云图。图4为Al2O3含量对尾焰温度与压强沿轴线变化情况的影响图。

图2 工况1下的尾焰温度及压强云图Fig.2 Contours of temperature and pressure under condition 1

图3 工况2下的尾焰温度及压强云图Fig.3 Contours of temperature and pressure under condition 2

图4 不同推进剂下尾焰温度与压强沿轴线变化情况Fig.4 Variation of temperature and pressure along the axis under different propellants

固体推进剂中加入铝粉可提高发动机比冲[17]，铝粉燃烧生成Al2O3，随着Al2O3的增加，尾焰温度明显升高，推进剂燃烧更充分，燃烧产物充分膨胀，喷出喷管后在外界大气压的作用下，形成锥形波与倒锥形波的反复出现，但Al2O3含量从10%增加到20%对尾焰压强的影响较小。

2.2 飞行高度对尾焰流场的影响

图5、图6是Al2O3含量10%，发动机处于不同飞行高度下，尾焰的温度云图与压强云图。图7为Al2O3含量10%时，不同飞行高度对尾焰温度与压强沿轴线变化情况的影响图。

图5 工况3下的尾焰温度及压强云图Fig.5 Contours of temperature and pressure under condition 3

图6 工况4下的尾焰温度及压强云图Fig.6 Contours of temperature and pressure under condition 4

图7 不同飞行高度下尾焰温度与压强沿轴线变化情况Fig.7 Variation of temperature and pressure along the axis under different heights

飞行高度从10 km增加到20 km，外界大气压降低，对过膨胀燃气的压缩作用降低，尾焰温度峰值明显后移，但其峰值几乎不变，这也说明了确实是Al2O3的含量对尾焰温度起主要影响作用；飞行高度对尾焰压强分布影响不大。

3 尾焰等离子体特征参数计算结果

尾焰等离子体浓度是表征尾焰等离子体属性的重要物理量；当电磁波频率小于尾焰等离子体频率时，电磁波将发生全反射而无法穿过尾焰[4，18-20]。计算以上等离子体特征参数，对于计算尾焰与电磁波的相互作用意义重大。

3.1 等离子体浓度

使用CFD后处理模块的自定义函数功能，获得尾焰等离子体浓度分布云图如图8所示。使用编程软件计算不同工况下等离子体浓度随轴线变化趋势如图9所示。二者吻合较好。

图8 四种工况下尾焰等离子体浓度云图Fig.8 Contours of plasma concentration under four conditions

图9 不同工况下尾焰等离子体浓度沿轴线变化情况Fig.9 Variation of plasma concentration along the axis under different working conditions

随着Al2O3含量的增加，尾焰温度升高，燃气获得更多的电离能，电离程度变高，等离子体浓度显著增高，高浓度区域范围增大，峰值位置后移；随着发动机工作高度增加，尾焰高温区域后移，等离子体浓度峰值位置后移。

3.2 等离子体频率

电磁波仅在频率高于等离子体频率时可以穿过等离子体[18-20]，等离子体频率对研究尾焰对电磁波的衰减有重要作用。使用CFD后处理模块的自定义函数功能，获得尾焰等离子体频率分布云图如图10所示。使用编程软件计算不同工况下等离子体频率随轴线变化趋势如图11所示。二者吻合较好。

图10 四种工况下尾焰等离子体频率云图Fig.10 Contours of plasma frequency under four conditions

图11 不同工况下尾焰等离子体频率沿轴线变化情况Fig.11 Variation of plasma frequency along the axis under different working conditions

根据简化后的等离子体频率公式[式(10)]可知，在此种简化条件下，等离子体频率为等离子体浓度的函数，二者正相关，因此趋势变化相似。

4 结论

研究了不同推进剂组成以及不同飞行高度下火箭尾焰的流场分布，计算了火箭尾焰等离子体浓度、等离子体频率随尾焰轴线方向的变化情况，给出了火箭尾焰等离子体特征参数在尾焰对称面上的分布情况，得出以下结论。

(1)在喷管扩张段，热能向动能转化，燃气速度升高的同时，温度下降，等离子体浓度降低；燃气喷出喷管后，在激波位置，温度迅速升高，等离子体浓度也随之增大。

(2)尾焰等离子体特征参数高度依赖尾焰温度，推进剂中铝粉含量增加，尾焰温度显著升高，等离子体特征参数(等离子体浓度、等离子体频率)峰值升高；火箭飞行高度增加，尾焰高温区域明显后移，等离子体特征参数(等离子体浓度、等离子体频率)峰值位置后移。

以往相关领域的研究，多以火箭尾焰流场经验值为基础计算尾焰等离子体参数，在此基础上针对尾焰等离子体对电磁波的折射以及衰减的研究的准确性难以得到保证。得到了尾焰等离子体特征参数在尾焰对称面上的分布情况，对后续计算尾焰对电磁波的折射及衰减提供了重要参考。