摘要：根据国外相关文献整理得出了包含83 种组分 550 步基元化学反应的肼类燃料/红烟硝酸(RFNA)详细化学反应动力学模型，并针对常规隐式和显式多步数值积分方法求解刚性化学反应动力学方程出现的质量分数为负值问题，提出了一种新型变步长隐式差分迭代算法；然后针对发动机开机启动过程的低压环境，采用该化学反应动力学模型及算法对自燃推进剂化学着火延迟进行了数值模拟，得出了 RFNA 组分含量、O/F 混合比、环境压力及温度对着火延迟时间的影响。结果表明：当RFNA 中的二氧化氮含量较低时，适当增加二氧化氮含量或提高O/F 混合比都可缩短化学着火延迟时间，但前者的效果更为明显；当环境压力和温度较低时，提高温度或压力能明显缩短化学着火延迟时间，但随着温度和压力的增大，这两者对化学着火延迟时间的影响会逐渐削弱。
关键词：自燃推进剂；化学反应动力学；差分算法；化学着火延迟
Numerical Simulation of Chemical Delay Time  of Liquid Rocket Hypergolic Propellant
Abstract: A detail hydrazine derivatives/red fuming nitric  acid (RFNA) chemical reaction kinetic model consisting of 550 elementary reactions among 83 species was constituted based on several foreign relevant literatures, and  a newly implicit interactive difference algorithm with variable time steps was proposed to tackle the problem of mass fraction non-positive in the integration procedure of stiff chemical kinetic equations using common implicit and explicit integral algorithm. The chemical delay time of hypergolic propellant was simulated using the chemical kinetic model and the algorithm on the condition of rocket engine low pressure start, and the impact of RFNA component content, O/F mix ratio, temperature and pressure on the chemical delay tim e was analyzed. The results show that increasing nitrogen dioxide and O/F mix ratio both could reduce the chemical delay time  at low concentration of nitrogen dioxide, but the former factor is much more  distinct. Increasing temperature and pressure could also improve chemical delay time at low temperature  and low pressure environment, but the effect would vanish gradually with much more increase of temperature and pressure.
Key words: hypergolic propellant; chemical reaction kinetic; difference algorithm; chemical delay time
　　引言
　　液体火箭自燃推进剂具有常温常压下正常贮存，不需                                                                  要任何点火源自燃着火燃烧的优点，被广泛应用于运载火箭发动机及空间姿轨控发动机，还被应用于新型凝胶推进剂的开发研制[1]。为了提高发动机开机启动可靠性及推进剂燃烧效率，缩短自燃推进剂着火延迟时间通常是发动机设计及推进剂研制所考虑的重要因素。
　　自燃推进剂在液体火箭发动机的开机启动过程中主要经历了推进剂雾化、喷雾液滴蒸发、气相混合及自燃着火燃烧四个阶段[2]，其中前三个阶段属于点火延迟中的物理延迟，而化学点火延迟指混合后的气相组分发生化学反应达到着火温度的时间。国内外对物理着火延迟进行过大量研究[3-5]，但针对自燃推进剂气相化学着火延迟方面的数值研究非常有限[6]，一是由于肼类自燃推进剂存在复杂化学反应，使得国内外现阶段还没有可完整描述肼类燃料/RFNA 化学反应过程的详细化学反应动力学模型；二是由于现有的常规隐式和显式数值积分方法[7,8](如LSODE软件包和VODE 软件包所采用的Gear 、Adams 方法)都不能确保多步复杂刚性化学反应动力学的质量分数积分结果为正值，导致化学动力学方程出现错误的数值积分结果。
　　针对自燃推进剂化学着火延迟研究存在的以上两方面问题，本文首先根据国外相关文献整理[9~11]得出了包含83种组分550 步基元化学反应的肼类燃料/ 红烟硝酸(RFNA)化学反应动力学模型，然后提出了一种处理多步复杂刚性化学动力学方程的新型变步长隐式差分迭代算法，并通过理论和仿真数据验证了模型及算法合理性；最后针对发动机开机启动过程的低压环境，通过求解等容绝热条件下的气相化学反应动力学方程，仿真分析了RFNA 组分含量、O/F 混合比、环境温度、环境压力对甲基肼(MMH)/ 红烟硝酸(RFNA)化学着火延迟时间的影响。
　　1  自燃推进剂化学反应动力学模型
　　肼类燃料与RFNA 的化学反应动力学包括肼类燃料及氧化剂的分解反应、肼类燃料与其他分子间的双分子反应、小分子及离子间的复合反应等。文献[9]给出了包含204 步基元反应的MMH高温分解反应动力学模型，文献[10] 给出了包含403 步基元化学反应的 MMH/NTO 反应动力学模型，文献[11] 的GRI_MECH3.1化学反应动力学模型包含了自燃推进剂小分子分解产物之间的基元化学反应。通过对以上文献中的化学反应动力学模型进行整理，得出了包含83种组分550 步基元化学反应的肼类燃料/ 红烟硝酸化学反应动力学模型，其中既包含可逆反应和不可逆反应，也包含三体反应及压力关联反应，因此适用于低压及高压环境下的化学反应动力学仿真。 化学反应动力学模型中的基元反应通式可表示为： 11 1, 2, ,ss NNij i ij i riiXM XM j N              (1) 其中： s N 、 r N 分别为化学反应动力学模型的组分数量和基元反应数量；i X 代表每个基元反应中的反应组分；M 代表三体反应以及压力关联反应中的活性组分；ij   、ij   分别为各基元反应中的反应组分的化学计量数。 根据反应动力学通式可得出各组分的净摩尔生成率： 11 1 ()ij ijrs s NN Nii FR iijijj jii ji iYY KKMW MW                      (2) 其中：FkK 、RkK 分别为前向和逆向化学反应常数，可通过Arrhenius定律及平衡常数进行计算； 为气相混合物的密度(3kg m )；i Y 为各组分质量分数； i M W 为各组分摩尔质量。 根据式(2)可得出各组分质量浓度变化率： iiidYMWdt                                   (3) 根据热力学第二定律，化学反应体系在等容绝热条件下的内能保持不变，即 01()consts Niiinh RT U                       (4) 其中：in 为组分摩尔质量；R 为理想气体常数；T 为反应体系温度；hi 为各组分在温度为T 时的理想气体焓值；0 U为反应体系初始总内能。 对式(4)求微分，可得出化学反应体系温度对时间的导数： 1,1()()ssNiiiNipiiihRTdTY dtCRMW                     (5) 其中： , p i C 为各组分定压摩尔比热容(-1 -1Jkmol K )。 式(3)和式(5)联立组成了化学反应动力学模型的等容绝热化学反应动力学方程组，通过求解该方程组可以得出化学反应体系温度及各组分浓度随时间的变化，并认为温度梯度最大时刻为自燃推进剂的化学着火延迟时间。
　　2  变步长隐式差分迭代算法
　　LSODE软件包和VODE 软件包中的Gear 线性多步隐式积分方法通常用于刚性常微分方程初值问题的求解[12]，但该方法没有对质量分数进行边界条件限制，因此不能保证多步复杂刚性化学动力学方程的积分结果始终保持为正值。文献[7]在质量守恒及质量分数为正的约束条件下，采用数学规划方法对常规积分方法所得的负值进行修正；虽然该修正算法确保了质量分数为正，但该方法所采用的数学规划算法比较复杂，且计算效率也会受修正算法的影响。文献[8]提出了一种分离单步积分算法(Split Single Reaction Integrator)，该算法的基本思想是首先得出每个基元反应的理论积分结果，然后采用分离算法思想分步得出完整化学动力学模型的数值结果，但该算法需要首先确定各基元反应理论解，并根据各基元反应的化学反应速率进行排序，这限制了计算程序在不同化学反应动力学模型间的通用性。 本文通过把化学反应净摩尔生成率分解为生成率和消耗率两部分，并采用变步长隐式差分迭代算法求解差分方程，在确保质量分数始终为正值的基础上提高了大时间步长的计算稳定性。
　　首先，净摩尔生成率i  分解为生成率Ci  和消耗率Di 两部分之差：CD ii i                                   (6) 其中： 11 1ij ijrs s NN N ii CF R iijj ijj ii ji iYY KKMW MW                   (7) 11 1ij ijrs s NN N ii DF R iijj ijj ii ji iYY KKMW MW                     (8) 消耗率Di  移项到式(3)左边，进行整理可得： Dii Ciiii idYMW Y MWdt Y                   (9) 然后，对式(5)及式(9)的时间导数进行一阶隐式差分，整理可得组分及温度的差分方程： 10 1()()() nDinCn iiiii itMW Y tMW YY      (10) 1 10,1{[( ) ( ) ] ( )}()( )ssNCn Dn n niii i nN nnipiiithRT TT YC RMW        (11) 其中：t  表示差分计算步长；上标 0、n 、 1 n  分别代表上一时刻，第n 迭代步及第 1 n  迭代步。
　　根据式(7)和式(8)对生成率和消耗率的定义可知，组分差分方程式(10) 左右两边的系数在各时刻及各迭代步都始终保持为正值，因此能保证各时刻的质量分数为正值。另一方面，该方法通过分离和移项处理，把质量生成率中的质量分数进行了隐式处理，并在计算过程中对每步迭代结果进行欠松弛处理，因此有效地提高了大时间步长的计算稳定性。
　　3  化学动力学模型及算法验证
　　在初始压力00.242bar P  ，初始温度0298K T  ，各反应物初始质量分数MMH0.221 Y  、NTO0.454 Y  、2 NO0.325 Y 条件下，图1 给出了采用所整理的肼类燃料/RFNA 多步化学反应动力学模型及变步长隐式迭代差分算法对其在绝热等容环境下的化学着火延迟仿真结果。 由图1(a) 给出的仿真步长和每步迭代次数可知，该变步长隐式迭代算法在迭代次数超过某限定值时自动对松弛因子及仿真步长进行调整，从而成功地为着火诱导区、热释放区及化学平衡区选取了合理的仿真步长，提高了计算稳定性和计算效率。 由图1(b)~(d)给出的反应物、主要生成产物及反应温度在热释放区附近随时间的变化可知，MMH/NTO 反应体系的反应温度及各组分质量浓度在热释放区内的变化非常剧烈，最终的平衡温度为 3290.47K ，主要平衡生成产物为2O 、2N 、OH、2HO、CO、2CO ，自燃着火延迟时间为3.65ms ，这与文献[10] 中的着火延迟时间非常吻合。 0 1 2 3 4 510-410-310-210-1Time (ms)   t   ( ms)0 1 2 3 4 50100200300400500Interative NumberIgnition Induction DomainEquilibrium DomainHeat Realse   DomianIgnition Induction DomainHeat ReleaseDomainEquilibriumDomain0 1 2 3 4 510-410-310-210-1Time (ms)   t   ( ms)0 1 2 3 4 50100200300400500Interative NumberIgnition Induction DomainEquilibrium DomainHeat Realse   DomianIgnition Induction DomainHeat ReleaseDomainEquilibriumDomain                    3.5 3.6 3.7 3.8 3.9 400.10.20.30.40.50.60.7Time(ms)Mass Fraction3.3 3.4 3.5 3.6 3.7 3.8 3.9 40500100015002000250030003500Tem p erature( K ) MMH NTO NO2Temperature Temperature3.5 3.6 3.7 3.8 3.9 400.10.20.30.40.50.60.7Time(ms)Mass Fraction3.3 3.4 3.5 3.6 3.7 3.8 3.9 40500100015002000250030003500Tem p erature( K ) MMH NTO NO2Temperature Temperature (a)  仿真步长及迭代次数                                            (b)  反应物质量分数及温度 3.5 3.6 3.7 3.8 3.9 400.10.20.30.40.5Time(ms)Mass Fraction3.3 3.4 3.5 3.6 3.7 3.8 3.9 40500100015002000250030003500Tem p erature( K ) O2 N2 OHTemperatureTemperature3.5 3.6 3.7 3.8 3.9 400.10.20.30.40.5Time(ms)Mass Fraction3.3 3.4 3.5 3.6 3.7 3.8 3.9 40500100015002000250030003500Tem p erature( K ) O2 N2 OHTemperatureTemperature                     3.5 3.6 3.7 3.8 3.9 400.10.20.3Time(ms)Mass Fraction3.3 3.4 3.5 3.6 3.7 3.8 3.9 40500100015002000250030003500Tem p erature( K ) H2O CO2TemperatureTemperature3.5 3.6 3.7 3.8 3.9 400.10.20.3Time(ms)Mass Fraction3.3 3.4 3.5 3.6 3.7 3.8 3.9 40500100015002000250030003500Tem p erature( K ) H2O CO2TemperatureTemperature (c)  主要产物质量分数及温度                                         (d)  主要产物质量分数及温度 图1  MMH/NTO 自燃着火延迟仿真结果 5004003002001000 0       1        2       3       4      5 Time/ms 10-1∆ t/ms  10-210-310-4Interative Number 3500300025002000150010005003.3   3.4     3.5    3.6     3.7   3.8   3.9     4 Time/ms Mass Fraction  0.50.40.30.20.1Temperature/K O2 N2 OH 3500300025002000150010005003.3   3.4     3.5   3.6     3.7   3.8   3.9     4 Time/ms Mass Fraction  0.30.20.1Temperature/K H2 OCO23500300025002000150010005003.3   3.4     3.5   3.6     3.7   3.8   3.9     4 Time/ms Mass Fraction  0.70.60.50.40.30.20.1Temperature/K 第25卷第4 期  Vol. 25 No. 4 2013年4 月  何博,  等:  液体火箭自燃推进剂化学着火延迟数值模拟 Apr., 2013  http: ∥www.china-simulation.com• 615 • 表1 在相同初始条件下，分别给出了图 1 中采用化学动力学模型求得的平衡温度和主要平衡产物质量分数，及采用Cantera软件包中的化学平衡程序[13]所得理论结果。由表可知，本文仿真结果与理论结果的平衡温度非常吻合，大多数主要平衡产物的质量分数也比较吻合，仅平衡产物氧气的误差达到19.7% ；这可能是因为化学动力学模型中包含大量组分及基元反应，造成计算过程中出现了舍入误差逐渐累积。 表1  平衡温度及主要平衡产物质量分数   T E Q / K O2  N2 OH H2 O CO CO2本文结果 3290.47 0.12 0.36 0.059 0.200 0.083 0.113Cantera 3293.78 0.15 0.35 0.062 0.192 0.068 0.1054  仿真结果及分析 图2 以液体火箭发动机开机启动过程的低压环境00.242bar P  、0298K T  为例，给出了 MMH/RFNA化学着火延迟时间随O/F混合比及红烟硝酸(RFNA)中二氧化氮(NO2)质量含量的变化。由图可知，当混合比和NO2 含量较低时，NO2 含量对化学着火延迟时间的影响非常明显，此时的NO2 含量几乎与着火延迟时间的对数成线性关系，即略微增加NO2 含量就可大幅降低化学着火延迟时间；但随着混合比及NO2 含量进一步增大，着火延迟时间下降到毫 0.15 0.17 0.19 0.21 0.23 0.25100101102103104Mass Fraction (w%)Autoignition Delay Time  (ms)   O/F =4   O/F =5   O/F =6   O/F =7P  =0.242bar   T  =298K 图2  着火延迟与混合比和NO2 含量的关系 290 300 310 320 330 340050100150200250Temperature  (K)Autoignition Delay Time (ms)   P  =0.242bar   P =0.4bar   P =0.6bar   P =0.8barw%=17.5%    O/F=5 图3  着火延迟与压力和温度的关系 秒量级，并逐渐趋于稳定。 图3 在NO2 质量含量w% 0.175  、混合比O/F 5  条件下，给出了发动机开机启动阶段化学着火延迟时间随燃烧室压力及温度的变化。由图可知，当环境压力较低时，提高环境温度可缩短着火延迟时间；低温环境下提高环境压力也能缩短着火延迟时间，但随着环境压力和温度的进一步升高，这两者对着火延迟时间的影响将逐渐消弱。
　　5  结论
　　本文根据国外文献整理得出了包含83 种组分550 步基元化学反应的肼类燃料/RFNA 化学反应动力学模型，并提出了一种能确保质量分数始终为正值的变步长隐式差分迭代算法。然后以火箭发动机开机启动过程的低压环境为例，仿真分析了RFNA 的NO2 含量、O/F 混合比、环境温度和压力对自燃推进剂化学着火延迟时间的影响。所得结果可为缩短自燃推进剂化学着火延迟时间，提高发动机开机启动可靠性提供以下参考。 (1)  在保证环境安全及不影响比冲性能的条件下，适当提高推进剂中NO2 含量可明显缩短化学着火延迟。 (2)  推进剂可在开机启动阶段提前进行预热，一方面可提高推进剂蒸汽温度，缩短化学着火延迟时间；另一方面还可提高蒸发速率，改善燃烧效率。 (3)  氧化剂可在开机启动阶段提前进入燃烧室，其蒸汽可提高燃烧室压力，并增大O/F 混合，达到缩短开机启动化学着火延迟时间的目的。
　　参考文献: [1]  丰松江,  何博,  聂万胜.  液体火箭凝胶推进剂燃烧特性研究进展[J].  火箭推进, 2009, 35(4):1-7. [2]  庄逢辰.  液体火箭发动机喷雾燃烧的理论、模型和应用[M]. 长沙:  国防科技大学出版社, 1995. [3]  Kushari A. Study of an Internally Mixed Liquid Injector for Active Control of Atomization Process [D]. USA: Georgia Institute of Technology, 2000. [4]  Pope D N. Numerical Simulation of Convective Fuel Droplet Vaporization and Combustion in a Low Pressure Zero-Gravity Environment [D]. Lincoln, USA: The University of Nebraska, 2001. [5]  Schmehl R, Rosskamp H, Willmann M. CFD Analysis of Spray Propagation and Evaporation Including Wall Film Formation and Spray/Film Interactions [J]. International Journal of Heat and Fluid Flow (S0142-727X), 1999, 20(5): 520-529. [6]  Dasarathy R. Kinetic Studies Using Chemical Delay Times of Anhydrous Hydrazine and Other Hypergolic Fuels [D]. Huntsville, USA: The University of Alabama, 2006. [7]  Sandu A. Positive Numerical Integration Methods for Chemical Kinetic Systems [J]. Journal of Computational Physics   (S0021-9991), 2001, 170(2): 589-602. 0.15    0.17    0.19    0.21    0.23    0.25Mass Fraction/ w% Autoignition Delay Time  /ms 104103102101100290     300     310     320     330    340Temperature/K  Autoignition Delay Time  /ms 250200150100500w%=17.5%  O/F=5  P=0.242barP=0.4bar P=0.6bar P=0.8bar P=0.242bar    T =298K O/P=4O/P=5O/P=6O/P=7 (a)  [8]  Nguyen K, Caboussat A, Dabdub D. Mass conservative, positive definite integrator for atmospheric chemical dynamics [J]. Atmospheric Environment (S1352-2310), 2009, 43(40): 6287-6295. [9]  Sun H Y, Catoire L, Law C K. Thermal Decomposition of Monomethyldrazine: Shock Tube Experiments and Kinetic Modeling [J]. International Journal of Chemical Kinetics (S0538-8066), 2009, 41(3): 176-186. [10]   Catoire L, Chaumeix N, Paillard C. Chemical Kinetic Model for Monomethylhydrazine/Nitrogen Tetroxide Gas-Phase Combustion and Hypergolic Ignition [J]. Journal of Propulsion and Power (S0748-4658), 2004, 20(1): 87-92. [11]   Smith G P, Golden D M, Frenklach M, Moriarty N W, Eiteneer B. GRI-Mech 3.0. [EB/OL]. (2000-2) [2010-10].   http://www.me.berkeley.edu/gri_mech. [12]   Radhakrishnan K, Hindmarsh A C, Description and Use of LSODE, the Livermore Solver for Ordinary Differential Equations [R]// LLNL report UCRL-ID-113855, 1993.   [13]   Goodwin D G. Installing Cantera 1.7 for Windows [EB/OL]. (2007-3) [2010-12].