-- 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

Attracteur de Lorenz

L’attracteur de Lorenz est une structure fractale correspondant au comportement à long terme de l'oscillateur de Lorenz. L'attracteur montre comment les différentes variables du système dynamique évoluent dans le temps en une trajectoire non périodique.

L'oscillateur est régi par le système d'EDO suivant:

$$ \left\{ \begin{array}{rcl} \dfrac{\partial x}{\partial t}&=&\sigma(y-x)\\ \dfrac{\partial y}{\partial t}&=&x(\rho-z)-y\\ \dfrac{\partial z}{\partial t}&=&xy-\beta z\\ \end{array} \right. $$

Resolution avec scipy

La librairie scipy permet l'utilisation directe des algorithmes d'intégration d'EDO utilisés précedement. On peut alors résoudre rapidement le système précédent:

In [10]:
rho,sigma,beta = 28.0,10.,8./3.

def K(state, t):
    x, y, z = state  # unpack the state vector
    return sigma * (y - x), x * (rho - z) - y, x * y - beta * z  # derivatives

state0 = [1.0, 1.0, 1.0]
t = np.arange(0.0, 40.0, 0.01)
states = odeint(K, state0, t)

Pour visualiser la trajectoire dans l'espace 3D, on peut par exemple utiliser la librairie Mayavi permettant une visualisation interactive:

In [14]:
s = mlab.plot3d(states[:,0],states[:,1],states[:,2],t,tube_radius=0.25, colormap='hot')
s
Out[14]:

La vidéo suivante montre plusieurs oscillateurs issue de différentes conditions initiales:

In [1]:
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)
Out[1]:

Retour en haut de la page
In [ ]: