Pendule en référentiel non galiléen

1.a. Système mécanique

Le pendule pesant est en liaison pivot avec un chariot, lequel glisse en translation rectiligne sur un bâti. On suppose que le bâti constitue un référentiel galiléen. Le moteur lié au bâti entraîne en rotation une roue. Le chariot est lié à cette roue par une bielle. On a ainsi un système de type manivelle-bielle qui permet de convertir le mouvement de rotation du moteur en mouvement de translation du chariot.

Figure pleine page

On note L la distance entre le centre de gravité du pendule et son axe de rotation, θ son angle avec la direction verticale. Soit un point O sur l'axe de rotation du moteur. On note a=OA la longueur de la manivelle et b=AB la longueur de la bielle. L'extrémité B de celle-ci (liaison pivot avec le chariot) a pour coordonnées (x,h).

L'angle que fait la manivelle avec l'axe horizontal Ox est α=ωt, où ω est la vitesse angulaire de rotation du moteur. L'angle que fait la bielle avec cet axe est β. En décomposant les vecteurs $\overrightarrow{OA}$ et $\overrightarrow{AB}$ sur la base orthonormée, on obtient :

$$\begin{array}{ccc}\hfill & x=acos\alpha +bcos\beta \hfill & \hfill \text{(1)}\\ \hfill & h=asin\alpha +bsin\beta \hfill & \hfill \text{(2)}\end{array}$$

Remarque : sur la figure ci-dessus, h est négatif.

En éliminant l'angle β de ces équations on obtient :

$$x=acos\alpha +\sqrt{b^2-{(h-asin\alpha )}^2}\text{(3)}$$

Voici deux exemples de courbe x(t) :

import numpy
from matplotlib.pyplot import *

def position(a,b,h,t):
    alpha = 2*numpy.pi*t
    return a*numpy.cos(alpha)+numpy.sqrt(b**2-numpy.power(h-a*numpy.sin(alpha),2))
    
figure()
N = 500
t = numpy.arange(N)*2.0/N
a=10
b=40
plot(t,position(a,b,0,t),label="h=0")
plot(t,position(a,b,-10,t),label="h=-10")
xlabel("t") 
ylabel("x")
legend(loc="lower right")
grid()
    
                    

Figure pleine page

Lorsque h=0 et $b^2\gg a^2$, la courbe est très proche d'une sinusoïde. On a en effet l'expression approchée suivante :

$$x=b+acos\alpha =b+acos\left(\omega t\right)\text{(4)}$$

On utilisera cette expression pour la modélisation du pendule.

1.b. Équation différentielle du mouvement

Le référentiel lié au chariot est non galiléen. L'accélération d'entraînement dans ce référentiel est :

$$a_e=\frac{d^{2}x}{d{t}^2}=-a{\omega }^{2}cos\left(\omega t\right)\text{(5)}$$

Les forces de pesanteur s'appliquent au centre de gravité G du pendule. Il en est de même des forces d'inertie d'entraînement (pour un référentiel en translation). Soit m la masse du pendule et J son moment d'inertie par rapport à son axe de rotation. Le théorème du moment cinétique appliqué, dans le référentiel du chariot, par rapport à l'axe de rotation, conduit à l'équation différentielle suivante :

$$J\frac{d^2\theta }{d{t}^2}+\nu \frac{d\theta }{dt}+mgLsin\theta =Lma{\omega }^{2}cos\theta cos\left(\omega t\right)\text{(6)}$$

On a introduit un frottement proportionnel à la vitesse angulaire, avec un coefficient ν.

Il faut ajouter à cette équation deux conditions initiales : $\theta \left(0\right)$ et $\dot{\theta }\left(0\right)$.

L'équation différentielle est non linéaire. Il est toutefois possible d'obtenir une équation linéaire en supposant que l'angle reste petit :

$$J\frac{d^2\theta }{d{t}^2}+\nu \frac{d\theta }{dt}+mgL\theta =Lma{\omega }^{2}cos\left(\omega t\right)\text{(7)}$$

La solution analytique exacte de cette équation est connue. En revanche, l'équation complète non linéaire ne peut être résolue que par des calculs numériques, par exemple avec la méthode d'Euler.

Pour simplifier l'étude on divise l'équation par J et on pose :

$${\omega }_0=\sqrt{\frac{mgL}{J}}\text{(8)}$$

Il s'agit de la pulsation propre du pendule linéaire (petits angles), c'est-à-dire sa pulsation d'oscillation lorsque le chariot est immobile (ω=0). Soit T la période correspondante. Il est intéressant de faire le changement de variable suivant pour le temps :

