Modelado balístico con Cx(M) y ρ(h) en Python

Modelado Numérico de la Balística Externa con Densidad del Aire Variable

Continuando con el análisis de problemas de balística externa, se examinan casos que involucran la dependencia del coeficiente de arrastre Cx con la velocidad y la densidad del aire con la altitud. Los casos cuarto y quinto requieren métodos numéricos, ya que una solución analítica es imposible. Se simula el vuelo de un proyectil con parámetros iniciales: velocidad 600 m/s, ángulo 55°, masa 50 kg.

La ley de resistencia del aire de 1943 define Cx como una función del número de Mach M = v/a. La densidad ρ(h) se aproxima mediante una dependencia exponencial ρ(h) = ρ₀ × exp(-h/H), donde H = 8000 m es la altura de escala.

Dependencia de Cx con el Número de Mach

El coeficiente Cx es una función por tramos:

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

La fuerza de arrastre F = 0.5 ρ(h) v² S Cx (v/|v|), donde S es el área transversal (0.01 m²). Ecuaciones de movimiento:

ẋ = vx

ẏ = vy

v̇x = -Fx/m

v̇y = -g - Fy/m

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

Implementación de la Integración Numérica

Para resolver el sistema de EDOs, se aplica el método de Runge-Kutta de cuarto orden (RK4). Paso de integración dt = 0.01 s, tiempo final t_end = 98.4 s.

Código completo de simulación en Python:

import numpy as np
import matplotlib.pyplot as plt

# Constantes
g = 9.81  # m/s²
m = 50.0  # masa del cuerpo, kg
S = 0.01  # área transversal, m²
rho0 = 1.225  # densidad del aire a nivel del mar, kg/m³
H = 8000  # altura de escala, m
a = 343  # velocidad del sonido, m/s
v0 = 600  # velocidad inicial, m/s
theta = np.radians(55)  # ángulo inicial, rad
dt = 0.01  # paso de integración, s
t_end = 98.4  # tiempo final, s

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

# Arreglos para almacenar resultados
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)


# Función para calcular 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


# Función de densidad del aire
def rho(h):
    return rho0 * np.exp(-h / H)


# Sistema de EDOs
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])


# Método de Runge-Kutta de cuarto orden
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


# Integración
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)

# Visualización
plt.figure(figsize=(10, 6))
plt.plot(X, Y, label='Trayectoria')
plt.xlabel('X, m')
plt.ylabel('Y, m')
plt.title('Trayectoria de Vuelo del Cuerpo')
plt.grid(True)
plt.legend()
plt.show()

# Resultados
print(f"Alcance de vuelo: {X[-1]:.2f} m")
print(f"Altitud máxima: {np.max(Y):.2f} m")

Resultados de la Simulación

Alcance de vuelo: 29 km. Altitud máxima: 11.5 km. La trayectoria considera la transición a velocidades supersónicas y el enrarecimiento del aire a altitud.

Supuestos clave:

• Modelo de densidad exponencial (simplificado; la atmósfera real requiere perfiles más precisos).
• Arrastre cuadrático sin considerar otros efectos (giro, viento).

| Parámetro | Valor |

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

| Alcance | 29,000 m |

| Altitud máx. | 11,500 m |

| Tiempo de vuelo | 98.4 s |

| Velocidad inicial | 600 m/s |

Conclusiones Clave

• Los métodos numéricos son esenciales para modelos realistas con Cx(M) y ρ(h) variables.
• RK4 garantiza precisión con dt=0.01 s para tareas balísticas.
• La ley de 1943 modela con precisión el régimen transónico (M≈0.8–1.0).
• ρ(h) exponencial proporciona una aproximación razonable hasta 15 km.

— Editorial Team