1 引言

金属燃料能显著提高推进剂的爆热和密度，其燃烧产物能抑制推进剂震荡燃烧，因此在推进剂的发展过程中发挥了重要作用^{[1, 2]}。以往金属燃料主要作为燃烧剂用于复合固体推进剂中，在推进剂配方中，其含量约为15%~25%^[3]，但近些年来由于纳米技术、燃料包覆技术的应用，也被作为主要燃料用在铝冰发动机、金属燃料水冲压发动机等新型推进剂发动机中^{[4, 5]}。在铝冰发动机中，金属燃料的含量提高到40%~50%^[6]，在金属燃料水冲压发动机推进剂中，金属含量达73%甚至更高，加入二次入水后，等效于固体火箭发动机的推进剂金属含量约为13%~ 37%^[5]。可见随着固体推进剂发动机的发展和变型，金属的质量含量不断提高。

金属燃料的使用也带来了一些问题。实验发现随着金属含量的增高，发动机的工作效率下降（η_Isp，η_C^*降低）^{[7, 8]}。其中特征速度效率η_C^*下降会导致发动机燃烧室压强估算不准确，影响实验正常进行，这是在发动机预先研究过程中就需要解决的问题^[9]。在HTPB推进剂火箭发动机中，推进剂金属含量在15%~30%变化时，η_C^*由98%下降为93%^[7]，变化不是特别明显，通常发动机实验能够正常进行，用经验公式对C^*进行修正后甚至能够得到准确的压强预估值^[10]。但这一情况仅在常规复合推进剂发动机中适用。在铝冰发动机中，推进剂金属含量为50%时，η_C^*仅为50%左右^[9]，而目前又没有合适的经验公式能够预测这一结果，我们只能诉诸于理论研究。研究发现导致η_C^*下降的主要因素包括推进剂燃烧效率、两相流、发动机热损失等，主流观点是燃烧效率为主导^[10~12]。但是研究也发现单从含金属推进剂的燃烧效率而言，HTPB复合推进剂燃烧效率在96%以上^[13]，铝冰推进剂也能达到85%以上^[14]，显然还有其他因素在影响η_C^*。从目前已有的理论研究方法来看，考虑两相流对η_C^*和燃烧室压强的影响是有价值且必要的^[4]。

发动机内流场中两相流包括燃烧室内两相流和喷管两相流，燃烧室内流速低两相流损失小，在喷管中，只有喉部以前的两相流动会影响燃烧室压强^[11]。这里需要说明的是，金属氧化物颗粒由于布朗运动对燃烧室压强的影响很小，通常可以被忽略^[15]，燃烧室压强是指发动机燃烧室内燃烧产物中气相的总压，而喷管喉部以前两相流动影响了气相的膨胀过程，因此影响了发动机燃烧室压强。按照这一思路，我们只需求解喷管喉部前两相流场，得到喷管两相混合物流量与燃烧室压强的关系，就可在燃烧室压强的计算中考虑两相流的影响。下面使用这一思路进行具体的推导、计算和验证。

2 零维燃烧室压强计算方法

在传统零维内弹道理论中，首先通过推进剂热力计算得到其理论特征速度C^*，然后再用流量守恒推导燃烧室平衡压强计算式^[11]。

质量守恒

$ \frac{{{\rm{d}}\left( {{\rho _{\rm{c}}}{V_{\rm{c}}}} \right)}}{{{\rm{d}}t}} = {\rho _{\rm{p}}}{A_{\rm{b}}}ap_{\rm{c}}^n - {{\dot m}_{\rm{t}}} $

(1)

流量公式

$ {{\dot m}_{\rm{t}}} = {p_{\rm{c}}}{A_{\rm{t}}}/{C^*} $

(2)

稳态工作时

$ {\rm{d}}\left( {{\rho _{\rm{c}}}{V_{\rm{c}}}} \right)/{\rm{d}}t \approx 0 $

推得

$ {p_{\rm{c}}} = {\left( {{\rho _{\rm{p}}}{C^*}a{A_{\rm{b}}}/{A_{\rm{t}}}} \right)^{1/\left( {1 - n} \right)}} $

(3)

式中ρ_p为推进剂密度，A_b为装药燃面面积，A_t为喷管喉部面积，a为推进剂燃速系数，n为压力指数。式（3）即为燃烧室平衡压强的基本计算式，为了修正由于推进剂不完全燃烧、两相流等因素造成的燃烧室压强下降，通常在括号内理论特征速度C^*前乘以修正量η_C^*，以对理论压强进行修正^[12]。

