Modélisation balistique avec Cx(M) et ρ(h) en Python

Modélisation numérique de la balistique externe avec densité d'air variable

Poursuivant l'analyse des problèmes de balistique externe, nous examinons ici les cas impliquant la dépendance du coefficient de traînée Cx à la vitesse et de la densité de l'air à l'altitude. Les quatrième et cinquième cas nécessitent des méthodes numériques, car une solution analytique est impossible. Le vol d'un projectile est simulé avec les paramètres initiaux : vitesse 600 m/s, angle 55°, masse 50 kg.

La loi de résistance de l'air de 1943 définit Cx comme une fonction du nombre de Mach M = v/a. La densité ρ(h) est approximée par une dépendance exponentielle ρ(h) = ρ₀ × exp(-h/H), où H = 8000 m est l'échelle de hauteur.

Dépendance de Cx au nombre de Mach

Le coefficient Cx est une fonction par morceaux :

• 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 force de traînée F = 0,5 ρ(h) v² S Cx (v/|v|), où S est la section transversale (0,01 m²). Équations du mouvement :

ẋ = vx

ẏ = vy

v̇x = -Fx/m

v̇y = -g - Fy/m

g = 9,81 m/s², a = 343 m/s.

Mise en œuvre de l'intégration numérique

Pour résoudre le système d'équations différentielles ordinaires, la méthode de Runge-Kutta d'ordre 4 (RK4) est appliquée. Pas d'intégration dt = 0,01 s, temps final t_end = 98,4 s.

Code de simulation complet en Python :

import numpy as np
import matplotlib.pyplot as plt

# Constantes
g = 9.81  # m/s²
m = 50.0  # masse du corps, kg
S = 0.01  # section transversale, m²
rho0 = 1.225  # densité de l'air au niveau de la mer, kg/m³
H = 8000  # échelle de hauteur, m
a = 343  # vitesse du son, m/s
v0 = 600  # vitesse initiale, m/s
theta = np.radians(55)  # angle initial, rad
dt = 0.01  # pas d'intégration, s
t_end = 98.4  # temps final, s

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

# Tableaux pour stocker les résultats
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)


# Fonction pour calculer 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


# Fonction de densité de l'air
def rho(h):
    return rho0 * np.exp(-h / H)


# Système d'équations différentielles ordinaires
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éthode de Runge-Kutta d'ordre 4
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


# Intégration
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)

# Visualisation
plt.figure(figsize=(10, 6))
plt.plot(X, Y, label='Trajectoire')
plt.xlabel('X, m')
plt.ylabel('Y, m')
plt.title('Trajectoire de vol du corps')
plt.grid(True)
plt.legend()
plt.show()

# Résultats
print(f"Portée de vol : {X[-1]:.2f} m")
print(f"Altitude maximale : {np.max(Y):.2f} m")

Résultats de la simulation

Portée de vol : 29 km. Altitude maximale : 11,5 km. La trajectoire prend en compte la transition à travers les vitesses supersoniques et l'amincissement de l'air en altitude.

Hypothèses clés :

• Modèle de densité exponentielle (simplifié ; l'atmosphère réelle nécessite des profils plus précis).
• Traînée quadratique sans considérer d'autres effets (rotation, vent).

| Paramètre | Valeur |

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

| Portée | 29 000 m |

| Altitude max. | 11 500 m |

| Temps de vol | 98,4 s |

| Vitesse initiale | 600 m/s |

Points clés à retenir

• Les méthodes numériques sont essentielles pour des modèles réalistes avec Cx(M) et ρ(h) variables.
• RK4 garantit la précision avec dt=0,01 s pour les tâches balistiques.
• La loi de 1943 modélise avec précision le régime transsonique (M≈0,8–1,0).
• ρ(h) exponentielle fournit une approximation raisonnable jusqu'à 15 km.

— Editorial Team