In [1]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

Références bibliographiques

[1] Xavier Rogé et Paul Narchi - Le vol du boomerang - Rapport de TIPE pour le concours Inter ENS 2009, Filière MP

Repérage et systèmes de coordonnées

Dans cette étude, le boomerang est modélisé comme un solide indéformable. De par la cinématique du boomerang, nous considérons qu'il est suffisant pour repérer complètement le boomerang de le décrire par deux variables :

  • la position d'un point de référence (on prendra comme référence le centre de gravité du boomerang dans la suite)
  • le plan de rotation du boomerang

Autrement dit, nous utiliserons les variables suivantes dans notre étude :

  • le vecteur position $\vec{OG}$, soit trois coordonnées de l'espace $(x, y, z)$
  • le vecteur unitaire normal au plan de rotation du boomerang $\vec{n}$, soit trois coordonnées de l'espace $(n_x, n_y, n_z)$

Le système d'axe local au boomerang sera constitué des trois vecteurs suivants :

  • un vecteur unitaire colinéaire au vecteur vitesse : $\vec{e_x'}$
  • un vecteur unitaire qui définit le plan du boomerang et donc son axe de rotation : $\vec{e_z'}$
  • un vecteur unitaire obtenu par produit vectoriel des deux précédents (de $z'$ avec $x'$) : $\vec{e_y'}$

Équations utilisées

On utilise les équations de la mécanique du solide indéformable. Celles-ci se déclinent en deux équations :

  • l'équation des forces (semblable à la mécanique du point)
  • l'équation des moments

Nous nous inspirons de la mise en équation du TIPE.

L'équation des forces permet d'obtenir la trajectoire du boomerang en terme de position alors que l'équation des moments permet d'écrire l'évolution du plan de rotation du boomerang au cours de la propagation.

Équation des forces

$$ \frac{d \vec{p}}{dt} = \sum \vec{F_{ext}} $$

Équation du moment dynamique

$$ \frac{d \vec{\sigma_G}}{dt} = \sum \vec{ M}(\vec{F_{ext}) } $$

Liste des sous-problèmes à résoudre

On constate qu'il y a un certain nombre de sous-problèmes à résoudre quand on s'attaque à la modélisation du boomerang. En particulier, on distingue :

  • calcul de l'intertie de rotation du boomerang (dépend du profil, peut être simplifié en contributions d'objets élémentaires de type "pâle")
  • calcul exact de la portance et de la traînée avec la méthode des potentiels (Joukowski)
  • calcul du point d'application des forces aérodynamiques résultantes
  • estimation de paramètres réalistes pour la simulation (par exemple à l'aide de vidéos, comme dans le TIPE)

Algorithme de résolution du problème

A chaque pas de temps

    Résoudre l'équation des forces : en déduire l'accélération, puis la variation de vitesse, puis la variation de position (en translation)

    Résoudre l'équation du moment : en déduire la variation angulaire de la direction du boomerang

Note to self : on sent qu'on va avoir des problèmes à cause de la résolution couplée du système : il faut examiner si on peut, ou non, découpler les équations.

Forces en présence

La force de portance

En première approximation, la force de portance est définie par rapport

Résolution de l'équation du moment

On calcule tout d'abord le moment cinétique du boomerang :

$$ \vec{\sigma_G} = J \: \omega_z \: \vec{e_z'} $$

La dérivée de cette équation par rapport au temps s'écrit : $$ \frac{d}{dt} \vec{\sigma_G} = J \frac{d \omega_z}{dt} \vec{e_z'} + J \omega_z \frac{d \vec{e_z'}}{dt} $$

On suppose que la variation de vitesse de rotation est nulle. On obtient donc : $$ \frac{d}{dt} \vec{\sigma_G} = J \omega_z \frac{d \vec{e_z'}}{dt} $$

Or, le théorème du moment cinétique stipule que la quantité de l'équation précédente est égale au moment des forces extérieures. Soit :

$$ J \omega_z \frac{d \vec{e_z'}}{dt} = - L \vec{e_y'} \wedge F \vec{e_z'} = - L \: F \vec{e_x'} $$

En intégrant l'équation précédente, on peut donc facilement trouver la variation d'orientation du boomerang au cours du temps.


In [2]:
v = array([1., 0., 0.])
x = array([0., 0., 0.])
n = array([0., 1., 0.])
L = 0.1
F = 1
omega_z = 10 # tours / s
J = 0.2 * .2 ** 2 # m * L^2

dt = 0.001

In [3]:
def calc_axes(v, n):
    e_xp = v / norm(v)
    e_zp = n / norm(n)
    e_yp = cross(e_zp, e_xp)
    return (e_xp, e_yp, e_zp)

Avant de pouvoir continuer, nous devons être capable de dessiner des vecteurs en 3D. La section suivante s'occupe de résoudre ce problème.

Dessin de vecteurs

Exemple trouvé sur StackExchange

La classe Arrow3D définie dans l'exemple dérive de FancyArrowPatch. Elle est initialisée avec des coordonnées 3D. La fonction init instancie un objet FancyArrowPatch de coordonnées nulles puis remplace ses propriétés _verts3d par les coordonnées $(x, y, z)$ fournie par l'utilisateur.

Une fois que cela est fait, la méthode draw est redéfinie afin de coorectement dessiner dans l'espace le vecteur en question.


In [65]:
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d

class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        FancyArrowPatch.draw(self, renderer)

Pour utiliser cette classe, on peut se servir de l'exemple ci-dessous :


In [66]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
a = Arrow3D([0,1],[0,1],[0,1], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(a)
c = Arrow3D([0,0.5],[0,1],[0,1], mutation_scale=20, arrowstyle="-|>")
ax.add_artist(c)


Out[66]:
<__main__.Arrow3D at 0x9b46570>

Construction d'une classe pour dessiner des vecteurs centrés sur l'origine

Pour dessiner un vecteur centré sur l'origine, nous spécialisons la classe Arrow3D en construisant une fonction init différente : celle-ci est appelée avec 3 coordonnées au lieu de 3 tuples comme le constructeur standard.


In [75]:
class Arrow3D_centered(Arrow3D):
    def __init__(self, x, y, z, *args, **kwargs):
        """"3D arrow class for plotting vectors starting at (0, 0, 0).
        Expected input for the constructor is are three coordinates: x, y and z."""
        super(Arrow3D_centered, self).__init__((0, x), (0, y), (0, z), mutation_scale=20, arrowstyle="-|>", *args, **kwargs)

On peut maintenant écrire un exemple pour l'utilisation de cette classe :


In [90]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

a = Arrow3D_centered(1, 0.5, .5)
ax.add_artist(a)

b = Arrow3D_centered(0.5, 0.25, 0.75)
ax.add_artist(b)

c = Arrow3D_centered(2, 4, -10)
ax.add_artist(c)

ax.set_zlim3d(auto=True)
ax.set_xlim3d(auto=True)
ax.set_ylim3d(auto=True)


Out[90]:
(0.0, 1.0)

Il semble que nous ayons un problème d'axes automatiques...

Simulation


In [77]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

for i in range(1000):
    (e_xp, e_yp, e_zp) = calc_axes(v, n)
    dn = - L * F / J / omega_z * dt * e_xp
    n += dn
    if i % 100:
        ax.add_artist(Arrow3D_centered(n[0], n[1], n[2]))


Dessin de vecteurs


In [7]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

def plot_vect(ax, x, y, z):
    xs = array((0, x))
    ys = array((0, y))
    zs = array((0, z))
    ax.plot(xs, ys, zs)
    
plot_vect(ax, 1, 1, 1)



In [8]:
%pylab qt


Populating the interactive namespace from numpy and matplotlib

In [9]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

for i in range(1000):
    (e_xp, e_yp, e_zp) = calc_axes(v, n)
    dn = - L * F / J / omega_z * dt * e_xp
    n += dn
    if i % 30:
        plot_vect(ax, n[0], n[1], n[2])

In [ ]: