金属燃料由于热值大、成本低，能显著提高推进剂的燃烧温度，在固体推进剂火箭发动机的发展中占有重要地位^[1]。以往金属主要作为燃烧剂用于复合推进剂中，近年来由于纳米材料技术、颗粒包覆技术的成熟和应用，包覆的纳米金属燃料颗粒也作为主要成分用在固体推进剂中^[2]。这些金属燃料以铝为主，使用纳米铝为燃料的铝冰固体推进剂中金属燃料质量含量高达40%~50%^[3]。

推进剂中较高的金属含量使其燃烧产物中凝相颗粒含量很高，喷管内流动是高凝相浓度两相流，传统理论中两相平衡流假设不再适用^[4]。如果仍按该理论设计铝冰发动机，试验就会失败^[5]。由于喷管内凝相颗粒速度低于燃气、温度高于燃气，研究人员曾提出用“常滞后”模型(假设凝相以一定比例滞后于燃气)求解喷管两相流^[6-7]。但由于实际流动并非两相常滞后流，该方法计算误差较大^[8]。为了准确求解喷管两相非平衡流，研究人员做了大量工作，主要研究方法是数值模拟，该方法准确度高，但计算量大且不便于工程实用^[9-11]。陈军用数值计算和数据拟合的方法，得到了喷管两相流量公式^[12]和零维压强公式^[13]，在燃烧室压强的计算中，该方法准确度较高但计算量大。

两相常滞后模型得到了喷管两相流的解析解，据此能得到考虑两相流影响的、简洁的压强修正计算式，但误差较大。为了沿用这一思路并降低压强计算误差，研究人员在高凝相浓度两相流的研究中提出了“两相等温流”模型^[14-15]。该模型考虑到凝相含量高时相间换热量大、混合物温度变化量小的情况，假设燃气、凝相温度在喷管中不变，凝相速度常滞后于燃气^[4]。使用该模型能得到喷管两相流的解析解。但是在实际中，由于燃气在喷管中需要膨胀做功，其温度是变化的，因此只能假设颗粒的温度不变化。所以提出“颗粒定温”模型，拟对一维喷管两相流进行求解。本文在两相等温流模型和颗粒定温模型下，推导喷管两相流量公式和考虑两相不平衡的燃烧室压强修正公式。然后再用50%铝含量的铝冰发动机压强实验结果^[15]对两种模型下得到的燃烧室压强修正公式进行验证。

1 两相等温流模型

在适当的假设下，发动机一维喷管两相流可以通过如下控制方程来求解^{[7, 16]}：

