固体火箭发动机零维内弹道方程组

目的

对于给定装药，燃面面积随肉厚的变化$A_{\mathrm{b}}\left(e\right)$是已知的，需要根据该曲线计算得到燃烧室压力随时间的变化关系$p_{\mathrm{c}}\left(t\right)$

连续方程

1 装药燃烧产气的流量

$${\dot{m}}_1={\rho }_{\mathrm{p}}{A}_{\mathrm{b}}r=a{\rho }_{\mathrm{p}}{A}_{\mathrm{b}}{p}_{\mathrm{c}}^{\mathrm{n}}$$这里${\rho }_{\mathrm{p}}$是装药的密度，$r$是线燃速，可用$a{p}_{\mathrm{c}}^{\mathrm{n}}$经验公式计算

2 喷管流出燃气的流量

$${\dot{m}}_2=\frac{p_{\mathrm{c}}{A}_{\mathrm{t}}}{c^{\ast }}$$式中$c^{\ast }$是推进剂的特征速度，$A_{\mathrm{t}}$是喷管的喉部面积，需要注意的是，该公式的使用前提是喷管处于最大流量状态（喉部音速），所以对发动机启动过程和关机可能不适用

3 燃烧室内积存的燃气流量

$${\dot{m}}_3=\frac{\mathrm{d}\left({\rho }_{\mathrm{g}}{V}_{\mathrm{f}}\right)}{\mathrm{d}t}$$式中$V_{\mathrm{f}}$是燃烧室空腔的容积，${\rho }_{\mathrm{g}}$是燃烧室中燃气的密度，假设燃气温度和分子量不变，将完全气体状态方程带入，有
$${\dot{m}}_3=\frac{1}{RT}\frac{\mathrm{d}\left(p_{\mathrm{c}}{V}_{\mathrm{f}}\right)}{\mathrm{d}t}=\frac{1}{{\left(c^{\ast }\Gamma \right)}^2}\frac{\mathrm{d}\left(p_{\mathrm{c}}{V}_{\mathrm{f}}\right)}{\mathrm{d}t}$$式中$R$是完全气体常数，$T$是燃气温度。为了避免直接计算燃气温度（需要化学平衡计算，不便于使用），故引进了$\Gamma$函数，使用${\left(c^{\ast }\Gamma \right)}^2$替换$RT$，可以证明二者是相等的
$$\Gamma =\sqrt{\gamma }{\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma +1}{2\left(\gamma -1\right)}}$$式中$\gamma$是比热比，可以简单地取为1.2
所以使用乘积的求导公式，有
$${\dot{m}}_3=\frac{1}{{\left(c^{\ast }\Gamma \right)}^2}\frac{\mathrm{d}\left(p_{\mathrm{c}}{V}_{\mathrm{f}}\right)}{\mathrm{d}t}=\frac{1}{{\left(c^{\ast }\Gamma \right)}^2}\left(p_{\mathrm{c}}\frac{\mathrm{d}{V}_{\mathrm{f}}}{\mathrm{d}t}+V_{\mathrm{f}}\frac{\mathrm{d}{p}_{\mathrm{c}}}{\mathrm{d}t}\right)$$

4 流量平衡

$${\dot{m}}_1={\dot{m}}_2+{\dot{m}}_3$$故有
$$a{\rho }_{\mathrm{p}}{A}_{\mathrm{b}}{p}_{\mathrm{c}}^{\mathrm{n}}=\frac{p_{\mathrm{c}}{A}_{\mathrm{t}}}{c^{\ast }}+\frac{1}{{\left(c^{\ast }\Gamma \right)}^2}\left(p_{\mathrm{c}}\frac{\mathrm{d}{V}_{\mathrm{f}}}{\mathrm{d}t}+V_{\mathrm{f}}\frac{\mathrm{d}{p}_{\mathrm{c}}}{\mathrm{d}t}\right)\text{\hspace{1em}(1)}$$

燃面推移方程

燃烧会导致燃面推移，有
$$\frac{\mathrm{d}e}{\mathrm{d}t}=a{p}_{\mathrm{c}}^{\mathrm{n}}\text{\hspace{1em}(2)}$$燃烧还会导致燃烧室空腔的容积扩大，有
$$\frac{\mathrm{d}{V}_{\mathrm{f}}}{\mathrm{d}e}=A_{\mathrm{b}}$$结合这两式，不难得到
$$\frac{\mathrm{d}{V}_{\mathrm{f}}}{\mathrm{d}t}=A_{\mathrm{b}}a{p}_{\mathrm{c}}^{\mathrm{n}}\text{\hspace{1em}(3)}$$

最终的微分方程组

将式$(3)$代入$(1)$可以化简为
$$a{\rho }_{\mathrm{p}}{A}_{\mathrm{b}}{p}_{\mathrm{c}}^{\mathrm{n}}=\frac{p_{\mathrm{c}}{A}_{\mathrm{t}}}{c^{\ast }}+\frac{1}{{\left(c^{\ast }\Gamma \right)}^2}\left(A_{\mathrm{b}}a{p}_{\mathrm{c}}^{\mathrm{n}+1}+V_{\mathrm{f}}\frac{\mathrm{d}{p}_{\mathrm{c}}}{\mathrm{d}t}\right)\text{\hspace{1em}(4)}$$将式$(2)(3)(4)$整理后，写在一起即可得到最终的微分方程组
$$\frac{\mathrm{d}{p}_{\mathrm{c}}}{\mathrm{d}t}=\frac{{\left(c^{\ast }\Gamma \right)}^2}{V_{\mathrm{f}}}\left(a{\rho }_{\mathrm{p}}{A}_{\mathrm{b}}{p}_{\mathrm{c}}^{\mathrm{n}}-\frac{p_{\mathrm{c}}{A}_{\mathrm{t}}}{c^{\ast }}\right)-\frac{A_{\mathrm{b}}a{p}_{\mathrm{c}}^{\mathrm{n}+1}}{V_{\mathrm{f}}}$$$$\frac{\mathrm{d}{V}_{\mathrm{f}}}{\mathrm{d}t}=A_{\mathrm{b}}a{p}_{\mathrm{c}}^{\mathrm{n}}$$$$\frac{\mathrm{d}e}{\mathrm{d}t}=a{p}_{\mathrm{c}}^{\mathrm{n}}$$给定初始条件后，即可使用数值计算方法求解