Toupie pesante à point fixe

1. Équations du mouvement

On suppose que la pointe de la toupie (point O) reste fixe par rapport au référentiel inertiel. En pratique, l'expérience est réalisée avec une liaison de type Cardan, c'est-à-dire un gyroscope dont le centre de masse est non confondu avec le point fixe du mouvement.

La mise en équation en vue d'une intégration numérique est faite selon la méthode exposée dans équations du mouvement d'un solide.

Soit R_s=(0x_sy_sz_s) le repère lié au solide coïncidant avec les axes propres du tenseur d'inertie au point O. L'axe Oz_s est l'axe de symétrie de la toupie, de moment d'inertie I₃. Les moments d'inertie I₁ et I₂ sont égaux.

Considérons tout d'abord le moment des forces de pesanteur au point O. Le centre de masse étant sur l'axe Oz_s, on pose L=OG. Soit le repère d'espace R=(0xyz) (lié au référentiel) où Oz est la direction verticale. La matrice Q orthogonale définie l'orientation de la toupie. On a tout d'abord :

$$\overrightarrow{O}G=L(q_{13}{\overrightarrow{e}}_x+q_{23}{\overrightarrow{e}}_y+q_{33}{\overrightarrow{e}}_z)\text{(1)}$$

Le moment du poids s'écrit :

$$\overrightarrow{N}=mgL(-q_{23}{\overrightarrow{e}}_x+q_{13}{\overrightarrow{e}}_y)\text{(2)}$$

La matrice colonne correspondante est donc :

$$n=mgL\left(\begin{array}{c}\hfill -q_{23}\hfill \\ \hfill q_{13}\hfill \\ \hfill 0\hfill \end{array}\right)\text{(3)}$$

Pour écrire les équations d'Euler, nous avons besoin des composantes du moment dans le repère R_s du solide. Ceci est obtenu par la transformation suivante :

$$n_s=Q^{T}n\text{(4)}$$

Pour modéliser les frottements dus à la rotation de la toupie, on introduit un moment dans la direction de l'axe de symétrie de la toupie :

$$f_s=-k\left(\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill {\omega }_{s3}\hfill \end{array}\right)\text{(5)}$$

3 paramètres constants interviennent dans les équations d'Euler :

$$\begin{array}{ccc}\hfill & r=\frac{I_3}{I_1}=\frac{I_3}{I_2}\hfill & \hfill \text{(6)}\\ \hfill & a={\omega }_A^2=\frac{mgL}{I_3}\hfill & \hfill \text{(7)}\\ \hfill & b=\frac{k}{I_3}\hfill & \hfill \text{(8)}\end{array}$$

Soit ω₀ la vitesse angulaire initiale de la toupie, c'est-à-dire la dérivée par rapport au temps de ψ à t=0 (à ne pas confondre avec la valeur initiale de ω_s3). Cette vitesse angulaire servira d'unité pour les variables ω_k. Pour une toupie rapide, on a :

$${\omega }_A\ll {\omega }_0\text{(9)}$$

Écrivons les équations d'Euler avec Mathematica :