$ 气相:\left\{ \begin{array}{l} \left( {1 - \varepsilon } \right)\dot m = \left( {1 - \phi } \right){\rho _{mg}}{u_g}A\\ \left( {1 - \phi } \right){\rho _{mg}}{u_g}{\rm{d}}{u_g} + \phi {\rho _{mp}}{u_p}{\rm{d}}{u_p} + {\rm{d}}P = 0\\ \left( {1 - \phi } \right){\rho _{mg}}{u_g}\left( {{C_{p,g}}d{T_g} + {u_g}{\rm{d}}{u_g}} \right) + \\ \;\;\;\;\;\;\;\;\phi {\rho _{mp}}{u_p}\left( {Cd{T_p} + {u_p}{\rm{d}}{u_p}} \right) = 0\\ P = {\rho _{mg}}{R_g}{T_g} \end{array} \right. $

(1)

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

(2)

式中：$\dot m$为喷管两相流量，ε为凝相颗粒的质量分数，ϕ为凝相体积分数，u_g、u_p、T_g、T_p分别为燃气和凝相的速度和温度，C_D为两相阻力系数，P为气相的压强，其他各参数见文献[7]。

求解上述方程的入口边界条件为：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为喷管出口压强^[2]。

对两相等温流模型(T_p=T_g=T₀，u_p=Ku_g)，推导喷管两相流量公式和燃烧室压强修正式。假定流动过程中各相热力性质不变，颗粒无相变，燃烧室内凝相含量为ε_m，则喷管任意截面上凝相质量流量百分数ε_i和凝相体积分数ϕ_i分别为^[7]

$ {\varepsilon _i} = \frac{{K{\varepsilon _m}}}{{1 - {\varepsilon _m}\left( {1 - K} \right)}} $

(3)

$ {\phi _i} = \frac{{{\varepsilon _i}}}{{{\varepsilon _i}\left( {1 - K{\rho _{{\rm{mp}}}}/{\rho _{{\rm{mg}}}}} \right) + K{\rho _{{\rm{mp}}}}/{\rho _{{\rm{mg}}}}}} $

(4)

而由凝相和燃气的质量守恒方程知：

$ \frac{{{\varepsilon _i}}}{{1 - {\varepsilon _i}}} = \frac{{{\phi _i}}}{{1 - {\phi _i}}}\frac{{{\rho _{{\rm{mp}}}}}}{{{\rho _{{\rm{mg}}}}}}K $

(5)

把速度滞后数K代入动量守恒方程，得到

$ \left[ {\left( {1 - {\phi _i}} \right){\rho _{{\rm{mg}}}} + {\phi _i}{\rho _{{\rm{mp}}}}{K^2}} \right]u{\rm{d}}u + {\rm{d}}P = 0 $

(6)

结合式(4)~(6)和燃气状态方程，最终得到：

$ \begin{array}{*{20}{c}} {\frac{{1 - {\varepsilon _i} + {\varepsilon _i}K}}{{1 - {\varepsilon _i}}}\frac{1}{{{R_g}{T_0}}}u{\rm{d}}u + }\\ {\left[ {\frac{{{\varepsilon _i}}}{{\left( {1 - {\varepsilon _i}} \right)K\left( {{\rho _{mp}}{R_g}{T_0}} \right)}} + \frac{1}{P}} \right]{\rm{d}}P = 0} \end{array} $

(7)

对式(7)从燃烧室到喷管任意位置积分，得

$ u = \sqrt {\frac{2}{{1 - {\varepsilon _i} + {\varepsilon _i}K}}\left[ {\frac{{{\varepsilon _i}}}{{K{\rho _{{\rm{mp}}}}}}\left( {{P_{\rm{c}}} - P} \right) + \left( {1 - {\varepsilon _i}} \right){R_{\rm{g}}}{T_0}\ln \frac{{{P_{\rm{c}}}}}{P}} \right]} $

(8)

这是喷管中任意位置燃气的流速。

对燃气、凝相的质量方程两端分别进行微分，在喷管喉部(满足(dA/A)_t=0)得到

$ 燃气:{\left( { - \frac{{{\rm{d}}{\phi _i}}}{{1 - {\phi _i}}} + \frac{{{\rm{d}}P}}{P} + \frac{{{\rm{d}}u}}{u}} \right)_t} = 0 $

(9)

$ 凝相:{\left( {\frac{{{\rm{d}}{\phi _i}}}{{{\phi _i}}} + \frac{{{\rm{d}}u}}{u}} \right)_t} = 0 $

(10)

由式(9)和(10)可求出：

$ {\left( {\frac{{{\rm{d}}P}}{P} + \frac{1}{{1 - {\phi _i}}}\frac{{{\rm{d}}u}}{u}} \right)_t} = 0 $

(11)

把式(11)代入式(6)以消除(dP)_t，而由于(du)_t≠0，因此可求出u_t。用式(4)消去其中ϕ_{i, t}，最终得

$ {u_t} = \sqrt {\frac{{\left( {1 - {\varepsilon _i}} \right){R_g}{T_0}}}{{1 - {\varepsilon _i} + {\varepsilon _i}K}}} \left( {\frac{{{\varepsilon _i}{P_t}}}{{\left( {1 - {\varepsilon _i}} \right)K{\rho _{{\rm{mp}}}}{R_g}{T_0}}} + 1} \right) $

(12)

把式(8)写为喷管喉部的形式，并把式(12)代入得

$ \left( {\frac{{{\phi _{i,0}}}}{{\left( {1 - {\phi _{i,0}}} \right)}}\frac{{{P_t}}}{{{P_c}}} + 1} \right) = \sqrt {\frac{{2{\phi _{i,0}}}}{{\left( {1 - {\phi _{i,0}}} \right)}}\left( {1 - \frac{{{P_t}}}{{{P_c}}}} \right) - 2\ln \frac{{{P_t}}}{{{P_c}}}} $

(13)

式中：ϕ_{i, 0}为燃烧室中颗粒体积分数，用式(4)求得，而燃烧室中ε_i=ε_m，K=1。

由式(13)不能直接求解喷管喉部压力P_t，但由热力计算求得ε_m，并估算出合适的燃气密度ρ_{mg, c}后，可以求得喷管喉部压力比P_t/P_c的数值，由此可求得喷管流量。喷管两相流量$\dot m$可写为

$ {{\dot m}_t} = \frac{{\left( {1 - {\phi _{i,t}}} \right){\rho _{mg,t}}{u_t}{A_t}}}{{1 - {\varepsilon _m}}} $

(14)

把燃气状态方程和式(3)、(4)和(12)代入式(14)，化简后得

$ {{\dot m}_t} = \frac{{{P_c}{A_t}}}{{\sqrt {{R_g}{T_0}} }} \cdot \frac{{{P_t}/{P_0}}}{{\sqrt {\left( {1 - {\varepsilon _m}} \right)\left( {1 - {\varepsilon _m} + {K^2}{\varepsilon _m}} \right)} }} $

(15)

通常由燃气参数计算的理论特征速度公式形式为：$C_g^* = \sqrt {{R_g}{T_0}} /{\mathit{\Gamma }_g}$，把该式代入式(15)，得到

$ {{\dot m}_t} = \frac{{{P_c}{A_t}}}{{C_g^ * }} \cdot \frac{{{P_t}/{P_c}}}{{{\mathit{\Gamma }_g}\sqrt {\left( {1 - {\varepsilon _m}} \right)\left( {1 - {\varepsilon _m} + {K^2}{\varepsilon _m}} \right)} }} $

(16)

式中P_t/P_c由方程(13)求得。

得到喷管两相流量公式后，就可以求燃烧室平衡压强的计算式。从推进剂燃面到喷管出口，发动机内流场质量守恒方程为

$ {\rho _c}\frac{{{\rm{d}}{V_c}}}{{{\rm{d}}t}} + {V_c}\frac{{{\rm{d}}{\rho _c}}}{{{\rm{d}}t}} = {\rho _p}{A_b}aP_c^n - {{\dot m}_t} $

(17)

式中：ρ_p为推进剂密度，A_b为装药燃面面积，a为推进剂燃速系数，n为压力指数。

发动机稳态工作时，近似有dρ_c/dt≈0，另外再引入dV_c/dt=A_b·aP_cⁿ，把燃气状态方程及式(16)代入式(17)，并忽略燃气密度ρ_{mg, c}的影响，得

$ {P_c} = {\left[ {\frac{{{\mathit{\Gamma }_g}\sqrt {\left( {1 - {\varepsilon _m}} \right)\left( {1 - {\varepsilon _m} + {K^2}{\varepsilon _m}} \right)} }}{{{P_t}/{P_c}}} \cdot {\rho _p}C_g^ * a\frac{{{A_b}}}{{{A_t}}}} \right]^{\frac{1}{{1 - n}}}} $

(18)

式中：P_t/P_c视为已知参量，由方程(13)求得。

式(18)与传统的两相平衡流理论下燃烧室压强公式：P_c=(ρ_pC^*aA_b/A_t)^1/(1-n)形式类似，但考虑了两相流的影响。由式(18)即可求得两相等温流模型下发动机稳态工作时零维平衡压强。

2 颗粒定温模型

颗粒定温模型假设凝相颗粒在喷管流动的过程中温度不变，而燃气的温度变化，凝相速度常滞后于燃气(T_p=T₀，u_p=K·u_g)。

该模型的流动参数用下标“x”表示，可知各截面凝相含量ε_x和体积分数ϕ_x仍满足关系式(3)、(4)。研究发现凝相体积分数只在其粒径小于1 nm的情况下对燃烧室压强有1%以上的影响^[7]，这种情况很少见，因此在本模型中忽略ϕ_x的影响。把ε_x表达式和颗粒定温流假设代入燃气动量方程，得到

$ \left( {1 + \frac{{{\varepsilon _x}}}{{1 - {\varepsilon _x}}}K} \right)u{\rm{d}}u + {R_g}T\frac{{{\rm{d}}P}}{P} = 0 $

(19)

方程组(1)中燃气的能量守恒方程可写为

$ u{\rm{d}}u = - \frac{{1 - {\varepsilon _x}}}{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}} \cdot {C_{p,g}}{\rm{d}}T $

