RocketPy 作为一款专业的开源火箭轨迹仿真工具，其核心价值在于对火箭飞行物理过程的精确数学建模。本文将系统梳理 RocketPy 中涉及的所有关键公式，涵盖从基础运动方程到复杂气动计算的完整体系，为深入理解和使用该工具提供理论支撑。
一、坐标系定义与变换公式
1.1 参考坐标系定义
RocketPy 采用两套主要坐标系来描述火箭运动：
惯性坐标系 A（地球固定坐标系）：
-
基向量：{a1,a2,a3}
-
a3指向地心（垂直向下）
-
a1和 a2在水平面内正交
物体坐标系 B（火箭固连坐标系）：
-
基向量：{b1,b2,b3}
-
b3沿火箭纵轴指向头部
-
原点位于火箭干质量中心
1.2 位置与速度表示
火箭总质心 ∗（包含刚性体 R、临时固体相 T和气体相 G）的位置和速度在惯性系中表示为：
Ar∗=xa1+ya2+za3
Av∗=dtdxa1+dtdya2+dtdza3
1.3 坐标系变换
使用欧拉参数（四元数）进行坐标系变换，变换矩阵为：
CAB=e02+e12−e22−e322(e1e2+e0e3)2(e1e3−e0e2)2(e1e2−e0e3)e02−e12+e22−e322(e2e3+e0e1)2(e1e3+e0e2)2(e2e3−e0e1)e02−e12−e22+e32
其中欧拉参数满足约束条件：e02+e12+e22+e32=1
二、运动方程核心公式
2.1 平移运动方程
火箭质心加速度方程（导轨阶段）：
mv˙+rCM′′=T−2m˙rCM′+m¨(rnoz−rCM)+A−(mga^3⋅b^3)b^3
其中：
-
m：时变质量
-
v：质心速度
-
rCM：质心相对于干质量中心的位置向量
-
rnoz：喷管出口相对于干质量中心的位置向量
-
T：推力向量
-
A：轴向气动力
-
撇号(')表示在物体坐标系B中的导数
2.2 旋转运动方程
刚体旋转动力学方程：
I⋅ω˙+ω×(I⋅ω)=Mtotal
其中：
-
I：时变惯性张量
-
ω：角速度向量
-
Mtotal：总力矩（气动力矩+推力力矩）
欧拉参数运动学方程：
e˙0e˙1e˙2e˙3=21−e1e0e3−e2−e2−e3e0e1−e3e2−e1e0ω1ω2ω3
三、力与力矩模型公式
3.1 推力计算
固体发动机推力：
T(t)=m˙prop⋅ve+(pe−pa)⋅Ae
其中：
-
m˙prop：推进剂质量流率
-
ve：排气速度
-
pe：喷管出口压力
-
pa：环境压力
-
Ae：喷管出口面积
推力向量方向：
假设推力沿火箭纵轴方向：
T=T(t)⋅b^3
3.2 气动力模型
轴向力公式：
A=21CDρ∥V∞∥2S(−b^3)
其中：
-
CD：阻力系数（马赫数函数）
-
ρ：空气密度
-
V∞：火箭相对于风的相对速度
-
S：参考面积
相对速度计算：
V∞=W(r)−v
W(r)是位置相关的风速向量。
3.3 升力面法向力
对于第 i个升力面：
压力中心速度：
vi=v+ω×ri
相对风速：
V∞,i=W(r+ri)−v−ω×ri
攻角计算：
αi=cos−1(∥V∞,i∥−V∞,i⋅b^3)
法向力公式：
Ni=CL(αi)S21ρ∥V∞,i∥2∥b^3×V∞,i∥b^3×V∞,i
3.4 重力模型
地球重力加速度：
g(h)=g0(Re+hRe)2
其中：
-
g0=9.80665m/s2（海平面重力加速度）
-
Re=6378137m（地球赤道半径）
-
h：海拔高度
重力向量：
Fg=−mg(h)a^3
四、质量特性与时变参数
4.1 质量变化方程
总质量时变函数：
m(t)=mdry+mprop(t)
推进剂质量消耗：
对于固体发动机：
m˙prop(t)=ρprop⋅Ab(t)⋅rb(t)
其中：
-
ρprop：推进剂密度
-
Ab(t)：燃烧面积（时间函数）
-
rb(t)：燃速（压力函数）
4.2 质心位置变化
总质心位置：
rCM(t)=m(t)mdryrdry+mprop(t)rprop(t)
质心速度（物体坐标系）：
rCM′=dtdrCMB
4.3 惯性张量变化
瞬时惯性张量：
I(t)=Idry+Iprop(t)
对于圆柱形推进剂药柱：
Iprop,xx(t)=Iprop,yy(t)=121mprop(t)(3R2+h2)+mprop(t)d2
Iprop,zz(t)=21mprop(t)R2
五、空气动力学系数模型
5.1 阻力系数模型
马赫数相关阻力系数：
CD=CD,0+CD,2M2+CD,6+M4CD,4M4
或通过查表插值获得。
基于Barrowman稳定性导数（用于法向力系数）：
CNα=i=1∑n(CNα)i
对于圆锥形头部：
(CNα)nose=2
对于圆柱体：
(CNα)body=0
对于梯形翼面：
(CNα)fin=1+1+(cr+ct2l)24N(ds)2
5.2 升力系数曲线
线性区域（小攻角）：
CL(α)=CLα⋅α
非线性区域（大攻角）：
CL(α)=CL,maxsin(2α)(经验公式)
或使用表格插值方法。
六、环境模型公式
6.1 国际标准大气模型（1976）
温度剖面：
T(h)=⎩⎨⎧T0+L0hT11T11+L1(h−20)⋮0≤h≤11km11km<h≤20km20km<h≤32km⋮
压力计算：
对流层（0-11 km）：
p(h)=p0(T(h)T0)R∗L0g0M
平流层（11-20 km）：
p(h)=p11exp(−R∗T11g0M(h−11000))
6.2 风速模型
指数风剖面：
W(z)=Wref(zrefz)α
风切变模型：
W(x,y,z)=Wx(z)Wy(z)0+Wgust(t)
6.3 密度计算
理想气体状态方程：
ρ(h)=RspecificT(h)p(h)
其中 Rspecific=287.058J/(kg⋅K)（干空气气体常数）
七、多体动力学与分离事件
7.1 多级火箭分离
分离前总质量：
mtotal,before=mpayload+i=1∑nmstage,i
分离后质量：
mtotal,after=mpayload+i=2∑nmstage,i
分离冲量：
Δv=mstage1FsepΔt
7.2 开伞动力学
降落伞阻力公式：
Dparachute=21CD,paraρAparav2
降落伞展开过程：
dtdApara={kinflationt00≤t≤tfullt>tfull
八、数值积分方法
8.1 状态向量定义
RocketPy 使用13维状态向量：
y=[x,y,z,x˙,y˙,z˙,e0,e1,e2,e3,ω1,ω2,ω3]T
8.2 微分方程系统
dtdy=f(t,y)
其中：
f=x˙y˙z˙x¨y¨z¨e˙0e˙1e˙2e˙3ω˙1ω˙2ω˙3
8.3 积分器选择
RocketPy 支持多种ODE求解器：
-
LSODA：自动在Adams（非刚性）和BDF（刚性）方法间切换
-
RK45：显式Runge-Kutta方法（4阶）
-
DOP853：高精度8阶Runge-Kutta方法
误差控制：
局部截断误差满足：
error≤rtol⋅∣y∣+atol
其中 rtol 为相对容差，atol 为绝对容差。
九、特殊飞行阶段公式
9.1 导轨约束阶段
在导轨长度 Lrail内，火箭运动受约束：
v⋅b^3=vrail(沿导轨方向)
ω=0(姿态固定)
导轨反作用力：
Rrail=(mga^3⋅n^rail)n^rail−(T+A)⋅n^rail
其中 n^rail为导轨法向。
9.2 助推段与滑行段
发动机关机条件：
mprop(t)≤mprop,min
或
t≥tburn
滑行段动力学：
推力 T=0，仅受气动力和重力作用。
9.3 再入与回收
弹道系数：
β=CDSm
过载系数：
n=g0∥a∥
其中 a为总加速度。
十、验证与精度控制公式
10.1 能量守恒验证
机械能变化率：
dtdE=T⋅v+A⋅v−m˙prop(21ve2+ρproppe−pa)
10.2 动量守恒验证
系统总动量变化：
dtdP=Fext+m˙propve
10.3 数值稳定性条件
CFL条件（显式方法）：
Δt≤CvmaxΔx
刚度检测：
stiffness=min∣λi∣max∣λi∣
其中 λi为雅可比矩阵特征值。
总结与应用建议
RocketPy 的公式体系体现了现代火箭动力学的完整数学模型，具有以下特点：
物理完整性：涵盖6自由度刚体动力学、变质量效应、复杂气动耦合
数值稳健性：采用自适应步长积分和误差控制
模块化设计：各物理模型可独立验证和替换
可扩展性：支持用户自定义力和力矩模型
实际应用建议：
-
对于初步设计，可使用简化气动模型（Barrowman方法）
-
对于高精度仿真，应导入CFD计算的气动系数
-
蒙特卡洛分析时，注意关键参数（如 CD、燃速）的概率分布
-
验证仿真结果时，对比能量、动量守恒关系
这些公式构成了RocketPy仿真的数学基础，理解这些公式有助于用户正确配置参数、解释结果，并在需要时扩展功能。官方文档和学术论文提供了更详细的推导和验证数据，建议深入阅读以掌握完整理论体系。
