第八章  固体火箭发动机内弹道计算

 内弹道计算是固体火箭发动机设计及性能预估的重要内容，在装药计算完成以后，结合喷管参数、点火系统参数、推进剂数据进行内弹道计算。影响内弹道计算精度的因素有很多，本文在计算方法上采用常规的计算方法，但是综合考虑了燃烧室热损失、喷管效率、喉部烧蚀、压强上升、工作和拖尾段、单一燃速公式和多燃速公式等各种情况。对于短粗的发动机即长径比小于6的情况，可以认为燃烧室内参数分布为均匀，即按照零维内弹道计算可以满足精度要求，如果长径比较大，侵蚀效应显著，则必须按一维进行内弹道计算。对于有两种以上推进剂所组成的装药，则必须进行双燃速内弹道计算。内弹道计算可分为典型的三段，即点火过程的压力上升阶段、平衡阶段和拖尾阶段。本文只介绍点火阶段和工作阶段，拖尾段由方程统一处理，不做阶段划分。本文针对工程中遇到的问题，介绍零维内弹道计算、双燃速零维内弹道计算、一维内弹道计算、双燃速一维内弹道计算和燃气发生器内弹道计算等。

8.1  零维内弹道计算

对于短粗的发动机或长细比小于6的情况，可以认为燃烧室内参数是均匀的，即燃烧室头部和尾部燃速、压强、温度分布等是一致的。此种内弹道计算称为零维内弹道（nondimensional ballistics）计算，零维内弹道计算假定装药燃面肉厚\({{A}_{bj}}\tilde{\ }{{e}_{j}}\)（j=1,…,n_g）已知（由装药设计软件来完成），推进剂数据已知，并考虑喷管烧蚀和效率。本文首先介绍计算模型，然后介绍计算方法和相关参数。其基础知识可参阅文献【1】。

8.1.1  点火计算模型

不考虑点火器点火药流率点火过程计算。

压强建立过程从点火开始，到喷管堵盖破裂压强P_ig最后到平衡压强P_c,eq，采用了两种计算模型：一种是利用公式（8-1）计算：

                 \(t=\frac{1}{1-n}\frac{{{V}_{c}}}{{{\Gamma }^{2}}{{C}^{*}}{{A}_{t}}}\ln \left[ \frac{{{\rho }_{P}}{{C}^{*}}aK-P_{ig}^{1-n}}{{{\rho }_{p}}{{C}^{*}}aK-P_{c}^{1-n}} \right]\) （8-1）

\(\Gamma \)－与气体参数相关的常数

V_c－为燃烧室空腔自由容积m³；

C*－推进剂特征速度m/s；

ρ_P－推进剂密度；

a－为推进剂的燃速系数；

K－发动机面喉比；

P_ig－点火压强；

P_c－燃烧室压强MPa；

n－推进剂压力指数。

8.1.2  平衡段及拖尾段

平衡段及拖尾段采用公式（8-2）和（8-3）计算：

\(\dot{m}={{K}_{a}}{{a}_{0}}P_{C}^{n}{{A}_{b}}{{\rho }_{P}}\exp \left[ {{\sigma }_{P}}\left( T-{{T}_{0}} \right) \right]\)                     （8-2）

\(\frac{d{{P}_{c}}}{dt}=\frac{{{\eta }_{c}}^{2}{{R}_{g}}{{T}_{g}}}{{{V}_{c}}}\left[ \dot{m}-\frac{{{P}_{c}}{{\eta }_{P}}{{A}_{t}}\varphi }{{{C}^{*}}{{\eta }_{c}}} \right]\times {{10}^{-6}}\)                   （8-3）

式中：

\({{\eta }_{P}}\)－为喷管流动总压恢复系数，对于长尾管\({{\eta }_{P}}\)<1，其它喷管取\({{\eta }_{P}}\)＝1；

\(\dot{m}\)－燃气的质量生成率kg/s；

T－计算温度；

T₀－参考温度，一般是设计状态温度，可选20℃；

\({{\eta }_{c}}\)－推进剂的特征速度修正系数，\({{\eta }_{c}}\le 1\)；

φ－喷管的流量修正系数；

R_g T_g－燃气的气体常数和燃烧温度。

A_b－装药燃烧面积，由装药设计软件计算；

\({{a}_{0}}\)－为推进剂在T₀时的燃速系数\(\left[ {\text{m}}/{\left( \text{MP}_{\text{a}}^{\text{n}} \right)\bullet \text{s}}\; \right]\)；

公式中，当压强小于P_af（工作终点压强）时全部工作过程结束。

（8-2）、（8-3）式可采用四阶RongKutta求解。

如果采用瞬时平衡压强法，（8-3）式简化为：

\(\frac{{{\eta }_{c}}^{2}{{R}_{g}}{{T}_{g}}}{{{V}_{c}}}\left[ \dot{m}-\frac{{{P}_{c}}{{\eta }_{P}}{{A}_{t}}\varphi }{{{C}^{*}}{{\eta }_{c}}} \right]\times {{10}^{-6}}=0.0\)

8.1.3  计算方法

已知的关系

1. 肉厚和燃面的对应关系\({{A}_{b}}_{j}\tilde{\ }{{W}_{j}}\)

由装药设计模块计算出的肉厚和装药各参量间的关系，并分成j（j=0,1,…,steps）分，通过给定肉厚的值可以插值来获得燃面的数据。

1. 喉衬的烧蚀关系

根据模拟发动机或条件相近、喉衬材料相同的有关试验绘制喉部烧蚀量和时间的关系曲线\(\Delta {{R}_{tj}}\tilde{\ }{{t}_{j}}\)，本文给出了喉衬的线烧蚀率E_k（1/m），并有如下关系：

\({{R}_{t}}={{R}_{t0}}\)      \(t\le {{s}_{t0}}\) 

   \({{R}_{t}}={{R}_{t0}}+{{E}_{k}}\left( t-{{s}_{t0}} \right)\times {{10}^{3}}\)  \( t>{{s}_{t0}}\)                  （8-4）

式中： \({{s}_{t0}}\)－喉部开始烧蚀时间（s）；

\({{R}_{t0}}\)－喉部初始半径（mm）。

如果不是线性关系，如烧蚀率是压强的0.8次方关系，则喉径要采用迭代的方式加以确定。

1. 燃速系数的关系