(20)

联立以上两个方程，得到

$ \frac{{1 - {\varepsilon _x} + {\varepsilon _x}K}}{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}}\frac{{{\gamma _g}}}{{{\gamma _g} - 1}} \cdot \frac{{{\rm{d}}T}}{T} = \frac{{{\rm{d}}P}}{P} $

(21)

从燃烧室到喷管任意位置对其进行积分，得

$ {\left( {\frac{T}{{{T_0}}}} \right)^{\frac{{1 - {\varepsilon _x} + {\varepsilon _x}K{\gamma _g}}}{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}{\gamma _g} - 1}}}} = \frac{P}{{{P_c}}} $

(22)

另一方面，假设燃烧室流速u₀≈0，对式(20)积分，并把式(22)代入，得到喷管任意截面流速：

$ u = \sqrt {\frac{{2\left( {1 - {\varepsilon _x}} \right)}}{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}}\frac{{{\gamma _g}}}{{{\gamma _g} - 1}}{R_g}{T_0}\left( {1 - {{\left( {\frac{P}{{{P_c}}}} \right)}^{\frac{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}{\gamma _g} - 1}}{{1 - {\varepsilon _x} + {\varepsilon _x}K~~{\gamma _g}}}}}} \right)} $

(23)

由式(23)可写出喷管中燃气流量的计算式，最终得到混合物总流量计算式：

