3、导弹制导模型基础：导弹六自由度运动方程、气动参数建模、推力与质量模型
各位同学，欢迎来到第三章。这一章，我们要啃一块硬骨头——导弹的数学模型。
说实话，做硬件在环仿真，最核心的就是模型。模型不准，你仿真做得再花哨，也是白搭。我见过太多团队，花大价钱买了仿真设备，结果模型里气动系数都是拍脑袋估的，最后半实物跑起来跟实弹数据对不上，那叫一个尴尬。
所以这一章，咱们把基础打牢。我会从六自由度运动方程讲起，再到气动参数怎么建模，最后说说推力和质量的变化。嗯，内容不少，但都是干货。
3.1 导弹六自由度运动方程
六自由度，说白了就是三个平动自由度（前后、左右、上下）加上三个转动自由度（俯仰、偏航、滚转）。你想想看，导弹在空中飞，不就是这六个量在变吗？
我个人习惯把方程分成两组：
- 质心运动方程——描述导弹飞到哪里去
- 绕质心转动方程——描述导弹姿态怎么变
3.1.1 坐标系选择
做仿真第一步，先定坐标系。我建议用下面这套：
| 坐标系 | 原点 | 说明 |
|---|---|---|
| 地面坐标系 | 发射点 | 惯性系，忽略地球自转 |
| 弹体坐标系 | 导弹质心 | 固连在弹体上，x轴指向弹头 |
| 速度坐标系 | 导弹质心 | x轴沿速度方向 |
这里有个坑——坐标系转换矩阵千万别搞反了。我曾经在项目里因为欧拉角旋转顺序写错，仿真出来的导弹直接往反方向飞，查了三天bug才发现。
3.1.2 质心运动方程
在地面坐标系下，质心运动方程长这样：
m * dV/dt = P + A + G其中：
- m 是导弹质量（随时间变化）
- V 是速度矢量
- P 是推力矢量
- A 是气动力矢量
- G 是重力矢量
展开成标量形式，就是三个方向的力平衡方程。我一般写成：
m * du/dt = P_x + A_x - m * g * sin(θ)
m * dv/dt = P_y + A_y + m * g * cos(θ) * sin(φ)
m * dw/dt = P_z + A_z + m * g * cos(θ) * cos(φ)这里 u、v、w 是速度在弹体系下的三个分量。θ 是俯仰角，φ 是滚转角。
3.1.3 绕质心转动方程
转动方程用的是动量矩定理：
J * dω/dt + ω × (J * ω) = MJ 是转动惯量矩阵，ω 是角速度矢量，M 是外力矩。
展开后：
J_x * dp/dt + (J_z - J_y) * q * r = M_x
J_y * dq/dt + (J_x - J_z) * p * r = M_y
J_z * dr/dt + (J_y - J_x) * p * q = M_zp、q、r 分别是滚转、俯仰、偏航角速度。
3.2 气动参数建模
气动参数，说白了就是导弹在空气中飞行时受到的力和力矩。这部分最头疼，因为气动数据通常来自风洞实验或CFD计算，而我们做仿真时往往拿不到完整的数据表。
怎么办？我的做法是——用系数法。
3.2.1 气动力系数
气动力可以写成：
A_x = C_x * q * S
A_y = C_y * q * S
A_z = C_z * q * S其中 q = 0.5 * ρ * V² 是动压，S 是参考面积。
C_x、C_y、C_z 是气动力系数，它们通常是攻角 α、侧滑角 β 和马赫数 Ma 的函数。
我常用的建模方式有两种：
- 查表法——把风洞数据做成二维/三维表格，仿真时插值
- 解析法——用多项式拟合气动系数
我个人更推荐查表法。为什么呢？因为解析法容易在边界处跑飞。我记得有一次，用二次多项式拟合升力系数，结果在大攻角时系数变成了负值，仿真出来的导弹直接失速了。
3.2.2 气动力矩系数
力矩系数类似：
M_x = C_mx * q * S * L
M_y = C_my * q * S * L
M_z = C_mz * q * S * LL 是参考长度，通常取弹体直径。
这里有个关键点——阻尼力矩。没有阻尼力矩，导弹的角运动会一直振荡，永远不收敛。阻尼力矩系数通常与角速度有关：
C_mz = C_mz_α(α, Ma) + C_mz_q * (q * L / V)C_mz_q 就是俯仰阻尼导数，一般是个负值。
3.3 推力与质量模型
推力模型相对简单，但容易出错的地方也不少。
3.3.1 推力模型
固体火箭发动机的推力通常用推力-时间曲线描述：
P(t) = P_profile(t)这个曲线一般由发动机厂家提供。如果没有，可以用简化模型：
P(t) = P_max * exp(-t / τ)τ 是时间常数，P_max 是最大推力。
但要注意，推力方向必须沿着弹体轴线。很多新手直接把推力当成标量加在速度方向上，这是错的。
3.3.2 质量模型
质量变化来自燃料消耗：
dm/dt = -m_dot(t)m_dot 是质量流率，通常与推力成正比：
m_dot = P(t) / (I_sp * g_0)I_sp 是比冲，g_0 是海平面重力加速度。
质心位置也会变。燃料烧掉后，质心会前移。这个变化对稳定性影响很大。我一般用线性插值：
x_cg(t) = x_cg_0 + (x_cg_f - x_cg_0) * (m_0 - m(t)) / (m_0 - m_f)3.4 模型集成与仿真流程
好了，现在我们把所有方程串起来。仿真流程是这样的：
- 初始化状态量（位置、速度、姿态、角速度、质量）
- 计算当前时刻的推力、质量流率
- 根据当前状态计算气动力和力矩
- 代入六自由度方程，求导
- 用数值积分（我常用四阶龙格-库塔）更新状态
- 更新时间步，回到第2步
代码框架大概这样：
def missile_sim_step(t, state):
# 解包状态
x, y, z, u, v, w, phi, theta, psi, p, q, r, m = state
# 计算推力
P = get_thrust(t)
m_dot = P / (I_sp * g_0)
# 计算气动力
V = sqrt(u**2 + v**2 + w**2)
alpha = atan2(w, u)
beta = asin(v / V)
Ma = V / speed_of_sound
Cx, Cy, Cz = get_aero_coeffs(alpha, beta, Ma)
q_dyn = 0.5 * rho * V**2
Ax = Cx * q_dyn * S
Ay = Cy * q_dyn * S
Az = Cz * q_dyn * S
# 计算力矩（略）
# 六自由度方程求导
du_dt = (P + Ax) / m - g * sin(theta)
# ... 其他导数
return [dx_dt, dy_dt, dz_dt, du_dt, ...]好了，这一章的内容就到这里。六自由度方程是基础，气动参数是灵魂，推力质量是动力源。三者缺一不可。
下一章，我们会把这些模型放到硬件在环仿真平台上跑起来。到时候你会发现，模型再精确，跟硬件一对接，还是会有各种意想不到的问题。嗯，那才是真正考验功力的时候。
公众号:蓝海资料掘金营，微信deep3321