Python 中的 Cx(M) 和 ρ(h) 弹道建模

考虑空气密度变化的外弹道数值模拟

继续分析外弹道问题，本文探讨了阻力系数Cx随速度变化、空气密度随高度变化的情况。第四和第五种情况需要数值方法，因为无法获得解析解。模拟弹丸飞行，初始参数为：速度600米/秒，角度55°，质量50千克。

1943年空气阻力定律将Cx定义为马赫数M = v/a的函数。密度ρ(h)采用指数近似：ρ(h) = ρ₀ × exp(-h/H)，其中H = 8000米为标高。

Cx随马赫数的变化

系数Cx为分段函数：

• M ≤ 0.73：Cx = 0.157
• 0.73 < M < 0.82：Cx = 0.033M + 0.133
• 0.82 ≤ M < 0.91：Cx = 0.161 + 3.9(M - 0.823)²
• 0.91 ≤ M < 1.00：Cx = 1.5M - 1.176
• M ≥ 1.00：Cx = 0.384 - 1.6(M - 1.176)²

阻力F = 0.5 ρ(h) v² S Cx (v/|v|)，其中S为横截面积（0.01平方米）。运动方程：

ẋ = vx

ẏ = vy

v̇x = -Fx/m

v̇y = -g - Fy/m

g = 9.81米/秒²，a = 343米/秒。

数值积分实现

为求解常微分方程组，采用四阶龙格-库塔法（RK4）。积分步长dt = 0.01秒，终止时间t_end = 98.4秒。

完整Python模拟代码：

import numpy as np
import matplotlib.pyplot as plt

# 常量
g = 9.81  # 重力加速度，米/秒²
m = 50.0  # 物体质量，千克
S = 0.01  # 横截面积，平方米
rho0 = 1.225  # 海平面空气密度，千克/立方米
H = 8000  # 标高，米
a = 343  # 声速，米/秒
v0 = 600  # 初始速度，米/秒
theta = np.radians(55)  # 初始角度，弧度
dt = 0.01  # 积分步长，秒
t_end = 98.4  # 终止时间，秒

# 初始条件
x, y = 0, 0
vx, vy = v0 * np.cos(theta), v0 * np.sin(theta)

# 存储结果的数组
t = np.arange(0, t_end, dt)
X, Y = np.zeros_like(t), np.zeros_like(t)
VX, VY = np.zeros_like(t), np.zeros_like(t)


# 计算cx(M)的函数
def cx(M):
    if M <= 0.73:
        return 0.157
    elif M < 0.82:
        return 0.033 * M + 0.133
    elif M < 0.91:
        return 0.161 + 3.9 * (M - 0.823) ** 2
    elif M < 1.00:
        return 1.5 * M - 1.176
    else:
        return 0.384 - 1.6 * (M - 1.176) ** 2


# 空气密度函数
def rho(h):
    return rho0 * np.exp(-h / H)


# 常微分方程组
def f(state, t):
    x, y, vx, vy = state
    v = np.sqrt(vx ** 2 + vy ** 2)
    M = v / a
    rho_h = rho(y)
    cx_val = cx(M)

    Fx = 0.5 * rho_h * v ** 2 * S * cx_val * vx / v
    Fy = 0.5 * rho_h * v ** 2 * S * cx_val * vy / v

    dxdt = vx
    dydt = vy
    dvxdt = -Fx / m
    dvydt = -g - Fy / m

    return np.array([dxdt, dydt, dvxdt, dvydt])


# 四阶龙格-库塔法
def rk4_step(state, t, dt):
    k1 = dt * f(state, t)
    k2 = dt * f(state + 0.5 * k1, t + 0.5 * dt)
    k3 = dt * f(state + 0.5 * k2, t + 0.5 * dt)
    k4 = dt * f(state + k3, t + dt)
    return state + (k1 + 2 * k2 + 2 * k3 + k4) / 6


# 积分
state = np.array([x, y, vx, vy])
for i, t_val in enumerate(t):
    X[i], Y[i], VX[i], VY[i] = state
    state = rk4_step(state, t_val, dt)

# 可视化
plt.figure(figsize=(10, 6))
plt.plot(X, Y, label='弹道轨迹')
plt.xlabel('X, 米')
plt.ylabel('Y, 米')
plt.title('物体飞行轨迹')
plt.grid(True)
plt.legend()
plt.show()

# 输出结果
print(f"飞行距离：{X[-1]:.2f} 米")
print(f"最大高度：{np.max(Y):.2f} 米")

模拟结果

飞行距离：29公里。最大高度：11.5公里。该轨迹考虑了超音速过渡和高空空气稀薄的影响。

关键假设：

• 指数密度模型（简化；真实大气需要更精确的剖面）。
• 二次阻力，未考虑其他效应（旋转、风）。

| 参数 | 数值 |

|-----------|-------|

| 射程 | 29,000 米 |

| 最大高度 | 11,500 米 |

| 飞行时间 | 98.4 秒 |

| 初始速度 | 600 米/秒 |

关键要点

• 数值方法至关重要：对于包含变量Cx(M)和ρ(h)的逼真模型。
• RK4确保精度：在弹道任务中，dt=0.01秒的步长下。
• 1943年定律：准确模拟跨音速区域（M≈0.8–1.0）。
• 指数ρ(h)：在15公里以下提供合理近似。

— Editorial Team