$ \begin{array}{*{20}{c}} {\dot m = \frac{{{\rho _{mg}}uA}}{{1 - {\varepsilon _m}}} = \frac{{{\rho _{mg}}A}}{{1 - {\varepsilon _m}}} \cdot }\\ {\sqrt {\frac{{2\left( {1 - {\varepsilon _x}} \right)}}{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}}\frac{{{\gamma _g}}}{{{\gamma _g} - 1}}{R_g}{T_0}\left( {1 - {{\left( {\frac{P}{{{P_c}}}} \right)}^{\frac{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}{\gamma _g} - 1}}{{1 - {\varepsilon _x} + {\varepsilon _x}K\;\;{\gamma _g}}}}}} \right)} } \end{array} $

(24)

把燃气状态方程代入，得到：

$ \begin{array}{*{20}{c}} {\dot m = \frac{{{P_c}A}}{{\left( {1 - {\varepsilon _m}} \right)\sqrt {{R_g}{T_0}} }}\sqrt {\frac{{2\left( {1 - {\varepsilon _x}} \right)}}{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}}\frac{{{\gamma _g}}}{{{\gamma _g} - 1}}} \cdot }\\ {\sqrt {\left( {{{\left( {\frac{P}{{{P_c}}}} \right)}^{2 - 2\left( {\frac{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}{\gamma _g} - 1}}{{1 - {\varepsilon _x} + {\varepsilon _x}K\;\;{\gamma _g}}}} \right)}} - {{\left( {\frac{P}{{{P_c}}}} \right)}^{2 - \frac{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}{\gamma _g} - 1}}{{1 - {\varepsilon _x} + {\varepsilon _x}K\;\;{\gamma _g}}}}}} \right)} } \end{array} $

(25)

把$\dot m$/A对P/P_c求导，在喷喉处导数为零，得到：

$ \frac{{{P_t}}}{{{P_c}}} = {\left( {\frac{{2 - 2\left( {\frac{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}}{{1 - {\varepsilon _x} + {\varepsilon _x}K}}\frac{{{\gamma _g} - 1}}{{{\gamma _g}}}} \right)}}{{2 - \frac{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}}{{1 - {\varepsilon _x} + {\varepsilon _x}K}}\frac{{{\gamma _g} - 1}}{{{\gamma _g}}}}}} \right)^{\frac{{1 - {\varepsilon _x} + {\varepsilon _x}K\;{\gamma _g}}}{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}\;{\gamma _g} - 1}}}} $

(26)

令：

