Table des matières Mathematica

Équations du mouvement d'un solide

1. Introduction

L'objectif de cette page est d'établir les équations différentielles du mouvement d'un solide indéformable, dans le but de faire leur intégration numérique. On se limite au cas où les forces agissant sur le solide sont connues a priori. Ces forces peuvent être non conservatives.

La méthode exposée ici consiste à utilise la matrice rotation qui définit l'orientation du solide dans l'espace. Elle peut s'appliquer aux problèmes de gyroscopes, éventuellement avec frottements, et aux mouvements dans un fluide. La méthode peut être utilisée en dynamique moléculaire, pour la simulation des molécules rigides.

2. Mouvement du centre de masse

Dans le cas d'un mouvement comportant des degrés de liberté de translation, le théorème de la résultante cinétique permet d'écrire l'équation du mouvement du centre de masse :

mdVGdt=Fext(1)

Lorsque la résultante des forces extérieures est connue, l'intégration de ces 3 équations différentielles du premier ordre se fait en ajoutant l'équation suivante :

dXGdt=VG(2)

En projection sur une base orthonormée, nous avons donc un système différentiel de la forme y'=f(y,t) comportant 6 équations.

3. Mouvement de rotation

3.a. Théorème du moment cinétique

Le mouvement d'un solide ayant un point fixe ou le mouvement dans le référentiel du centre de masse se fait avec 3 degrés de liberté de rotation. Dans les deux cas, l'équation du mouvement de la rotation est obtenue avec le théorème du moment cinétique (en un point fixe ou au centre de masse) :

dσdt=Next(3)

3.b. Repère d'espace et repère lié au solide

Soit O le point fixe du mouvement. On note Oxyz un repère cartésien lié au référentiel dans lequel le mouvement du solide est étudié; on l'appellera repère d'espace. On définit également un repère lié au solide Oxsyszs.

Les composantes d'un vecteur forment une matrice colonne (3x1). On désignera par une lettre minuscule ce type de matrice. Par exemple, σ désigne la matrice colonne des composantes du moment cinétique dans le repère d'espace et σs désigne les composantes du même vecteur dans le repère du solide.

Les composantes d'un tenseur de rang 2 forment une matrice (3x3) qui sera désignée par une lettre majuscule. Par exemple T désigne les composantes du tenseur d'inertie dans le repère d'espace et Ts les composantes du même tenseur dans le repère du solide.

3.c. Matrice rotation

Soit xi la matrice (3x1) des coordonnées d'un point particulier du solide (indice i) dans le repère d'espace et xis ses coordonnées dans le repère du solide. La rotation du solide est donnée par une matrice (3x3) Q orthogonale [1] :

xi=Qxis(4)

La matrice Q est orthogonale; elle vérifie :

QQT=1(5)

Comme il s'agit d'une rotation, le déterminant de Q est +1.

La matrice Q sera choisie pour représenter la position à tout instant du solide. Cela correspond à 9 paramètres alors que 3 sont en principe suffisants pour décrire l'orientation du solide dans l'espace. Par exemple, les 3 angles d'Euler définis plus loins définissent cette orientation. Cependant, les angles d'Euler conduisent à des équations différentielles difficiles à intégrer numériquement. Par ailleurs, la matrice Q permet d'obtenir rapidement les coordonnées d'un point quelconque du solide. Elle peut être fournie directement à un programme de visualisation pour obtenir une représentation graphique de l'objet.

3.d. Vitesse angulaire

En dérivant l'équation (4) par rapport au temps :

vi=dQdtxis=dQdt QTxi=Ωxi(6)

Ω est la matrice vitesse angulaire. Pour une rotation infinitésimale du solide pendant l'intervalle de temps dt, le déplacement d'un point est :

xi(t+dt)=(1+Ωdt)xi(t)(7)

La rotation infinitésimale est donc donnée par la matrice orthogonale 1+Ωdt . La condition d'orthogonalité s'écrit :

(1+Ωdt)(1+ΩTdt)=1(8)

En inversant la rotation infinitésimale, on obtient :

(1+Ωdt)(1-Ωdt)=1(9)

La comparaison des deux équations précédentes montre que la matrice Ω est antisymétrique :

ΩT=-Ω(10)

On posera donc :

Ω=(0-ω3ω2ω30-ω1-ω2ω10)(11)