$$t\text{'}=\frac{t}{T}=\frac{{\omega }_0}{2\pi }t\text{(9)}$$

Le temps t' est sans dimensions. Pour faire le changement de variable, on utilise les identités suivantes :

$$\begin{array}{ccc}\hfill & \frac{1}{dt}=\frac{{\omega }_0}{2\pi }\frac{1}{dt\text{'}}\hfill & \hfill \text{(10)}\\ \hfill & \frac{1}{d{t}^2}=(\frac{{\omega }_0}{2\pi }{)}^2\frac{1}{dt{\text{'}}^2}\hfill & \hfill \text{(11)}\end{array}$$

On obtient ainsi l'équation différentielle

$$\frac{d^2\theta }{dt{\text{'}}^2}+\frac{2\pi \nu }{J{\omega }_0}\frac{d\theta }{dt\text{'}}+{\left(2\pi \right)}^{2}sin\theta =\frac{Lma}{J}\left(\frac{2\pi \omega }{{\omega }_0}{)}^{2}cos\theta cos\right(\frac{\omega }{{\omega }_0}2\pi t\text{'})\text{(12)}$$

Cette équation est adimensionnée, car les variables θ et t' sont sans dimensions. Lorsqu'on résout numériquement une équation du mouvement, on a toujours intérêt à la mettre sous forme adimensionnée. La période des petites oscillations libres est T'=1. On voit apparaître trois paramètres sans dimensions :

$$\begin{array}{ccc}\hfill & \gamma =\frac{2\pi \nu }{J{\omega }_0}\hfill & \hfill \text{(13)}\\ \hfill & A=\frac{Lma}{J}\hfill & \hfill \text{(14)}\\ \hfill & f=\frac{\omega }{{\omega }_0}\hfill & \hfill \text{(15)}\end{array}$$

Voici finalement l'équation différentielle à étudier :

$$\frac{d^2\theta }{dt{\text{'}}^2}+\gamma \frac{d\theta }{dt\text{'}}+{\left(2\pi \right)}^{2}sin\theta =A{\left(2\pi \right)}^2{f}^{2}cos\theta cos\left(2\pi ft\text{'}\right)\text{(16)}$$

1.c. Intégration numérique

Pour résoudre numériquement cette équation avec les deux conditions initiales, il faut la mettre sous la forme d'un système différentiel du premier ordre, en introduisant la vitesse angule du pendule, que l'on notera v_a :

$$\begin{array}{ccc}\hfill & \frac{d\theta }{dt\text{'}}=v_a\hfill & \hfill \text{(17)}\\ \hfill & \frac{d{v}_a}{dt\text{'}}=-\gamma v_a-{\left(2\pi \right)}^{2}sin\theta -A{\left(2\pi \right)}^2{f}^{2}cos\theta cos\left(2\pi ft\right)\hfill & \hfill \text{(18)}\end{array}$$

Le système est défini par une fonction (on étudie tout d'abord un système sans frottement) :

import math
import scipy.integrate

gamma = 0.0
A = 0.3
f = 1.0
w02 = (2*math.pi)**2

def systeme(Y,t):
    return [Y[1],-w02*math.sin(Y[0])-gamma*Y[1]-A*w02*f**2*math.cos(Y[0])*math.cos(2*math.pi*f*t)]
                 

On fixe une condition initiale :

Yi = [math.pi*0.5,0.0]
                 

On calcule les instants pour lesquels on veut les valeurs approchées :

T = 100.0
t = numpy.arange(start=0,stop=T,step=0.01)
                 

On fixe la tolérance absolue et la tolérance relative :

atol=1e-8
rtol=1e-8
                 

On fait l'intégration numérique :

tab_y = scipy.integrate.odeint(systeme,Yi,t,rtol=rtol,atol=atol)
yt = tab_y.transpose()
theta = yt[0]
va = yt[1]
figure()
plot(t,theta)
xlabel("t")
ylabel("theta")
axis([0,20,-2,2])
grid()
                 

figB.pdf

Pour analyser le mouvement, on peut faire une analyse spectrale :

import numpy.fft

tfd = numpy.fft.fft(theta)
N = theta.size
spectre = numpy.absolute(tfd)/N
freq = numpy.arange(N)*1.0/T
figure()
plot(freq,spectre)
xlabel("freq")
ylabel("theta")
axis([0,5,0,spectre.max()])
                 

figC.pdf

On peut aussi représenter la trajectoire dans l'espace des phases. La dimension de l'espace des phases est égale au nombre de degrés de libertés du système (au sens de système du premier ordre). Ici, il y a trois degrés de liberté (θ,v_a,α).

On peut écrire le système en faisant apparaître explicitement le degré de liberté correspondant à la rotation de la manivelle :

$$\begin{array}{ccc}\hfill & \frac{d\theta }{dt\text{'}}=v_a\hfill & \hfill \text{(19)}\\ \hfill & \frac{d\alpha }{dt\text{'}}=2\pi f\hfill & \hfill \text{(20)}\\ \hfill & \frac{d{v}_a}{dt\text{'}}=-\gamma v_a-{\left(2\pi \right)}^{2}sin\theta -A{\left(2\pi \right)}^2{f}^{2}cos\theta cos\left(\alpha \right)\hfill & \hfill \text{(21)}\end{array}$$

La représentation graphique (à deux dimensions) de la trajectoire dans l'espace des phases à 3 dimensions n'est pas faisable. Une première solution est de représenter la projection sur le plan (θ,v_a) :

figure()
plot(theta,va)
xlabel("theta")
ylabel("va")
grid()
                  

figD.pdf

Une représentation simplifiée est possible en remarquant que la variable α est de période 2π. En trace alors dans le plan de projection seulement les points pour lesquels α est multiple de 2π. Cette représentation est appelée section de Poincaré.

t = numpy.arange(start=0,stop=T,step=1.0/f)
tab_y = scipy.integrate.odeint(systeme,Yi,t,rtol=rtol,atol=atol)
yt = tab_y.transpose()
theta = yt[0]
va = yt[1]
figure()
plot(theta,va,".")
xlabel('theta')
ylabel('va')
grid()
                   

figE.pdf

Les points se répartissent sur une courbe fermée, et remplissent cette courbe lorsque t tend vers l'infini. Il s'agit d'un mouvement quasipériodique.

Voyons ce qu'il se passe si l'on augmente l'amplitude A. Pour cela, il faut augmenter la longueur de la manivelle.

A = 2.0
T = 100.0
t = numpy.arange(start=0,stop=T,step=0.01)
Yi = [math.pi*0.5,0.0]
tab_y = scipy.integrate.odeint(systeme,Yi,t,rtol=rtol,atol=atol)
yt = tab_y.transpose()
theta = yt[0]
va = yt[1]
figure()
plot(t,theta)
xlabel("t")
ylabel("theta")
grid()
                     

figF.pdf

On obtient un mouvement chaotique, au cours duquel le pendule fait des tours complets, dans un sens puis dans l'autre. Une caractéristique importante d'un mouvement chaotique est la très grande sensibilité aux conditions initiales. Refaisons le calcul avec un angle initial légèrement différent :

A = 2.0
T = 100.0
t = numpy.arange(start=0,stop=T,step=0.01)
Yi = [math.pi*0.5+1e-3,0.0]
tab_y = scipy.integrate.odeint(systeme,Yi,t,rtol=rtol,atol=atol)
yt = tab_y.transpose()
theta = yt[0]
va = yt[1]
plot(t,theta)
xlabel("t")
ylabel("theta")
grid()
                     

figG.pdf

Au début, les deux mouvements sont très proches, mais l'écart grandit très vite. On obient finalement deux trajectoires complètement différentes. En ce sens, l'évolution d'un système chaotique n'est pas prévisible bien que son équation différentielle soit déterministe. Il faudrait une précision infinie sur la condition initiale pour prévoir l'état du système à un instant t quelconque, ce qui est impossible en pratique.

Voyons l'analyse spectrale de ce mouvement chaotique :

tfd = numpy.fft.fft(theta)
N = theta.size
spectre = numpy.absolute(tfd)/N
freq = numpy.arange(N)*1.0/T
figure()
plot(freq,spectre)
xlabel("freq")
ylabel("theta")
axis([0,1,0,spectre.max()])
                 

figH.pdf

Un mouvement chaotique a un spectre continu. Voyons la section de Poincaré :

t = numpy.arange(start=0,stop=T,step=1.0/f)
tab_y = scipy.integrate.odeint(systeme,Yi,t,rtol=rtol,atol=atol)
yt = tab_y.transpose()
theta = yt[0]%(2*math.pi)-math.pi
va = yt[1]
figure()
plot(theta,va,".")
xlabel('theta')
ylabel('va')
grid()
                   

figI.pdf

Pour un mouvement chaotique, les points de la section de Poincaré se répartissent sur une surface, et non pas sur une courbe.

1.d. Animation

f=1.0

A=0.3

gamma=0.0