\([r={{K}_{a}}a{{P}^{n}}\)

\(a={{a}_{0}}\exp \left[ {{\sigma }_{P}}\left( T-{{T}_{0}} \right) \right]\)                   （8-5）

式中：K_a－为燃速修正系数，缺省K_a=1，如果所给的燃速数据为燃速仪中的数据，在应用到发动机时，必须进行修正，此时\({{K}_{a}}\ne 1\)。

8.1.4  计算公式

1. 上升段计算

如果采用公式（8-1）计算时，将压强建立过程从喷管堵盖破裂压强P_ig到平衡压强P_c,eq0，共分为N段

\({{P}_{c,i}}={{P}_{ig}}+\frac{i}{N}\left( {{P}_{c,eq0}}-{{P}_{ig}} \right)\)                             （i=0,1,…,N）

利用（8-1）式可求出对应的时间。

如果采用固相点火模型，参考8.5节燃气发生器内弹道计算。

1. 平衡段及拖尾段压强计算

将（8-2）式代入（8-3）式利用Ronge Kutta求解（8-6）式可以求得P_ci~t_i       

\(\frac{d{{P}_{c}}}{dt}=f({{P}_{c}},t)\)                          （8-6）

关于Ronge Kutta法的详细知识，请参阅清华大学出版社《C语言算法程序集锦》，在此不作详细介绍。

1. 其它参数计算（下标j为离散时间步j节点）

喉部面积

\({{A}_{tj}}=\pi R_{tj}^{2}\)                                                         （8-7）

相邻两点平均压强（Mpa）：

\({{\bar{P}}_{cj}}=\frac{1}{2}\left( {{P}_{cj}}+{{P}_{c\left( j-1 \right)}} \right)\)                                              （8-8）

相邻两点的时间间隔（s）:

\(\Delta t={{t}_{j}}-{{t}_{j-1}}=\frac{{{W}_{j}}-{{W}_{j-1}}}{a\bar{P}_{cj}^{n}}\)                                （8-9）

推力（kN）：

\({{F}_{j}}={{\xi }_{F,j}}{{C}_{F,j}}{{A}_{tj}}{{P}_{cj}}{{\sigma }_{f}}{{\eta }_{P}}\times {{10}^{3}}\)                               （8-10）

相邻两点平均推力

\({{\bar{F}}_{j}}=\frac{1}{2}\left( {{F}_{j}}+{{F}_{j-1}} \right)\)                                                  （8-11）

下面介绍几个概念：

发动机燃烧过程结束时的压强，为燃烧终点压强P_tb，在这之后是排气过程直到工作终点压强\({{P}_{af}}=0.3\text{MPa}\)（一般情况如此，如果是高空工作，则还要低）为止。但是大部分的装药燃面都有拖尾现象，因此难以判断燃烧过程的终止，有些固体发动机，已知有效工作终

点压强P_eff，常取\({{P}_{eff}}=1.73\text{MPa}\)。本文建议采用有效工作终点压强，代替燃烧工作终点压

强，即\({{P}_{tb}}={{P}_{eff}}\)燃烧过程终止。一些参考书中取拖尾段和平衡段压强曲线的切线的角平分线

与压强曲线交点作为燃烧室过程结束点。

推力冲量和压强冲量\({{I}_{Fj}},{{I}_{Pj}}\)

\({{I}_{Fj}}=\sum\limits_{j=1}^{j}{{{F}_{j}}\left( {{t}_{j}}-{{t}_{j-1}} \right)}\)                                       （8-12）

\({{I}_{Pj}}=\sum\limits_{j=1}^{j}{{{P}_{j}}\left( {{t}_{j}}-{{t}_{j-1}} \right)}\)                                        （8-13）

j=bf燃烧时间推力冲量和压强冲量\({{I}_{tb}},{{I}_{Pb}}\)

\({{I}_{tb}}=\sum\limits_{j=1}^{bf}{{{F}_{j}}\left( {{t}_{j}}-{{t}_{j-1}} \right)}\)                                        （8-14）

\({{I}_{Pb}}=\sum\limits_{j=1}^{bf}{{{P}_{j}}\left( {{t}_{j}}-{{t}_{j-1}} \right)}\)                                        （8-15）

j=ef有效推力冲量和压强冲量\({{I}_{ef}},{{I}_{Pef}}\)

\({{I}_{ef}}=\sum\limits_{j=1}^{ef}{{{F}_{j}}\left( {{t}_{j}}-{{t}_{j-1}} \right)}\)                                       （8-16）

\({{I}_{Pef}}=\sum\limits_{j=1}^{ef}{{{P}_{j}}\left( {{t}_{j}}-{{t}_{j-1}} \right)}\)                              （8-17）

j=af推力冲量和压强冲量\({{I}_{af}},{{I}_{Paf}}\)

\(I=\sum\limits_{j=1}^{af}{{{F}_{j}}\left( {{t}_{j}}-{{t}_{j-1}} \right)}\)                                      （8-18）

\({{I}_{P}}=\sum\limits_{j=1}^{j}{{{P}_{j}}\left( {{t}_{j}}-{{t}_{j-1}} \right)}\)                                       （8-19）

求出以上参数，可以求出平均推力和平均压强：

燃烧时间平均推力和平均压强：

\({{\bar{F}}_{tb}}={{{I}_{tb}}}/{{{t}_{b}}}\;\)                                                       （8-20）

\({{\bar{P}}_{tb}}={{{I}_{Pb}}}/{{{t}_{b}}}\;\)                                                  （8-21）

有效工作时间平均推力和平均压强：

\({{\bar{F}}_{ef}}={{{I}_{ef}}}/{{{t}_{ef}}}\;\)                                                       （8-22）

\({{\bar{P}}_{ef}}={{{I}_{Pef}}}/{{{t}_{ef}}}\;\)                                                    （8-23）

工作时间平均推力和平均压强：

\({{\bar{F}}_{ta}}={{{I}_{{}}}}/{{{t}_{a}}}\;\)                                                               （8-24）

\({{\bar{P}}_{ta}}={{{I}_{P}}}/{{{t}_{a}}}\;\)                                                                 （8-25）

排出燃气的质量流率（kg/s）：

\({{\bar{M}}_{gr}}={{10}^{6}}\frac{1}{2{{C}^{*}}}\left[ {{P}_{Cj}}{{A}_{tj}}+{{P}_{Cj-1}}{{A}_{tj-1}} \right]\)                                   （8-26）

8.1.5  分析讨论

对于常规发动机，上述方法可以满足要求，对于特殊情况，必须特殊处理，以下将讨论喷管过膨胀、长尾喷管、斜切喷管及有推力终止情况下的内弹道计算。

（1）非设计状态下的内弹道计算

喷管处于设计状态下工作时（P_e=P_a），能获得最大推力，但是火箭发动机工作时工作高度、环境压强P_a都在不断的变化中，因此绝大多数都是在非设计状态下工作。在欠膨胀状态下P_e>P_a，膨胀是连续的，在喷管扩张段下游形成一系列的激波。在喷管出口面外的自由射流，形成了连续的扩散－收敛的自由射流边界。此时推力系数可以按标准的公式进行计算。

\({{C}_{F}}=\Gamma \sqrt{\frac{2k}{k-1}\left[ 1-{{\left( \frac{{{P}_{e}}}{{{P}_{a}}} \right)}^{\frac{k-1}{k}}} \right]}+\frac{{{A}_{e}}}{{{A}_{t}}}\left( \frac{{{P}_{e}}}{{{P}_{c}}}-\frac{{{P}_{a}}}{{{P}_{c}}} \right)\)              （8-27）

在过膨胀状态（P_e<P_a）工作时，出现的物理现象比较复杂。根据P_e<P_a的程度不同，燃气

在喷管的流动状态就不同。由实验结果可知当\({{\text{P}}_{\text{e}}}\text{/}{{\text{P}}_{\text{a}}}\ge \text{0}\text{.3}\)时喷管的流动正常；当P_e/P_a<0.3~0.4

时喷管中出现激波，激波强度与压强比P_e/P_a有关，激波下游流动就分离使燃气压强P_e恢复到环境压强P_a。

过膨胀喷管流动是否出现分离可以用下式判断：

\(\frac{{{P}_{i}}}{{{P}_{a}}}=\frac{2}{3}{{\left( \frac{{{P}_{a}}}{{{P}_{c}}} \right)}^{{}^{1}/{}_{5}}}\)                      （8-28）

式中：P_i——分离压强；

如果Pe/Pa<Pi/Pa则流动分离。流动分离后，喷管修正性能就要损失。下面介绍推力系数法。

当Pe/Pa>Pi/Pa时，喷管流动是正常的，可按标准推力系数计算公式计算。

当Pe/Pa<Pi/Pa时，喷管流动出现分离，Pi是起始分离的压强，Ps是出现分离涡流的起始压强，P_0.95是压强等于环境压强的0.95倍，此点实际上是说明分离涡流的范围。P_0.95后压强等于环境压强。

燃气分离后喷管推力系数由两项合成。即分离点上游部分的推力系数及分离点下游部分的推力系数。分离点上游的推力系数可用（8-27）式计算。分离点下游的推力系数\({{C}_{{{F}_{u}}}}\)可使用Kalt Badal经验关系式：

\({{C}_{{{F}_{s}}}}=0.55\left[ \frac{{{P}_{i}}+{{P}_{0.95}}}{{{P}_{c}}} \right]\left( {{\varepsilon }_{0.95}}-{{\varepsilon }_{i}} \right)+0.975\frac{{{P}_{a}}}{{{P}_{c}}}\left( {{\varepsilon }_{e}}-{{\varepsilon }_{0.95}} \right)\)        （8-29）

若\({{\varepsilon }_{i}}\le \frac{{{\varepsilon }_{e}}}{1.604}+0.377\)，则\({{\varepsilon }_{0.95}}-{{\varepsilon }_{i}}=\frac{{{\varepsilon }_{i}}-1}{2.4}\)

若\({{\varepsilon }_{i}}>\frac{{{\varepsilon }_{e}}}{1.604}+0.377\)，则\({{\varepsilon }_{0.95}}-{{\varepsilon }_{i}}=\frac{{{\varepsilon }_{e}}-{{\varepsilon }_{i}}}{1.45}\)

式中：\({{\varepsilon }_{i}}\)——起始分离点扩张比；

\({{\varepsilon }_{0.95}}\)环境压强0.95倍处的压强。最后总的推力系数为：

\({{C}_{F}}={{C}_{{{F}_{u}}}}+\Delta {{C}_{{{F}_{s}}}}-\left( \frac{{{P}_{a}}}{{{P}_{c}}} \right){{\varepsilon }_{e}}\)                  （8-30）

\({{\varepsilon }_{i}}\)可以按下式计算：

\({{\varepsilon }_{i}}={{(\frac{{{r}_{i}}}{{{r}_{t}}})}^{2}}=\frac{{{(\frac{2}{k+1})}^{\frac{1}{k-1}}}{{(\frac{k-1}{k+1})}^{\frac{1}{2}}}}{{{\left[ {{(\frac{{{P}_{i}}}{{{P}_{C}}})}^{\frac{2}{k}}}-{{(\frac{{{P}_{i}}}{{{P}_{C}}})}^{\frac{k+1}{k}}} \right]}^{\frac{1}{2}}}}\)                   （8-31）

\({{P}_{i}}={{P}_{a}}\frac{2}{3}{{\left( \frac{{{P}_{a}}}{{{P}_{c}}} \right)}^{{}^{1}/{}_{5}}}\)                       （8-32）

更为详细内容可参考喷管设计中关于分离流的计算。

（2）亚声速长尾喷管内弹道计算

亚声速长尾喷管是为了将发动机的重心放在全弹重心位置，而在喷管收敛段加一管段的加长喷管，由于管段的摩擦作用，使发动机喷管入口的总压降低，长尾管内弹道计算主要是计算管段的总压损失系数\({{\eta }_{P}}\)。长尾管参数如图（4-4）所示，

计算采用（8-33）式：

\(\frac{1}{\lambda _{1}^{2}}+\ln \lambda _{1}^{2}-\frac{1}{\lambda _{2}^{2}}-\ln \lambda _{2}^{2}=\frac{k}{k+1}f\frac{L}{D}\)                 （8-33）

式中： D为管段的直径，单位m；

L为管段的长度，单位m；

f为管段的摩擦系数，可取f=0.00925，根据粗糙度不同，可选适当的值。

入口段的速度系数\({{\lambda }_{2}}\)可由喷管收缩比\({{\varepsilon }_{i}}\)，采用下式计算：

\(q\left( {{\lambda }_{2}} \right)={{\varepsilon }_{i}}\)                           （8-34）

将\({{\lambda }_{2}}\)代入（8-33）式，利用迭代计算，可以计算出\({{\lambda }_{1}}\)，总压恢复系数采用（8-35）式计算。

\({{\eta }_{P}}=\frac{q\left( {{\lambda }_{1}} \right)}{q\left( {{\lambda }_{2}} \right)}\)                          （8-35）

将\({{\eta }_{P}}\)代入（8-3）式，即可计算长尾喷管发动机内弹道。

（3）推力终止装置内弹道计算

一般弹道导弹固体火箭发动机，要求具有推力终止功能。目前得到广泛使用的推力终止方案多采用反向喷管推力终止装置。这种方案主要是在发动机壳体上，对称地开若干个孔，在孔的位置上安装与发动机轴线成某一角度的反向喷管。其功能是：当推力终止时，由打开机构打开反向喷管，燃气流经反向喷管排出时建立反推力，来抵消发动机主推力。燃烧室压强按质量守恒方程，考虑反向喷管的喉部面积，反向推力按多个斜切喷管计算，发动机总推力为主喷管推力减去反向推力。详细可参见第四章4.10节。

这里定义了一个反向喷管打开时间t_r，图8-1（b）、（c）显示了喷管打开后内弹道特征：

（a）单纯主喷管工作内弹道

（b）反向喷管打开后压强曲线               

（c）反向喷管打开后推力曲线

图8-1  反向喷管打开时压强和推力曲线

由计算结果可以看出，反向喷管打开时反向推力瞬间较大，逐渐向平衡过度。此时压强曲线和推力曲线并不一致。

（4）推进剂数据

内弹道计算首先必须核实推进剂数据，燃速数据、质量特性数据和热力学参数，对内弹道性能有不同程度的影响。

1）燃速数据

燃速数据常见的有燃速系数a、压强指数n、燃速温度敏感系数\({{\sigma }_{P}}\)。如果燃烧室压强变

化较大，a和n通常不是一个固定的参数，这里分为三种情况：

A）a和n是固定，由试验测定并提供给设计。

B）指定压强指数n，然后给出推进剂在标准压强P_c下的燃速，此时，a可由下式计算得出：\(a={{r}_{b}}/P_{C}^{n}\)。

C）如果压强变化较大，有时提供一系列燃速压强数据表\(\left( {{P}_{i}},{{r}_{i}} \right),i=1,…,{{n}_{b}}\)，共有n_b组，计算之前首先拟合出n_b-1组燃速系数和压力指数，供内弹道计算用。可采用下式：

\({{n}_{i}}=\frac{\ln \left( {{{r}_{i}}}/{{{r}_{i+1}}}\; \right)}{\ln \left( {{{P}_{i}}}/{{{P}_{i+1}}}\; \right)}\).  

\({{a}_{i}}=\frac{{{r}_{i}}}{{{P}^{{{n}_{i}}}}}\)    \(i=1,…,{{n}_{b}}-1\)                             （8-36）

在内弹道计算时，根据计算的压力选择燃速系数，然后根据选定的燃速系数再计算压力，直至相邻两次的计算压强满足给定的误差为止。

                     

D）有时直接给出不同压强段的燃速公式，首先指定压强，选择燃速公式，计算工作压强，与指定值比较然后迭代计算来使用燃速公式。

E）燃速温度敏感系数\({{\sigma }_{P}}\)，一般发动机都是在常温下设计，然后验算高低温。推进剂所给出的燃速通常也是常温的T₀，高低温下的燃速数据使用下式计算：\({{r}_{b}}={{r}_{b0}}{{e}^{{{\sigma }_{P}}\left( T-{{T}_{0}} \right)}}\)，对于战术武器和战略武器高低温范围有所不同。推进剂燃速温度敏感系数通常是在0.15mm/K数量级。

F）由于工作条件的差异可能会导致实际燃速与实验发动机燃速不一致使用时，对推进剂燃速必须加以修正。

2）热力学参数

热力学参数是由推进剂生产厂家根据配方进行计算和试验获得的燃烧气体发的特性参数，主要是：燃温T_g，气体常数R_g，比热比k，特征速度C^*和凝相含量，热力学计算给出的是定压下的热力学数据，随着燃烧室压强的变化，这些参数应当是不固定的，但是，计算时一般都作常数处理。另外对于高低温，推进剂的焓值是不同的，特征速度也是变化的，准确的计算应当考虑温度的影响。

3）装药数据，装药数据为装药设计计算出的燃面肉厚数据表，计算时采用插值计算。

（5）喷管数据

喷管数据分为三类：常规喷管、长尾喷管和斜切喷管。

• 常规喷管：
• 喷管数量n：单位 无；
• 初始喉径 R_t0：单位 m；
• 开始烧蚀时间 S_t0：单位 s；
• 烧蚀斜率 T_k：单位 m/K，如果喉部不烧蚀T_k=0；
• 喷管效率 h_F：考虑到各种损失对推力系数加以修正，h_F<1，效率的问题在总体设计一章中已经讨论过
• 喷管面积膨胀比e_A；
• 喷管安装角\(\gamma \)：喷管轴线与发动机轴线的夹角；
• 长尾喷管：

与常规喷管相比喷管多了如下三个数据：

• 尾管长度Lc： 长尾喷管圆柱段长度，单位mm；
• 尾管半径Rc： 长尾喷管圆柱段半径，单位mm；
• 尾管摩擦系数f： 长尾喷管圆柱段摩擦系数。

对于非长尾管喷管入口段总压恢复系数\({{\eta }_{P}}=1\)，总压恢复系数是由喷管参数计算出来的。长尾管越长总压恢复系数越小。

• 斜切喷管：
• 喉面收缩系数Cd： 由于径缩引起实际流通面积减小；
• 斜切角\(\varphi \)： 喷管出口平面与喷管轴线的夹角；
• 喷管轴线倾角： 与喷管安装角同。
• 推力终止喷管数据：

对于大型发动机有推力终止装置的发动机，反向喷管一般是斜切喷管：

• 反向喷管打开时间t_r： 在t_r之后，主喷管与终止喷管同时工作，推力方向相反。

5）其它数据

• 工作终点压强 P_af 单位  MPa，一般可取 P_af£3Mpa；
• 有效工作压强 P_ef 单位  MPa, 一般可取 P_ef£7Mpa；
• 参考温度 T₀ 给出推进剂参数的温度，通常是设计温度，单位K；
• 工作温度 T 装药的计算初温，单位K；
• 点火压强 P_ig 喷管堵盖破裂压强 MPa；
• 燃烧效率 h_c  考虑到燃烧室热损失所引起的压强降低 h_c£0；
• 如果P_c<P_af，拖尾段结束；
• 工作高度 H 单位 m, 通过调用大气计算程序，计算环境压强 P_a，对于导弹飞行过程，工作高度是变化的，因此以地面试验和计算的内弹道数据必须加以修正；
• 体积装填分数V_l，根据装药的初始体积计算燃烧室初始自由容积。

计算结果应为以下参量随时间的变化：

装药肉厚（m）；压强（MPa）；推力（kN）；工作时间（s）；压强冲量（Mpa.s）；推力冲量（kN·s）；推力系数；喉部面积（m²）；质量流率（kg/s）。

图8-2  典型内弹道曲线

图8-2示为某发动机在-40℃，+20℃，+50℃三个工作温度下的典型内弹道。

8.2.1  双燃速计算内弹道数学模型

（1）简化热力计算模型

1）符号表

A_t                 喷喉面积

H_TN          燃料的焓

a_i                 线性燃速常数

k               绝热指数

b_i              线性燃速系数

L              装药长度

B（k）=\(\sqrt{k}{{\left( \frac{2}{k+1} \right)}^{\frac{k+1}{2(k-1)}}}\)

M             燃气分子量

C_D                质量流率系数

\(\overset{.}{\mathop{m}}\,\)          每秒质量排出量

C_q               流量损失系数

n                 压力指数

C_p                  定压比热

p                 燃烧室压强

C_v                   定容比热

Q                热量

D                装药外径

R                燃气气体常数

D_k                  燃烧室内径

A_b               燃面面积

d                 装药内径

T                燃气温度

d_t                    喉部直径

T_p               定压燃烧温度

E                温度T时的内能

T_v                   定容燃烧温度

E_N              燃料的比内能

t                  时间

e                  燃掉肉厚

r                  燃速

f₀                折算火药力

a₁               指数燃速系数

f₁                定容火药力

V                燃烧室自由容积

H_T              温度T时燃气的焓

\({{\rho }_{p}}\)               火药密度

\(\rho \)              燃气密度

i=1,2         表示两种火药或两级参数

eq               下标平衡

b        下标表示燃面

p        下标表示定压

2）假设和基本关系式

1. 燃烧室 2. 点火装置 3. 第一种装药 4. 第二种装药 5. 喷管

图8-3表示出了单室双药双推力发动机的基本结构示图。两种火药燃烧后，形成的各自的燃烧产物，还可能进行某些化学热力反应，形成动态平衡的燃气，由喷管排出。这个过程能用热力计算来完成。为了简化且不影响问题的解，提出以下几点假设，避免热力计算过程，得到问题的解。

• 假设两种火药燃烧后生成的燃气，在燃烧室内不再进行化学热力反应而是以各自的燃气独立存在，组成混合燃气，由喷管排除。
• 燃气在喷管中流动呈冻结平衡状态。
• 混合燃气满足理想气体状态方程。
• 系统是绝热的。

根据假设（1），利用混合气体平衡参数的计算方法，可以用两种火药燃气的参数，求得混合燃气的有关参数。

这些参数是：

M=\(\frac{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{M}_{i}}{{A}_{bi}}}}{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{A}_{bi}}}}\)                       （8-39）