$ \mathit{\Omega } = \frac{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}}{{1 - {\varepsilon _x} + {\varepsilon _x}K}}\frac{{{\gamma _g} - 1}}{{{\gamma _g}}} $

(27)

则式(26)可写为

$ \frac{{{P_t}}}{{{P_c}}} = {\left( {\frac{{2 - 2\mathit{\Omega }}}{{2 - \mathit{\Omega }}}} \right)^{\frac{1}{\mathit{\Omega }}}} $

(28)

把该式代入喷管流量式(25)中，得到

$ \begin{array}{*{20}{c}} {{{\dot m}_t} = \frac{{{P_c}{A_t}}}{{\left( {1 - {\varepsilon _m}} \right)\sqrt {{R_g}{T_0}} }} \cdot }\\ {\sqrt {\frac{{2\left( {1 - {\varepsilon _x}} \right)}}{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}}\frac{{{\gamma _g}}}{{{\gamma _g} - 1}}\left( {{{\left( {\frac{{2 - 2\mathit{\Omega }}}{{2 - \mathit{\Omega }}}} \right)}^{\frac{{2 - 2\mathit{\Omega }}}{\mathit{\Omega }}}} - {{\left( {\frac{{2 - 2\mathit{\Omega }}}{{2 - \mathit{\Omega }}}} \right)}^{\frac{{2 - \mathit{\Omega }}}{\mathit{\Omega }}}}} \right)} } \end{array} $

(29)

令：

$ {\mathit{\Gamma }_x} = \sqrt {\frac{{2\left( {1 - {\varepsilon _x}} \right)}}{{1 - {\varepsilon _x} + {\varepsilon _x}{K^2}}}\frac{{{\gamma _g}}}{{{\gamma _g} - 1}}\left( {{{\left( {\frac{{2 - 2\mathit{\Omega }}}{{2 - \mathit{\Omega }}}} \right)}^{\frac{{2 - 2\mathit{\Omega }}}{\mathit{\Omega }}}} - {{\left( {\frac{{2 - 2\mathit{\Omega }}}{{2 - \mathit{\Omega }}}} \right)}^{\frac{{2 - \mathit{\Omega }}}{\mathit{\Omega }}}}} \right)} $

(30)

式中ε_x用式(3)求得。

最终喷管两相流量公式可写为

$ {{\dot m}_t} = \frac{{{P_c}{A_t}}}{{\left( {1 - {\varepsilon _m}} \right)\sqrt {{R_g}{T_0}} /{\mathit{\Gamma }_x}}} $

(31)

同样令：$C_g^* = \sqrt {{R_g}{T_0}} /{\mathit{\Gamma }_g}$，式(31)可写为

$ {{\dot m}_t} = \frac{{{P_c}{A_t}}}{{\left( {1 - {\varepsilon _m}} \right){\mathit{\Gamma }_g}/{\mathit{\Gamma }_x} \cdot C_g^ * }} $

(32)

式(32)形式上与式(16)相同。根据质量守恒的原则，采用与第1节相同的推导方法，最终得到颗粒定温模型下燃烧室平衡压强的计算公式为

$ {P_{\rm{c}}} = {\left[ {\left( {1 - {\varepsilon _m}} \right)\frac{{{\mathit{\Gamma }_g}}}{{{\mathit{\Gamma }_x}}} \cdot {\rho _{\rm{p}}}C_g^ * a{K_{\rm{b}}}} \right]^{\frac{1}{{1 - n}}}} $

(33)

式中Γ_x由式(30)求得。

3 计算和验证 3.1 模型本身的验证

两相等温和颗粒定温两种模型通过简化假设把两相流引入到燃烧室压强的计算中，但这些假设是否合理还需要验证。两相流数值计算方法往往能得到与实验接近的结果^{[8, 17]}，可用于验证本文的两相流模型。选用Fluent软件中Euler模型进行数值计算和验证，计算网格和边界条件如图 1所示。

该喷管来自于Bates-A发动机，该发动机稳态工作时喷管两相流量近似为9.153 kg/s^[18]，在该流量入口条件下进行计算。

1) 两相等温流模型的验证。