Le vecteur vitesse angulaire est le vecteur dont les composantes sur le repère d'espace sont :

ω=(ω1ω2ω3)(12)

Les composantes de la vitesse s'expriment alors comme un produit vectoriel de deux matrices colonnes :

vi=Ωxi=ωxi(13)

La matrice vitesse angulaire est en fait la matrice des composantes d'un tenseur antisymétrique de rang 2. Ses composantes dans le repère du solide sont données par la matrice :

Ωs=QTΩQ=QTdQdt(14)

Avec l'équation (6), on en déduit :

dQdt=QΩs(15)

3.e. Équations d'Euler

Considérons les composantes du moment cinétique (au centre de masse ou en un point fixe) sur le repère d'espace et sur le repère du solide :

σ=Qσs(16)

En dérivant par rapport au temps :

dσdt=dQdtσs+Qdσsdt(17)

Le théorème du moment cinétique s'écrit en projection sur le repère d'espace :

dσdt=n(18)

n désigne la matrice colonne (3x1) des composantes du moment des forces extérieures dans le repère d'espace.

Les composantes du moment des forces dans le repère du solide sont données par :

ns=QTn(19)

Le théorème du moment cinétique s'écrit donc :

QTdQdtσs+dσsdt=ns(20)

D'autre part :

QTdQdt=QTQΩs=Ωs(21)

On obtient finalement l'équation suivante :

dσsdt+Ωsσs=ns(22)

qui s'écrit aussi avec les composantes du vecteur vitesse angulaire dans le repère du solide :

dσsdt+ωsσs=ns(23)

Le moment cinétique s'exprime en fonction du vecteur vitesse angulaire par le tenseur d'inertie. Dans le repère lié au solide, on a :

σs=Tsωs(24)

Ts est la matrice des composantes du tenseur d'inertie au point O considéré pour le moment cinétique (centre de masse ou point fixe). On choisit à présent le repère (Oxsyszs) constitué par les axes propres du tenseur d'inertie. On désigne ces trois axes par les indices 1,2,3. Les moments d'inertie par rapport aux axes propres sont notés Ik. L'équation (23) conduit alors aux équations d'Euler :

I1dωs1dt=ωs2ωs3(I2-I3)+ns1(25)I2dωs2dt=ωs3ωs1(I3-I1)+ns2(26)I3dωs3dt=ωs1ωs2(I1-I2)+ns3(27)

3.f. Équations différentielles du mouvement de rotation

Les équations d'Euler associées à l'équation (15) permettent de définir un système différentiel du premier ordre. Notons qij les éléments de la matrice Q. L'équation (15) conduit aux 9 équations différentielles suivantes :

dq11dt=ωs3q12-ωs2q13(28)dq12dt=-ωs3q11+ωs1q13(29)dq13dt=ωs2q11-ωs1q12(30)dq21dt=ωs3q22-ωs2q23(31)dq22dt=-ωs3q21+ωs1q23(32)dq23dt=ωs2q21-ωs1q22(33)dq31dt=ωs3q32-ωs2q33(34)dq32dt=-ωs3q31+ωs1q33(35)dq33dt=ωs2q31-ωs1q32(36)

Les 12 équations (25) à (36) constituent un système différentiel du premier ordre pour les 12 variables qij et ωk. Les composantes du moment des forces extérieures sont soit évaluées directement sur le repère du solide, soit d'abord calculées sur le repère d'espace puis converties avec la relation

ns=QTn(37)

Remarquons que les 3 degrés de libertés nécessitent en principe 6 équations du premier ordre. Nous avons donc ici 6 équations redondantes. Cependant, la méthode offre l'avantage de fournir directement la matrice de la transformation orthogonale. De plus, la condition d'orthogonalité de Q pourra être vérifiée pendant l'intégration numérique, ce qui donne un moyen de contrôler la précision du calcul.

Voyons à présent comment définir les conditions initiales avec les angles d'Euler.

3.g. Angles d'Euler

Les angles d'Euler permettent de décomposer la transformation orthogonale Q en trois rotations. Considérons les repères R=(Oxyz) et Rs=(Oxsyszs) initialement confondus. On fait tout d'abord tourner le repère Rs d'un angle φ autour de l'axe Oz. La matrice de rotation correspondante est :

Qφ=(cosφ-sinφ0sinφcosφ0001)(38)

La seconde transformation est une rotation d'un angle θ autour de l'axe Oxs :

