对火箭进行数学建模 | Aaron212的Github博客

Aaron212的Github博客

对火箭进行数学建模

Jan 21, 2023

引入

21世纪是航天的世纪，我们的科学和技术的发展已经使得人类能够离开地球，探索宇宙。航天学是一门重要的学科，它研究如何在宇宙中运行航天，而火箭技术是航天学的核心部分，它是人类离开地球的基础。

提出问题

什么样的火箭才能带着有效载荷离开地球、飞向深空呢？

建立模型

地球上的火箭在起飞时会受到三个力的作用：重力、阻力与推力。可以用如下物理公式表达

[\Sigma F = F_t - F_f -G]

转换为数学公式是

[m \frac{d^{2} z}{d t^{2}}=F{t}(t) - F{f}\left(z, \dot{z}^{2}\right) - F_{G}(z, t)]

其中$t$代表时间，$z$代表火箭的高度，$\dot{z}$表示火箭高度关于时间的求导。

推力

大多数进入轨道的火箭都是液体燃料火箭，而不是固体燃料火箭。对于液体燃料火箭，推力将在上升过程中由油门控制，其可以根据火箭的状态，随时间调整推力以控制轨迹。我们假设模型火箭的推力在上升过程中保持不变，因此推力的数学表示可以是

[F{t}(t)= \left{ \begin{aligned} &C &, t<t{b} \ &0 &, t \geq t_{b} \end{aligned} \right.]

其中$C$代表某个常数，$t_b$代表了火箭耗尽燃料的时间。换句话说，推力在火箭耗尽燃料前一直保持不变，之后变为零。

阻力

作用在火箭上的阻力可以分为两部分：大气阻力和激波阻力。

大气阻力

大气阻力是寄生阻力的一种形式，是由于周围空气分子与火箭表面之间的摩擦而产生的与火箭运动相反的阻力。由大气表面摩擦引起的阻力方程为：

[F_f(t)= \frac{1}{2} C_D \rho S v^2]

其中$\rho$为空气密度$\mathrm{(kg/m^3)}$，$v$为火箭的速度，$C_D$为阻力系数，$S$为火箭的横截面积。每个火箭的直径都是给定的，因此$S$是已知的。如果不对火箭进行物理测试，则很难确定阻力系数，但我们将假设$C_D= 0.75$，这代表介于圆柱体和圆锥体之间的飞行形状。这是一个比较高的风阻系数；相比之下，机翼的$C_D$ 值远低于 0.05。

空气密度的计算更复杂，因为随着火箭上升穿过大气层，它会呈指数下降。对于86 公里以下的高度(对流层、平流层和中间层)，空气密度可由理想气体定律给出

[\rho(z)= \frac{M \cdot T(z)}{R \cdot P(z)}]

其中$R$是气体常数$\mathrm{(8.314 J \cdot K^{-1}\cdot mol^{-1})}$，$M$是空气的摩尔质量$\mathrm{(0.029 kg/mol)}$。大气温度$T(z)$和压力$P(z)$也取决于高度，并由下面的气压方程式给出。

[\begin{gathered} T(z)=T{0}-L z \ P(z)=P{0}\left(\frac{T{0}}{T{0}-L z}\right)^{\frac{g_{0} M}{R L}} \end{gathered}]

其中$T_0$和$P_0$是特定海拔区域的基线温度和压力，$L$是该区域的温度递减率，$g_0$是海平面重力加速度，这里取$\mathrm{9.81 m/s^2}$。

86公里以上是热层，温度上升得更快，因为大气层不再为空间遮挡太阳光线。压力和空气密度现在非常小，二者随着高度的四次关系下降，公式如下

[\begin{gathered} P(z)=A z^{4}+B z^{3}+C z^{2}+D z+E \ \rho(z)=A z^{4}+B z^{3}+C z^{2}+D z+E \end{gathered}]

每一项的系数在不同的高海拔地区是不同的。有关这些系数的表格以及$T_0$、$P_0$和$L$的值，请参阅地球大气的综合数学模型(http://www.braeunig.us/space/atmmodel.htm)。

激波阻力

一般来说，一旦火箭超过其临界马赫数(一般是0.8马赫)并接近音障，激波阻力就会发挥作用。当一个物体接近音速时，它开始比它前面的空气波移动得更快，导致阻力迅速增加。1马赫的阻力峰值会产生音障，这是火箭试图离开地球大气层的主要障碍。当在模型中不考虑激波阻力时，火箭能够毫无困难地到达1,000,000 米远的太空。

阻力和马赫数之间关系的精确计算很困难，并且在很大程度上依赖于火箭的具体形状。此外，大多数关于计算波浪阻力的可用信息来自空气动力学工程领域，该领域涉及飞机，而不是火箭。然而，普朗特-格劳厄脱法则有一个方程式可以帮助在估算跨音速阻力的影响。

但是在这种情况下使用普朗特-格劳厄脱法则需要进一步说明。首先，目前普遍认为普朗特-格劳厄脱法则对跨音速区的精确计算无效，仅对马赫数在0.8以下和1.2-5之间有效。我们这里选择忽略这个条件，因为这个区域下该法则非常接近音障的效果，足以满足模型。

普朗特-格劳厄脱法则指出，在处理可压缩流体时，空气动力学系数(例如压力、升力、阻力等)必须除以$\beta$，即普朗特-格劳厄脱因子，它取决于自由流马赫数。马赫数定义为速度与声速之比，而声速又取决于空气温度。普朗特-格劳厄脱法则的公式形式如下

[\begin{gathered} C{d}=\frac{c{d 0}}{\beta} \ \beta=\sqrt{1-M{\infty}^{2}} \ M{\infty}(v, z)=\frac{v{\infty}}{a{\infty}} \ a{\infty}(z)=\sqrt{\gamma R{\mathrm{sp}} T(z)} \end{gathered}]

在空气动力学中，通常使用随飞机移动的参考系进行计算。马赫数、速度和声速$(a)$下面的下标$\infty$指的是每一个的”自由流”值，意思是在这个参照系中火箭周围的移动空气的速度。在(4)中，$\gamma$是比热比(平均约为 1.4)，$R_{sp}$是比气体常数(干燥空气为$\mathrm{287 J \cdot kg^{-1} \cdot K^{-1}}$)。对于高于 1 的马赫数(超音速)，(2) 中平方根下的项是相反的。

下图演示了10km高空，速度与阻力系数$C_d$的关系

重力

重力的方程非常简单：$G = mg$。然而，在这种情况下，计算重力加速度和火箭质量会引入新的问题。

火箭的质量总是会随着燃料的排出而减少。继续假设推力恒定，因此火箭流出的质量流量恒定，我们可以对火箭质量建模，公式如下所示

[m(t)=m_{0} - \dot{m} t]

火箭质量与时间呈线性关系。距地球表面高度为 z 的重力加速度的公式如下所示

[g(z)=\frac{G M{E}}{\left(z+R{E}\right)^{2}}]

其中 $R_E$ 是地球的半径 $(6,371,000 \mathrm{m})$，$M_E$ 是地球的质量$(5.972 \times 1024 \mathrm{kg})$。于是火箭受到的重力可以如下表示

[F_{G}=m(t) \cdot g(z)]

给定火箭的推力和比冲，可以从齐奥尔科夫斯基火箭方程式计算质量流率。比冲以秒为单位，是衡量火箭发动机效率的指标。

[\dot{m}=\frac{F{t}}{g{0} I_{s p}}]

因为质量随着时间的推移而减小，而引力随着高度的增加而减小，所以重力取决于时间和高度。我们现在已经将运动方程简化为仅依赖于时间，高度作为时间的函数，速度(高度的导数)作为时间的函数。我们现在可以应用欧拉-克罗墨方法求解高度和速度并确定火箭的一维轨迹。

所有的方程如下 $$ \begin{gathered}

F_{G}(t,z)=m(t) \cdot g(z) \

m(t)=m_{0}-\dot{m} t \

g(z)=\frac{G M{E}}{\left(z+R{E}\right)^{2}}\

\dot{m}=\frac{F{t}}{g{0} I_{s p}}\

T(z)=\left{

\begin{aligned}

&T_{0}-L z &, z \leq 86 \mathrm{km}\

&1000 - A e^{\frac{B(x-86)}{C +x}} &, z > 86 \mathrm{km}

\end{aligned}

\right.\

P(z)=\left{

\begin{aligned}

&P{0}\left(\frac{T{0}}{T{0}-L z}\right)^{\frac{g{0} M}{R L}} &, z \leq 86 \mathrm{km}\

&A’ z^{4}+B’ z^{3}+C’ z^{2}+D’ z+E’ &, z > 86 \mathrm{km}

\end{aligned}

\right.\

\rho(z)=\left{

\begin{aligned}

&\frac{M \cdot T(z)}{R \cdot P(z)} &, z \leq 86 \mathrm{km}\

&A’’ z^{4}+B’’ z^{3}+C’’ z^{2}+D’’ z+E’’ &, z > 86 \mathrm{km}

\end{aligned}

\right.\

F_f(t)= \frac{1}{2} C_d \rho S v^2 \

C{d}=\frac{c{d 0}}{\beta} \

\beta=\sqrt{1-M_{\infty}^{2}} \

M{\infty}(v, z)=\frac{v{\infty}}{a_{\infty}} \

a{\infty}(z)=\sqrt{\gamma R{sp} T(z)}\

F_{t}(t)=

\left{

\begin{aligned}

&C &, t<t_{b} \

&0 &, t \geq t_{b}

\end{aligned}

\right.\

%FI

\frac{dv}{dt} = \frac{F_f(t) - F_f(t) - F_G(t,z)}{m(t)}\

%if v < 0: # flip drag force vector if rocket falls

%drag = drag*-1

\frac{dz}{dt} = v

\end{gathered}$$

其中1-4式计算了火箭所受的重力，5-12式计算了火箭所受的阻力，13式计算了火箭的推力，14-15式通过前13个式子计算了火箭关于时间的速度与高度。

需要特别注意的是，5-8式只是一个近似，对于模拟火箭升空，应选择更加复杂的分段函数来拟合大气不同高度的温度、密度与压力。

求解模型

我们以SpaceX的单级猎鹰1e为例，其数据如下表所示

火箭高度 (m)直径 (m)推力 (kN)载油质量 (kg)I_{sp}$ (s)猎鹰1e27.41.756046760304

将数据代入，得到结果如图1所示

速度曲线中65-95秒左右的平坦线代表火箭接近音障。相应的阻力突然增加代表了压缩性和激波阻力的影响。注意在95秒左右直接穿过音障的阻力图似乎表现出更加不稳定，可能是因为此时的推力和阻力正在争夺对合力的主导权。一旦突破音障，阻力很快就会变得可以忽略不计，因为空气波在火箭后面移动，空气密度变得非常稀薄。因此，跃升至超音速对于离开地球大气层至关重要。尽管跃升需要很大的推力和坚固的航天器，但由于火箭的阻力较大，选择保持亚音速飞行反而总体上效果会比较差。

我们可以看到，单级猎鹰1e飞行了大约169秒，然后燃料耗尽，推力降为零。到那时它已经达到了大约 125公里的高度，低于近地轨道所需的200公里。这与我们的预期是一致的：历史上没有任何一枚单级火箭从地球进入轨道。毕竟单级入轨航天器(SSTO) 的主题是当前研究的一个热门领域。

为了进入太空，我们的模拟火箭需要两级。如图2所示，我们可以看到第一级完成了大部分的累活。穿过音障后，火箭只需很少的额外推力或速度增加即可进入轨道，这正是因为高空和超音速下的阻力可以忽略不计。

总结

由上面两个模拟可以看出，激波阻力是火箭发射中最需要克服的难题，因此现在的火箭至少需要两级才能进入轨道，且第一级的推力需要足够大，来跨越音障造成的激波阻力。而跨越音障后，升空之路畅通无阻，只需要一个推力小的发动机就可以把火箭送入轨道。

但是需要特别注意的是，真实的火箭并不是直上直下，而是发射后会倾斜一定角度，将一小部分加速度应用到提升轨道速度，这样可以进一步节省燃油。这个模型只研究了火箭发射后垂直上升的情况，没有考虑火箭的倾角。不过可以确定是给火箭增加倾角后，空气阻力的的公式将会更加复杂。

通过建立并求解这个模型，我们能更好地理解航天学的核心——火箭，如何离开大地、飞向深空。

boyahub / Master-Dialogue

对火箭进行数学建模 | Aaron212的Github博客 #212