R=\(\frac{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{R}_{i}}{{A}_{bi}}}}{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{A}_{bi}}}}\)                       （8-40）

\({{C}_{p}}=\frac{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{C}_{pi}}{{A}_{bi}}}}{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{A}_{bi}}}}\)                       （8-41）

\({{C}_{v}}=\frac{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{C}_{vi}}{{A}_{bi}}}}{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{A}_{bi}}}}\)                       （8-42）

k=\({}^{{{C}_{p}}}/{}_{{{C}_{v}}}\)                           （8-43）

T_p=\(\frac{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{T}_{pi}}{{A}_{bi}}}}{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{A}_{bi}}}}\)                       （8-44）

\({{T}_{v}}=\frac{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}}{{T}_{v}}{{A}_{bi}}}{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{A}_{bi}}}}\)                       （8-45）

\({{f}_{0}}=\frac{\sum\limits_{i=0}^{2}{{{\rho }_{pi}}{{r}_{i}}{{f}_{0i}}{{A}_{bi}}}}{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{A}_{bi}}}}\)                       （8-46）

\({{f}_{1}}=\frac{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}}{{f}_{1i}}{{A}_{bi}}}{\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{A}_{bi}}}}\)                        （8-47）

E=C_vT                             （8-48）

E_N=C_vT_v                                        （8-49）

H_T=C_pT=E_T+RT                        （8-50）

H_TN=C_pT_p=E_N+RT_p                      （8-51）

3）内弹道计算方程

在得到燃烧室中混合燃气的热力参数后，根据质量平衡和能量关系，可以得到下列基本方程式。

连续方程：

            \(\frac{d}{dt}\left( \frac{PV}{RT} \right)=\sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}}-\overset{.}{\mathop{m}}\,\overset{\bullet }{\mathop m}\,\)                 （8-52）

       \(\overset{.}{\mathop{m}}\,=\frac{{{C}_{D}}\cdot B(k)}{\sqrt{RT}}{{A}_{t}}\cdot P\)                   （8-53）

展开式（8-52）左端并引用\(\frac{dV}{dt}=\sum\limits_{i=1}^{2}{{{A}_{bi}}}{{r}_{i}}\)关系，经推导整理得到：

V\(\cdot \frac{dP}{dt}=\sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}\cdot RT-\frac{{{C}_{q}}B(k)}{\sqrt{RT}}{{A}_{t}}\cdot P\cdot RT}-P\sum\limits_{i=1}^{2}{{{A}_{bi}}{{r}_{i}}+\frac{PV}{(RT)}\cdot \frac{d}{dt}}(RT)\)

V\(\cdot \frac{dP}{dt}=\sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}\cdot RT-\frac{{{C}_{q}}B(k)}{\sqrt{RT}}{{A}_{t}}\cdot P\cdot RT}-P\sum\limits_{i=1}^{2}{{{A}_{bi}}{{r}_{i}}+\frac{PV}{(RT)}\cdot \frac{d}{dt}}(RT)\)                                    （8-54）

             \(\frac{d(RT)}{dt}=T\cdot \frac{dR}{dt}+R\cdot \frac{dT}{dt}\)                   （8-55）

由式（8-52）可以得到：

\(\frac{dR}{dt}=\frac{{{\rho }_{pi}}\cdot {{a}_{1.1}}\cdot {{A}_{b1}}{{\rho }_{p2}}\cdot {{a}_{1.2}}\cdot {{A}_{b2}}({{R}_{2}}-{{R}_{1}})({{n}_{2}}-{{n}_{1}})\cdot {{P}^{{{n}_{1}}+{{n}_{2}}-1}}}{{{\left( \sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{A}_{bi}}{{r}_{i}}} \right)}^{2}}}\times \frac{dp}{dt}={{F}_{R}}\cdot {{P}^{{{n}_{1}}+{{n}_{2}}-1}}\cdot \frac{dP}{dt}\)

                                                  （8-56）

式中

\({{F}_{R}}=\frac{{{\rho }_{p1}}\cdot {{a}_{1.1}}{{A}_{b1}}\cdot {{\rho }_{p2}}\cdot {{a}_{1.2}}\cdot {{A}_{b2}}\cdot ({{R}_{2}}-{{R}_{1}})\cdot ({{n}_{2}}-{{n}_{1}})}{{{\left( \sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{A}_{bi}}{{r}_{i}}} \right)}^{2}}}\)

      （8-57）

式（8-45）和（8-56）代入（8-54）式中得到：

      \(V\left[ 1-\frac{{{F}_{R}}\cdot {{P}^{{{n}_{1}}+{{n}_{2}}}}}{R} \right]\cdot \frac{dp}{dt}=\)

 \(\sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}\cdot RT-\frac{{{C}_{q}}\cdot B(k)}{\sqrt{RT}}}{{A}_{t}}\cdot P\cdot RT-P\sum\limits_{i=1}^{2}{{{A}_{bi}}{{r}_{i}}+\frac{PV}{T}}\cdot \frac{dT}{dt}\)

                 （8-58）

式中\(\frac{dT}{dt}\)项还须由能量方程得到。

能量平衡方程：

         \(\frac{d}{dt}\left( \frac{VP}{RT}E \right) ={{E}_{N}}\cdot \sum\limits_{i=1}^{2}{{{A}_{bi}}\cdot {{\rho }_{pi}}{{r}_{i}}}-\overset{\bullet }{\mathop m}\,\overset{.}{\mathop{m}}\,{{H}_{T}}-\frac{dQ}{dt}\)             （8-59）

展开式（8-59）左端得到：

         \(\frac{d}{dt}\left( \frac{VP}{RT}E \right)＝E\cdot \frac{d}{dt}\left( \frac{VP}{RT} \right)+\frac{VP}{RT}\cdot \frac{dE}{dt}\)               （8-60）

将式（8-52）两端乘以E得到：

         E\(\cdot \frac{d}{dt}\left( \frac{VP}{RT} \right)＝E\cdot \sum\limits_{i=1}^{2}{{{A}_{bi}}\cdot {{\rho }_{pi}}{{r}_{i}}}－E\cdot m\)             （8-61）

综合（8-61）（8-60）（8-59）三式得到：

         \(\frac{VP}{RT}\frac{dE}{dt}=\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{A}_{bi}}{{r}_{i}}\cdot {{C}_{v}}\cdot ({{T}_{v}}-T)-\overset{.}{\mathop{m}}\,RT-\frac{dQ}{dt}}\overset{\bullet }{\mathop m}\,\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{A}_{bi}}{{r}_{i}}\cdot {{C}_{v}}\cdot ({{T}_{v}}-T)-\overset{.}{\mathop{m}}\,RT-\frac{dQ}{dt}}\)         （8-62）

式中

         \(\frac{dE}{dt}={{C}_{v}}\cdot \frac{dT}{dt}+T\cdot \frac{d{{C}_{v}}}{dt}\)                     （8-63）

由式（8-54）可以得到：

         \(\frac{d{{C}_{v}}}{dt}=\frac{{{P}^{{{n}_{1}}+{{n}_{2}}-1}}}{(\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{A}_{bi}}{{r}_{i}}{{)}^{2}}}}\cdot {{A}_{bi}}\cdot {{\rho }_{p1}}{{a}_{1.1}}{{\rho }_{p2}}\cdot {{A}_{b2}}\cdot {{a}_{1.2}}({{n}_{1}}-{{n}_{2}})\)

                \(\cdot ({{C}_{v1}}-{{C}_{v2}})\cdot \frac{dp}{dt}={{F}_{c}}\cdot {{P}^{{{n}_{1}}+{{n}_{2}}-1}}\cdot \frac{dp}{dt}\)                  （8-64）

式中

           \({{F}_{c}}=\frac{{{A}_{b1}}{{\rho }_{p1}}\cdot {{a}_{1.1}}\cdot {{A}_{b2}}{{\rho }_{p2}}{{a}_{1.2}}({{n}_{1}}-{{n}_{2}})({{C}_{v1}}-{{C}_{v2}})}{{{\left( \sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{A}_{bi}}{{r}_{i}}} \right)}^{2}}}\)          （8-65）

将式（8-63）、（8-64）代入到（8-62）式中，两端乘以R/C_v，代入R/C_v=k-1关系，整理得到：

\(\frac{VP}{T}\cdot \frac{dT}{dt}+\frac{V{{P}^{{{n}_{1}}+{{n}_{2}}}}}{{{C}_{v}}}\cdot {{F}_{c}}\cdot \frac{dP}{dt}\)

\(=R\cdot ({{T}_{v}}-T)\sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}-(k-1)\overset{.}{\mathop{m}}\,RT-\frac{R}{{{C}_{v}}}\cdot \frac{dQ}{dt}}\frac{VP}{T}\cdot \frac{dT}{dt}+\frac{V{{P}^{{{n}_{1}}+{{n}_{2}}}}}{{{C}_{v}}}\cdot {{F}_{c}}\cdot \frac{dP}{dt}\)

\(=R\cdot ({{T}_{v}}-T)\sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}-(k-1)\overset{.}{\mathop{m}}\,RT-\frac{R}{{{C}_{v}}}\cdot \frac{dQ}{dt}}\) （8-66）

将（8-66）式中的\(\frac{dT}{dt}\)代入式（8-58）中并考虑取热损失为常量、热损失系数为常数x。

\(\frac{R}{{{C}_{v}}}\cdot \frac{dQ}{dt}=(1-x)R{{T}_{v}}\cdot \sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}}\)                （8-67）

最后得到连续方程和能量方程：

\(V\left( 1\frac{{{P}^{{{n}_{1}}+{{n}_{2}}}}}{{{C}_{v}}}\cdot {{F}_{c}}-\frac{{{P}^{{{n}_{1}}+{{n}_{2}}}}}{R}\cdot {{F}_{R}} \right)\cdot \frac{dP}{dt}\)

\(=x\cdot R{{T}_{v}}\sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}-kRT\cdot \frac{{{C}_{q}}B(k)}{\sqrt{RT}}{{A}_{t}}\cdot P-P\sum\limits_{i=1}^{2}{{{A}_{bi}}{{r}_{i}}}}\)

                   （8-68）

\(\frac{VP}{T}\cdot \frac{dT}{dt}=R(x{{T}_{v}}-T)\sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}-(k-1)\frac{{{C}_{q}}B(k)}{\sqrt{RT}}\cdot {{A}_{t}}\cdot PRT}\)

\(-\frac{{{P}^{{{n}_{1}}+{{n}_{2}}}}}{{{C}_{v}}\left( 1+\frac{{{P}^{{{n}_{1}}+{{n}_{2}}}}}{{{C}_{v}}}\cdot {{F}_{c}}-\frac{{{P}^{{{n}_{1}}+{{n}_{2}}}}}{R}{{F}_{R}} \right)}(xR{{T}_{v}}\sum\limits_{i=1}^{2}{{{A}_{bi}}{{\rho }_{pi}}{{r}_{i}}-kRT\frac{{{C}_{q}}B(k)}{\sqrt{RT}}{{A}_{t}}P-P\sum\limits_{i=1}^{2}{{{A}_{bi}}{{r}_{i}}}}\)

                   （8-69）

式（8-68）、（8-69）就是单室双药双推力发动机的内弹道主要方程。

燃速方程：

         \(\frac{d{{e}_{i}}}{dt}={{a}_{1i}}{{P}^{ni}}\)                           （8-70）

或

         \(\frac{d{{e}_{i}}}{dt}={{a}_{i}}+{{b}_{i}}P\)                          （8-71）

燃面的变化，根据具体药型可以表示为肉厚变化的函数。较普遍的是等面变化规律。

4）等面燃烧情况下的平衡压力计算公式

等面燃烧情况具有普遍的意义，在这个条件下能得到两级平衡压力的计算公式。第二级，即剩余单一火药装药燃烧时的平衡压力，计算公式是大家熟悉的：

\({{P}_{eq1}}={{\left( \frac{{{\rho }_{p1}}{{a}_{1.1}}{{A}_{b1}}\sqrt{x{{R}_{1}}\cdot {{T}_{pi}}}}{B(k){{C}_{q}}{{A}_{t}}} \right)}^{\frac{1}{1-{{n}_{1}}}}}\)                   （8-72）

第一级，即两种火药同时燃烧时的平衡压力计算，可以由\(\frac{dp}{dt}=\frac{dT}{dt}=0\)

及T=xT_p条件建立的质量和能量平衡的关系得到：

                \(\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}{{A}_{bi}}=\frac{{{C}_{q}}\cdot B(k){{A}_{t}}\cdot P}{\sqrt{x\cdot R{{T}_{p}}}}={{C}_{D}}\cdot {{A}_{t}}P}\)            （8-73）

和

\(\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{r}_{i}}R{{T}_{p}}=\sum\limits_{i=1}^{2}{{{\rho }_{pi}}{{A}_{bi}}{{r}_{i}}{{R}_{i}}{{T}_{pi}}}}\)                （8-74）

利用式（8-71）线性燃烧经验关系，并取k=\[\frac{{{k}_{1}}+{{k}_{2}}}{2}\]，由（8-73）和（8-74）两式联立求解得到平衡压力的解析式。

  \({{P}_{eq2}}=\frac{\left( \overline{b}\cdot \overline{c}+\overline{a}\cdot \overline{d} \right)+\sqrt{{{\left( \overline{b}\cdot \overline{c}+\overline{a}\cdot \overline{d} \right)}^{2}}+4\overline{a}\overline{c}\left[ {{\left( \frac{{{C}_{D1}}}{{{K}_{N1}}} \right)}^{2}}-\overline{b}\cdot \overline{d} \right]}}{2\times \left[ {{\left( \frac{{{C}_{D1}}}{{{K}_{N1}}} \right)}^{2}}-\overline{b}\overline{d} \right]}\)      （8-75）

式中：

      \(\overline{a}=\overline{{{a}_{2}}}\left( \frac{\overline{{{a}_{1}}}}{\overline{{{a}_{2}}}}+\frac{{{A}_{b2}}}{{{A}_{b1}}} \right)\)                 

\(\overline{b}=\overline{{{b}_{2}}}\left( \frac{\overline{{{b}_{1}}}}{\overline{{{b}_{2}}}}+\frac{{{A}_{b2}}}{{{A}_{b1}}} \right)\)

  \(\overline{c}=\overline{{{a}_{2}}}\left[ \frac{\overline{{{a}_{1}}}}{\overline{{{a}_{2}}}}+\frac{{{A}_{b2}}}{{{A}_{b1}}}{{\left( \frac{{{C}_{D1}}\cdot B({{k}_{2}})}{{{C}_{D2}}B({{k}_{1}})} \right)}^{2}} \right]\)

  \(\overline{d}=\overline{{{b}_{2}}}\left[ \frac{\overline{{{b}_{1}}}}{\overline{{{b}_{2}}}}+\frac{{{A}_{b2}}}{{{A}_{b1}}}{{\left( \frac{{{C}_{D1}}\cdot B({{k}_{2}})}{{{C}_{D2}}B({{k}_{1}})} \right)}^{2}} \right]\)

  \(\overline{{{a}_{i}}}={{a}_{i}}\cdot {{\rho }_{pi}}\)            

  \(\overline{{{b}_{i}}}={{b}_{i}}{{\rho }_{pi}}\)

  \({{C}_{Di}}=\frac{{{C}_{q}}B(k)}{\sqrt{x{{R}_{i}}{{T}_{pi}}}}\)                  \({{K}_{N1}}=\frac{{{A}_{b1}}}{{{A}_{t}}}\)

由式（8-71）和（8-75）看到平衡压力是燃面比\({{{A}_{b2}}}/{{{A}_{b1}}}\;\)，第一级面喉比K_N1的函数。调节这两个参数就可以得到所需要的P_eq2与P_eq1。 

（2）重量成分法计算模型

1）固体火箭发动机实现双推力的方案

根据固体发动机的推力公式：

F=CFA_tP_c                            （8-76）

其中: 

\({{P}_{c}}={{({{C}^{*}}\cdot \rho \cdot a\cdot {{{A}_{b}}}/{{{A}_{t}}}\;)}^{\frac{1}{1-n}}}\)                      （8-77）

A_t喷喉面积。

       \({{C}_{f}}=\sqrt{\frac{2k}{k-1}\left[ {{\left( 1-{{{p}_{e}}}/{{{p}_{c}}}\; \right)}^{\frac{k-1}{k}}} \right]}+\frac{{{A}_{e}}}{{{A}_{t}}}\cdot \frac{{{p}_{e}}-{{p}_{a}}}{{{p}_{c}}}\)           （8-78）

将（8-77）式代入（8-76）

\(F={{C}_{f}}{{({{C}^{*}}\cdot \rho \cdot a\cdot {{A}_{b}})}^{\frac{1}{1-n}}}\cdot {{A}_{t}}\cdot {{\left( \frac{1}{{{A}_{t}}} \right)}^{\frac{1}{1-n}}}\)

 \(\overset{{}}{\mathop {}}\,={{C}_{f}}{{({{C}^{*}}\cdot \rho \cdot a\cdot {{A}_{b}})}^{\frac{1}{1-n}}}\cdot A_{t}^{\frac{-n}{1-n}} \)

            （8-79）

从（8-78）可知推力系数C_f主要取决于推进剂的比热比和喷管膨胀比，各种推进剂的比热比变化不大，所以对推力系数的影响也不大。而主要影响推力系数的参数是膨胀比，要改变膨胀比的主要办法是改变A_t，要在发动机工作过程中改变A_t，喷管设计很复杂，而也大大增加了喷管的设计质量，因而大多数小型固体火箭发动机，要改变推力，不采用这种办法。这时推力是下列因素的函数：

         \(F={{C}_{f}}\left( {{C}^{*}},\rho ,a,{{A}_{b}} \right)\)                       （8-80）

从（8-80）式可知，只要改变式中任何一个因素，就可以得到改变推力的目的。我们的发动机是采用两种不同燃速的推进剂和改变装药燃面的办法，获得双推力的。

推进剂的燃速公式一般采用r=apⁿ，要改变燃速就意味着要改变燃速系数a和压力指数n。一般希望n取小值，对发动机工作有利，因此改变燃速主要就是改变a。对于一种推进剂而言在同一工作压力下，使a发生变化是不可能的。因而简捷的办法是采用两种推进剂，他们有各自的性能。再辅之以变燃面设计，就可以获得双推力。

2）重量成分法实现单室双推力发动机的内弹道计算

由于在发动机中装入两种不同燃速的推进剂，因而在内弹道计算方面就不象单推力发动机那样简单。在内弹道计算之前要进行一些假设。对两种推进剂的热力参数要进行处理。主要有下列几条假设：

1. 两种推进剂的燃烧产物是物理混合的理想气体，他们之间存在热交换，但系统没有质量变化，总焓不变，是一个绝热等熵系统。
2. 燃烧产物在燃烧室中是一维定常流动。
3. 在燃烧过程中，热损失为一常数，用损失系数X考虑。
4. 装药全面点燃是压力达到点火压力P_ig时实现的。

• 上升段压强时间计算

   根据教科书的推导结果，单推力发动机上升段的计算公式为：

         \({dP}/{dt}\;=\frac{XRT}{{{V}_{0}}}\left( {{A}_{b}}\cdot r\cdot \rho -{{C}_{D}}{{A}_{t}}{{P}_{c}} \right)\)                 （8-81）

因为我们设计的双推力发动机，使用了两种推进剂，而且在助推段要同时燃烧，所以上

式中的R，T，r，ρ，C等都各不相同。那么怎样来确定这些参数呢？我们采用了重量成分

法，根据各自燃烧重量来确定他们在混合气体中的比例。按比例来确定混合气体的R，T，r，

ρ，C是推进剂的燃气生成量，如果是两种推进剂同时燃烧，则燃气生成量应是\({{A}_{b1}}{{r}_{1}}{{\rho }_{1}}+{{A}_{b2}}{{r}_{2}}{{\rho }_{2}}\)。这时（8-81）式修改为：

\({dP}/{dt}\;=\frac{XRT}{{{V}_{0}}}\left( {{A}_{b1}}\cdot {{r}_{1}}\cdot {{\rho }_{1}}{{A}_{b2}}{{r}_{2}}{{\rho }_{2}}-{{C}_{D}}{{A}_{t}}{{P}_{c}} \right)\)           （8-82）

要求的混合气体的R、T、C_D，首先得求两种推进剂燃气生成量的重量成分和比容。

第一种推进剂燃气的重量成分：

      \({{q}_{1}}={{{A}_{b1}}\cdot {{r}_{1}}\cdot {{\rho }_{1}}}/{\left( {{A}_{b1}}\cdot {{r}_{1}}\cdot {{\rho }_{1}}+{{A}_{b2}}\cdot {{r}_{2}}\cdot {{\rho }_{2}} \right)}\;\)               （8-83）

第二种推进剂燃气的重量成分：

     \({{q}_{2}}={{{A}_{b2}}\cdot {{r}_{2}}\cdot {{\rho }_{2}}}/{\left( {{A}_{b1}}\cdot {{r}_{1}}\cdot {{\rho }_{1}}+{{A}_{b2}}\cdot {{r}_{2}}\cdot {{\rho }_{2}} \right)}\;\)              （8-84）

其中：   A_b1，A_b2为两种推进剂的燃面，由两种推进剂的药型而定。

        r₁,r₂为两种推进剂的燃速。

        ρ₂,ρ₂为两种推进剂的密度。

        \({{r}_{1}}={{a}_{1}}p_{c}^{{{n}_{1}}}\)                           （8-85）

        \({{r}_{2}}={{a}_{2}}p_{c}^{{{n}_{2}}}\)                           （8-86）

其中：a₁,a₂两种推进剂的燃速常数。

n₁,n₂两种推进剂的燃速指数。

第一种推进剂燃气的比容

     \({{v}_{1}}=\frac{{{{q}_{1}}}/{{{M}_{1}}}\;}{{{{q}_{1}}}/{{{M}_{1}}}\;+{{{q}_{2}}}/{{{M}_{2}}}\;}\)                      （8-87）

第二种推进剂燃气的比容

     \({{v}_{2}}=\frac{{{{q}_{2}}}/{{{M}_{2}}}\;}{{{{q}_{1}}}/{{{M}_{1}}}\;+{{{q}_{2}}}/{{{M}_{2}}}\;}\)                      （8-88）

其中M₁、M₂两种推进剂燃气的平均分子量，两种推进剂燃气的折合分子量M

     \(M={{M}_{1}}{{\upsilon }_{1}}+{{M}_{2}}{{\upsilon }_{2}}\)                        （8-89）

两种推进剂燃烧的混合气体的气体常数R

     \(R={{{R}_{0}}}/{M}\;\)                           （8-90）

式中R₀理想气体常数，混合燃气的折合分子量。

两种推进剂燃气的混合温度T

     \(T=\frac{{{T}_{1}}{{C}_{p1}}{{q}_{1}}+{{T}_{2}}{{C}_{p2}}{{q}_{2}}}{{{C}_{p1}}{{q}_{1}}+{{C}_{p2}}{{q}_{2}}}\)                      （8-91）

其中T₁,T₂两种推进剂燃气的定压燃温。C_p1,C_p2两种推进剂燃气的定压比热。

V₀   发动机点火之前的自由容积，当常量处理。

A_t   为喷管喉部面积。

P_c   发动机燃烧室压力。

C_D   两种推进剂燃烧的混合气体流量系数。

        \({{C}_{D}}=\frac{\Gamma }{\sqrt{R\cdot T}}\)                         （8-92）

R，T可按照（8-91），（8-92）求得。

        \(\Gamma =\sqrt{k}{{\left( \frac{2}{k+1} \right)}^{\frac{\left( k+1 \right)}{2\left( k-1 \right)}}}\)                      （8-93）

