前几天看到了，不过这几天很忙没时间回，今天有点时间我仔细看了一下。

首先我觉得楼主在推力与总冲的计算上是没有问题的，用了一些近似，认为喷嘴面积远小于内部截面积。如果要考虑喷嘴面积，那么火箭加速度会带来一阶修正，而流体力学会带来二阶修正。我对水火箭不是很了解，不过如果补上和喷嘴面积有关的量应该会好一些，毕竟这对计算机来说计算量不算什么，当然，这不一定能减小误差，因为更复杂的流体力学效应与火箭的偏转影响会更大。

其次是通过总冲判断飞行高度和初速度的，这里我觉得问题可能挺大，只能用作估算。首先由于内部气体是绝热膨胀，火箭质量变化不是线性的，此时用平均质量不是很合理，另外由实验视频可以看出，喷水时间和飞行时间是同一数量级，用瞬间工作的模型来估计是不合理的。论坛常见的固体燃料火箭工作时间相比飞行时间小了至少1-2个数量级，但是水火箭不同，这一点倒是非常有必要，能大大增加精度。

我也用了ai辅助，让GPT帮我推导，我自己检查了一遍。

以下可能有一些量的字母或角标与一楼不同，V为总体积，Vg为气体体积，V0为初始气体体积，P0为大气压，P为初始气压。（emm……的确楼主的标注比较正常……但是我是已经复制粘贴过来并排版了，所以将就看看？）

如果不考虑一阶修正（推力与总冲与一楼一致，略去）：

$P_g(t)=P\left(\frac{V_0}{V_g(t)}\right)^\gamma$

$M(t)=m+\rho\left[V-V_g(t)\right]$

$\frac{dV_g}{dt} = A\sqrt{ \frac{2\left[p_g(t)-P_0\right]}{\rho} }$

$\frac{dv}{dt} = \frac{ 2A\left[p_g(t)-P_0\right] } { m+\rho\left[V-V_g(t)\right] } -g$

程序用龙格库塔法求解这个微分方程组即可获得更加精确的解，当然还需要结合喷水结束的条件，这个在后续写。

如果不忽略喷嘴面积，简略的推导如下：

$p_g(t) = P\left(\frac{V_0}{V_g(t)}\right)^\gamma$

$M(t) = m+\rho\left[V-V_g(t)\right]$

$L(t) = \frac{V-V_g(t)}{S}$

$u(t) = \sqrt{ \frac{ 2\left[p_g(t)-P_0\right] } { \rho\left[ 1-\left(\frac{A}{S}\right)^2 - \frac{2\rho A L(t)}{M(t)} \right] } }$

代入得，喷流速度为：

$u(t) = \sqrt{ \frac{ 2\left[ P\left(\frac{V_0}{V_g(t)}\right)^\gamma-P_0 \right] } { \rho\left[ 1-\left(\frac{A}{S}\right)^2 - \frac{ 2\rho A\left[V-V_g(t)\right] } { S\left\{ m+\rho\left[V-V_g(t)\right] \right\} } \right] } }$

又由于$\frac{dV_g}{dt} = A u(t)$，可以求得完整的微分方程组。

$\begin{cases} \displaystyle \frac{dV_g}{dt} = A \sqrt{ \frac{ 2\left[ P\left(\frac{V_0}{V_g}\right)^\gamma-P_0 \right] } { \rho\left[ 1-\left(\frac{A}{S}\right)^2 - \frac{ 2\rho A\left(V-V_g\right) } { S\left[m+\rho\left(V-V_g\right)\right] } \right] } }, \\[2.2em] \displaystyle \frac{dv}{dt} = \frac{ 2A\left[ P\left(\frac{V_0}{V_g}\right)^\gamma-P_0 \right] } { \left[m+\rho\left(V-V_g\right)\right] \left[ 1-\left(\frac{A}{S}\right)^2 - \frac{ 2\rho A\left(V-V_g\right) } { S\left[m+\rho\left(V-V_g\right)\right] } \right] } -g, \\[2.2em] \displaystyle \frac{dh}{dt}=v. \end{cases}$

初始条件略

接下来通过已经求解出的量求推力与总冲：

$T(t) = \rho A u^2(t)$

代入得，推力为：

$T(t) = \frac{ 2A\left[ P\left(\frac{V_0}{V_g(t)}\right)^\gamma-P_0 \right] } { 1-\left(\frac{A}{S}\right)^2 - \frac{ 2\rho A\left[V-V_g(t)\right] } { S\left\{ m+\rho\left[V-V_g(t)\right] \right\} } }$

总冲积分求解：

$I = \int_0^{t_b} \frac{ 2A\left[ P\left(\frac{V_0}{V_g(t)}\right)^\gamma-P_0 \right] } { 1-\left(\frac{A}{S}\right)^2 - \frac{ 2\rho A\left[V-V_g(t)\right] } { S\left\{ m+\rho\left[V-V_g(t)\right] \right\} } } \,dt$

这里需要注意结束条件为水全部喷完或者内外压差相等。

$V_{g,b} = \min\left[ V,\, V_0\left(\frac{P}{P_0}\right)^{1/\gamma} \right]$

既然有要用计算机数值求解，就最好直接求解微分方程组。

我记得我高中的时候做过一到物竞复赛模拟题，好像是苏国珍的，好像是再什么假设下可以求出解析解，我想找但是找不到了，干脆不找了，反正工程上只需要数值解。

我让GPT生成了一个计算程序，比较简陋，也没考虑开伞啥的，因为我现在没啥时间，你可以拿到程序后继续改改。

water_rocket.py19.83KBPY2次下载

water_rocket.exe39.53MBEXE0次下载

我使用你的数据计算了一下，发现理论高度大概25m，总冲17.6N·s，很接近，我开始还以为会有很大的差异。我最近没什么时间，现在也不太确定是计算器的问题还是就是这样，如果是就是这样的话，那么你的模型其实工程上完全可用，但是要考虑的是发射倾斜，内部流体阻力等等问题，这个就不确定了，楼主可以在此基础上进一步研究。我觉得空气阻力反而不是很大，因为速度很低。

这个程序是理想化的，而且没有考虑气体种类选择，牙管阻力等等，楼主可以试试只移植计算部分来修改，大概有10%到20%的精度提升，这样你的程序可能会更好。还可以给个建议，就是用多文件程序，这样结构更清晰，不过我这个也是单文件的程序。