燃烧室压强估算在固体火箭发动机的设计中是非常重要的，压强预估不准确很可能导致实验失败^[1]。零维燃烧室平衡压强公式目前最常用于压强工程计算，在实际使用中发现，随着推进剂中金属含量的提高，发动机工作效率逐渐下降^[2]，其中速度效率的下降证明了压强理论预估结果的准确度在降低^[3]。而金属燃料由于密度高、反应放热量大，在固体推进剂的发展过程中占有重要地位^[4]。近年来随着低温固体推进剂火箭发动机、金属燃料冲压发动机的发展^[5-6]，固体推进剂中金属含量逐渐增大，这些发动机推进剂中金属含量大于15%，因此被称为高金属含量发动机。金属含量增大导致压强预估准确度(η_C^*)下降。

通常η_C^*由推进剂燃烧效率、燃烧室及喷管喉部前热损失和两相流共同决定，由从实验获得的经验公式计算^[7]。经验公式有特定适用范围，在固体火箭发动机发动机中，当金属含量从15%变化为30%时，η_C^*由98.2%下降到94%^[2]；在金属镁水冲压发动机中，推进剂金属含量为50%，水燃比为0.3~3时，η_C^*在75%~85%变化^[6]；在铝冰低温固体火箭发动机中，推进剂中铝含量在40%~ 50%变化时，η_C^*在40%~60%^[8]。研究表明，含铝17%的复合推进剂的燃烧效率大于96%^[9]，而镁水冲压发动机中金属燃烧效率约为94%^[10]，而铝冰推进剂的实测铝燃烧效率大于85%^[11]，除燃烧效率外还有两相流、热损失等在影响η_C^*，在燃烧室压强理论计算中应予以考虑。

数值计算能得到包括燃烧室压强在内的多种流场参数，还能考虑材料烧蚀、推进剂侵蚀燃烧等情况对压强的影响^[12-13]，但是涉及到由两相流、热损失等因素导致η_C^*变化进而影响压强计算结果准确度的研究却很少^[14-15]。有学者在两相流数值计算的基础上，用数据拟合的方法得到喷管两相流量公式^[16]，并把其用于零维两相内弹道的计算，得到了与实验吻合的结果^[17-18]，但是其喷管流量公式在颗粒含量及粒径选取、燃气属性等方面没有明确规律，因此其结果仅适用于文献中的发动机，用于铝冰等高金属含量发动机时结果不准确。本文对文献[16]中计算方法进行完善和改进，分别对铝含量为16%~18%的复合推进剂发动机、铝含量为40%~50%的铝冰发动机进行两相流数值计算、数据拟合和燃烧室压强公式的推导和修正(总称为拟合公式法)，得到燃烧室压强的修正计算式，并对其计算结果进行验证。

1 两相流数值计算模型

当推进剂燃烧产物为两相流时，数值计算方法包括欧拉-欧拉方法和欧拉-拉格朗日方法^[19-20]。其中欧拉-拉格朗日方法的控制方程如下所示^[20]：

