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})\).

%matplotlib inline

import numpy as N
import scipy.integrate as SI
import matplotlib.pyplot as P

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.*N.pi*rad**3 * rho            # Masse du boulet [kg]
alpha = 0.5*cx*rhoAir*N.pi*rad**2 / mass  # Coefficient de frottement par unité de masse
print("Masse du boulet: {:.2f} kg".format(mass))
print("Coefficient de frottement par unité de masse: {} S.I.".format(alpha))

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 *= N.pi / 180.  # Inclinaison [rad]
z0 = (0., 0., v0 * N.cos(alt), v0 * N.sin(alt)) # (x0, y0, vx0, vy0)

Temps caractéristique du système: \(t = \sqrt{\frac{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 = N.sqrt(mass / (g * alpha))
print("Temps caractéristique: {:.1f} s".format(tc))
t = N.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 * N.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 = SI.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("temps de coll. t(y~0) = {:.0f} s".format(t[ypos][-1])) # Dernier instant pour lequel y>0
print("portée x(y~0) = {:.0f} m".format(zs[ypos, 0][-1])) # Portée approximative du canon
#print "y(y~0) = {:.0f} m".format(zs[ypos, 1][-1]) # ~0
print("vitesse(y~0): {:.0f} m/s".format(N.hypot(zs[ypos, 2][-1], zs[ypos, 3][-1])))

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

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

This page was generated from Exercices/canon.ipynb.