两相等温流假设是在发动机喷管中凝相含量较高的情况下提出的^[13]，因此该模型特别适用于铝冰发动机的两相流求解。但是为了比较各种金属含量下的流动情况，假设使用图 1中喷管的发动机装药分别为：铝含量为15%、27%的HTPB推进剂(配方为Al/AP/HTPB，凝相含量ε分别为28.3%和47.6%)，铝含量为40%、50%的铝冰推进剂(ALICE，配方为Al/H₂O，凝相含量分别为75.6%和93.7%)，其中凝相含量由热力计算^[19]得到。再假设凝相粒径为5 μm，对上述四种情况下喷管两相流进行数值计算，得到各情况下凝相和燃气温度随空间变化。把温度值用热力计算求得的燃烧室温度进行无量纲化，最终得到四种配方下凝相和燃气的无量纲温度随推进剂铝含量和轴线坐标的变化如图 2所示。从图 2中可以看出，随着推进剂金属铝含量的增高，喷管中气相和凝相的温度变化逐渐变小。推进剂铝含量为15%时，从喷管入口到出口气相和凝相的温度约变化了35%，而推进剂铝含量增高到50%时，温度变化约14%。由此证明推进剂金属含量较高，其产物中凝相含量较高时，两相等温流模型的相关假设在喷管中是成立的。

2) 颗粒定温模型的验证。

颗粒定温模型假设燃气的温度在喷管中变化，而凝相的温度则不变。为了验证颗粒定温假设，在不同凝相粒径下对50%铝含量的铝冰推进剂发动机在图 1所示喷管中的两相流进行计算，得到凝相温度随喷管轴向坐标的变化如图 3所示。

从图 3中可以看出，随着凝相颗粒粒径增大，其温度在喷管中的变化逐渐减小，当粒径从5 μm变化到200 μm时，在整个喷管中凝相温度变化由14%下降为4%左右，这说明颗粒温度不变的假设在金属含量较高的发动机喷管中是适用的。另一方面，颗粒温度变化小说明凝相和燃气间换热量减少，对燃气而言就是减少了热量补充，因此燃气在喷管中温度下降会更多一些。这在喷管喉部以前体现得很明显，从数值计算结果可以得到燃气和凝相温度在喷管喉部前变化如图 4所示。

从图 4中可以看出，随着凝相颗粒粒径的增大，燃气的温度变化得越来越大，而凝相的温度变化得越来越小。当凝相粒径为5 μm时，燃气和凝相的温度变化过程基本一致，但粒径增大到200 μm后，从喷管入口到喷管喉部，凝相温度变化约2%，而燃气温度变化了14%。由此可见，凝相含量高且粒径较大时，简化假设中不能忽略燃气温度的变化。在喷管的扩张段，由于惯性作用，较大凝相颗粒集聚在喷管轴线附近对燃气集中换热，使燃气的温度回升，同时喷管壁面处会出现无颗粒区^[9-10]。但是由于已经过了喉部，这一现象对燃烧室压强的影响不大，不是本文研究范畴。从上述分析可以发现，在喷管喉部以前，颗粒定温流假设是符合铝冰发动机的喷管实际流动情况的。

3.2 凝相速度滞后系数的选择

在使用两相等温流模型或颗粒定温模型进行计算前，需要假定一个合适的颗粒速度滞后系数K。该系数是喷管两相流的特征参数，因此与凝相含量ε、凝相颗粒粒径d_p有关。使用数值计算的方法，可以确定在不同ε、d_p时，速度滞后系数K=u_p/u_g的数值，进而得到其计算公式，最终用于压强修正式(18)和(33)的工程计算。

选用图 1中喷管计算，假设喷管两相流量为9.153 kg/s，首先选定凝相粒径为50 μm，在铝冰发动机装药铝含量分别为40%、43.5%、47.6%、50%四种情况下分别进行计算。最终得到喷管喉部前凝相速度滞后系数K随凝相含量ε和空间位置的变化如图 5所示。

从图 5中可以看出，K值在喷管喉部前基本保持不变，因此喷管收敛段内近似是两相速度常滞后流动。另外，不同凝相含量下K值不同。

凝相粒径对K值的影响通过在凝相粒径为1、5、25、50、100、200 μm的情况下计算得到。下面首先对铝含量为50%的铝冰发动机喷管两相流进行计算，得到喷管喉部前轴线上不同位置处K随粒径d_p的变化如图 6所示。

