📚 课程目录第4章导弹质点弹道仿真

4. 导弹质点弹道仿真

任务定义

建立导弹质点运动模型，在给定初始速度、弹道倾角和发射高度条件下，对无控飞行弹道进行数值仿真，输出射程与高度随时间的变化曲线，并提取最大射程、最大高度及对应时刻作为弹道特征参数。

质点运动模型

将导弹视为质量不变的质点，忽略姿态动力学、气动控制及地球曲率影响，在二维铅垂平面内建立运动方程。采用地面坐标系，x为水平射程，y为高度，v为速度，θ为弹道倾角（速度与水平面夹角）。作用力包括重力mg和与速度平方成正比的阻力D = 0.5ρv²SC_D，其中ρ为空气密度（取常数1.225 kg/m³），S为参考面积，C_D为阻力系数。运动方程组如下：

dx/dt = v·cosθ
dy/dt = v·sinθ
dv/dt = -(D/m) - g·sinθ
dθ/dt = -(g·cosθ)/v

仿真终止条件为y ≤ 0（落地）。采用四阶Runge-Kutta法进行数值积分，时间步长取0.01 s。

仿真代码实现

import numpy as np
import matplotlib.pyplot as plt

# 物理参数
m = 500.0          # 质量 kg
S = 0.25           # 参考面积 m^2
CD = 0.35          # 阻力系数
rho = 1.225        # 空气密度 kg/m^3
g = 9.81           # 重力加速度 m/s^2

# 初始条件
v0 = 300.0         # 初速 m/s
theta0 = np.radians(45.0)   # 弹道倾角 rad
x0 = 0.0
y0 = 1000.0        # 发射高度 m

# 仿真参数
dt = 0.01          # 时间步长 s
t_max = 200.0      # 最大仿真时间 s

def derivatives(state):
    """计算状态导数 [dx/dt, dy/dt, dv/dt, dtheta/dt]"""
    x, y, v, theta = state
    D = 0.5 * rho * v**2 * S * CD
    dxdt = v * np.cos(theta)
    dydt = v * np.sin(theta)
    dvdt = -(D / m) - g * np.sin(theta)
    dthetadt = -(g * np.cos(theta)) / v
    return np.array([dxdt, dydt, dvdt, dthetadt])

def rk4_step(state, dt):
    """四阶Runge-Kutta单步积分"""
    k1 = derivatives(state)
    k2 = derivatives(state + 0.5 * dt * k1)
    k3 = derivatives(state + 0.5 * dt * k2)
    k4 = derivatives(state + dt * k3)
    return state + (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4)

# 数值积分主循环
state = np.array([x0, y0, v0, theta0])
t = 0.0
trajectory = [(t, state.copy())]

while t < t_max and state[1] > 0.0:
    state = rk4_step(state, dt)
    t += dt
    trajectory.append((t, state.copy()))

# 提取结果
t_arr = np.array([item[0] for item in trajectory])
x_arr = np.array([item[1][0] for item in trajectory])
y_arr = np.array([item[1][1] for item in trajectory])
v_arr = np.array([item[1][2] for item in trajectory])
theta_arr = np.array([item[1][3] for item in trajectory])

# 弹道特征参数
max_range = x_arr[-1]                    # 落地时射程
max_height = np.max(y_arr)               # 最大高度
t_max_height = t_arr[np.argmax(y_arr)]   # 达到最大高度时刻
impact_time = t_arr[-1]                  # 落地时刻

print(f"落地时刻: {impact_time:.2f} s")
print(f"最大射程: {max_range:.2f} m")
print(f"最大高度: {max_height:.2f} m (t={t_max_height:.2f} s)")
print(f"落地速度: {v_arr[-1]:.2f} m/s")

# 绘制弹道曲线
plt.figure(figsize=(10, 6))
plt.plot(x_arr, y_arr, 'b-', linewidth=1.5)
plt.xlabel('射程 x (m)')
plt.ylabel('高度 y (m)')
plt.title('导弹无控飞行弹道')
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.savefig('trajectory.png', dpi=150)
plt.show()

运行与验证

运行方式：将上述代码保存为 missile_trajectory.py，在已安装NumPy和Matplotlib的Python环境中执行 python missile_trajectory.py。

验证方法：程序输出弹道特征参数并显示弹道曲线图。判定标准如下：

弹道曲线从(0, 1000)出发，先上升后下降，最终落地时y=0，曲线连续光滑无跳变；
最大高度出现在弹道上升段结束处，对应时刻弹道倾角θ过零；
落地时射程为正值，落地速度大于初始速度的50%（因重力加速作用）；
将阻力系数CD设为0时，射程应显著增大，验证阻力模型生效。

预期输出（基于给定初始条件）：落地时刻约33.5 s，最大射程约4200 m，最大高度约2400 m（出现在约12 s处）。该结果由上述代码与输入参数直接唯一确定，可作为代码正确性的参考判据。