-- Mécanique du Vol --

TP n°1

Introduction au calcul de trajectoire

simon.marie@lecnam.net
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import scipy
#Option pour afficher les figures dans le notebook et eviter le plt.show():
%matplotlib inline
plt.rc('xtick', labelsize=20)
plt.rc('ytick', labelsize=20)

On se propose dans ce TP de se familiariser avec l'environement de Jupyter et de Python en calculant numériquement la trajectoire d'un projectile. Ce TP est une introduction au calcul de trajectoire et servira de base pour l'étude des phases importantes de la mécanique du vol.

1 - Position du problème

Considérons un projectile légèrement planant catapulté avec une force $\mathbf{C_0}$ initiale orientée dans la direction $\theta_0$. Le projectile est soumis à son poid $\mathbf{P}=m\mathbf{g}$ et aux forces aérodynamique de trainée et de portance $\mathbf{F_D}$ et $\mathbf{F_L}$ dont on donne les coefficients. On considérera que la force exercée par la catapulte agit sur le projectile pendant une durée $t_0=1$s.

Le principe fondamentale de la dynamique appliqué à ce projectile donne:

$$ m\dfrac{d\mathbf{U}}{dt}=\mathbf{C_0}+\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$.

En projetant cette équation sur les deux axes du repère $(Oxy)$ on obtient:

$$ m\dfrac{du_x}{dt}=C_0\cos(\theta_0)-F_D\cos(\theta)-F_L\sin(\theta) $$$$ m\dfrac{du_y}{dt}=C_0\sin(\theta_0)-F_D\sin(\theta)+F_L\cos(\theta)-mg $$

Avec: $$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)$$

Ainsi, pour connaitre à chaque instant la position du projectile, il est nécessaire de connaitre avec précision sa vitesse. Pour connaitre sa vitesse il est nécessaire de connaitre son accélération, et donc les forces qui s'exercent sur le projectile. Or ces forces dépendent de la vitesse, il est donc impossible de résoudre analytiquement ce problème.

On utilise donc des approximations successives en décomposant le temps en itérations. On se propose d'étudier dans ce TP les méthodes pour calculer une trajectoire.

2 - Calcul de trajectoire

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

In [ ]:
m=1.     # Masse de 1kg
D=0.1    # longueur charactéristique du projectile
g=9.81   # Acceleration de la pesanteur
S=D**2   # Surface de référence
rho=1.22 # Masse volumique standard de l'air

# Aérodynamique
Cd=0.05  # Coeff de trainée
Cl=0.05   # Coeff de portance

# Propulsion:
C0=80                # Force du catapultage
t0=2                 # Temps du catapultage
theta0=np.pi/6       # Angle du catapultage

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

avec $$ K(\mathbf{U},t)=\dfrac{\mathbf{C_0}(t)+\mathbf{P}+\mathbf{F_D}(\mathbf{U})+\mathbf{F_L}(\mathbf{U})}{m} $$

Le principe de l'intégration temporelle par itération est basé sur la procédure suivante:

  1. Je connais la vitesse du projectile à l'itération $n$
  2. J'en déduis le bilan des forces qui s'exercent sur celui-ci: Calcul de K(U,t) => Je connais donc l'accélération.
  3. Je calcul donc la vitesse de l'itération $n+1$
  4. Je calcul la nouvelle position
  5. Je recommence à l'étape 1

Pour réaliser l'étape 2. il faut donc créer une fonction qui calcul le bilan des forces et renvoi l'accélération sur les deux axes $x$ et $y$:

In [ ]:
#### Fonction qui calcul le bilan des forces: ####
def K(U,t):
    # !!! Ici U est un vecteur à deux composantes U=[ux,uy]
    Fd=...   # Force de trainee
    Fl=...   # Force de portance
    theta=np.arctan2(U[1],U[0])
    Kx=(1./m)*(...)
    Ky=(1./m)*(...)
    return np.array([Kx,Ky])

