ENUNCIADO DEL EJEMPLO 31

Sea el sistema material formado por un disco de radio R y masa m que rueda sin deslizar por el plano horizontal describiendo una circunferencia de radio d y estando vertical en todo momento, y por una partícula de masa M que está ensartada en el borde del círculo, pudiendo deslizar sobre dicho borde como si estuviera ensartada en un aro liso. Se pide obtener las ecuaciones del movimiento.

Cargamos los paquetes de Maple que vamos a emplear.

> restart:

> with(linalg):with(plots):with(plottools):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

Cargamos también el paquete mecapac3d, para lo cuál es necesario indicar previamente donde se encuentra la librería en el disco duro.

> libname:="C:\",libname:

> with(mecapac3d):

El sistema formado por el disco y la partícula tiene dos grados de libertad, que quedan perfectamente definidos con la coordenadas generalizadas psi y theta. Psi es el ángulo girado por el disco al rodar sin deslizar y theta es el ángulo que nos indica la posición de la partícula en el disco.

> cg:=[psi,theta];

cg := [psi, theta]

Para definir el sistema emplearé un subsistema de referencia que acompañe al sólido en su movimiento pero que no acompañe al disco en su giro alrededor de su eje perpendicular, que se traslade y describa la circunferencia que describe el sistema (lo que supone girar alrededor del eje Z del sistema absoluto). Definiré, por este orden, los siguientes elementos: la posición del centro del sistema de referencia móvil, la matriz de rotación del subsistema, las matrices de rotación cuyo producto por la derecha (giros relativos) nos da la matriz rotación del disco en el subsistema y dicha matriz de rotación.

> censub:=[d*cos(psi),d*sin(psi),R]:

> rotsub:=rota(psi,3):

> rotprimera:=rota(Pi/2,2):

> rotsegunda:=rota(d/R*psi,1):

> rotdis:=evalm(rotsegunda&*rotprimera):

Defino el disco y la partícula, así como unos elementos gráficos para representar los ejes del subsistema.

> dis:=[disco,[0,0,0],rotdis,m,R]:

> par:=[punto,0,R*cos(theta),R*sin(theta),M]:

> ejexx:=[vector,[0,0,0],[3,0,0],blue]:ejeyy:=[vector,[0,0,0],[0,3,0],blue]:ejezz:=[vector,[0,0,0],[0,0,3],blue]:

Defino el subsistema con todos los elementos que van referidos a él.

> subsis:=[subsistema2,censub,rotsub,[dis,par,ejexx,ejeyy,ejezz]]:

Defino ahora unos elementos gráficos para representar los ejes del sistema de referencia fijo.

> ejeY:=[vector,[0,0,0],[0,5,0],green]:

> ejeX:=[vector,[0,0,0],[5,0,0],red]:

> ejeZ:=[vector,[0,0,0],[0,0,5],blue]:

> TO := [texto,[0,0,-1],"O"]:

> TY := [texto,[0,5,1],"Y"]:

> TZ := [texto,[0,0,6],"Z"]:

> TX := [texto,[6,0,0],"X"]:

Defino el sistema general con todos los elementos referidos a él, incluido el subsistema.

> sistema:=[subsis,ejeX,ejeY,ejeZ,TO,TX,TY,TZ];

sistema := [[subsistema2, [d*cos(psi), d*sin(psi), R], rotsub, [[disco, [0, 0, 0], rotdis, m, R], [punto, 0, R*cos(theta), R*sin(theta), M], [vector, [0, 0, 0], [3, 0, 0], blue], [vector, [0, 0, 0], [...sistema := [[subsistema2, [d*cos(psi), d*sin(psi), R], rotsub, [[disco, [0, 0, 0], rotdis, m, R], [punto, 0, R*cos(theta), R*sin(theta), M], [vector, [0, 0, 0], [3, 0, 0], blue], [vector, [0, 0, 0], [...sistema := [[subsistema2, [d*cos(psi), d*sin(psi), R], rotsub, [[disco, [0, 0, 0], rotdis, m, R], [punto, 0, R*cos(theta), R*sin(theta), M], [vector, [0, 0, 0], [3, 0, 0], blue], [vector, [0, 0, 0], [...

Obtenemos ahora las energías del sistema y la Lagrangiana.

> T:=fT(sistema):

> V:=fV(sistema):

> L:=simplify(T-V);

L := 1/2*M*R^2*theta1^2-M*g*R*sin(theta)-m*g*R+1/8*psi1^2*m*R^2+1/2*m*d^2*psi1^2-M*d*psi1*R*sin(theta)*theta1+1/2*M*psi1^2*R^2*cos(theta)^2+1/2*M*d^2*psi1^2-M*g*R

Ecuaciones del movimiento

> ecua:=ec_lag();

ecua := [1/4*diff(psi(t), `$`(t, 2))*m*R^2+m*d^2*diff(psi(t), `$`(t, 2))-M*d*R*cos(theta(t))*diff(theta(t), t)^2-M*d*R*sin(theta(t))*diff(theta(t), `$`(t, 2))+M*diff(psi(t), `$`(t, 2))*R^2*cos(theta(t...ecua := [1/4*diff(psi(t), `$`(t, 2))*m*R^2+m*d^2*diff(psi(t), `$`(t, 2))-M*d*R*cos(theta(t))*diff(theta(t), t)^2-M*d*R*sin(theta(t))*diff(theta(t), `$`(t, 2))+M*diff(psi(t), `$`(t, 2))*R^2*cos(theta(t...

Damos valores a los parámetros para poder representar el sistema en una situación determinada y para poder realizar después la integración numérica.

> g:=9.8:m:=1:M:=1:R:=1:d:=3:

> fG([evalf(Pi/4),evalf(Pi/2)]);

[Plot]

Integración numérica, indicando valores iniciales de coordenadas y velocidades generalizadas.

> res:=fint([0,1,0,0]):

Gráficas del movimiento

> odeplot(res,[t,theta(t)],0..3);

[Plot]

> odeplot(res,[t,psi(t)],0..3);

[Plot]

Animación del movimiento

> dibu3(5.7,35);

[Plot]

>