Ballistics modeling with Cx(M) and ρ(h) in Python

Numerical Modeling of External Ballistics with Variable Air Density

Continuing the analysis of external ballistics problems, cases involving the dependence of the drag coefficient Cx on velocity and air density on altitude are examined. The fourth and fifth cases require numerical methods, as an analytical solution is impossible. The flight of a projectile is simulated with initial parameters: velocity 600 m/s, angle 55°, mass 50 kg.

The 1943 air resistance law defines Cx as a function of the Mach number M = v/a. Density ρ(h) is approximated by an exponential dependence ρ(h) = ρ₀ × exp(-h/H), where H = 8000 m is the scale height.

Dependence of Cx on Mach Number

Coefficient Cx is a piecewise function:

• 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)²

Drag force F = 0.5 ρ(h) v² S Cx (v/|v|), where S is the cross-sectional area (0.01 m²). Equations of motion:

ẋ = vx

ẏ = vy

v̇x = -Fx/m

v̇y = -g - Fy/m

g = 9.81 m/s², a = 343 m/s.

Implementation of Numerical Integration

To solve the system of ODEs, the 4th-order Runge-Kutta method (RK4) is applied. Integration step dt = 0.01 s, final time t_end = 98.4 s.

Full simulation code in Python:

import numpy as np
import matplotlib.pyplot as plt

# Constants
g = 9.81  # m/s²
m = 50.0  # body mass, kg
S = 0.01  # cross-sectional area, m²
rho0 = 1.225  # air density at sea level, kg/m³
H = 8000  # scale height, m
a = 343  # speed of sound, m/s
v0 = 600  # initial velocity, m/s
theta = np.radians(55)  # initial angle, rad
dt = 0.01  # integration step, s
t_end = 98.4  # final time, s

# Initial conditions
x, y = 0, 0
vx, vy = v0 * np.cos(theta), v0 * np.sin(theta)

# Arrays to store results
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)


# Function to compute 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


# Air density function
def rho(h):
    return rho0 * np.exp(-h / H)


# System of ODEs
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])


# 4th-order Runge-Kutta method
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


# Integration
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)

# Visualization
plt.figure(figsize=(10, 6))
plt.plot(X, Y, label='Trajectory')
plt.xlabel('X, m')
plt.ylabel('Y, m')
plt.title('Body Flight Trajectory')
plt.grid(True)
plt.legend()
plt.show()

# Output results
print(f"Flight range: {X[-1]:.2f} m")
print(f"Maximum altitude: {np.max(Y):.2f} m")

Simulation Results

Flight range: 29 km. Maximum altitude: 11.5 km. The trajectory accounts for the transition through supersonic speeds and air thinning at altitudes.

Key assumptions:

• Exponential density model (simplified; real atmosphere requires more accurate profiles).
• Quadratic drag without considering other effects (spin, wind).

| Parameter | Value |

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

| Range | 29,000 m |

| Max. altitude | 11,500 m |

| Flight time | 98.4 s |

| Initial velocity | 600 m/s |

Key Takeaways

• Numerical methods are essential for realistic models with variable Cx(M) and ρ(h).
• RK4 ensures accuracy with dt=0.01 s for ballistic tasks.
• The 1943 law accurately models the transonic regime (M≈0.8–1.0).
• Exponential ρ(h) provides a reasonable approximation up to 15 km.

— Editorial Team