$ 燃气:\left\{ \begin{array}{l} \nabla \cdot \left( {{\rho _g}{\mathit{\boldsymbol{u}}_g}} \right) = 0\\ \nabla \cdot \left( {{\rho _g}{\mathit{\boldsymbol{u}}_g}{\mathit{\boldsymbol{u}}_g}} \right) = - \nabla P + \nabla \cdot \tilde \tau + {X_g}\\ \begin{array}{*{20}{c}} {\nabla \cdot \left( {{\rho _g}{\mathit{\boldsymbol{u}}_g}{C_{p,g}}{T_g}} \right) = \nabla \cdot \left( {\nabla {T_g}} \right) + }\\ {{\tau _{rT}} \cdot \left( {{T_p} - {T_g}} \right) + {{\bar S}_g}} \end{array} \end{array} \right. $

(1)

$ 颗粒:\left\{ \begin{array}{l} \int\limits_A {{n_p}{\mathit{\boldsymbol{u}}_p} \cdot {\rm{d}}A} = 常数\\ \frac{{{\rm{d}}{\mathit{\boldsymbol{u}}_p}}}{{{\rm{d}}t}} = - \frac{{{\mathit{\boldsymbol{u}}_p} - {\mathit{\boldsymbol{u}}_g}}}{{{\tau _{rv}}}} + {X_p}\\ \frac{{{\rm{d}}{T_p}}}{{{\rm{d}}t}} = - {\tau _{rT}} \cdot \left( {{T_p} - {T_g}} \right) + {{\bar S}_p} \end{array} \right. $

(2)

式中：u_g、u_p、T_g、T_p分别为气相和颗粒的速度和温度，P为气相的压强，$ \boldsymbol{\widetilde \tau} $为粘性应力张量，τ_rv为两相阻力因子，τ_rv为两相传热因子，其计算式参见文献[20]。X_g和X_p分别为燃气或颗粒受力，包括传质量力、Magnus力等，S_g和S_p分别为燃气或颗粒的热源相，包括相变热、辐射传热等^[19]。方程中计算所需的燃气属性参数、颗粒百分含量等通过推进剂热力计算得到^[21]。

控制方程考虑了燃气粘性对流动的影响，实际上该方程是湍流时均方程，而对于湍流本身，本文选用涡粘性模型计算，该模型中湍流的动力粘性系数μ_t则用k-ε双方程模型计算(其中ε为湍流耗散率，k为湍流脉动动能)，其控制方程如下^[22]：

$ \left\{ \begin{array}{l} \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho k\mathit{\boldsymbol{u}}} \right) = \nabla \cdot \left[ {\left( {\mu + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _k}}}} \right)\nabla k} \right] + {G_{\rm{k}}} - \rho \varepsilon \\ \begin{array}{*{20}{c}} {\frac{{\partial \left( {\rho \varepsilon } \right)}}{{\partial t}} + \nabla \cdot \left( {\rho \varepsilon \mathit{\boldsymbol{u}}} \right) = \nabla \cdot \left[ {\left( {\mu + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _\text{k}}}}} \right)\nabla \varepsilon } \right] + }\\ {{C_{1\varepsilon }}\frac{\varepsilon }{k}{G_k} - {C_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k}} \end{array} \end{array} \right. $

(3)

式中：$ {{G}_{k}}={{\mu }_{t}}\left( \frac{\partial {{u}_{\rm{i}}}}{\partial {{x}_{\rm{j}}}}+\frac{\partial {{u}_{\rm{j}}}}{\partial {{x}_{\rm{i}}}} \right)\frac{\partial {{u}_{\rm{i}}}}{\partial {{x}_{\rm{j}}}} $, $ {{\mu }_{t}}=\rho {{C}_{\mu }}\frac{{{k}^{2}}}{\varepsilon } $, C_μ、C_1ε和C_2ε为经验常数。

在湍流的影响下，欧拉-拉格朗日方法的计算模型为随机轨道模型。在方程(2)的基础上，颗粒动量方程变为：