从图 6中可以看出，凝相粒径越大，K越小，其变化量也逐渐变小，说明粒径越大两相流越趋近于速度常滞后流动。

在不同凝相含量ε、凝相粒径d_p时分别进行上述计算，并对计算结果进行统计和分析。考虑到喷管入口效应的影响，取喷管轴向坐标为0.03~0.14 m凝相速度滞后系数的平均值为常滞后系数K，得到K随ε和d_p的变化如表 1所示。

用最小二乘法拟合K值随凝相含量ε和粒径d_p的变化关系式，结果如下(R²＞0.97)

$ K = 1.058\;8d_{\rm{p}}^{ - 0.9363 \cdot \varepsilon + 0.3499} $

(35)

使用式(35)可以求得喷管中凝相速度常滞后系数的近似值。该式适用于铝含量为40%~50%，凝相含量为75.6%~93.7%的铝冰推进剂，其产物在喷管中凝相粒径为1~200 μm的铝冰发动机。

3.3 压强修正公式的验证

对两相等温流模型和颗粒定温流模型下得到的燃烧室压强修正公式分别进行计算，并与传统压强计算公式、数值计算方法得到的结果进行比较，用压强的实验结果验证模型的准确性。

铝冰发动机装药为锥柱形，工作过程中燃面近似恒定，实验测得的燃烧室平衡压强约为90 Pa^[15]。装药配方Al/H₂O比为1：1，密度ρ_p=1 620 kg/m³，燃速系数a=4.85×10^-4 m/s·Pa^-n，压力指数n=0.229 9。通过热力计算^[19]得到的计算参数如表 2所示。

当发动机装药的燃层厚web=10 mm时，燃喉比K_b=416.75。计算假定凝相全部为液态氧化铝，在喷管中无相变，C=1 886.9 J/(kg·K)，ρ_mp=3 020 kg/m³。另外，根据上述参数和燃气状态方程可求得燃烧室内燃气密度ρ_{mg, c}=0.8 kg/m³。喷管中凝相的粒径是比较难确定的，且通常呈一定规律分布，但可用质量平均粒径d₄₃近似表示。在文献[20]中，推荐使用如下公式计算凝相颗粒的平均粒径：

$ \begin{array}{*{20}{c}} {{d_{43}}\left[ {{\rm{ \mathsf{ μ} m}}} \right] = 3.630\;4D_{\rm{t}}^{0.2932} \cdot }\\ {\left[ {1 - \exp \left( { - 0.000\;816\;3{\xi _{\rm{c}}}{P_{\rm{c}}}\tau } \right)} \right]} \end{array} $

(36)

式中：D_t为喷管喉径，in; ξ_c为凝相浓度，mol/100g;P_c为燃烧室压强，psi; τ为滞留时间，ms, τ=ρ_cV_c/$\dot m$。

由式(36)可求得凝相颗粒粒径d_p≈3.92 μm。根据已知的凝相含量和凝相粒径，由式(35)求得凝相速度常滞后系数K=0.515。在两相等温流模型、颗粒定温流模型下分别对发动机稳态燃烧室压强进行计算，再采用传统压强公式、数值方法进行计算，并与实验比较，结果如表 3所示。

从表 3中计算结果可以看出，传统压强公式估算的铝冰发动机燃烧室压强与实验相差很大，而考虑两相不平衡效应的影响后，压强估算的精度大大提高。其中，两相等温流模型下压强修正公式的计算结果与实验误差约为13%，而颗粒定温模型下压强修正公式结果与实验误差约为6.6%，该误差值已接近工程误差的范围。虽然这些模型下压强结果与实验压强相差仍较大，但相比传统压强公式，上述两种模型的计算精度都有明显的提高。

4 结论

1) 在铝冰发动机燃烧室压强的计算中，应该考虑两相不平衡效应的影响，对传统压强计算公式进行修正，以降低压强估算的误差；

2) 颗粒定温模型相比两相等温流模型考虑了喷管中燃气的温度变化，因此更接近实际情况，精度更高。在该模型下，压强修正公式的计算结果相比实验误差低于7%，这说明该模型可用于指导铝冰发动机试验。