上述公式中，喷管流量公式（2）是通过求解一维等熵定常喷管流动得到的^[11]。当推进剂含有金属燃料，燃烧产物为两相流时，燃气流动不等熵，而颗粒不参与膨胀，因此喷管中混合物的流量会发生变化，这一变化通过求解一维喷管两相流得到。通常一维定常喷管两相流用如下方程组求解^{[15, 16]}。

燃气

$ \left\{ \begin{array}{l} \left( {1 - \varepsilon } \right)\dot m = \left( {1 - \phi } \right){\rho _{{\rm{mg}}}}{u_{\rm{g}}}A\\ \left( {1 - \phi } \right){\rho _{{\rm{mg}}}}{u_{\rm{g}}}{\rm{d}}{u_{\rm{g}}} + \phi {\rho _{{\rm{mp}}}}{u_{\rm{p}}}{\rm{d}}{u_{\rm{p}}} + {\rm{d}}P = 0\\ \left( {1 - \phi } \right){\rho _{{\rm{mg}}}}{u_{\rm{g}}}\left( {{C_{{\rm{p,g}}}}{\rm{d}}{T_{\rm{g}}} + {u_{\rm{g}}}{\rm{d}}{u_{\rm{g}}}} \right) + \\ \phi {\rho _{{\rm{mp}}}}{u_{\rm{p}}}\left( {Cd{T_{\rm{p}}} + {u_{\rm{p}}}{\rm{d}}{u_{\rm{p}}}} \right) = 0\\ P = {\rho _{{\rm{mg}}}}{R_{\rm{g}}}{T_{\rm{g}}} \end{array} \right. $

(4)

颗粒

$ \left\{ \begin{array}{l} \varepsilon \dot m = \phi {\rho _{{\rm{mp}}}}{u_{\rm{p}}}A\\ {u_{\rm{p}}}\frac{{{\rm{d}}{u_{\rm{p}}}}}{{{\rm{d}}x}} = \frac{3}{4}\frac{{\left( {1 - \phi } \right){\rho _{{\rm{mg}}}}{C_{\rm{D}}}}}{{{\rho _{{\rm{mp}}}}{d_{\rm{p}}}}}\left( {{u_{\rm{g}}} - {u_{\rm{p}}}} \right)\left| {{u_{\rm{g}}} - {u_{\rm{p}}}} \right|\\ {u_{\rm{p}}}C\frac{{{\rm{d}}{T_{\rm{p}}}}}{{{\rm{d}}x}} = \frac{{6Nu}}{{{\rho _{{\rm{mp}}}}d_{\rm{p}}^2}}\frac{{{\mu _{\rm{g}}}{C_{{\rm{p,g}}}}}}{{\mathit{Pr}}}\left( {{T_{\rm{g}}} - {T_{\rm{p}}}} \right) \end{array} \right. $

(5)

式中ε为两相流中颗粒质量分数，ϕ为颗粒体积分数，u_g，u_p，T_g，T_p分别为气相和颗粒的速度和温度，p为气相的压强，其他参数见文献[15]。

求解上述方程的入口边界条件为u_g = u_p = u₀ ≈0，p = p_c，T_g = T_p = T₀ ≈ T_f，出口边界条件为p = p_e。其中T₀，u₀为燃烧室温度和流速，T_f为推进剂燃烧温度，p_c为燃烧室压强，p_e为喷管出口压强^[11]。

直接求解上述方程的解析解很困难，但若通过适当的简化，给出燃气和颗粒参数的变化规律，则可由这些规律推导喷管中混合物流量的计算式。为了便于分析，定义 K = u_p/u_g为两相流速度系数，定义L=(T₀ -T_p)/(T₀ -T_g)为两相流温度系数，假定合适的K和L值，根据其值的不同，在不同情况下把两相流简化模型分为以下几种^[17]，见表 1。

表中C₁，C₂表示K或L为常数，T_g或T_p等于T₀表示燃气或颗粒的温度保持为燃烧室温度且不变，无下标的参数表示燃气和颗粒的参数值相等。下面在各个喷管两相流模型下推导对应的喷管流量公式，把它们写为与流量公式（2）类似的形式，然后再推导零维燃烧室压强的计算式。

（1）两相冻结流模型

两相冻结流模型假设颗粒留在燃烧室中不进入喷管，也不影响喷管中燃气的流动^[18]，因此喷管中燃气质量流量的计算式为

$ {{\dot m}_{\rm{t}}} = {{\dot m}_{{\rm{g,t}}}} = {p_{\rm{c}}}{A_{\rm{t}}}/\left[ {\sqrt {{R_{\rm{g}}}{T_0}} /{\mathit{\Gamma }_{\rm{g}}}} \right] $

(6)

式中Γ_g由燃气等熵指数γ_g计算得到

$ {\mathit{\Gamma }_{\rm{g}}} = \sqrt {{\gamma _{\rm{g}}}{{\left[ {2/\left( {{\gamma _{\rm{g}}} + 1} \right)} \right]}^{\frac{{{\gamma _{\rm{g}}} + 1}}{{{\gamma _{\rm{g}}} - 1}}}}} $

(7)

带下标“g”的表示参数两相流中气相的参数。

（2）两相平衡流模型

两相平衡流模型假设颗粒与燃气的参数处处相等^[18]，在该模型下，两相流混合物可以简化为具有等效参数的一种流体，而该流体在喷管中等熵膨胀时具有与纯气相具有类似的热力学规律^[15]。经过推导^[19]，可得喷管中两相混合物流量计算式为

$ {{\dot m}_{\rm{t}}} = {p_{\rm{c}}}{A_{\rm{t}}}/\left[ {\sqrt {1 - \varepsilon } \sqrt {{R_{\rm{g}}}{T_0}} /{\mathit{\Gamma }_{{\rm{mix,e}}}}} \right] $

(8)

式中Γ_{mix, e}的计算式与式（7）类似，但其中γ_g被替换成为γ_{mix, e}，后者为两相平衡混合物的等熵指数，计算式为

$ {\gamma _{{\rm{mix,e}}}} = {\gamma _{\rm{g}}} \cdot \frac{{1 + \varepsilon /\left( {1 - \varepsilon } \right) \cdot C/{C_{{\rm{p,g}}}}}}{{1 + \varepsilon /\left( {1 - \varepsilon } \right) \cdot C/{C_{{\rm{p,g}}}} \cdot {\gamma _{\rm{g}}}}} $

(9)

式中C，C_{p, g}分别为颗粒和燃气的比热。

（3）两相常滞后模型

两相常滞后模型假设颗粒的速度、温度始终以一定的程度落后于燃气，表现为速度系数K，温度系数L均为常数。而经过推导可知，温度系数L 可以通过K求得^[20]

$ L = \frac{{K{\lambda _{\rm{g}}}}}{{K{\lambda _{\rm{g}}} + 3{\mu _{\rm{g}}}C \cdot \left( {1 - K} \right)}} $

(10)

式中C为颗粒的比热容，μ_g，λ_g分别为燃气的导热系数和粘性系数。

常滞后混合物也具有与纯气相、两相平衡混合物相类似的热力学关系，通过推导^[19]得喷管中常滞后混合物的流量公式为

$ {{\dot m}_{\rm{t}}} = {p_{\rm{c}}}{A_{\rm{t}}}/\left[ {\sqrt {1 - \varepsilon } \sqrt {{R_{\rm{g}}}{T_0}} /{\mathit{\Gamma }_{{\rm{mix,c}}}}} \right] $

(11)

式中Γ_{mix, c}的计算式为

$ {\mathit{\Gamma }_{{\rm{mix,c}}}} = \sqrt {\frac{{1 + \varepsilon /\left( {1 - \varepsilon } \right) \cdot LC/{C_{{\rm{p,g}}}}}}{{1 + \varepsilon /\left( {1 - \varepsilon } \right) \cdot {K^2}}}\frac{{2{\gamma _{\rm{g}}}}}{{{\gamma _{\rm{g}}} - 1}} \cdot {\gamma _{{\rm{mix,c}}}}} $

(12)

$ {\gamma _{{\rm{mix,c}}}} = {\left( {\frac{2}{{{\gamma _{{\rm{mix,c}}}} + 1}}} \right)^{\frac{2}{{{\gamma _{{\rm{mix,c}}}} - 1}}}} - {\left( {\frac{2}{{{\gamma _{{\rm{mix,c}}}} + 1}}} \right)^{\frac{{{\gamma _{{\rm{mix,c}}}} + 1}}{{{\gamma _{{\rm{mix,c}}}} - 1}}}} $

(13)

$ {\gamma _{{\rm{mix,c}}}} = {\gamma _{\rm{g}}} \cdot \frac{{1 + \varepsilon /\left( {1 - \varepsilon } \right) \cdot LC/{C_{{\rm{p,g}}}}}}{{1 + \varepsilon /\left( {1 - \varepsilon } \right) \cdot L{\gamma _{\rm{g}}}C/{C_{{\rm{p,g}}}}}} $

(14)

式中γ_{mix, c}为两相常滞后混合物的等熵指数，K和L为颗粒的速度、温度滞后系数。

（4）两相等温流模型

两相等温流模型假设在喷管流动中燃气和颗粒的温度均保持不变，适用于大负载比（颗粒质量分数很大）的喷管两相流，已被实验验证^[21~23]。通过质量守恒方程和动量方程，最终可以推导得喷管喉部压强与燃烧室压强之比p_t/p_c的确定式

$ \left( {\frac{{{\phi _{i,0}}}}{{\left( {1 - {\phi _{i,0}}} \right)}}\frac{{{p_{\rm{t}}}}}{{{p_{\rm{c}}}}} + 1} \right) = \sqrt {\frac{{2{\phi _{i,0}}}}{{\left( {1 - {\phi _{i,0}}} \right)}}\left( {1 - \frac{{{p_{\rm{t}}}}}{{{p_{\rm{c}}}}}} \right) - 2\ln \frac{{{p_{\rm{t}}}}}{{{p_{\rm{c}}}}}} $

(15)

式中ϕ_{i, 0}为燃烧室中颗粒的体积分数，由下式计算

$ {\phi _{{\rm{i}},0}} = \frac{\varepsilon }{{\varepsilon + \left( {1 - \varepsilon } \right) \cdot {\rho _{{\rm{mp}}}}/{\rho _{{\rm{mg,c}}}}}} $

(16)

式中ρ_{mg, c}为燃烧室中燃气密度，ρ_mp为颗粒密度。

这样就可以求得压力比p_t/p_c，由此值可以推导得到喷管两相流量的计算公式

$ {{\dot m}_{\rm{t}}} = \frac{{{p_{\rm{c}}}{A_{\rm{t}}}}}{{\sqrt {{R_{\rm{g}}}{T_0}} }} \cdot \frac{{{p_{\rm{t}}}/{p_{\rm{c}}}}}{{\sqrt {\left( {1 - \varepsilon } \right)\left( {1 - \varepsilon + {K^2}\varepsilon } \right)} }} $

(17)

式中K为颗粒的速度常滞后系数，K=1时为速度平衡的等温流（两相等温流模型-1）。

（5）速度平衡等温流模型

在研究两相等温流时发现^[24]，由于喷管流动过程中温度不变（dT _g = 0），因此当两相流负载比很大时，两相等温混合物的比热比γ_{mix, i} ≈ 1。进一步假设颗粒与燃气的速度平衡（K=1），可以得到更简单的喷管流场解，最终可得喷管流量计算式为

$ {{\dot m}_{\rm{t}}} = \frac{{{p_{\rm{c}}}{A_{\rm{t}}}}}{{\sqrt {{R_{\rm{g}}}{T_0}} }}\frac{1}{{\sqrt {{\rm{e}}\left( {1 - \varepsilon } \right)} }} $

(18)

式中e为自然底数。

（6）颗粒定温模型

颗粒定温模型是两相等温流模型的改进，该模型假设在喷管中燃气的温度下降以提供其混合物热能转化为动能的推动力。此时从两相流三大守恒方程出发，推导发现颗粒定温的混合物也有类似于纯气相、两相平衡、两相常滞后模型的等熵关系，最终得到喷管中混合物流量公式

$ {{\dot m}_{\rm{t}}} = {p_{\rm{c}}}{A_{\rm{t}}}/\left[ {\left({1 -\varepsilon } \right)\sqrt {{R_{\rm{g}}}{T_0}} /{\mathit{\Gamma }_{{\rm{mix, x}}}}} \right] $

(19)

式中Γ_{mix, x}的计算式为

$ {\mathit{\Gamma }_{{\rm{mix,x}}}} = \sqrt {\frac{{2\left( {1 - \varepsilon } \right)}}{{1 - \varepsilon - \varepsilon {K^2}}}} \frac{{{\gamma _{\rm{g}}}}}{{{\gamma _{\rm{g}}} - 1}}\left( {{{\left( {\frac{{2 - 2\mathit{\Omega }}}{{2 - \mathit{\Omega }}}} \right)}^{\frac{{2 - \mathit{\Omega }}}{\mathit{\Omega }}}} - {{\left( {\frac{{2 - 2\mathit{\Omega }}}{{2 - \mathit{\Omega }}}} \right)}^{\frac{{2 - \mathit{\Omega }}}{\mathit{\Omega }}}}} \right) $

(20)

$ \mathit{\Omega } = \frac{{1 - \varepsilon - \varepsilon {K^2}}}{{1 - \varepsilon + \varepsilon K}}\frac{{{\gamma _{\rm{g}}} - 1}}{{{\gamma _{\rm{g}}}}} $

(21)

式中K为颗粒的速度常滞后系数，K=1时为速度平衡的颗粒定温模型（颗粒定温模型-1）。

上述在两相流简化模型下得到的喷管流量公式都可写为式（2）中C^*的修正式

$ {{\dot m}_{\rm{t}}} = {p_{\rm{c}}}{A_{\rm{t}}}/\left( {{\chi _{{C^*}}} \cdot {C^*}} \right) $

(22)

式中χ_C^* 在不同两相流模型下有不同的计算式，可分别列写出如下形式

冻结流

$ {\chi _{{C^*}}} = {C^*} \cdot {\mathit{\Gamma }_{\rm{g}}}/\sqrt {{R_{\rm{g}}}{T_0}} $

平衡流

$ {\chi _{{C^*}}} = {C^*} \cdot {\mathit{\Gamma }_{{\rm{mix,e}}}}/\sqrt {\left( {1 - \varepsilon } \right){R_{\rm{g}}}{T_0}} $

常滞后流

$ {\chi _{{C^*}}} = {C^*} \cdot {\mathit{\Gamma }_{{\rm{mix,c}}}}/\sqrt {\left( {1 - \varepsilon } \right){R_{\rm{g}}}{T_0}} $

(23)

等温流

$ {\chi _{{C^*}}} = {C^*} \cdot \left( {{p_{\rm{t}}}/{p_{\rm{c}}}} \right)/\sqrt {\left( {1 - \varepsilon + {K^2}\varepsilon } \right){R_{\rm{g}}}{T_0}} $

速度平衡等温

$ {\chi _{{C^*}}} = {C^*}/\sqrt {{\rm{e}}\left( {1 - \varepsilon } \right){R_{\rm{g}}}{T_0}} $

颗粒定温流

$ {\chi _{{C^*}}} = {C^*} \cdot {\mathit{\Gamma }_{{\rm{mix,x}}}}/\left[ {\left( {1 - \varepsilon } \right) \cdot \sqrt {{R_{\rm{g}}}{T_0}} } \right] $

上述格式中C^*为由热力计算得到的特征速度，其他各个参量由热力计算结果及制定的颗粒速度系数K最终都可得到，说明考虑两相流影响的喷管流量修正量是可以计算的。

另外还需指明的是，推进剂两相燃烧产物中气相的参数（带下标“g”）是通过热力计算^[25]得到的，热力计算过程中，得到燃烧产物组分之后，气相参数的计算公式为

$ \left\{ \begin{array}{l} \varepsilon = \sum\limits_{j = 1}^L {\left( {{M_j}{n_j}} \right)/1000} \\ {M_{\rm{g}}} = 1000\left( {1 - \varepsilon } \right)/\sum\limits_{j = L + 1}^N {{n_j}} \\ {C_{{\rm{p,g}}}} = \sum\limits_{j = L + 1}^N {\left( {{n_j}{C_{p,j}}} \right)/\left( {1 - \varepsilon } \right)} \\ {\gamma _{\rm{g}}} = {C_{{\rm{p,g}}}}/\left[ {{C_{{\rm{p,g}}}} - {R_0}/{M_{\rm{g}}}} \right]\\ {R_{\rm{g}}} = {R_0}/{M_{\rm{g}}} \end{array} \right. $

(24)

式中j=1，2，…，N为组分代号，其中1~L表示凝相组分，L+1~N为燃气组分，R₀ = 8.314J/(mol ·K)。

以喷管流量公式（22）为基础，通过流量守恒，最终得到发动机燃烧室压强的计算式为

$ {p_{\rm{c}}} = {\left[ {{\chi _{{C^*}}} \cdot {\rho _{\rm{p}}}{C^*}a{A_{\rm{b}}}/{A_{\rm{t}}}} \right]^{1/\left( {1 - n} \right)}} $

(25)

由该公式可以计算考虑两相流影响的发动机燃烧室压强。

3 多维燃烧室压强计算方法

除了用零维燃烧室压强公式估算燃烧室压强外，还可以用发动机内流场数值计算的方法得到燃烧室压强的理论计算结果。由于数值计算通常对发动机二维轴对称/三维内流场，或者喷管流场进行，因此我们可称之为多维燃烧室压强计算方法。推进剂燃烧产物为两相流时，求解方法包括欧拉-欧拉方法和欧拉拉格朗日方法^[26]，其中欧拉-欧拉方法的模型控制方程可写为如下形式。

燃气

$ \left\{ \begin{array}{l} \nabla \cdot \left[ {\left( {1 - \phi } \right){\rho _{{\rm{mg}}}}{\mathit{\boldsymbol{u}}_{\rm{g}}}} \right] = - {{\dot m}_{{\rm{p,g}}}}\\ \nabla \cdot \left[ {\left( {1 - \phi } \right){\rho _{{\rm{mg}}}}{\mathit{\boldsymbol{u}}_{\rm{g}}}{\mathit{\boldsymbol{u}}_{\rm{g}}}} \right] = - \nabla P + \nabla \cdot \tilde \tau + \frac{{{\mathit{\boldsymbol{u}}_{\rm{p}}} - {\mathit{\boldsymbol{u}}_{\rm{g}}}}}{{{\tau _{{\rm{rv}}}}}} + {\mathit{\boldsymbol{X}}_{\rm{g}}}\\ \nabla \cdot \left[ {\left( {1 - \phi } \right){\rho _{{\rm{mg}}}}{\mathit{\boldsymbol{u}}_{\rm{g}}}{C_{{\rm{p,g}}}}{T_{\rm{g}}}} \right] = \nabla \cdot \left( {\nabla {T_{\rm{g}}}} \right) + {\tau _{{\rm{rT}}}} \cdot \left( {{T_{\rm{p}}} - {T_{\rm{g}}}} \right) + {{\bar S}_{\rm{g}}} \end{array} \right. $

(26)

颗粒

$ \left\{ \begin{array}{l} \nabla \cdot \left[ {\phi {\rho _{{\rm{mp}}}}{\mathit{\boldsymbol{u}}_{\rm{p}}}} \right] = {{\dot m}_{{\rm{p,g}}}}\\ \nabla \cdot \left[ {\phi {\rho _{{\rm{mp}}}}{\mathit{\boldsymbol{u}}_{\rm{p}}}{\mathit{\boldsymbol{u}}_{\rm{p}}}} \right] = - \frac{{{\mathit{\boldsymbol{u}}_{\rm{p}}} - {\mathit{\boldsymbol{u}}_{\rm{g}}}}}{{{\tau _{{\rm{rv}}}}}} + {\mathit{\boldsymbol{X}}_{\rm{p}}}\\ \nabla \cdot \left[ {\phi {\rho _{{\rm{mp}}}}{\mathit{\boldsymbol{u}}_{\rm{p}}}C{T_{\rm{p}}}} \right] = - {\tau _{{\rm{rT}}}} \cdot \left( {{T_{\rm{p}}} - {T_{\rm{g}}}} \right) + {{\bar S}_{\rm{p}}} \end{array} \right. $

(27)

式中ṁ_{p, g}为颗粒与燃气的质量交换量，$\tilde \tau $表示燃气粘性应力张量。τ_rv是两相阻力因子，τ_rT是两相传热因子，其计算式参见文献[15, 26]。X_g和X_p分别为燃气和颗粒所受其他力，包括传质量力、Magnus力等，S_g和S_p为燃气和颗粒所受的其他热量源相，包括相变热、辐射传热等。

欧拉-拉格朗日模型的控制方程此处不再列出，可参见文献[26]。在计算燃气的流动时，这两个模型不仅考虑燃气对颗粒的粘性阻力，还考虑燃气自身的粘性，并且还能考虑湍流的影响。相比于表 1中各个一维喷管两相流简化模型，多维两相流场数值计算的方法完整地考虑了两相滞后效应，更接近实际情况，因此通常相比零维公式计算精度也更高一些^{[27, 28]}。

另外，由于多维两相流控制方程需要通过数值迭代法求解，因此也被称为两相流数值解法。本文使用Fluent软件对发动机多维流场和燃烧室压强进行计算，在Fluent中，欧拉-欧拉模型可以分为Euler模型和Mixture模型，其中Euler模型更接近实际，Mixture模型不考虑温度变化；欧拉-拉格朗日模型可以分为颗粒固定轨道模型和随机轨道模型，其中随机轨道模型能考虑湍流对颗粒运动的影响。

4 结果与讨论

在HTPB复合推进剂发动机和铝冰发动机中，分别使用上述燃烧室压强计算方法进行压强估算，与实验结果进行比较，以验证方法可行性。

（1）HTPB复合推进剂发动机

BATES-A^{[7, 29]}发动机是美国空军研制的标准实验发动机，它曾用于研究不同推进剂配方对发动机性能的影响^[7]。该发动机装药形式为内端燃管形，设计的燃面在稳态工作过程中近似恒定。当推进剂配方为AP/Al/HTPB=69/21/10时，通过实验测得发动机稳态工作时燃烧室压强为6.52MPa^[29]。我们首先对这种情况的发动机压强进行计算。

在上述推进剂配方下，使用NASA CEA^[25]软件进行热力计算，热力计算的条件是该发动机的设计压强6.89MPa，燃烧产物参数化学冻结在燃烧室。最终得到燃烧产物中燃气参数、推进剂特征速度、颗粒含量。将其与发动机、装药尺寸列如表 2所示。

根据表中参数及燃气状态方程可求得燃烧室内燃气密度ρ_{mg, c} =4.13kg/m³。而推进剂中金属铝的燃烧产物为凝相氧化铝颗粒，在液态情况下，其比热容C= 1886.9J/kg·K^-1。颗粒粒径通常呈一定规律分布，但可用质量平均粒径d₄₃近似表示^[30]。在文献[31]中，推荐使用如下公式计算颗粒平均粒径

$ {\mathit{d}_{{\rm{43}}}}\left[ {\mu m} \right] = 3.6304D_{\rm{t}}^{0.2932}\left( {1 - \exp \left( { - 0.0008163{\xi _{\rm{c}}}{p_{\rm{c}}}\tau } \right)} \right) $

(28)

式中D_t[in]为喷管喉径，ξ_c[mol /100g]为凝相摩尔浓度，p_c[psi]为燃烧室压强，τ[ms]为滞留时间，τ = ρ_cV_c/ṁ。

根据表 2中数据及发动机参数^[7]可以求得，颗粒的平均粒径d_p ≈4.37μm。另外，在常滞后模型、等温流模型及颗粒定温模型下的零维燃烧室压强公式中，需事先已知颗粒速度系数K。但是在流动情况未知时，颗粒速度系数不好确定。借助数值计算工具，对BATES-A发动机在文献[7]中不同配方下分别进行计算^[19]，可以得到不同凝相含量ε下速度系数K在喷管喉部的变化规律如图 1。

从图中可以看出，不同凝相含量下，颗粒速度系数变化不大，在压强计算中，速度系数K可近似取为K=0.79，而由公式（10）可以求得，L=0.739。在上述各个参数下，使用零维燃烧室压强计算方法、多维（二维轴对称）燃烧室压强计算法，分别进行计算，得到结果如表 3所示。

从表中计算结果可以看出，两相冻结流模型、速度平衡的等温流模型（包括两相等温流模型-1）、速度平衡的颗粒定温模型的计算结果与实验的误差大于15%，说明这几个模型偏离了实际。基本计算式和两相平衡流模型都基于两相平衡流假设，最终压强结果与实验误差大于或接近10%，说明不考虑两相流滞后效应时压强计算结果与实际有一定偏差。这一偏差虽然不会影响实验的进行，通过经验公式也可以修正^[10]，但没有根本解决问题。

使用两相常滞后模型、两相等温流模型、颗粒定温模型计算能把压强计算结果的误差降低到6%以内，说明考虑两相滞后效应影响后，模型的计算精度相比两相平衡流模型有较大的提高。其中两相等温流模型不考虑气相和颗粒的温度变化，相比较于颗粒定温模型，误差明显偏高，说明实际流动中气相的温度是变化的。而两相常滞后模型由于同时考虑了颗粒温度变化，相较于颗粒定温模型，误差仍偏大。另一方面，从表中多维流场数值计算法的结果可以看出，多维燃烧室压强计算方法的计算精度很高，计算结果与实验的误差能下降到2.5%以内。这主要是由于多维内流场数值计算的方法完整地、未简化地考虑了两相不平衡效应，因此更接近实际情况。而多维流场计算由于需要通过数值方法求解控制方程，求解过程复杂，相比零维公式计算量大大增加。

（2）铝冰发动机

铝冰发动机的一个显著特点是燃烧产物中凝相颗粒质量含量达到90%以上。现有一铝冰推进剂配方为Al/H₂O=50/50，把它用在某一锥柱内燃装药的固体火箭发动机中，实验测得该发动机稳定工作段的燃烧室压强约为9.12MPa^[32]。下面取装药燃层厚为10mm的情况讨论，首先用热力计算^[25]得到燃烧产物中燃气的参数、颗粒含量和特征速度，把它们同推进剂、发动机参数一起列如表 4所示。

根据表中参数及燃气状态方程可求得燃烧室内燃气密度ρ_{mg, c}=0.8kg/m³。由公式（28）以及文献[19]中发动机的相关参数，最终可以求得铝冰发动机凝相燃烧产物的质量平均粒径 d_p≈3.28μm。

当铝冰推进剂中金属含量为40%~50%时，使用数值计算的手段，可以得到凝相颗粒在喷管喉部的速度系数K=0.35~0.55变化^[19]，当推进剂铝含量为50%时，可近似取K=0.4，由公式（10）可以求得，温度系数L=0.701。

在上述已知参数下，分别使用零维燃烧室压强计算方法、多维（二维轴对称）压强计算法进行计算。并根据表 3的计算结果及分析，舍去两相冻结流模型、速度平衡的等温流模型（包括两相等温流模型-1）和速度平衡的颗粒定温模型。另外，多维压强计算方法中，颗粒固定轨道模型和随机轨道模型的计算结果相差不大，Mixture模型不考虑温度的变化，偏离实际。因此下面仅用剩余的零维模型压强公式、Euler模型和颗粒随机轨道模型进行计算。最终压强计算结果如表 5所示。

从表 5中结果可以看出，如果不考虑两相流不平衡效应的影响，铝冰发动机的燃烧室压强预估值远远高于实验值，误差达到57.4%（基本计算式和平衡流模型），以此值设计的发动机，必然导致实验失败。考虑两相流不平衡效应后，用零维燃烧室压强计算方法和一维两相流模型，能把燃烧室压强计算结果与实验的误差降低到17.5%以下，其中使用两相等温流模型和颗粒定温模型估算的燃烧室压强结果误差仅为11.5%左右，这说明当两相流负载比很大时，颗粒温度不变的假设接近实际，可以用此假设做理论计算。另一方面，用多维流场数值计算完整地考虑了两相不平衡效应，燃烧室压强的计算结果与实验的误差下降为2.3%以下。

5 结论

通过考虑喷管中两相流对燃烧室压强的影响，本文根据不同的两相流简化模型提出了几种燃烧室压强计算方法。通过比较各种方法的计算结果与实验结果，得出以下结论：

（1）在金属含量较高的固体火箭发动机，如本文计算所用的21%铝含量的HTPB推进剂发动机、50%铝含量的铝冰发动机的压强计算中，考虑两相流的影响是必要的。在铝冰发动机中，使用基本压强计算方法预估的结果与实验误差大于57%，而用两相流修正后，误差下降到11%左右。

（2）在考虑两相流影响的零维压强计算方法中，两相常滞后模型、两相等温流模型、颗粒定温模型都能得到较好的预测结果，在HTPB推进剂发动机中，这几个模型的计算结果相比实验误差低于6%，在铝冰发动机中，两相等温流模型和颗粒定温模型相较于常滞后模型计算结果更准确。

（3）多维两相内流场数值计算的方法计算精度高，在本文算例中，误差低于2.5%，但是由于需要通过迭代求解多维控制方程组，计算过程比零维方法更复杂。在实际计算中，应酌情选择不同的计算方法。