-- Introduction au Python scientifique --

Mécanique numérique

Intégration temporelle
Application au calcul de trajectoire

simon.marie@lecnam.net
In [1]:
import numpy as np
import vtk
import scipy

#Integration scipy
from scipy.integrate import odeint
#Pour les animations
from JSAnimation.IPython_display import display_animation
In [2]:
from mayavi import mlab
mlab.init_notebook('x3d',900,500)
Notebook initialized with x3d backend.
In [3]:
import matplotlib
matplotlib.style.use('fivethirtyeight')
import matplotlib.pyplot as plt

#Option pour afficher les figures dans le notebook et eviter le plt.show():
%matplotlib inline   

On se propose dans ce TP de calculer numériquement la trajectoire d'une balle de tennis dans l'air lancée à une vitesse initiale $U_0$ dans la direction $\theta_0$ et soumis à son poid $\mathbf{P}=m\mathbf{g}$ et aux forces aérodynamiques de trainée et de portance $\mathbf{F_D}$ et $\mathbf{F_L}$.

Le principe fondamentale de la dynamique appliqué à la balle donne:

$$ m\dfrac{d\mathbf{U}}{dt}=\mathbf{P}+\mathbf{F_D}+\mathbf{F_L} $$

Avec $\mathbf{U}=\left(\begin{array}{c} u_x \\ u_y \end{array}\right)$ et $U=\lVert{\mathbf{U}}\lVert$.

$$F_D=\dfrac{1}{2}\rho U^2 S C_d$$$$F_L=\dfrac{1}{2}\rho U^2 S C_l$$$$\theta=arctan\left(\dfrac{u_y}{u_x}\right)$$

On défini les paramètres suivants qui serviront de référence pendant le TP:

In [4]:
#Atmosphere
rho=1.24 # Masse volumique standard de l'air
g=9.81 # Acceleration de la pesanteur

#Balle
balle=dict()
balle['masse']=0.58 # Masse d'une balle de tennis
balle['diam']=0.066 # Diametre de la balle
balle['u']=[0,0]
balle['pos']=[0,0]
balle['surf']=np.pi*balle['diam']**2/4. # Surface de référence
balle['Cd']=0.4
balle['Cl']=0.1

#Terrain
L=23.77

# Etat initial:
def init(balle,U0=20,th0=0.):
    balle['U0']=U0
    balle['theta0']=th0*np.pi/180.
    balle['u'][0]=balle['U0']*np.cos(balle['theta0'])
    balle['u'][1]=balle['U0']*np.sin(balle['theta0'])
    balle['pos'][0]=0.
    balle['pos'][1]=2.8

Le calcul de trajectoire (position $\mathbf{X}=(x,y)$ à chaque instant) se ramène donc à un problème du type: $$ \dfrac{d\mathbf{X}}{dt}=Kp(\mathbf{U},t) $$

et

$$ \dfrac{d\mathbf{U}}{dt}=Ku(\mathbf{U},t) $$
In [5]:
def Ku(balle,u):
    Fd=0.5*rho*balle['surf']*balle['Cd']*(u[0]**2+u[1]**2)
    Fl=0.5*rho*balle['surf']*balle['Cl']*(u[0]**2+u[1]**2)
    theta=np.arctan2(u[1],u[0])
    Kx=(1./balle['masse'])*(-Fd*np.cos(theta)-Fl*np.sin(theta))
    Ky=(1./balle['masse'])*(-Fd*np.sin(theta)+Fl*np.cos(theta)-balle['masse']*g)
    return np.array([Kx,Ky])

def Kp(u):
    Kx=u[0]
    Ky=u[1]
    return np.array([Kx,Ky])

Schéma d'Euler $\mathcal{O}(\Delta t)$

Le schéma d'euler explicite s'écrit de la façon suivante (Df cours):

$$ U^{n+1}=U^{n}+\Delta t K(U^{n},t) $$

Implémenter le schéma d'euler pour le calcul de la vitesse à chaque instant, puis en déduire le calcul de la position à chaque instant.

Comparer vos résultats avec les valeurs expérimentales fourni dans le fichier Trajectoire.dat.

In [6]:
def step_eul(balle,dt):
    balle['u']=balle['u']+dt*Ku(balle,balle['u'])
    balle['pos']= balle['pos']+dt*Kp(balle['u'])

Schéma de Runge-Kutta $\mathcal{O}(\Delta t^2)$

Le schéma de Runge_Kutta permet d'augmenter l'ordre d'intégration et consiste à calculer des étapes intermédiaires entre 2 itérations. (On parle de sous-iterations). Ainsi on peut écrire l'ordre 2 sous la forme:

$$ k_1=K(U^{n},t)\\ k_2=K(U^{n}+\dfrac{k_1\Delta t}{2},t+\dfrac{\Delta t}{2})\\ U^{n+1}=U^{n}+\Delta t k_2 $$

Implémenter le schéma de Runge-Kutta d'ordre 2 pour le calcul de la vitesse. On gardera le schéma d'Euler pour le calcul de la position. On veillera à bien réinitialiser les grandeurs initiales:

In [7]:
def step_rk(balle,dt):
    k1=Ku(balle,balle['u'])
    k2=Ku(balle,balle['u']+dt*k1/2.)
    balle['u']=balle['u']+dt*k2
    l1=Kp(balle['u'])
    l2=Kp(balle['u']+dt*l1/2.)
    balle['pos']= balle['pos']+dt*l2

Trajectoire

In [8]:
t,Nt,dt=0,300,0.01
init(balle,U0=210./3.6,th0=-7.)
Xt=np.zeros((Nt,2))
In [9]:
fig=plt.figure(figsize=(15,4),dpi=80, facecolor='w', edgecolor='k')
ax=fig.add_subplot(111)
rect = matplotlib.patches.Rectangle([11.885-0.05,0], 0.1, 0.914,facecolor="black")
ax.add_patch(rect);
ax.set_facecolor((0.88, 0.42, 0.02))
for i in np.arange(Nt):
    t+=dt
    step_rk(balle,dt)    
    if balle['pos'][1]<=0:
        balle['u'][1]=-0.95*balle['u'][1]
        balle['pos'][1]=0
    if balle['pos'][0]<=0:
        balle['u'][0]=-balle['u'][0]
        balle['pos'][0]=0
    if balle['pos'][0]>=(L+2):
        [x0,y0]=balle['pos']
        init(balle,U0=125./3.6,th0=3.)
        balle['u'][0]=-balle['u'][0]
        balle['pos'][0]=x0
        balle['pos'][1]=y0        
    Xt[i,:]= balle['pos']
    plt.plot(balle['pos'][0],balle['pos'][1],'.',color=(0.9,0.9,0.1))
    plt.ylim(0,3.);plt.xlim(0,L+2)
In [29]:
fig = plt.figure(figsize=(15,2))
ax = plt.axes(xlim=(0,L+2),ylim=(0,3.))
rect = matplotlib.patches.Rectangle([11.885-0.05,0], 0.1, 0.914,facecolor="white")
ax.add_patch(rect)
line, = ax.plot([],[],'o',color=(0.9,0.9,0.1),Markersize=12)
ax.set_facecolor((0.88, 0.42, 0.02))

def animate(data):
    line.set_data(data[0],data[1])
    return line

anim = matplotlib.animation.FuncAnimation(fig, animate, frames=Xt, interval=dt*1000)
display_animation(anim, default_mode='loop')
Out[29]:


Once Loop Reflect