k为两种推进剂燃烧的混合气体比热比。

\(k=\frac{{{C}_{D}}}{{{C}_{v}}}=\frac{1}{1-{1.99}/{M{{C}_{p}}}\;}\)                    （8-94）

M为两种推进剂燃烧的混合气体的折合分子量。

C_p为两种推进剂燃烧的混合气体的定压比热。

        \({{C}_{p}}={{C}_{p1}}\cdot {{q}_{1}}+{{C}_{p2}}\cdot {{q}_{2}}\)                     （8-95）

把上述推导出的各参数代入（8-82）式，就可以求得上升段的压力和时间。计算时可按变燃面计算。

• 第一次（助推级）平衡段压强计算。

一种推进剂的计算公式为

\({{A}_{b}}\cdot r\cdot \rho ={{C}_{D}}\cdot {{A}_{t}}\cdot {{P}_{c}}\)                        （8-96）

两种推进剂同时燃烧时，要对（2-56）式进行修改，主要是改变燃气生成量公式，流量公式不变，现变化如下：

        \({{A}_{b1}}\cdot {{r}_{1}}\cdot {{\rho }_{1}}+{{A}_{b2}}\cdot {{r}_{2}}\cdot {{\rho }_{2}}={{C}_{D}}\cdot {{A}_{t}}\cdot {{P}_{c}}\)                 （8-97）

式中燃面A_b1,A_b2随肉厚变化，密度ρ₂,ρ₂是不变的，燃速r₁,r₂随计算压力变化。C_D是R,T,k的函数，按重量成分计算，计算方法前边已推导出公式。

• 第一次压强下降段的计算

当第一次压强下降时，我们认为助推级装药已工作完毕，并不留余药。此时的自由容积即点火初期的自由容积与已燃烧掉的装药所占容积之和。为了计算方便，认为此时的自由容积是一个常数，故可利用下式计算

\({dp}/{dt}\;=\frac{XRT}{V}\left( {{A}_{b2}}\cdot {{r}_{2}}\cdot {{\rho }_{2}}-{{C}_{D}}{{A}_{t}}{{P}_{c}} \right)\)                 （8-98）

式中R,T,C_D都采用助推级燃烧完毕时续航级的瞬时混合值。和当时的重量成分有关。V取当时的瞬时值.

• 第二次压力平衡段计算

此时的计算公式与单推力发动机的平衡计算公式一样

         \({{A}_{b2}}\cdot {{r}_{2}}\cdot {{\rho }_{2}}={{C}_{D}}\cdot {{A}_{t}}\cdot {{P}_{c}}\)                      （8-99）

\({{P}_{c}}={{\left( {{{A}_{b2}}{{\rho }_{2}}\cdot {{a}_{2}}C_{2}^{*}}/{{{A}_{t}}}\; \right)}^{\frac{1}{1-{{n}_{2}}}}}\)

燃面A_b2是第一级燃烧结束时，第二级的燃面。

• 二级下降段压强－时间曲线的计算

假设设计的第二级药型是锥管，则燃烧结束时不再存在余药。燃气在燃烧室中是纯膨胀过程。推导结果是：

\(P{{P}_{ce}}-\frac{R{{T}_{0}}\times {{C}_{D}}{{A}_{t}}}{{{V}_{c}}}\left( t-{{t}_{c}} \right)\)                  （8-100）

R,T,C_D是二级火药的参数。

V_c是装药烧完后的发动机自由容积。

t_c是二级装药烧完时的时间。

P_ce第二次平衡压力。

P,t是第二次下降段的压强与对应的工作时间。

8.2.3  内弹道计算方法

（1）计算公式

上面给出了双燃速混合气体的热力参数处理的数学物理模型，本节将要介绍程序编制过程的处理方法和计算方法。尽管从实际上整个燃烧过程可以分为一级压强建立过程、一级平衡阶段、一级压强下降过程、二级压强平衡过程和二级压强下降过程这5个阶段。但是，本文分为两个阶段处理，即：一级压强建立过程和以后的整个工作过程。

压强建立过程从喷管堵盖破裂压强P_ig到平衡压强P_c,eq，采用公式（8-11）计算。

后续的整个工作过程采用公式（8-101）（8-102）（8-103）计算：

\({{\dot{m}}_{1}}={{K}_{a1}}{{a}_{1}}P_{C}^{n1}{{A}_{b1}}{{\rho }_{P1}}\exp \left[ {{\sigma }_{P1}}\left( T-T0 \right) \right]\)               （8-101）

\({{\dot{m}}_{2}}={{K}_{a2}}{{a}_{2}}P_{C}^{n2}{{A}_{b2}}{{\rho }_{P2}}\exp \left[ {{\sigma }_{P2}}\left( T-T0 \right) \right]\)              （8-102）

\(\frac{dPc}{dt}=\frac{{{\eta }_{c}}^{2}{{R}_{g}}{{T}_{g}}}{{{V}_{c}}}\left[ \dot{m}-\frac{{{P}_{c}}{{\eta }_{P}}{{A}_{t}}\varphi }{{{C}^{*}}{{\eta }_{c}}} \right]\times {{10}^{-6}}\)              （8-103）

式中：热力参数中带下标1的为缓燃推进剂，带下标2的为速燃推进剂，不带下标的为混合燃气参数

混合燃气参数在给定燃烧室压强时，采用（8-104）各式计算

   \({{{\dot{m}}}_{1}}={{K}_{a1}}{{a}_{1}}{{P}_{c}}^{n1}{{\rho }_{P1}}{{A}_{b1}}\exp \left[ {{\sigma }_{P1}}\left( T-{{T}_{0}} \right) \right]  \)

    \({{{\dot{m}}}_{2}}={{K}_{a2}}{{a}_{2}}{{P}_{c}}^{n2}{{\rho }_{P2}}{{A}_{b2}}\exp \left[ {{\sigma }_{P2}}\left( T-{{T}_{0}} \right) \right] \) 

    \({{{\dot{m}}}_{1}}+{{{\dot{m}}}_{2}}=\frac{{{\eta }_{P}}{{P}_{c}}{{A}_{t}}\varphi }{C*} \) 

    \(C*=\sqrt{\frac{C*_{1}^{2}{{{\dot{m}}}_{1}}+C*_{2}^{2}{{{\dot{m}}}_{2}}}{{{{\dot{m}}}_{1}}+{{{\dot{m}}}_{2}}}} \) 

    \({{C}_{P}}=\frac{{{C}_{P1}}{{{\dot{m}}}_{1}}+{{C}_{p2}}{{{\dot{m}}}_{2}}}{{{{\dot{m}}}_{1}}+{{{\dot{m}}}_{2}}} \) 

    \({{C}_{V}}=\frac{{{C}_{V1}}{{{\dot{m}}}_{1}}+{{C}_{V2}}{{{\dot{m}}}_{2}}}{{{{\dot{m}}}_{1}}+{{{\dot{m}}}_{2}}} \) 

    \(k={{{C}_{P}}}/{{{C}_{V}}}\; \) 

   \( {{R}_{g}}=\frac{{{R}_{g1}}{{{\dot{m}}}_{1}}+{{R}_{g2}}{{{\dot{m}}}_{2}}}{{{{\dot{m}}}_{1}}+{{{\dot{m}}}_{2}}} \) 

    \({{T}_{g}}=\frac{{{C}_{P1}}{{T}_{g2}}{{{\dot{m}}}_{1}}+{{C}_{P2}}{{T}_{g2}}{{{\dot{m}}}_{2}}}{{{C}_{P1}}{{{\dot{m}}}_{1}}+{{C}_{P2}}{{{\dot{m}}}_{2}}} \) 

          （8-104）

公式中，当压强小于P_af时全部工作过程结束。

燃烧室压强通过上式迭代求得。为了叙述问题方便，设求压强用一个函数来表示：

\({{P}_{c}}=f\left( {{K}_{a1}},{{K}_{a2}},{{a}_{1}},{{a}_{2}},{{n}_{1}},{{n}_{2}},C{{*}_{1}},C{{*}_{2}},{{\sigma }_{P1}},{{\sigma }_{P2}},{{A}_{b1}},{{A}_{b2}},V,T,{{T}_{0}},{{A}_{t}},\varphi ,{{\eta }_{P}} \right)\)

                                                                                                    （8-105）

在第一级，由于Ab₁>0,A_b2>0，因此\({{\dot{m}}_{1}}>0,{{\dot{m}}_{2}}>0\)，两种推进剂都参与燃烧，在第一阶段拖尾中，由于A_b2逐渐减小，速燃推进剂所占的分量就越来越小，直到第二级A_b2=0,\({{\dot{m}}_{2}}=0\)。（8-105）对一种推进剂也适用。由于（8-105）式是根据通用的连续方程导出的，因此也适用于排气过程。在整个工作过程中喉部面积At是增大的，燃烧室空腔容积V是增大的，由装药数据可以求出在任意肉厚时的空腔容积。

（2）装药数据

装药可分为缓燃和速燃两部分。装药的燃面数据由专门的装药设计软件来完成，对于串装装药分为两种情况，一种是先后燃，另一种是同时燃。

在同时燃的情况下，如果存在一个共同的交界面，在交界面上认为速燃装药按照原有的规律燃烧，而缓燃装药则生成一个圆台面（图8-4）。

图8-4  缓燃装药和速燃装药交界面处理

图中粗线为燃面位置，柱面的参数取决于两种推进剂的燃速比

\(\alpha {{\sin }^{-1}}\left( {{{r}_{1}}}/{{{r}_{2}}}\; \right)\)                       （8-106）

如果 r₁=r₂，则\(\alpha 90{}^\circ \)，可以认为无中间的圆台面。在装药设计软件里，在缓燃装药可以加一圆台体，且使\(\alpha \)按（8-118）式求得。缓燃药和速燃药分别进行设计，然后调入该内弹道程序进行计算。

如果两种装药之间是较为复杂的搭接面需要特殊处理，例如星形之间、翼形之间，把附加燃面找出来，然后各种进行计算。

如果不存在共同交界面，则各自燃烧。

在先后燃的情况下，需要给出二级药开始燃烧时第一级药所处肉厚，为了使两级过度圆滑，在装药计算时也应考虑内腔形体的构造，如图8-5所示。

由abcd构成的回转体为二级装药的初始形状，弧2为弧1当第一级结束时的位置，图8-5中D/d=r1/r2，所以\(d=D\frac{{{r}_{2}}}{{{r}_{1}}}\)，其中r₁/r₂为一二级燃速比。

图8-5中W（单位：m）为第二级开始燃烧时第一级的肉厚，在第一阶段只有第一级参与燃烧，在过渡阶段两级同时燃烧，在第二阶段只要第二级燃烧。

在过渡阶段型面是不规则的，我们可以用圆弧或直线近似逼近，只要不至于造成大的误差即可。要根据实际情况处理两级过渡。

图8-5  一级药与二级药交界面处理

一般我们把第一级作为速燃推进剂，也可能会有低燃速。但是一级推进剂总是在第一级燃烧，而二级推进剂则有可能跨过两个阶段。

对于套装装药一般很难计算，只有特殊情况设计的才能计算出来，如图8-6的星形装药，图中间距大者为速燃药，小者为缓燃药，交界面的形状取决于两种药的燃速比，可见采用两种药合理搭配，一可以增大初始燃气生成量，减少内孔燃烧初始面积小的缺点，二可以降低二级燃气生成量，降低星形装药带来的后期增面。同时也可以减小后期的拖尾，因此，套装装药合理设计可以获得理想的内弹道。但是这种装药需要两次浇注，而且设计难度大。

对于三燃速装药也可采用重量成分法处理混合燃气的参数，如果点火药燃气参数同时考虑，则为四种燃气参数，处理方法相同。

图8-6  星形套装装药

 

8.3  一维内弹道计算

对于长细比较大的发动机，不可避免地会出现侵蚀，关于侵蚀的机理可参见文献【2】有详细的介绍。考虑侵蚀对内弹道的影响必须进行一维内弹道计算。对于大型发动机尽量避免或减小侵蚀效应。侵蚀具有初始压力峰高和长拖尾等特点，侵蚀并不总是有害的，对于增面性较强的装药，有时利用侵蚀峰值来调平内弹道曲线。侵蚀一方面与推进剂有关，例如侵蚀燃烧燃速特性，另一方面与装药设计有关，如通气参量。推进剂侵蚀燃烧特性传统地采用中断燃烧试验，也可以根据发动机试验结果曲线采用参数辨识的方法确定侵蚀燃烧关系式。采用通气参量经验公式来修正侵蚀燃速关系式，则可以给出侵蚀的初始效应，但不能计算出内弹道的拖尾现象。本文在装药计算的基础上，结合喷管参数、推进剂数据来进行内弹道计算。

严格来讲，发动机内弹道是三维的，但是三维内弹道计算是一个复杂的数值计算，其燃烧的不均匀性，不能按照平行层燃烧来处理，作为工程计算，很不经济。这里采用简化的一维内弹道来进行计算，假定装药在任意截面燃速、参数分布是均匀的，不同截面可以不同，装药为内孔燃烧的加质流动，装药设计只需要提供每一截面不同肉厚下的燃烧周长和通气面积即可。