$ \frac{{{\rm{d}}{\mathit{\boldsymbol{u}}_p}}}{{{\rm{d}}t}} = \frac{{\left( {{{\mathit{\boldsymbol{\bar u}}}_g} + {{\mathit{\boldsymbol{u'}}}_g}} \right)}}{{{\tau _{rv}}}} + {\mathit{\boldsymbol{F}}_{i,p}} $

(4)

式中u_g和u′_g为燃气的时均速度和脉动速度。

计算的物理模型是喷管二维轴对称模型，以18%铝含量的复合推进剂发动机为例^[23]，喷管流场网格及边界条件图 1所示。

在以上物理模型的基础上，用Fluent软件求解两相流控制方程，用标准k-ε模型求解气相的湍流流动，用颗粒随机轨道模型求解颗粒受湍流影响的参数，最终可以得到燃气总压沿喷管的变化。一般喷管入口总压容易测量，因此取图 1监测线上燃气总压为燃烧室压强，上述计算方法为燃烧室压强的数值方法。

2 燃烧室压强拟合公式法

喷管两相流动对于气相不等熵，因此喷管各截面总压是变化的，由等熵流动理论得到的喷管流量公式：$ {{\dot{m}}_{\rm{t}}}={{P}_{c}}{{A}_{t}}/{{C}^{*}} $不再适用，此时流量和压强的关系应该用上一节所述的数值计算方法求得。不同颗粒含量ε、颗粒粒径d_p下，颗粒对燃气流动的影响不同，喷管各截面的总压变化也不一样。数值计算不能直接得到喷管流量$ {{\dot{m}}_{t}} $的计算式，但如果在发动机每一种可能的工作状态，即不同喷管流量、颗粒含量和颗粒粒径下进行计算，就可以得到每一种ε-d_p组合下$ {{\dot{m}}_{t}}-{{P}_{c}} $的对应关系。根据经验，$ {{\dot{m}}_{t}}-{{P}_{c}} $是线性关系^[24]，而每种ε、d_p下该关系的斜率均不同。使用数据拟合的方法，以ε和d_p为自变量，可得到$ {{\dot{m}}_{t}}-{{P}_{c}} $线性关系之系数的计算式为g_x=g(d_p, ε)。与等熵喷管流量公式比较，最终可把喷管流量计算式写为

$ {{\dot m}_t} = {P_c}{A_t}/\left( {{g_x} \cdot {C^ * }} \right) $

喷管流量公式是燃烧室压强计算式的基础，上述流量公式实质是对等熵流量公式的修正，而这一修正最终会体现为燃烧室压强的修正。经过推导^[24]，燃烧室压强的修正计算式可写为

$ {P_c} = {\left( {{g_x}{\rho _p}{C^ * }a{K_b}} \right)^{1/\left( {1 - n} \right)}} $

由于修正量g_x的计算式是由数据拟合得来的，而数据拟合以数值计算为基础，每一次数值计算都得到了压强P_c，因此燃烧室压强的修正公式相当于把压强的数值解表示成了简单、通用的公式。

拟合公式法的具体实施过程如图 2所示。为了保证上述计算不过于复杂，不考虑颗粒粒径分布的影响，用质量平均粒径d₄₃代替d_p，实际计算证明，这不会带来很大的误差^[27]。正常情况下，每一种配方对应的颗粒含量和燃气属性是一定的，因此在数值计算的过程中应根据热力计算结果设定相应的参数^[21]。对于不同配方、不同发动机，Δε、$ \Delta {{\dot{m}}_{b}} $和Δd_p都不是恒定不变的，而是根据需要使ε_j、$ {{\dot{m}}_{b, i}} $和d_{p, k}满足发动机的指定工作状态，以进行数值计算。

在该方法的执行过程中，如果计算对象只有一个发动机，则参数的变化范围是理论计算的误差。这一误差主要由热力计算造成，热力计算得到的参数有化学平衡/冻结之分，喷管不同截面参数也不一样，理论计算无法避免这一问题，而通过上述方法能显示出颗粒含量等参数对压强的影响规律。在目前的推进剂体系下，推进剂配方和颗粒粒径在不同发动机中通常变化不大，因此由该方法得到的压强公式可用于参数相近的新发动机。

3 发动机计算过程与结果

首先列出HTPB推进剂发动机的计算过程和结果，然后给出铝冰发动机的计算结果。

3.1 复合推进剂发动机

以Farrow的4个复合推进剂发动机为例进行计算^[23]，这些发动机的装药配方及喷管尺寸相同，装药形式为内燃管形，已知参数如表 1所示。

根据这些已知参数，可以发现当这4个发动机稳态工作时，推进剂燃烧产物入口流量变化$ {{\dot{m}}_{b}} $范围是0.025~0.121 kg/s。而考虑到热力计算的误差、热力计算中化学平衡/冻结假设、HTPB分子式等因素的影响，颗粒含量的变化ε范围可选为30%~34%，对应于铝含量ω约为16%~18%的推进剂。原文献没有给定金属燃料的粒径，但对于HTPB复合推进剂，有学者通过实验测得燃烧室及喷管中氧化铝颗粒粒径的变化范围是0.01~100 μm，在颗粒聚集的情况下，最大粒径可能达到84~275 μm，这些颗粒粒径可能呈单峰、双峰或三峰分布，主要分布在2~10 μm的区域^[25]。因此，最终选取的颗粒粒径d_p变化范围是1~10 μm。

计算时把喷管流量变化范围分成11份(M=10)，$ {{\dot{m}}_{b, j}}=0.025+i\Delta {{\dot{m}}_{b}} $，$ \Delta {{\dot{m}}_{b}}=0.00\ 958\ \rm{kg/s} $；把颗粒含量变化范围分成9份(N=8)，ε_j=30%+jΔε，Δε=0.5%；再把颗粒粒径变化范围分成10份(O=9)，d_{p, k}=1+k·Δd_p，Δd_p=1 μm，其中i、j、k分别为从0到M、N、O变化的整数。下面在各个流量、颗粒含量、粒径下分别用数值计算的方法求解内流场和燃烧室压强。

首先保持颗粒含量为ε₀=30%，颗粒粒径为d_{p, 0}=1 μm，在不同喷管流量$ {{\dot{m}}_{b, j}} $下对喷管二维轴对称流场进行计算。最终可得到喷管出口流量$ {{\dot{m}}_{t}} $与燃烧室压强(用燃气总压表示)P_c的关系如图 3所示。

从图 3可以看出，在d_p=30%，d_p=1 μm时，喷管流量$ {{\dot{m}}_{t}} $与燃烧室压强P_c呈线性关系。使用最小二乘法拟合这一关系可得：$ {{\dot{m}}_{t}}={{f}_{x}}\cdot {{P}_{c}} $，其中P_c的单位为atm，f_x=0.00 284 kg/(s·atm^-1)，拟合均方误差R² > 0.999。

在不同颗粒含量ε_j=30%+j·Δε(j > 1)下重新计算$ {{\dot{m}}_{t}} $与P_c的关系，分别对各颗粒含量下的$ {{\dot{m}}_{t}} $-P_c变化关系进行拟合，统计拟合公式中斜率f_x与ε_j的关系(拟合均方误差R² > 0.999)，据此求解g_x=A_t/(f_x·C^*)×101 325，其中A_t=4.56×10^-5 m²，C^*=1 588.2 m/s。最终得到的结果如表 2所示。

把表 2中ε_j和g_x的数据值绘制成散点图，并拟合成曲线，如图 4所示。

在不同颗粒粒径d_p下按照上述步骤分别进行计算，可以得到多个g_x-ε_j变化关系和关系式，把它们绘制在一张图中，如图 5所示。

依次对图中不同颗粒粒径d_p下的ε-g_x关系进行拟合，假设拟合公式为g_x=α·ε+β，确保线性公式中斜率α和截距β的最小二乘法拟合均方误差R² > 0.989。统计各次拟合的结果，分析颗粒粒径与拟合公式系数的关系，如表 3所示。

根据表 3中的数据，绘制散点图，再拟合d_p与α、β之间关系的计算式，其中α拟合的均方误差R²=0.982，β拟合标准误差小于0.01。结果如图 6所示。

可以得到g_x随ε和d_p变化的具体计算关系式为

$ {{g}_{x}}=-(0.76588d_{p}^{0.1964}\varepsilon +1.2575 $

(5)

该式适用于燃烧产物中颗粒含量ε为30%~34%，氧化铝颗粒粒径1~10 μm的情况。

由此可知，喷管流量公式为

$ {{\dot m}_t} = \frac{{{P_c}{A_t}}}{{{g_x} \cdot {C^ * }}} $

(6)

由喷管流量公式可以推得^[25]，发动机零维燃烧室压强的修正公式为

$ {P_c} = {\left[ {{g_x}{\rho _p}{C^ * }a{K_b}} \right]^{\frac{1}{{1 - n}}}} $

(7)

式中：ρ_p为推进剂密度，C^*为特征速度，a为燃速系数，n为压力指数，K_b为发动机燃喉比。相比于零维基本压强计算式，压强修正公式中在括号里多了考虑两相流影响的修正量g_x。由于该修正量是由数值计算得来的，因此式(7)相当于用简便的方法计算了燃烧室压强的数值解。

在工程计算中，推进剂配方中金属含量稍有变化时，用热力计算也很不方便，在式(5)中，把颗粒ε含量替换为金属含量ω。通过热力计算可以得到HTPB推进剂中ε与ω的关系^[21]。

从图 7中可以看出，当推进剂中铝含量小于28%时，均满足：ε=1.887 3ω。把其代入式(5)中，得到：

$ {g_x} = - 1.44544\omega d_p^{0.1964} + 1.2575 $

(8)

式(8)适用于铝含量ω为15.9%~18%的HTPB推进剂固体火箭发动机。其中粒径d_p因发动机及其工作状态不同而不一样，为便于计算用质量平均粒径(d₄₃)表示。平均粒径有多种预估方法，文献[26]特别指出了一种被推荐的模型，该模型的粒径公式为

$ \begin{array}{*{20}{c}} {{d_{43}}\left[ {{\rm{ \mathsf{ μ} m}}} \right] = 3.6304D_t^{0.2932} \cdot }\\ {\left( {1 - \exp \left( { - 0.0008163{\xi _c}{P_c}\tau } \right)} \right)} \end{array} $

(9)

式中：D_t为喷管喉径，ξ_c为凝相摩尔浓度，P_c为燃烧室压强，τ为滞留时间，$ \tau ={{\rho }_{c}}{{V}_{c}}/\rm{\dot{m}} $。

通过式(9)可预估凝相的粒径，粒径变化范围是1~10 μm。把这一粒径用在式(8)中可求得燃烧室压强的修正量g_x，该公式适用于铝含量为16%~18%的HTPB发动机。把修正量g_x用于式(7)，最终可以求解得到燃烧室压强。

3.2 铝冰推进剂发动机

当推进剂中金属含量更高时，可采用相同的方法对压强进行修正。在铝冰发动机中，推进剂的铝含量ω为40%~50%时^[5]，热力计算得ε为75.6%~93.7%。根据已知的发动机装药和推进剂参数^[5]，选择喷管流量$ {{\dot{m}}_{b}} $的变化范围为2~5 kg/s。而铝冰推进剂以纳米铝为燃料，研究表明纳米颗粒的凝相燃烧产物粒径比微米铝的产物粒径小^[28]，因此计算选择颗粒粒径d_p的变化范围为0.05~50 μm。使用同样的方法，最终得到喷管流量和燃烧室压强修正量g_x的计算式为(计算过程参见文献[27])：

$ \begin{array}{*{20}{c}} {{g_{x,铝冰}} \ = \left( {0.09196 - 0.4047\omega } \right) \cdot \ln {d_p} - }\\ {1.4976\omega + 1.4229} \end{array} $

(10)

用上述公式可计算铝含量为40%~50%的铝冰发动机燃烧室压强的修正量。上述公式中凝相粒径d_p本应用最新的测量结果表示，但这一测量工作目前仍未展开，相关研究表明，纳米铝颗粒燃烧产物为微米氧化铝。因此为了简化分析，取式(9)的计算结果作为铝冰推进剂凝相燃烧产物粒径d_p的近似值。

4 燃烧室压强计算方法的验证

在复合推进剂发动机中，金属含量为16%~18%的情况最常见。下面这几个发动机是美国实用的导弹/航天发动机^[29]，使用零维压强修正计算式和基本计算式得到的结果如表 4所示。

从表中计算结果可以看出，压强修正公式估算的燃烧室压强相比未修正的基本公式的估算结果误差明显下降，由10%~13%下降到2%~5.8%。这证明了压强修正计算式可行且计算结果准确度提高。

压强估算不准确的问题在铝冰发动机中表现得很严重，现有铝冰发动机实验包括美国普渡大学的实验和国内哈工程的实验^[30]。2个发动机的参数及基本压强公式、压强修正计算式的计算结果如表 5所示。

从表中数据可以看出，使用数值解修正的压强计算公式计算的燃烧室平衡压强与实验的误差相比基本压强公式明显降低，在国内的铝冰发动机上，相对误差由57%下降到6%以内，在普渡大学的铝冰发动机上，燃烧室压强的相对误差由50%下降到15%左右。虽然最终理论压强与实验值的绝对误差仍然较大，但相比基本计算公式，相对误差已明显降低，这在铝冰发动机设计和实验中已经意味着发动机实验失败的可能性下降。而普渡大学的铝冰发动机实验数据并未详细给定，上述计算的已知参数是总结他们多篇文献的结果得到的^[27]，因此压强的绝对误差偏大也有可能是实验参数的选择造成的。另一方面，发动机燃烧室压强的理论结果的准确度受多种因素的影响，包括推进剂燃烧效率、发动机热损失等。总而言之，上述计算结果的误差在工程中是可以接受的。

由上述计算和验证可以看出，用拟合公式法得到的零维压强修正公式在实际使用中是适用的。

5 结论

1) 对铝含量为16~18%的复合推进剂来说，应用压强修正公式(7)求解得到的燃烧室压强，计算结果与实验的相对误差小于5.6%；

2) 对铝含量达到40%~50%的铝冰发动机来说，用数值解修正的压强计算公式(7)计算的燃烧室平衡压强与实验的误差相比基本压强公式降低，计算精度明显提高；

3) 通过数据拟合的方法，把燃烧室压强的数值解表示为压强的修正计算式，这既保留了数值计算相比一维纯气相理论的准确性，又保留了压强公式在工程计算中的便捷性。对本文计算指定参数范围外的发动机，或其他种类推进剂发动机，也可采用拟合公式法，在指定参数范围内，使用热力计算和流场数值计算，得到燃烧室压强修正公式。