Qθ=(1000cosθ-sinθ0sinθcosθ)(39)

Pour finir, on fait tourner Rs d'un angle ψ autour de l'axe Ozs :

Qψ=(cosψ-sinψ0sinψcosψ0001)(40)

La matrice Q, qui par définition permet de calculer les coordonnées d'un point dans le repère R en fonction de ses coordonnées dans le repère du solide Rs est :

Q=QφQθQψ(41)

Calculons cette matrice avec Mathematica :

Q=RotationMatrix[phi,{0,0,1}].RotationMatrix[theta,{1,0,0}].RotationMatrix[psi,{0,0,1}]
				
( cos(φ) cos(ψ)-cos(θ) sin(φ) sin(ψ) -cos(ψ) cos(θ) sin(φ)-cos(φ) sin(ψ) sin(φ) sin(θ) cos(ψ) sin(φ)+cos(φ) cos(θ) sin(ψ) cos(φ) cos(ψ) cos(θ)-sin(φ) sin(ψ) -cos(φ) sin(θ) sin(ψ) sin(θ) cos(ψ) sin(θ) cos(θ))

Les angles d'Euler peuvent être utilisés pour définir les valeurs initiales des éléments qij de la matrice Q. D'autre part, ces angles peuvent être calculés au cours de l'intégration en fonction des qij, ce qui peut être utile pour calculer certains moments de force ou pour suivre le mouvement du solide.

Les valeurs initiales du vecteur vitesse angulaire dans le repère du solide sont obtenues avec la relation :

Ωs=QTdQdt(42)
Q = Q/.{phi->phi[t],theta->theta[t],psi->psi[t]};
OmegaS = Simplify[Transpose[Q].D[Q,t]]
				
( 0 -cos(θ(t)) φ'(t)-ψ'(t) cos(ψ(t)) sin(θ(t)) φ'(t)-sin(ψ(t)) θ'(t) cos(θ(t)) φ'(t)+ψ'(t) 0 -sin(ψ(t)) sin(θ(t)) φ'(t)-cos(ψ(t)) θ'(t) sin(ψ(t)) θ'(t)-cos(ψ(t)) sin(θ(t)) φ'(t) sin(ψ(t)) sin(θ(t)) φ'(t)+cos(ψ(t)) θ'(t) 0)

On constate que les composantes du vecteur vitesse angulaire dans le repère du solide s'expriment en fonction des dérivées des angles d'Euler par la relation :

ωs=A(φ'(t)θ'(t)ψ'(t))(43)

où la matrice A est :

A={{Sin[theta[t]]*Sin[psi[t]],Cos[psi[t]],0},{Cos[psi[t]]*Sin[theta[t]],-Sin[psi[t]],0},{Cos[theta[t]],0,1}}
( sin(ψ(t)) sin(θ(t)) cos(ψ(t)) 0 cos(ψ(t)) sin(θ(t)) -sin(ψ(t)) 0 cos(θ(t)) 0 1)

La donnée des angles d'Euler initiaux et de leur vitesse angulaire initiale permet donc de calculer les valeurs initiales des variables ωsk.

Inversement, on peut exprimer les dérivées des angles d'Euler en fonction des composantes de la vitesse angulaire dans le repère du solide :

(φ'(t)θ'(t)ψ'(t))=A-1ωs(44)

La matrice inverse de A est :

invA=Simplify[Inverse[A]]
( csc(θ(t)) sin(ψ(t)) cos(ψ(t)) csc(θ(t)) 0 cos(ψ(t)) -sin(ψ(t)) 0 -cot(θ(t)) sin(ψ(t)) -cos(ψ(t)) cot(θ(t)) 1)

Les 3 équations (44) et les équations d'Euler (25,26,27) forment un système différentiel de 6 équations du premier ordre pour les angles d'Euler et les composantes de la vitesse angulaire dans le repère du solide. La matrice ci-dessus a toutefois l'inconvénient de ne pas être définie pour θ=0 et θ=π. L'intégration numérique de ces équations différentielles pose des problèmes au voisinage de ces valeurs particulières de θ. C'est pourquoi il est préférable d'utiliser la matrice rotation comme expliquée précédemment. Les angles d'Euler sont utilisés seulement pour obtenir la condition initiale.

Références
[1]  H. Goldstein,  Classical mechanics,  (Addison-Wesley, 1980)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.