8.3.1  内弹道计算模型

（1）控制体控制方程法

压强建立过程从喷管堵盖破裂压强P_ig到平衡压强P_c,eq，采用公式（8-1）计算.

图8-7  控制体示意图

平衡段采用燃气在通道内流动的动量守恒和质量守恒方程计算。对于图8-7所示内通道，在1-1和2-2截面间的燃气，其质量守恒和动量守恒方程刻分别写成如下形式：

\({{\dot{m}}_{2}}={{\dot{m}}_{1}}+{{\left( {{{\dot{m}}}_{b}} \right)}_{1-2}}\)                      （8-107）

和

\({{\dot{m}}_{2}}{{v}_{2}}+{{p}_{2}}{{A}_{2}}={{\dot{m}}_{1}}{{v}_{1}}+{{p}_{1}}{{A}_{1}}+{{p}_{1}}\left( {{A}_{2}}-{{A}_{1}} \right)\)            （8-108）

其中\({{\left( {{{\dot{m}}}_{b}} \right)}_{1-2}}\)为1-1和2-2截面间的燃气生成量，A₁和A₂分别为通道两端的横截面积。由于：

\(\dot{m}v+pA=\dot{m}\left( v+\frac{p}{\rho v} \right)\)                    （8-109）

将公式\(v={{a}^{*}}\lambda \)

 

   \({{a}^{*2}}=\frac{2\gamma }{\gamma +1}R{{T}_{0}}  \)

   \(\frac{p}{\rho }=R{{T}_{0}}\tau \left( \lambda  \right)  \)

   \(\tau \left( \lambda  \right)=1-\frac{\gamma -1}{\gamma +1}{{\lambda }^{2}}  \)

代入（8-109）式，得：

\(\dot{m}v+pA=\frac{\gamma +1}{2\gamma }\dot{m}{{a}^{*}}\left( \lambda +\frac{1}{\lambda } \right)=\frac{\gamma +1}{2\gamma }\dot{m}{{a}^{*}}z\left( \lambda  \right)\)            （8-110）

将（8-110）式代入（8-108）式，整理后得：

\(\frac{{{{\dot{m}}}_{2}}}{{{{\dot{m}}}_{1}}}z\left( {{\lambda }_{2}} \right)=z\left( {{\lambda }_{1}} \right)+\frac{2{{p}_{1}}{{A}_{1}}\gamma \left( {{{A}_{2}}}/{{{A}_{1}}-1}\; \right)}{\left( \gamma +1 \right){{{\dot{m}}}_{1}}{{a}^{*}}}\)              （8-111）

燃气质量流率可写成：

     \(\dot{m}=\frac{2\gamma pA\lambda }{\left( \gamma +1 \right){{a}^{*}}\tau \left( \lambda  \right)}\)                         （8-112）

将（112）式代入（111）式得：

\(\frac{{{{\dot{m}}}_{2}}}{{{{\dot{m}}}_{1}}}z\left( {{\lambda }_{2}} \right)=z\left( {{\lambda }_{1}} \right)+\frac{\tau \left( {{\lambda }_{1}} \right)}{{{\lambda }_{1}}}\left( \frac{{{A}_{2}}}{{{A}_{1}}}-1 \right)\)               （8-113）

如果\({{\lambda }_{1}}>0\)且通道内有加质直接采用（8-113）式计算\({{\lambda }_{2}}\)。如果通道内有加质且\({{\lambda }_{1}}=0\)则采用下式计算\({{\lambda }_{2}}\)：

\(z\left( {{\lambda }_{2}} \right)=\frac{{{p}_{2}}{{A}_{2}}\gamma }{\left( \gamma +1 \right){{\left( {{{\dot{m}}}_{b}} \right)}_{1-2}}{{a}^{*}}}\)                   （8-114）

对于装药如果沿轴向划分成许多段，可以近似认为\({{p}_{2}}={{p}_{1}}\)

如果通道内无加质即\({{\dot{m}}_{2}}={{\dot{m}}_{1}}\)，且\({{\lambda }_{1}}>0\)，采用（8-115）式计算：

\(z\left( {{\lambda }_{2}} \right)=z\left( {{\lambda }_{1}} \right)+\frac{\tau \left( {{\lambda }_{1}} \right)}{{{\lambda }_{1}}}\left( \frac{{{A}_{2}}}{{{A}_{1}}}-1 \right)\)                 （8-115）

对于内通道，沿轴向各点的压强、温度、流速和燃速是不同的，基础燃速：

\({{r}_{0}}=a{{p}^{n}}\)                          （8-116）

如果考虑侵蚀燃烧，并采用格林公式，则任意点的燃速为：

     \(r={{r}_{0}}\varepsilon =a{{p}^{n}}\left[ 1+{{K}_{cr}}\left( \lambda -{{\lambda }_{cr}} \right) \right]\)                  （8-117）

\({{\lambda }_{cr}}\)为临界速度系数，如果\(\lambda \le {{\lambda }_{cr}}\)则ε=1.0。

对于装药一维流动，取装药末端气流的出口速度等于喷管的入口速度，由\(q\left( {{\lambda }_{2-2}} \right)=\frac{{{A}_{t}}}{{{A}_{i}}}=\frac{{{A}_{t}}}{{{A}_{p}}}\)计算装药的出口速度系数\({{\lambda }_{2-2}}\)，式中A_t为喉部面积，A_p为装药末端通气面积。

计算的方法是选定头部压强p_1-1，计算获得装药末端出口截面的速度系数\({{{\lambda }’}_{2-2}}\)，如果\({{{\lambda }’}_{2-2}}\)足够接近\({{\lambda }_{2-2}}\)，满足迭代的精度要求，本步计算结束，否则修正头部压强重新计算：

\({{p}_{1-1}}={{{p}’}_{1-1}}\frac{{{{{\lambda }’}}_{2-2}}}{{{\lambda }_{2-2}}}\)                        （8-118）

式中\({{{p}’}_{1-1}}\)为上一步头部压强。

如果全部无质量加入，则进入到拖尾段。拖尾段从有效工作压强P_ef下降到最小工作压强P_af。这一过程是余药的燃烧和发动机的排气过程并存，采用公式（8-119）计算：

\({{\left( \frac{{{P}_{c}}}{{{P}_{eq}}} \right)}^{1-n}}=1-\left[ 1-{{\left( \frac{{{P}_{ef}}}{{{P}_{eq}}} \right)}^{1-n}} \right]{{e}^{-\frac{t}{\tau }}}\)                  （8-119）

公式中，P_eq为余药产生的压强，\(\tau =\frac{1}{1-n}\frac{{{V}_{c}}}{{{\Gamma }^{2}}{{C}^{*}}{{A}_{t}}}\)，当压强小于P_af时全部工作过程结束。

（2）一维准定常内弹道计算法

平衡段压强计算也可以采用（8-120）式计算^【2】：

\(\frac{dP}{dx}=\frac{\rho {{V}^{2}}}{{{A}_{p}}}\frac{{{a}^{2}}}{{{a}^{2}}-{{V}^{2}}}\frac{d{{A}_{p}}}{dx}-\frac{{{\rho }_{p}}rV\Pi }{{{A}_{p}}}\left[ \frac{2{{a}^{2}}+\left( k-1 \right){{V}^{2}}}{{{a}^{2}}-{{V}^{2}}} \right] \)

\(\frac{dV}{dx}=\frac{{{\rho }_{p}}r\Pi }{\rho {{A}_{p}}}\frac{{{a}^{2}}+k{{V}^{2}}}{{{a}^{2}}-{{V}^{2}}}-\frac{{{a}^{2}}V}{{{A}_{p}}\left( {{a}^{2}}-{{V}^{2}} \right)}\frac{d{{A}_{P}}}{dx} \)

\({{C}_{p}}T+\frac{{{V}^{2}}}{2}={{C}_{p}}{{T}_{0}} \)

\(P=\rho RT \)

   （8-120）

 

8.3.2  内弹道计算方法

平恒段压强计算是求解一维变截面通道流动的压强和速度的分布方程组，似乎各点的气流参数仅是位置x的函数，与时间无关。但应该看到，方程组所包含的A_b、A_p和Π都是随时间在变化的。只是在解上述方程时，对应于某一时刻，将A_b、A_p和Π看成常数。到另一时刻，将A_b、A_p和Π却是另一组常数。由此可得到下一时刻下，相应的气流参数沿轴向分布。由此可见通道内各点的参数仍是x和t的函数。即：\(P=P\left( x,t \right),V=V\left( x,t \right),{{A}_{p}}={{A}_{p}}\left( x,t \right)\) \(P=P\left( x,t \right),V=V\left( x,t \right),{{A}_{p}}={{A}_{p}}\left( x,t \right)\)

\(P=P\left( x,t \right),V=V\left( x,t \right),{{A}_{p}}={{A}_{p}}\left( x,t \right)\)等。

当t=0时，这些参数的起始值\({{A}_{p}}={{A}_{p0}},\Pi ={{\Pi }_{0}},{{A}_{b}}={{A}_{b0}}\)均可由装药的起始形状尺寸给定。随着燃烧的进行，每一时间间隔Δt，根据燃去的装药肉厚e=rΔt,按各种装药形状导得的\({{A}_{b}}(e),{{A}_{p}}(e),\Pi (e)\)可求得不同时刻的\({{A}_{b}}(t),{{A}_{p}}(t),\Pi (t)\),进而求得不同时刻气流参数沿轴向的分布。

下面介绍本系统使用的符号及其含义：

• 轴向截面编号，取值范围（1→Secs）

j—          肉厚方向编号，取值范围（1→Steps）

x_i—        i截面x坐标；

e_j— j步骤的肉厚

П_i,j—     i截面j步骤的燃烧周界；

A_pi,j—     i截面j步骤的通气面积；

在t时刻或j步骤

V_i—        i截面速度（m/s）;

P_i—        i截面压强（Pa）;

W_i—       i截面的肉厚（m）；

r_i― i截面的燃速（m/s）

在装药前端i=1或x₀=0

\(V\left[ 1 \right]=0.0  \)

   \(P\left[ 1 \right]={{p}_{1}}  \)

   \(T\left[ 1 \right]={{T}_{1}}={{T}_{01}}={{T}_{0}}  \)

   \(\rho \left[ 1 \right]={{\rho }_{1}}=\frac{{{p}_{1}}}{R{{T}_{0}}}  \)

  （8-121）

在装药末端i=Secs或X=L处燃气流参数是由前面的截面逐个计算出来的，并作为喷管的进口参数：

随着时间的推移，当燃烧室的压力小于有效工作压强时，计算即可进入拖尾阶段。

计算涉及的函数很多,下面将主要介绍燃面的求法及燃烧线和通气面积的计算。由装药设计软件计算出燃面。燃面为理想的平行层推移。然而实际的燃面由于侵蚀的作用，推移如图中的控制体所示。

取计算的一控制体，左端面为1，右端面为2，高为ΔX_i端面的形状与装药的结构有关。端面的通气面积和燃烧周长为：

当 \({{e}_{j-1}}<{{W}_{i}}\le {{e}_{j}}\)

   \({{\Pi }_{i}}={{\Pi }_{i,j-1}}+\frac{{{W}_{i}}-{{e}_{j-1}}}{{{e}_{j}}-{{e}_{j-1}}}\left( {{\Pi }_{i,j}}-{{\Pi }_{i,j-1}} \right)  \)

   \({{A}_{p}}_{i}={{A}_{pi,j-1}}+\frac{{{W}_{i}}-{{e}_{j-1}}}{{{e}_{j}}-{{e}_{j-1}}}\left( {{A}_{pi,j}}-{{A}_{pi,j-1}} \right)  \)

         （8-122）

体积和燃烧面积为

   \({{V}_{i}}=\frac{{{A}_{p1}}+{{A}_{p2}}}{2}\Delta {{x}_{i}}  \)

   \({{A}_{bi}}=\pi \left( {{R}_{1}}+{{R}_{2}} \right)\sqrt{\Delta x_{i}^{2}+{{\left( {{R}_{2}}-{{R}_{1}} \right)}^{2}}}  \)

       （8-123）

式中：

     \({{R}_{1}}=\frac{{{\Pi }_{1}}}{2\pi }\)，\({{R}_{2}}=\frac{{{\Pi }_{2}}}{2\pi }\)                    （8-124）

（8-120）（8-123）（8-124）式是近似的数值计算结果。

   \({{V}_{2}}=V\left[ Secs \right]  \)

   \({{P}_{2}}=P\left[ Secs \right] \)

   \({{T}_{2}}=T\left[ Secs \right] \)

   \({{\rho }_{2}}=\frac{{{P}_{2}}}{R{{T}_{2}}}=\rho \left[ Secs \right] \)

                     （8-125）

在装药末端，装药的生成量应等于喷管的排出量

\({{\rho }_{2}}{{V}_{2}}{{A}_{p}}\left[ Secs \right]=\frac{\Gamma }{\sqrt{R{{T}_{0}}}}{{P}_{02}}{{A}_{t}}\)                  （8-126）

