Précession elliptique du pendule sphérique
|
> |
restart; |
|
> |
with(plots):with(linalg): |
Nous allons vous exposer une méthode
permettant de calculer l’angle de précession pour le pendule sphérique.
|
|
|
Etablissement des équations différentielles du mouvement:
On
part de la RFD appliquée au pendule sphérique :
Bilan
des forces
: Poids
:
P = mgêz
Tension du fils : T=Têr
On obtient donc l'équation suivante:
|
> |
m*diff(r(t)êr,t$2)=m*g*êz-Têr; |
|
|
|
On exprime ensuite les vecteurs x, y et z en coordonnées sphériques pour pouvoir sélectionner l’équation en Thêta qui nous intéresse pour le calcul de la précession elliptique. Donc par projection dans la base sphérique, on obtient :
|
> |
êx:=êr*sin(theta)*cos(phi)+ê(theta)*cos(theta)*cos(phi)-ê(phi)*sin(phi);
|
On obtient donc la RFD en coordonnées sphériques :
|
> |
eq1:=T=m*g*cos(theta(t))+m*l*((diff(theta(t),t))^2+(diff(phi(t),t)^2*(sin(theta(t)))^2));
|
A ce stade du calcul, il est préférable de passer en coordonnées cartésiennes dans l'équation 2 qui correspond à notre équation suivant ê(theta).En effet, l'avantage de ces coordonnées est que pour des déplacements infinitésimaux, on pourra négliger certains termes selon x et y alors qu'en coordonnées sphériques, seul theta pourra être négligé. Il est donc nécessaire d'établir quelques relations nous permettant de procéder à ce changement de coordonnées.
Utilisons donc les relations de x, y et z en fonction de theta et phi:
|
> |
x(t):=l*sin(theta(t))*cos(phi(t));
|
|
> |
aa:=simplify(x(t)^2+y(t)^2);
|
De plus, en simplifiant l'expression ci-dessous par sin(theta), on retrouve la même expression que celle de eq3. On en déduit donc que l'expression ci-dessous est nulle:
|
> |
b:=diff(diff(phi(t),t)*sin(theta(t))^2,t); |
L'expression précédente peut s'interpréter comme la conservation de moment cinétique. On applique le moment cinétique selon (Oz) pour établir une relation très importante qui nous permettra de passer d'une dépendance en x et y à une dépendance en x ou y uniquement.
|
> |
OP := vector(['x','y','z']); |
|
> |
v := vector([diff('x(t)',t),diff('y(t)',t),diff('z(t)',t)]); |
|
> |
L[0]:=crossprod(OP,v); |
|
> |
L[Delta]:=dotprod(L[0],vector([0,0,1])); |
Or, d'après "ae" et "b", on sait que le dérivée de cette expression en fonction de t est nulle:
|
> |
c:=diff('x'*(diff('y'(t), t))-'y'*(diff('x'(t),t)),t); |
Ainsi en remplaçant dans eq2, à l'aide les différentes relations entre x, y,z et phi, theta, on obtient ainsi l'équation eq2 en coordonnées cartésiennes:
|
> |
eq22:='x(t)'*diff('y(t)',t$2)+'y(t)'*diff('y(t)',t$2)+1/l^2*('x(t)'*diff('x(t)',t)+'y(t)'*diff('y(t)',t))^2*1/(1-r^2/l^2)+1/l^2*('x(t)'*diff('y(t)',t)-'y(t)'*diff('x(t)',t))+w0^2*r^2*sqrt(1-r^2/l^2)=0; |
On va ensuite simplifier notre équation eq22 et surtout la partie:" 1/l^2*(x(t)*diff(x(t),t)+y(t)*diff(y(t),t))^2*1/(1-r^2/l^2) +1/l^2*(x(t)*diff(y(t),t) - y(t)*diff(x(t),t)) ", toujours à l'aide des équations aa, ab, ac... mais aussi de la relation sin(tetha)=r/l:
|
> |
eq23:='x(t)'*diff('y(t)',t$2)+'y(t)'*diff('y(t)',t$2)+r^2/l^2*(diff('x(t)',t)^2+diff('y(t)',t)^2)+r^2/(l^4*(1-r^2/l^2))*('x(t)'*diff('x(t)',t)+'y(t)'*diff('y(t)',t))^2+w0^2*r^2*sqrt(1-r^2/l^2)=0; |
Ensuite en utilisant la relation xÿ=¨xy, on pourra supprimer dans l'équation eq22 le terme en ÿ ou en ¨x, et finalement en considérant de petites oscillations, c'est à dire que x,y << l, on obtient les deux équations suivantes:
|
> |
da:=diff('x(t)',t$2)+omega0^2*'x(t)'-(1/2)*omega0^2*('x(t)'/l^2)*('x(t)'^2+'y(t)'^2)+('x(t)'/l^2)*(diff('x(t)',t)^2+diff('y(t)',t)^2)=0;
|
|
> |
restart; |
Déplacements infinitésimaux :
Attardons nous sur les déplacements infinitésimaux:
En négligeant les termes en x ou y d'ordre supérieures à 2 , on obtient ces deux équations en posant les conditions initiales du pendule sphérique:
|
> |
eq3:=diff(x(t),t$2)+omega0^2*x(t)=0;
|
|
> |
solx:=dsolve({eq3,x(0)=A,D(x)(0)=0},x(t));
|
|
> |
A:=a; |
On obtient un système d'équations paramétrées caractéristique d'une ellipse, parcourue à la vitesse angulaire omega0 et dont les axes principaux sont fixes.
Déplacements non infinitésimaux :
Cependant pour les déplacements non infinitésimaux, on observe bien une précession elliptique : c'est-à-dire qu'au cours du déplacement du pendule, l'axe principal de l'ellipse décrite par la trajectoire subit une rotation. Comme on est dans un référentiel galiléen et que l'on a négligé les frottements, il parait normal que la vitesse de rotation Omega soit constante. Nous posons donc un nouveau repère (O,X,Y,z) qui va suivre la rotation de l'axe principal de l'ellipse. Ainsi, en négligeant certains termes, nous considèrerons que le mouvement dans ce nouveau repère est quasiment identique à celui obtenu lors de déplacements infinitésimaux. Soit alpha l'angle de précession. Comme Omega est constant et en considérant que alpha est nul à l'instant t=0, on a : alpha=Omega*t. Par la projection des vecteurs ou par emploi d'une matrice rotation, on obtient:
|
> |
alpha:=Omega*t;
|
Or on considère, comme indiqué précedemment, que le mouvement dans le nouveau repère est elliptique de pulsation omega proche de omega0, on peut donc écrire :
|
> |
X:=a*cos(omega*t);
|
Ce qui donne finalement:
|
> |
x; |
En prenant en compte que w-w0<<w0 et W<<w0, on exprime grâce aux égalités précédentes les termes qui apparaissent dans l'équation da. Ainsi, nous pourrons les réinjecter dans celle-ci pour obtenir une nouvelle équation.
|
> |
ea:=simplify(diff(x,t$2)+omega0^2*x); |
|
> |
eaa:=factor(-a*cos(omega*t)*omega^2*cos(Omega*t)-a*cos(omega*t)*cos(Omega*t)*Omega^2-2*b*cos(omega*t)*omega*cos(Omega*t)*Omega+omega0^2*a*cos(omega*t)*cos(Omega*t)); |
|
> |
eab:=factor(2*a*sin(omega*t)*omega*sin(Omega*t)*Omega+b*sin(omega*t)*omega^2*sin(Omega*t)+b*sin(omega*t)*sin(Omega*t)*Omega^2-omega0^2*b*sin(omega*t)*sin(Omega*t)); |
or ea=eaa+eab:
|
> |
eaa+eab; |
Voici donc l'expression de "diff(x,t$2)+w0²*x". Donnons maintenant les expressions de "x²+y²" et de "diff(x,t)²+diff(y,t)²":
|
> |
eb:=simplify(x^2+y^2); |
|
> |
ec:=factor(simplify(diff(x,t)^2+diff(y,t)^2)); |
On exprime ensuite les expressions "x(x^2+y^2)" et "x(diff(x(t),t)^2+diff(y(t),t)^2)" ,de notre équation du pendule sphérique, en fonction de omega et Omega pour pouvoir ensuite regrouper les termes en cosWt.coswt et en sinWt.sinwt et les comparer avec ceux de l'équation ea. On obtient donc, en supprimant les plus grandes harmonique puisque leurs coefficients sont très inférieurs devant les coefficients des plus faibles harmoniques, les deux équations suivantes:
|
> |
'x(t)'*(diff('x(t)',t)^2+diff('y(t)',t)^2)=omega^2/4*((a^3+3*a*b^2)*cos(omega*t)*cos(Omega*t)-(3*b*a^2+b^3)*sin(omega*t)*sin(Omega*t)); |
|
> |
'x(t)'*('x(t)'^2+'y(t)'^2)=1/4*((3*a^3+a*b^2)*cos(omega*t)*cos(Omega*t)-(b*a^2+3*b^3)*sin(omega*t)*sin(Omega*t)); |
On obtient donc un système de deux équations à deux inconnues:
|
> |
solve({a*(omega0^2-omega^2)-2*b*omega0*Omega-(omega0^2*a/(8*l^2))*(a^2-5*b^2),b*(omega0^2-omega^2)-2*a*omega0*Omega-(omega0^2*b/(8*l^2))*(b^2-5*a^2)},{Omega,omega^2}); |
Nous obtenons ainsi les expressions de w² et W. Calculons maintenant l'angle a[n] formé au cours de n oscillations du pendule. On prendra w0 comme pulsation propre des oscillations donc T=2*p/w0 ( alors que T vaut en toute rigueur 2*p/w)
|
> |
Omega:=3/8*b*omega0*a/l^2;
|
|
> |
alpha; |
|
> |
alpha[n]:=n*subs(t=T,alpha);
|
Voici donc l'expression de l'angle de précession.
Calcul complémentaire:
En complément, nous pouvons vérifier la cohérence de notre résultat en étudiant un cas particulier du pendule sphérique : le pendule conique. On lance le pendule de telle manière que la trajectoire soit un cercle. Ce qui nous donne a=b=r et theta est une constante. L'équation eq2 devient donc:
|
> |
EQ2:=-(diff(phi(t),t))^2*sin(theta)*cos(theta)+omega0^2*sin(theta)=0; |
On exprime diff(phi(t),t) qui est la solution positive (sol[1]).
|
> |
sol:=solve(EQ2,diff(phi(t),t)); |
|
> |
sol[1]; |
Or:
|
> |
theta:=arccos(sqrt(1-(r/l)^2)); |
Effectuons un développement limité à l'ordre 1 en considérant r<<l. On pose par exemple u =r/l
|
> |
sol[1]; |
D'une autre manière, on remarque que diff(f(t),t)=w+W. En utilisant les expressions de W et w² que l'on a obtenues, il vient:
|
> |
dphi:=omega+Omega; |
|
> |
a:=r; |
|
> |
dphi; |
on pose de nouveau u=r/l et on fait un développement limité à l'ordre 2:
|
> |
assume(omega0>0):Dphi:=taylor(1/4*(-2*omega0^2*(-8+2*u^2))^(1/2)+3/8*u^2*omega0,u=0,3); |
On retrouve bien notre expression.