Notebook originel: canon.ipynb

Trajectoire d’un boulet de canon

Nous allons intégrer les équations du mouvement pour un boulet de canon soumis à des forces de frottement « turbulentes » (non-linéaires):

\[\ddot{\mathbf{r}} = \mathbf{g} - \frac{\alpha}{m}v\times\mathbf{v}.\]

Cette équation différentielle non linéaire du 2d ordre doit être réécrite sous la forme de deux équations différentielles couplées du 1er ordre:

\[\begin{split}\begin{cases} \dot{\mathbf{r}} &= \mathbf{v} \\ \dot{\mathbf{v}} &= \mathbf{g} - \frac{\alpha}{m}v\times\mathbf{v}. \end{cases}\end{split}\]

Il s’agit donc de résoudre une seule équation différentielle du 1er ordre en \(\mathbf{z} = (\mathbf{r},\mathbf{v})\).

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

%matplotlib inline

Valeurs numériques pour un boulet de canon de 36 livres:

g = 9.81        # Pesanteur [m/s2]
cx = 0.45       # Coefficient de frottement d'une sphère
rhoAir = 1.2    # Masse volumique de l'air [kg/m3] au niveau de la mer, T=20°C
rad = 0.1748/2  # Rayon du boulet [m]
rho = 6.23e3    # Masse volumique du boulet [kg/m3]
mass = 4./3.*np.pi*rad**3 * rho            # Masse du boulet [kg]
alpha = 0.5*cx*rhoAir*np.pi*rad**2 / mass  # Coefficient de frottement par unité de masse
print(f"Masse du boulet: {mass:.2f} kg")
print(f"Coefficient de frottement par unité de masse: {alpha} S.I.")

Masse du boulet: 17.42 kg
Coefficient de frottement par unité de masse: 0.0003718994604243878 S.I.

Conditions initiales:

v0 = 450.            # Vitesse initiale [m/s]
alt = 45.            # Inclinaison du canon [deg]
alt *= np.pi / 180.  # Inclinaison [rad]
z0 = (0., 0., v0 * np.cos(alt), v0 * np.sin(alt))  # (x0, y0, vx0, vy0)

Temps caractéristique du système: \(t = \sqrt{m/(g\alpha)}\) (durée du régime transitoire). L’intégration des équations se fera sur un temps caractéristique, avec des pas de temps significativement plus petits.

tc = np.sqrt(mass / (g * alpha))
print(f"Temps caractéristique: {tc:.1f} s")
t = np.linspace(0, tc, 100)

Temps caractéristique: 69.1 s

Définition de la fonction \(\dot{\mathbf{z}}\), avec \(\mathbf{z} = (\mathbf{r},\mathbf{v})\)

def zdot(z, t):
    """Calcul de la dérivée de z=(x, y, vx, vy) à l'instant t."""

    x, y, vx, vy = z
    alphav = alpha * np.hypot(vx, vy)

    return (vx, vy, -alphav * vx, -g - alphav * vy) # dz/dt = (vx,vy,x..,y..)

Intégration numérique des équations du mouvement à l’aide de la fonction scipy.integrate.odeint:

zs = odeint(zdot, z0, t)

Le tableau zs contient les valeurs de \(z\) à chaque instant \(t\): il est donc de taille (len(t),4).

ypos = zs[:,1]>=0 # y>0?
print(f"temps de coll. t(y~0) = {t[ypos][-1]:.0f} s") # Dernier instant pour lequel y>0
print(f"portée x(y~0) = {zs[ypos, 0][-1]:.0f} m") # Portée approximative du canon
#print(f"y(y~0) = {zs[ypos, 1][-1]:.0f} m") # ~0
print(f"vitesse(y~0): {np.hypot(zs[ypos, 2][-1], zs[ypos, 3][-1]):.0f} m/s")

temps de coll. t(y~0) = 36 s
portée x(y~0) = 3966 m
vitesse(y~0): 140 m/s

fig, ax = plt.subplots()
ax.plot(zs[ypos, 0], zs[ypos, 1])
ax.set(xlabel="x [m]", ylabel="y [m]", title="Trajectoire du boulet");

Cette page a été générée à partir de canon.ipynb.