2. 建立火箭动力学模型
任务说明
本案例实现一个二维平面内的火箭动力学模型，包含推力、重力和空气动力（升力与阻力）的作用。模型输出火箭在惯性坐标系下的加速度分量，并基于给定时间步长更新位置与速度状态。该模型可作为后续制导算法仿真的被控对象。
坐标系与受力定义
采用地面惯性坐标系（x水平向右，y竖直向上）。火箭状态由位置 (x, y)、速度 (vx, vy)、质量 m 和俯仰角 θ（速度方向与水平面的夹角）描述。受力包括：
- 推力 T：沿火箭纵轴方向，与速度方向夹角为攻角 α，因此推力方向角为 θ+α；
- 重力：竖直向下，大小为 mg；
- 空气动力：分解为阻力 D（与速度反向）和升力 L（垂直于速度方向，向上为正）。
加速度由牛顿第二定律计算，状态更新采用一阶欧拉积分。
代码实现
import numpy as np
class Rocket2D:
"""二维火箭动力学模型（推力、重力、空气动力）"""
def __init__(self, x0=0.0, y0=0.0, vx0=0.0, vy0=0.0,
m0=500.0, g=9.81, rho=1.225, S=0.25,
Cd=0.3, Cl_alpha=2.0):
self.x = x0
self.y = y0
self.vx = vx0
self.vy = vy0
self.m = m0
self.g = g # 重力加速度
self.rho = rho # 空气密度
self.S = S # 参考面积
self.Cd = Cd # 阻力系数
self.Cl_alpha = Cl_alpha # 升力线斜率
def speed(self):
return np.hypot(self.vx, self.vy)
def flight_path_angle(self):
"""速度方向角（与水平面夹角）"""
return np.arctan2(self.vy, self.vx)
def aerodynamic_forces(self):
"""计算阻力D和升力L（N）"""
v = self.speed()
q = 0.5 * self.rho * v**2 # 动压
D = q * self.S * self.Cd
# 升力系数与攻角成正比（小攻角假设）
alpha = 0.0 # 由外部传入，此处仅返回函数
return D, q * self.S * self.Cl_alpha * alpha
def derivatives(self, T, alpha):
"""计算加速度分量 ax, ay 及质量变化率 dmdt"""
v = self.speed()
theta = self.flight_path_angle()
# 空气动力
D = 0.5 * self.rho * v**2 * self.S * self.Cd
L = 0.5 * self.rho * v**2 * self.S * self.Cl_alpha * alpha
# 推力分量
thrust_angle = theta + alpha
Tx = T * np.cos(thrust_angle)
Ty = T * np.sin(thrust_angle)
# 阻力方向（与速度反向）
if v > 1e-6:
Dx = -D * self.vx / v
Dy = -D * self.vy / v
else:
Dx = Dy = 0.0
# 升力方向（垂直于速度，逆时针旋转90度）
if v > 1e-6:
Lx = -L * self.vy / v
Ly = L * self.vx / v
else:
Lx = Ly = 0.0
# 合力
Fx = Tx + Dx + Lx
Fy = Ty + Dy + Ly - self.m * self.g
ax = Fx / self.m
ay = Fy / self.m
# 质量变化率（假设质量流率恒定，由外部给定）
dmdt = 0.0 # 默认无质量变化
return ax, ay, dmdt
def step(self, T, alpha, dt, mdot=0.0):
"""使用欧拉法更新状态"""
ax, ay, dmdt = self.derivatives(T, alpha)
# 更新速度
self.vx += ax * dt
self.vy += ay * dt
# 更新位置
self.x += self.vx * dt
self.y += self.vy * dt
# 更新质量
self.m += mdot * dt
if self.m < 0:
self.m = 0.0
return self.x, self.y, self.vx, self.vy, self.m
# 示例：初始化火箭并执行一步仿真
if __name__ == "__main__":
# 初始状态：位于原点，静止，质量500kg
rocket = Rocket2D(x0=0.0, y0=0.0, vx0=0.0, vy0=0.0, m0=500.0)
# 推力10000N，攻角0度，时间步长0.1s
T = 10000.0
alpha = 0.0
dt = 0.1
# 执行一步
x, y, vx, vy, m = rocket.step(T, alpha, dt)
print(f"位置: ({x:.3f}, {y:.3f}) m")
print(f"速度: ({vx:.3f}, {vy:.3f}) m/s")
print(f"质量: {m:.1f} kg")关键说明
代码中 aerodynamic_forces 方法仅用于展示升力系数与攻角的关系，实际计算在 derivatives 中完成。攻角 α 由外部制导律给出，本模型不包含攻角动力学。质量变化率通过 mdot 参数传入 step 方法，默认无质量消耗。
验证方法
运行上述示例，预期输出为：
位置: (0.000, 0.000) m
速度: (2.000, 0.000) m/s
质量: 500.0 kg验证逻辑：初始静止，推力水平（α=0），无空气动力（v=0），因此水平加速度 a_x = T/m = 10000/500 = 20 m/s²，经过 0.1s 后 vx = 2.0 m/s，位移 = 0.5×20×0.01 = 0.1m，但欧拉法使用更新前速度计算位移，故 x = 0 + 0×0.1 = 0。竖直方向仅受重力，vy 初始为 0，ay = -g，因此 vy = -0.981 m/s，但示例输出显示 vy=0，这是因为代码中 step 先更新速度再更新位置，且示例中 vy0=0，经过一步后 vy = -0.981 m/s，与预期不符。修正示例验证：
# 正确验证代码
rocket = Rocket2D(x0=0.0, y0=0.0, vx0=0.0, vy0=0.0, m0=500.0)
T = 10000.0
alpha = 0.0
dt = 0.1
x, y, vx, vy, m = rocket.step(T, alpha, dt)
print(f"vx = {vx:.3f} m/s (预期 2.000)")
print(f"vy = {vy:.3f} m/s (预期 -0.981)")
print(f"x = {x:.3f} m (预期 0.000)")
print(f"y = {y:.3f} m (预期 -0.049)")判定标准：vx 接近 2.0 m/s，vy 接近 -0.981 m/s，x 保持 0，y 约为 -0.049 m（因欧拉法使用旧速度计算位移，y 位移 = vy0×dt = 0，但 vy 更新后为 -0.981，实际 y 更新使用更新后的 vy，因此 y = -0.981×0.1 = -0.0981 m。修正预期：y = -0.098 m）。
运行与验证
将代码保存为 rocket_dynamics.py，在 Python 3.8+ 环境中运行（需安装 NumPy）。执行 python rocket_dynamics.py 后，观察输出是否符合上述预期值。若 vx 与 vy 的计算结果与手算一致（误差小于 0.001），则模型实现正确。可通过修改初始速度、推力或攻角，进一步验证空气动力项是否按预期方向作用。