Q={{q11[t],q12[t],q13[t]},{q21[t],q22[t],q23[t]},{q31[t],q32[t],q33[t]}};
n={-q23[t],q13[t],0};
ns=Transpose[Q].n;
fs={0,0,omegaS3[t]};
euler={omegaS1'[t]==(1-r)*omegaS2[t]*omegaS3[t]+a*r*ns[[1]]-b*r*fs[[1]],
       omegaS2'[t]==(r-1)*omegaS3[t]*omegaS1[t]+a*r*ns[[2]]-b*r*fs[[2]],
	   omegaS3'[t]==a*ns[[3]]-b*fs[[3]]}
			 

$$\left(\begin{array}{c}\hfill \text{omegaS1}\text{'}\left(t\right)=ar(\text{q13}\left(t\right)\text{q21}\left(t\right)-\text{q11}\left(t\right)\text{q23}\left(t\right))+(1-r)\text{omegaS2}\left(t\right)\text{omegaS3}\left(t\right)\hfill \\ \hfill \text{omegaS2}\text{'}\left(t\right)=ar(\text{q13}\left(t\right)\text{q22}\left(t\right)-\text{q12}\left(t\right)\text{q23}\left(t\right))+(r-1)\text{omegaS1}\left(t\right)\text{omegaS3}\left(t\right)\hfill \\ \hfill \text{omegaS3}\text{'}\left(t\right)=-b\text{omegaS3}\left(t\right)\hfill \end{array}\right)$$

Bien que la dernière équation soit directement intégrable, on la laissera dans le système pour plus de généralité.

Les équations différentielles vérifiées par les composantes de la transformation orthogonale :

OmegaS={{0,-omegaS3[t],omegaS2[t]},{omegaS3[t],0,-omegaS1[t]},{-omegaS2[t],omegaS1[t],0}};
d = Q.OmegaS;
dQ = {q11'[t]==d[[1,1]],q12'[t]==d[[1,2]],q13'[t]==d[[1,3]],
      q21'[t]==d[[2,1]],q22'[t]==d[[2,2]],q23'[t]==d[[2,3]],
	  q31'[t]==d[[3,1]],q32'[t]==d[[3,2]],q33'[t]==d[[3,3]]}
			 

$$\left(\begin{array}{c}\hfill \text{q11}\text{'}\left(t\right)=\text{omegaS3}\left(t\right)\text{q12}\left(t\right)-\text{omegaS2}\left(t\right)\text{q13}\left(t\right)\hfill \\ \hfill \text{q12}\text{'}\left(t\right)=\text{omegaS1}\left(t\right)\text{q13}\left(t\right)-\text{omegaS3}\left(t\right)\text{q11}\left(t\right)\hfill \\ \hfill \text{q13}\text{'}\left(t\right)=\text{omegaS2}\left(t\right)\text{q11}\left(t\right)-\text{omegaS1}\left(t\right)\text{q12}\left(t\right)\hfill \\ \hfill \text{q21}\text{'}\left(t\right)=\text{omegaS3}\left(t\right)\text{q22}\left(t\right)-\text{omegaS2}\left(t\right)\text{q23}\left(t\right)\hfill \\ \hfill \text{q22}\text{'}\left(t\right)=\text{omegaS1}\left(t\right)\text{q23}\left(t\right)-\text{omegaS3}\left(t\right)\text{q21}\left(t\right)\hfill \\ \hfill \text{q23}\text{'}\left(t\right)=\text{omegaS2}\left(t\right)\text{q21}\left(t\right)-\text{omegaS1}\left(t\right)\text{q22}\left(t\right)\hfill \\ \hfill \text{q31}\text{'}\left(t\right)=\text{omegaS3}\left(t\right)\text{q32}\left(t\right)-\text{omegaS2}\left(t\right)\text{q33}\left(t\right)\hfill \\ \hfill \text{q32}\text{'}\left(t\right)=\text{omegaS1}\left(t\right)\text{q33}\left(t\right)-\text{omegaS3}\left(t\right)\text{q31}\left(t\right)\hfill \\ \hfill \text{q33}\text{'}\left(t\right)=\text{omegaS2}\left(t\right)\text{q31}\left(t\right)-\text{omegaS1}\left(t\right)\text{q32}\left(t\right)\hfill \end{array}\right)$$

Les conditions initiales sont données à l'aide des angles d'Euler :

$$\begin{array}{ccc}\hfill & \psi \left(0\right)=0,\psi \text{'}\left(0\right)={\omega }_0\hfill & \hfill \text{(10)}\\ \hfill & \theta \left(0\right)={\theta }_0,\theta \text{'}\left(0\right)=0\hfill & \hfill \text{(11)}\\ \hfill & \phi \left(0\right)=0,\phi \text{'}\left(0\right)={\omega }_p\hfill & \hfill \text{(12)}\end{array}$$

ω_p est la vitesse de précession initiale. Les conditions initiales doivent être exprimées pour les composantes du vecteur vitesse angulaire sur les axes de R_s et pour les éléments de la matrice Q :

Q0=RotationMatrix[phi,{0,0,1}].RotationMatrix[theta,{1,0,0}].RotationMatrix[psi,{0,0,1}];
dQ0 = D[Q0/.{phi->phi[t],theta->theta[t],psi->psi[t]},t];
Qinit = Q0/.{phi->0,theta->theta0,psi->0};
dQinit = dQ0/.{phi'[t]->omegaP,psi'[t]->omega0,theta'[t]->0,phi[t]->0,
               theta[t]->theta0,psi[t]->0};
OmegaS = Transpose[Qinit].dQinit;
init = {omegaS1[0]==OmegaS[[3,2]],omegaS2[0]==OmegaS[[1,3]],omegaS3[0]==OmegaS[[2,1]],
        q11[0]==Qinit[[1,1]],q12[0]==Qinit[[1,2]],q13[0]==Qinit[[1,3]],
		q21[0]==Qinit[[2,1]],q22[0]==Qinit[[2,2]],q23[0]==Qinit[[2,3]],
		q31[0]==Qinit[[3,1]],q32[0]==Qinit[[3,2]],q33[0]==Qinit[[3,3]]}
			 

$$\left(\begin{array}{c}\hfill \text{omegaS1}\left(0\right)=0\hfill \\ \hfill \text{omegaS2}\left(0\right)=\text{omegaP}sin\left(\text{theta0}\right)\hfill \\ \hfill \text{omegaS3}\left(0\right)=cos\left(\text{theta0}\right)(\text{omega0}cos\left(\text{theta0}\right)+\text{omegaP})+\text{omega0}{sin}^2\left(\text{theta0}\right)\hfill \\ \hfill \text{q11}\left(0\right)=1\hfill \\ \hfill \text{q12}\left(0\right)=0\hfill \\ \hfill \text{q13}\left(0\right)=0\hfill \\ \hfill \text{q21}\left(0\right)=0\hfill \\ \hfill \text{q22}\left(0\right)=cos\left(\text{theta0}\right)\hfill \\ \hfill \text{q23}\left(0\right)=-sin\left(\text{theta0}\right)\hfill \\ \hfill \text{q31}\left(0\right)=0\hfill \\ \hfill \text{q32}\left(0\right)=sin\left(\text{theta0}\right)\hfill \\ \hfill \text{q33}\left(0\right)=cos\left(\text{theta0}\right)\hfill \end{array}\right)$$

Le système de 12 équations différentielles à intégrer est finalement :

equations=Join[euler,dQ]
			 

2. Intégration avec Mathematica

3. Intégration avec Scilab

La fonction suivante effectue l'intégration numérique pour des paramètres et des conditions initiales donnés :

function [t,y]=solution(r,a,b,omegaP,omega0,theta0,tmax,te),
 c=cos(theta0),
 s=sin(theta0),
 function deriv=systeme(t,y),
  deriv(1)=(1-r)*y(2)*y(3)+a*r*(y(6)*y(7)-y(4)*y(9)), // omegaS1
  deriv(2)=(r-1)*y(1)*y(3)+a*r*(y(6)*y(8)-y(5)*y(9)), // omegaS2
  deriv(3)=-b*y(3), // omegaS3
  deriv(4)=y(3)*y(5)-y(2)*y(6), // q11
  deriv(5)=y(1)*y(6)-y(3)*y(4), // q12
  deriv(6)=y(2)*y(4)-y(1)*y(5), // q13
  deriv(7)=y(3)*y(8)-y(2)*y(9), // q21
  deriv(8)=y(1)*y(9)-y(3)*y(7), // q22
  deriv(9)=y(2)*y(7)-y(1)*y(8), // q23
  deriv(10)=y(3)*y(11)-y(2)*y(12), // q31
  deriv(11)=y(1)*y(12)-y(3)*y(10), // q32
  deriv(12)=y(2)*y(10)-y(1)*y(11), // q33
 endfunction
 init=[0;omegaP*s;omega0+omegaP*c;1;0;0;0;c;-s;0;s;c];
 t=[0:te:tmax];
    tolA=1d-7;
    tolR=1d-10;
 y=ode(init,0,t,tolR,tolA,systeme);
endfunction
   

La fonction suivante calcule les différentes quantités à tracer :

function [a1,a2,a3]=calcul(y)
 a1 = y(1,:);
 a2 = a1
 a3 = a1;
 s=size(a1);
 for i=1:s(2),
  a1(i)=y(4,i)*y(4,i)+y(7,i)*y(7,i)+y(10,i)*y(10,i),
  a2(i)=y(12,i),
  a3(i)=y(6,i)/sqrt(1-y(12,i)*y(12,i)),
 end;
endfunction
   

Voici par exemple la toupie sans frottements déjà considérée plus haut :

tmax=500;
[t,y]=solution(0.2,2*%pi*0.1,0,0,2*%pi,30*%pi/180,tmax,0.1);
[a1,a2,a3]=calcul(y);
   

plotH=scf();
plot2d(t,a1,style=5,rect=[0,0.99,tmax,1.01]);
xtitle('q11^2+q21^2+q22^2','t','');
   

plotH.pdf

plotI=scf();
plot2d(t,a2,style=5,rect=[0,0,tmax,1]);
xtitle('cos(theta)','t','');
   

plotI.pdf

plotJ=scf();
plot2d(t,a3,style=5,rect=[0,-1,tmax,1]);
xtitle('sin(phi)','t','');
   

plotJ.pdf

plotK=scf();
param3d(y(6,:),y(9,:),y(12,:),35,45,'x@y@z',[3,4],[-1,1,-1,1,-1,1]);
   

plotK.pdf