Schéma d'Euler

La façon la plus simple de mettre à jour la vitesse est d'utiliser une approximation d'ordre 1 de la dérivé temporelle. C'est le schéma d'Euler, il s'écrit de la façon suivante (cours AER211):

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

Ou de la même façon pour la position:

$$ X^{n+1}=X^{n}+\Delta t 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. Faites varier le pas de temps $dt$. Que remarquez vous ?

En choisissant les paramètres proposés, votre projectile doit retomber à 1870m de son point de départ.

In [ ]:
# Vitesse initiale:
U0=0 
ux=U0*np.cos(theta0)
uy=U0*np.sin(theta0)

# Position Initiale:
x=0
y=1.

# Integration
dt=0.1      # Pas de temps
t=0         # On initialise le temps à 0
C=C0        # La force de la catapulte 
fin=60      # On se fixe un temps maximum d'intégration
col='r'     
fig=plt.figure(figsize=(15,12))
while y>0 and t<fin:     # On répète les étapes tant que le projectile est en l'air et que le temps est plus petit que la limite
    t+=dt   # On incrémente le compteur du temps
    if t>t0: # si le temps est plus grand que t0 la catapulte n'a plus d'effet
        C=0
        col='k'
    [ux,uy]=[ux,uy]+...               ##### CALCUL de la NOUVELLE VITESSE
    [x,y]=[x,y]+dt*np.array([ux,uy])  ##### CALCUL de la NOUVELLE POSITION
    # Représentation
    plt.subplot(211)
    plt.plot(x,y,'.'+col)
    plt.subplot(212)
    plt.plot(ux,uy,'.'+col)

# On affine le graphique
plt.subplot(211)
plt.grid()
plt.title('t= '+str(round(t,3))+' s, Xmax= '+str(round(x,3))+' m',fontsize=20.)
plt.ylabel('Altitude',fontsize=20.)
plt.xlabel('Position',fontsize=20.)
plt.subplot(212)
plt.grid()
plt.ylabel('Uy',fontsize=20.)
plt.xlabel('Ux',fontsize=20.)

Discussion:

Schéma de Runge-Kutta

Si les accélération sont très importantes, le schéma d'Euler peut devenir très peu précis. Ainsi, 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 [ ]:
# Vitesse initiale:
U0=0 
ux=U0*np.cos(theta0)
uy=U0*np.sin(theta0)

# Position Initiale:
x=0
y=1.

# Integration
dt=0.001
t=0
C=C0
col='r'
fin=60
fig=plt.figure(figsize=(15,12))
while y>0 and t<fin:
    t+=dt
    if t>t0:
        C=0
        col='k'
    k1=...
    k2=...
    [ux,uy]=[ux,uy]+...
    [x,y]=[x,y]+dt*np.array([ux,uy])
    plt.subplot(211)
    plt.plot(x,y,'.'+col)
    plt.subplot(212)
    plt.plot(ux,uy,'.'+col)

plt.subplot(211)
plt.grid()
plt.title('t= '+str(round(t,3))+' s, Xmax= '+str(round(x,3))+' m',fontsize=20.)
plt.ylabel('Altitude',fontsize=20.)
plt.xlabel('Position',fontsize=20.)
plt.subplot(212)
plt.grid()
plt.ylabel('Uy',fontsize=20.)
plt.xlabel('Ux',fontsize=20.)

Discussions:

3 - Influence des paramètres

Lorsque vous avez bien compris le fonctionement des éléments précédent:

  • Essayez de faire varier les paramètres (Propulsion, portance, trainée) en essayant de prédire le résultat.
  • Augmentez la valeur de $C_0$ et observez les écarts entre les deux méthodes.

Dans les TP suivants, l'intégration temporelle sera intégrée dans un "simulateur" qui servira à simuler et à étudier les différentes phase de vol d'un avion.

Conclusion

Présentez une synthèse de vos résultats et des étapes de calcul de votre TP.

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