又由气动函数

\({{P}_{02}}={}^{{{P}_{2}}}/{}_{\pi \left( {{\lambda }_{2}} \right)}\)                      （8-127）

\({{\rho }_{2}}={{{P}_{2}}}/{R{{T}_{2}}}\;\)                      （8-128）

\({{\lambda }_{2}}{{{V}_{2}}}/{{{a}_{*}}}\;\)                       （8-129）

\(\pi \left( {{\lambda }_{2}} \right)=\left( 1-\frac{k-1}{k+1}\lambda _{2}^{2} \right)\)                     （8-130）

\({{a}_{*}}=\sqrt{\frac{2k}{k+1}R{{T}_{0}}}\)                        （8-131）

\(\Gamma ={{\left( \frac{2}{k+1} \right)}^{\frac{k+1}{2\left( k-1 \right)}}}\sqrt{k}\)                     （8-132）

在求解方程组计算气流参数沿轴向分布之前，必须首先确定装药头部的压强P1，作为一个边界条件，它本身也是由计算获得。在t=0时可采用平衡压强

\({{P}_{1}}={{10}^{6}}{{\left[ {{\rho }_{p}}{{C}^{*}}a{{K}_{N}} \right]}^{\frac{1}{1-n}}}\)

在\(t\ne 0\)时可取上一时刻的P₁值作为下一时刻的初值。

• 取前一截面参数计算 j 截面参数

\(\lambda =\frac{V}{\sqrt{\frac{2k}{k+1}R{{T}_{0}}}}\)

• 计算侵蚀速度系数

\(\varepsilon =  1.0\) 当\(\lambda <{{\lambda }_{th}}  \)

   \(\varepsilon =  1+{{K}_{th}}\left( \lambda -{{\lambda }_{th}} \right) \)    当\(\lambda >{{\lambda }_{th}}  \)

• 计算静温

\(T={\left( {{C}_{p}}{{T}_{0}}-\frac{{{V}^{2}}}{2} \right)}/{{{C}_{p}}}\;\)

• 计算燃气的密度

\(\rho \frac{P}{RT}\)

• 计算当地音速

\(a=\sqrt{kTR}\)

• 计算燃速

\({{r}_{b}}=\varepsilon a{{P}^{n}}\left[ 1+{{\sigma }_{P}}\left( {{T}_{i}}-{{T}_{0}}^{‘} \right) \right]\)

• 计算当前肉厚，该截面前一时刻肉厚＋Δt时间内烧去的肉厚

\(W={{W}_{j}}+{{r}_{b}}t\)

• 计算装药在该截面的通气面积和燃烧周长

通过上述计算，由i截面的数据可计算出i+1截面的气流参数。对于给定的P1参数从头计算到尾可以获得装药末端的参数，并得到装药末端静压P₂和速度系数\({{\lambda }_{2}}\)。由装药末端通气面积和喷管喉部面积可以计算收敛段的\(q\left( {{\lambda }_{2}}^{‘} \right)\)，可以计算出喷管的入口速度系数\(\lambda _{2}^{‘}\)，如果计算收敛应满足\(\left| {{\lambda }_{2}}-\lambda _{2}^{‘} \right|\le \varepsilon \)，如果不收敛采用\(\frac{{{\lambda }_{2}}}{\lambda _{2}^{‘}}{{P}_{1}}\)修正头部压强，重新计算直到满足误差要求为止。

8.3.3  算例

下面以一个具体的例子介绍计算过程。

[例] 已知：某火箭发动机推进剂的密度\({{\rho }_{P}}=1.69\times {{10}^{3}}kg/{{m}^{3}}\)，燃气的比热比k=1.222，特征速度C*＝1575m/s，燃速系数\(a=3.317\times {{10}^{-3}}\text{m/s}{{\left( \text{MPa} \right)}^{\text{n}}}\)，压强指数n=0.4365，装药星孔数N=7，角分数\(\varepsilon 0.9\)，星边半角\({\theta }/{235{}^\circ }\;\)，星孔过渡圆弧半径r=3mm，星根半径r₁=6mm，特征长度l=31mm，装药肉厚e₀=23.45mm，外径D=114.9mm,装药长度L=1560mm，燃烧室内径D_I=117.9mm，喷喉半径R_t＝18.5mm，膨胀比为14。假定发动机工作在10000m高度。假定喷管的工作效率为\({{\eta }_{F}}0.96\)，燃烧室的热效率为\({{\eta }_{C}}0.99\)。推进剂的气体常数为Rg=287N.m/kg.K，燃气温度为Tg=2900K。推进剂温度敏感系数\({{\sigma }_{P}}=0.002(1/K)\)，推进剂的初温和给出推进剂燃速特性的参考温度都是20℃。推进剂的侵蚀函数为：

\(\varepsilon \left( \lambda  \right)=1\) ， \(\left( \lambda \le 0.11 \right)  \)

   \(\varepsilon \left( \lambda  \right)=1+2.47\left( \lambda -0.11 \right) \),   \(\left( \lambda >0.11 \right)  \)

试计算一维内弹道。

（1）装药计算

在装药设计模块绘制一个星角，草图如图8-8（a），三维效果图如图8-8（b）：

（a）草图

（b）三维效果图

图8-8  星形装药设计

在装药计算时计算给出装药通气面积和燃烧周长数据，格式如下：

InternalBals_3D ←数据类型标识

• 面数

0．000000←截面1

…………………

1．5640←截面n=20处轴线坐标

0.000000 ←起始肉厚

0.279263,0.002363 ←1截面燃烧周长m，通气面积m²

……………

0279263,0.002363 ←20截面燃烧周长m，通气面积m²

……………

0.028000 ←结束肉厚

0.000000,0.010132 ←燃烧周长m，通气面积m²

 ……………

0.000000,0.010132 ←燃烧周长m，通气面积m²

该格式是在一种肉厚下，沿轴向计算20个截面的燃烧周长和通气面积。在下一个肉厚下，又是20个截面。直到所有燃面为零时才终止。

    输入喷管数据和用户选择数据后，可进行计算，计算结果曲线如图8-9，为了比较侵蚀的影响，将未考虑侵蚀时的内弹道示于图8-9中。

图8-9  考虑侵蚀与不考虑侵蚀内弹道计算结果曲线

为了体现不同温度对内弹道的影响，将计算结果列于同一图中（图8-10）。

图8-10  不同温度下一维内弹道计算结果

如果没有装药数据和截面划分段数数据则程序不能执行计算。如果出现了浮点溢出，可检查有关数据的合理性。如果出现浮点溢出，则有可能是装药设计得不合理，要么是燃面与通气面积之比太大，要么是喉通比太大，或者是侵蚀燃速数据不合理。

8.4  双燃速一维内弹道计算

双推力固体火箭发动机在进程战术导弹上得到普遍应用。其主要原因是：固体推进剂有了更新的发展，能为双推力固体火箭发动机提供满足使用性能的推进剂；导弹总体要求解决单推力固体火箭发动机的近界可用过载过大，远界可用过载不足的突出问题；双推力发动机使导弹在整个飞行过程中能量分配合理，结构简单，重量轻，特别是能使导弹进行主动段攻击，提高了导弹的机动性。

使用双燃速发动机，不仅可获得较大的两级推力比，而且可提高发动机的装填分数，因此，双燃速发动机是双推力发动机的主要形式。

双燃速发动机采用两级串联药柱装药时，往往选择两种不同的推进剂，其性能参数，如密度，燃速等均不同。发动机工作时，特征速度，比冲等也不相同。其后端装药的热力学参数随着截面的位置和两种燃气的质量生成率的不同而变化。另外，如果前端装药的通气面积比较小，则在前端装药部分就可能发生侵蚀效应,出现前、后端装药同时出现侵蚀的现象。

目前，双推力发动机已广泛用于中小型战术导弹，其中相当部分使用的是双推进剂发动机。

在前面双燃速内弹道和一维内弹道讨论的基础上综合利用两者的知识解决双燃速一维内弹道问题。

（1）基本假设

• 两种推进剂的燃烧产物是物理混合的理想气体，它们之间存在热交换，但系统没有质量变化，总焓不变，是一个绝热等熵系统。
• 两种燃气混合时，不发生任何化学反应，而是以各自的燃气独立存在，组成混合气体，由喷管排出。
• 混合燃气满足理想气体状态方程。
• 燃烧产物在燃烧室中是一维定型流动。
• 燃烧过程中热损失为一常数，用损失系数表示。
• 装药全面点燃是压力达到点火压力时实现的。
• 装药为串装，且速燃放在后端，缓燃放在前端。
• 装药设计由专门的装药设计软件来完成，提供给该系统的是：各肉厚下的燃烧周长和通气面积。

内弹道计算模型采用与8.3节相同的模型，混合气体参数采用与8.2节相同的模型。

（2）内弹道计算方法

对于串装装药，前段内弹道计算同一维内弹道。而对于后段则是一维与双燃速的结合。下面介绍计算过程：

1）计算初始压强作为迭代的起始压强P₀。当t=0时，参数的起始值\({{A}_{p}}={{A}_{p0}},\Pi ={{\Pi }_{0}},{{A}_{b}}={{A}_{b0}}\)均可由装药的起始形状尺寸给定。利用零维内弹道计算公式计算出压强P₀。

2）选取装药1（前段）中一控制体，以头部压强P₀气流速度V₀=0，从截面1到截面Secs分别计算。获得在装药一末端气体的质量流率、压强和流速等参数。

3）计算装药2入口气体参数。根据装药通气面积Ap计算入口气流的λ₂，采用下式：

\(Z\left( {{\lambda }_{2}} \right)=\frac{1}{{{\lambda }_{2}}}+{{\lambda }_{2}}={{P}_{1}}{{A}_{P2}}\frac{2k}{k+1}\frac{1}{{{{\dot{m}}}_{1}}{{a}_{cr}}}\)

4）选取装药2（后段）中一控制体，以前段末端压强P₂气流速度，从截面1到截面SSecs分别计算。获得在装药2末端气体的质量流率、压强和流速等参数。控制体内，m₁为前段装药流入质量流率，m₂₁为后段装药控制体入口端装药2的质量生成率，m₂₂为控制体内装药2的质量生成率。

图8-11  装药2控制体示意图

5）根据喷管收缩比计算喷管气流速度系数。喷管入口面积取装药末端出口通气面积。

6）根据内流场计算的末端气体流动速度系数与步骤5）的速度系数比较，如果精度满足要求，则迭代终止，否则修正头部压强，重新迭代计算，直到满足要求为止。

所用公式同一维内弹道计算，兹不赘述。参数的含义同零维内弹道。

图8-12  有侵蚀和无侵蚀内弹道曲线对比

从比较可以看出，侵蚀能提高初始压力和压力峰，且有较长的拖尾。

图8-13  不同温度下计算结果比较

（3）需要注意的问题

• 装药设计和计算应考虑交界面的特性。
• 对于等截面装药应注意适当地划分其轴向段数，轴向划分段数太少，则在计算过程中龙格－库塔算法会不收敛，得不到正确的计算结果；轴向划分段数太多，则计算耗时较多，另外，计算结果数据过于庞大，不便于分析。一般情况下，对于燃速较大，侵蚀效应比较明显的等截面装药，其轴向划分段数应相对较多；对于燃速较小，侵蚀不明显的等截面装药，其轴向划分段数可以相对较少。
• 采用前、后端装药的形式虽然能较好地实现两级推力，但是，当速燃药烧完后，燃烧室壁将直接暴露在高温燃气下，这就需要在壳体设计阶段注意燃烧室的热防护。
• 如果计算出现浮点溢出，有可能是装药设计不合理，如面通比太大或喉通比较大。
• 调节两级的内弹道曲线，最重要的是燃速数据和装药的结构数据。
• 这种方法只适用于串装装药，对于套装装药不适用。

8.5  燃气发生器内弹道计算

燃气发生器工作时间较短，对于一些特殊用途的燃气发生器，如弹射用燃气发生器内弹道曲线变化较为剧烈，内弹道计算要综合考虑点火因素、质量守恒和能量守恒等因素的影响。燃气发生器一般是为其它系统提供一定质量流率的燃气，而不在于产生多大的推力，因此喷管的主要参数是喉部面积。其设计过程是根据质量流率的要求，选择适当的药型，根据工作时间和所选择的推进剂，确定装药的几何参数，根据燃烧室自由容积和所选推进剂，选择点火要量和粒度。在计算出装药燃面肉厚数据后即可进行内弹道计算。本节要介绍其基本方程及程序设计。

8.5.1  内弹道基本假设

1. 假设发生器内燃气均匀分布，从而燃气的P、r、T沿筒轴向各处是均一的；
2. 假设燃气的成分、物理化学性质是固定不变的，从而爆温T₀、绝热指数k、比热C _v和C _p等均为常量；
3. 在发生器内燃气无流动，燃气在喷管中为一维绝热等熵流动；
4. 假设装药在表面温度达到T_s时，瞬时全面点燃，装药的燃烧服从几何燃烧定律，并且燃烧是均匀的、完全的；
5. 考虑流量损失、散热损失以及摩擦等，引入能量损失系数X _t。

8.5.2  基本方程

1. 固相点火理论模型

\({{\rho }_{s}}{{C}_{s}}\frac{\partial T\left( x,t \right)}{\partial t}={{\lambda }_{s}}\frac{{{\partial }^{2}}T\left( x,t \right)}{\partial {{x}^{2}}}+Q\),    \((x\ge 0,t>0)\)            （8-133）

式中：T（x，t）—固相温度，K；

r_s—     固相装药的密度，kg/m³；

C_s—     固相装药的比热，J/kg . K；

λ_s—   固相装药的热传导率，W/m . K；

Q—     单位体积装药化学反应放热速率，J/m³ . s。

1. 质量守恒方程

\(\frac{d}{dt}\left( {{\rho }_{b}}V \right)={{M}_{b}}+{{M}_{ig}}-{{G}_{c}}\).                    （8-134）

式中：r_b—     发生器燃气密度，kg/m³；

V—     发生器内自由容积，m³；

M_b—   主装药燃气生成率，kg/s；

M_ig—  点火药燃气生成率，kg/s；

G_c—    经喷管流出的燃气秒流量，kg/s；当发生器内的压强大于破膜压力P_I时，喷口堵盖破裂，燃气开始喷出。在此之前G_c=0。此后由下式求出：

\({{G}_{c}}={{\phi }_{2}}\frac{{{P}_{b}}{{A}_{t}}}{\sqrt{R{{T}_{b}}}}\sqrt{\frac{2k}{k-1}\left[ {{\left( \frac{{{P}_{c}}}{{{P}_{b}}} \right)}^{\frac{2}{k}}}-{{\left( \frac{{{P}_{c}}}{{{P}_{b}}} \right)}^{\frac{k+1}{k}}} \right]} \),       \(1\ge \frac{{{P}_{c}}}{{{P}_{b}}}\ge {{\left( \frac{2}{k+1} \right)}^{\frac{k}{k-1}}} \)

 

 \({{G}_{c}}={{\phi }_{2}}\sqrt{k}{{\left( \frac{2}{k+1} \right)}^{\frac{k+1}{2\left( k-1 \right)}}}{{P}_{b}}{{A}_{t}}/\sqrt{R{{T}_{b}}} \),    \(\frac{{{P}_{c}}}{{{P}_{b}}}<{{\left( \frac{2}{k+1} \right)}^{\frac{k}{k-1}}} \)

  （8-135）

式中：T _b—   发生器内燃气的平均温度，K；

P _b—    发生器内燃气的平均压强，Pa；

P _c—    弹射缸内燃气的平均压强，Pa；

R—      燃气的气体常数，J/kg · K；

k—      燃气的绝热指数；

A_t—     燃气发生器的喷喉面积，m²；

\[{{\phi }_{2}}\]—    喷管流量修正系数。

1. 能量守恒方程

\(\frac{d}{dt}\left( {{C}_{v}}{{\rho }_{b}}V{{T}_{b}} \right)={{M}_{b}}{{X}_{1}}{{C}_{pb}}{{T}_{pb}}+{{M}_{ig}}+{{C}_{pig}}{{T}_{pig}}-{{G}_{c}}{{C}_{p}}{{T}_{b}}\)       （8-136）

式中：C _v—   发生器内燃气的定容比热，J/kg · K；

C _{p b}— 装药燃气的定压比热，J/kg · K；

T _{p b}— 装药燃气的定压燃烧温度，K；

C _pig— 点火药燃气的定压比热，J/kg · K；

T _pig— 点火药燃气的定压燃烧温度，K；

C _p—   喷管出口燃气的定压比热，J/kg · K；

X _l—    发生器内燃气热损失修正系数。

1. 状态方程

\({{P}_{b}}={{\rho }_{b}}R{{T}_{b}}\)                           （8-137）

1. 主装药燃面变化规律，燃面变化规律由装药设计部分完成，并提供燃面肉厚数据表供内弹道计算调用。

A_b=f（e）                             （8-138）

将以上（8-133）~（8-138）式联立组成方程组，进行计算求解，即得出发生器内弹道参数。

8.5.3  计算流程

点火过程主要包括在密闭空腔建压和对主装药表面的加热，当压强建立达喷管堵盖破裂时即有气体流出，当装药表面达到推进剂发火点时，认为装药全表面瞬时燃烧。点火药被认为是球形颗粒，当其燃烧完后，点火药的工作结束。计算流程如图8-15所示。

图8-14  点火子程序流程图

当装药全表面燃烧时即进入主要工作阶段。如果点火药没有烧完，则混合气体为两种燃气混合物，否则只有一种燃气。当主装药烧完或压强低于用户规定的值，工作过程结束。

 

图8-15  主装药燃烧子程序流程图

8.5.4  算例

一燃气发生器装药采用双基药。装药为单孔柱型，内侧面、外侧面、端面全方位燃烧，装药数量为一根，药柱的长度为0.07m，外径0.036m，内径0.02m。燃速公式为：

\(u=0.01236{{P}^{0.367}}m/s\)（压强为MPa）

与装药有关的其它参数如下：

1. 火药力f _v=1164240J/kg；
2. 火药燃烧温度T_p=3152K；
3. 燃气绝热指数k=1.21；
4. 燃气气体常数R=305.26J/kg·K；
5. 燃速温度系数为0.2%。

点火药选用黑火药，点火药量为0.001kg。有关的参数如下：

1. 火药力f _vig=252840 J/kg；
2. 燃烧温度T _pig=2355K；
3. 绝热指数k _ig=1.1；
4. 气体常数R _ig=97.6J/kg·K。

图8-16  内弹道计算结果压强曲线

8.6  比冲效率计算

在内弹道计算过程中比冲效率的计算至关重要，因为它直接影响推力预示的准确性。在发动机设计过程中理论计算往往提供的是理论比冲，而实际发动机比冲是考虑各种损失后修正的结果，如果在理论预估能够较为准确地预估发动机实际比冲，可以大大节约研制成本。目前，理论上难以精确预示实际比冲，往往是根据现有的发动机根据回归分析或统计、拟合的处理方法给出这类发动机预估的经验公式，这些公式有一定的适用范围。本文首先介绍理论比冲的计算方法，然后再考虑比冲效率修正。

8.6.1  理论比冲计算

（1）由推进剂热力参数计算

由用户给定的推进剂的k、特征速度C*、环境压强P_a及喷管的膨胀比ε_A，采用下式，按一维冻结流进行计算。

\({{C}_{F}}=\Gamma \sqrt{\frac{2k}{k-1}\left[ 1-{{\left( \frac{{{P}_{e}}}{{{\eta }_{p}}{{P}_{c}}} \right)}^{\frac{k-1}{k}}} \right]}+{{\varepsilon }_{A}}\left( \frac{{{P}_{e}}}{{{\eta }_{p}}{{P}_{c}}}-\frac{{{P}_{a}}}{{{\eta }_{p}}{{P}_{c}}} \right)\)           （8-139）

其中η_P燃烧室总压恢复系数

显然当没有尾管时η_P=1。

\({{I}_{s}}_{D}={{C}_{F}}{{C}^{*}}\)                         （8-140）

用该式计算的是冻结流比冲。

（2）由热力计算程序直接计算出推进剂在工作条件下的平恒流理论比冲

直接调用热力计算程序费时，将热力计算结果制成表，计算时根据燃烧室压强和喷管膨胀比插值计算获得比冲，这样也较为麻烦。通常在发动机设计时，理论计算的是冻结流比冲，而比冲效率都是修正平衡流比冲，因此这里引出一个参数，冻结流与平衡流比冲之比Q，对于给定的推进剂Q基本是常数：

\({{I}_{s}}={{I}_{sD}}/Q\)                          （8-141）

IsD 采用（8-139）（8-140）式计算。

8.6.2  能量损失模型

发动机的能量损失包括有燃烧室热效率和喷管的效率。下面将介绍模型中有关处理各项损失的方法。

在初步设计时我们认为发动机内的流动为一维平衡等熵定常绝热流。依此我们计算出它的理论比冲，但是，实际上，这些比冲是有各种各样损失的，为此应对其进行修正。

\({{\xi }_{0}}={{I}_{s,ex}}/{{I}_{s,th}}\)                       （8-142）

其中：x₀–冲量系数；

I_s,ex–实际比冲；

I_s,th–理论比冲。

只要能给出冲量系数x₀，实际比冲便可求出。目前常用的两种方法，即解析法和多元回归分析法，本文都采用。

对于有些小型发动机，由于绝热层的烧蚀也将产生一部分附加的燃气，因此，引入一个由于惰性气体的流出而引起的比冲效率因子：

\({{\xi }_{b}}=\sqrt{1+{{\varphi }_{b}}}\)                         （8-143）

式中：\({{\varphi }_{b}}\)为惰性排出物与装药质量之比。对于一些大型的贴壁浇铸的发动机尽管绝热层占据一定的质量，但并非参与燃烧，这种情况下可取\({{\varphi }_{b}}=0.0\)。而对于那些绝热层完全暴露在燃气中的发动机，则该项不能忽略。综上所述，比冲效率可表示为：

\({{\xi }_{0}}={{\xi }_{c}}{{\xi }_{N}}{{\xi }_{b}}\)                        （8-144）

1.7.2节能量特性中介绍了SPP的经验公式及回归分析方法，本章补充一些经验公式共参考。

（1）Landsbaum 8变量回归分析经验公式

应用多元回归分析法确定的经验公式为：

\({{\xi }_{0}}=1.1328+0.0122\ln {{D}_{t}}-0.0465\ln \beta -0.230Al-0.00448\ln \varepsilon +0.0153\ln (1-{{\overline{L}}_{N}}) \)

 \( -0.00842\ln S-0.00548\ln {{P}_{c}}+0.00333\ln \rho\)     （8-145）

式中：Dt–喷管喉径（cm）;

b–反映喷管扩张损失参数（deg）；

\(\beta =\left( \overline{\alpha }+2{{\alpha }_{e}} \right)/3 \)

\(\overline{\alpha }=t{{g}^{-1}}\left[ {{D}_{t}}\left( {{\varepsilon }^{0.5}}-1 \right)/2{{L}_{N}} \right] \)                 （8-146）

其中：a_e为喷管出口扩张半角（deg）；

LN为喷管喉部至出口段的长度。

e –喷管扩张比；

Al–推进剂中铝粉含量的百分比；

\({{\bar{L}}_{N}}\)–喷管的嵌入分数\({{\bar{L}}_{N}}\)＝喷管的潜入长度/燃烧室空腔长度；

S –喷管进口面积与喉部面积之比；

Pc –燃烧室平均压强（kPa）;

r –喷喉前圆弧曲率半径与喉部半径之比；

（2）Landsbaum 4变量回归分析经验公式

应用多元回归分析法确定的经验公式为：

\({{\xi }_{0}}=1.0605+0.0111\ln {{D}_{t}}-0.0328\ln \beta -0.245Al-0.00617\ln \varepsilon \)    （8-147）

式中：Dt–喷管喉径（cm）;

b –反映喷管扩张损失参数（deg）；

\(\beta =\left( \overline{\alpha }+2{{\alpha }_{e}} \right)/3 \)

 \(\overline{\alpha }=t{{g}^{-1}}\left[ {{D}_{t}}\left( {{\varepsilon }^{0.5}}-1 \right)/2{{L}_{N}} \right] \)                 （8-148）

其中：a_e为喷管出口扩张半角（deg）。

e–喷管扩张比；

Al–推进剂中铝粉含量的百分比；

（3）国内大发动机8变量回归分析经验公式

应用多元回归分析法确定的经验公式为：

\({{\xi }_{0}}=1.094+0.01239\ln {{D}_{t}}-0.03553\ln \beta -0.2471Al-0.005448\ln \varepsilon +0.008054\ln (1-{{\overline{L}}_{N}}) \)

 \( -0.002976\ln S-0.004868\ln {{P}_{c}}+0.004393\ln \rho\)                                                                                  （8-149）

（4）国内大发动机4变量回归分析经验公式

应用多元回归分析法确定的经验公式为：

\({{\xi }_{0}}=1.04013+0.011111\ln {{D}_{t}}-0.03028\ln \beta -0.2181Al-0.006874\ln \varepsilon \)  （8-150）

（5）国内小发动机8变量回归分析经验公式

应用多元回归分析法确定的经验公式为：

\({{\xi }_{0}}=1.068+0.004818\ln {{D}_{t}}+0.00598\ln \beta +0.2336Al-0.00314\ln \varepsilon -0.00907\ln (1-{{\overline{L}}_{N}}) \)

\( -0.00973\ln S-0.0221\ln {{P}_{c}}-0.003613\ln \rho\)  

  （8-151）

（6）国内小发动机4变量回归分析经验公式

应用多元回归分析法确定的经验公式为：

\({{\xi }_{0}}=1.04013+0.011111\ln {{D}_{t}}-0.03028\ln \beta -0.2181Al-0.006874\ln \varepsilon \)      （8-152）

这些公式具有一定的适应范围，可以根据实际发动机选择验证使用。