ENUNCIADO DEL EJEMPLO 30

  Una varilla de masa m se halla ensartada en un aro circular de radio R por uno de sus extremos, pudiendo deslizar éste libremente sobre el aro, sometida a su propio peso. A su vez, al aro puede girar respecto de su diámetro vertical. Obtener las ecuaciones del movimiento.

> restart:

Cargamos los paquetes que vamos a emplear.

> 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

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

> with(mecapac3d):

El sistema tiene tres grados de libertad, que quedan definidos con las coordenadas generalizadas siguientes: theta, phi, beta. Theta representa el giro del aro alrededor del eje z. Phi es el ángulo que nos indica la situación del extremo de la varilla en contacto con el aro. Beta es el ángulo girado por la varilla.

> cg:=[theta,phi,beta];

cg := [theta, phi, beta]

Definimos a continuación el aro y la varilla, con  las coordenadas de sus centro de gravedad y sus matrices de rotación. Para calcular las matrices de rotaciñon se han descompuesto en giros absolutos, de forma que las matrices se han multiplicado por la izquierda.

> xg:=[0,0,0]:

> ar:=[aro,xg,rottot,m1,r]:

> rot1:=rota(Pi/2,1):

> rot2:=rota(theta,3):

> rottot:=evalm(rot2&*rot1):

> va:=[varilla,[r*sin(phi)*cos(theta)+h/2*sin(beta)*cos(theta),r*sin(phi)*sin(theta)+h/2*sin(beta)*sin(theta),-r*cos(phi)-h/2*cos(beta)],rotvar,m2,h]:

> rot3:=rota(-beta,2):

> rot4:=rota(theta,3):

> rotvar:=evalm(rot4&*rot3):

Definimos a continuación elementos gráficos, como los ejes del sistema de referencia fijo o los ángulos.

> a1:=[angulo,[0,0,-r],[0,0,0],[r*sin(phi)*cos(theta),r*sin(phi)*sin(theta),-r*cos(phi)],1]:

> s1:=[segmento,[0,0,0],[r*sin(phi)*cos(theta),r*sin(phi)*sin(theta),-r*cos(phi)],red]:

> ejeY:=[vector,[0,-r,0],[0,r,0],green]:

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

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

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

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

> TZ := [texto,[0,0,r+1],"Z"]:

> TX := [texto,[r+1,0,0],"X"]:

> ep:=[vector,[0,0,0],vector([7*cos(theta),7*sin(theta),0]),black]:

> v1:=[vector,[r*sin(phi)*cos(theta),r*sin(phi)*sin(theta),-r*cos(phi)],vector([2*cos(phi)*cos(theta),2*cos(phi)*sin(theta),2*sin(phi)]),green]:

> v2:=[vector,[r*sin(phi)*cos(theta),r*sin(phi)*sin(theta),-r*cos(phi)],vector([-2*sin(theta),2*cos(theta),0]),yellow]:

Finalmente definimos el sistema con todos sus elementos.

> sistema:=[ar,va,a1,ejeX,ejeY,ejeZ,TO,TX,TY,TZ,v1,v2,s1];

sistema := [[aro, [0, 0, 0], rottot, m1, r], [varilla, [r*sin(phi)*cos(theta)+1/2*h*sin(beta)*cos(theta), r*sin(phi)*sin(theta)+1/2*h*sin(beta)*sin(theta), -r*cos(phi)-1/2*h*cos(beta)], rotvar, m2, h]...sistema := [[aro, [0, 0, 0], rottot, m1, r], [varilla, [r*sin(phi)*cos(theta)+1/2*h*sin(beta)*cos(theta), r*sin(phi)*sin(theta)+1/2*h*sin(beta)*sin(theta), -r*cos(phi)-1/2*h*cos(beta)], rotvar, m2, h]...sistema := [[aro, [0, 0, 0], rottot, m1, r], [varilla, [r*sin(phi)*cos(theta)+1/2*h*sin(beta)*cos(theta), r*sin(phi)*sin(theta)+1/2*h*sin(beta)*sin(theta), -r*cos(phi)-1/2*h*cos(beta)], rotvar, m2, h]...sistema := [[aro, [0, 0, 0], rottot, m1, r], [varilla, [r*sin(phi)*cos(theta)+1/2*h*sin(beta)*cos(theta), r*sin(phi)*sin(theta)+1/2*h*sin(beta)*sin(theta), -r*cos(phi)-1/2*h*cos(beta)], rotvar, m2, h]...

Calculamos ahora la energía cinética y la potencial para obtener así la Lagrangiana.

> T:=fT(sistema):

> V:=fV(sistema):

> L:=simplify(T-V);

L := 1/6*m2*h^2*theta1^2+1/6*beta1^2*m2*h^2+1/4*theta1^2*m1*r^2+m2*g*r*cos(phi)+1/2*m2*g*h*cos(beta)+1/2*m2*r*sin(phi)*phi1*h*sin(beta)*beta1-1/2*m2*r^2*theta1^2*cos(phi)^2-1/6*m2*h^2*theta1^2*cos(bet...

Ecuaciones del movimiento.

> ecua:=ec_lag();

ecua := [1/3*m2*h^2*diff(theta(t), `$`(t, 2))+1/2*diff(theta(t), `$`(t, 2))*m1*r^2-m2*r^2*diff(theta(t), `$`(t, 2))*cos(phi(t))^2+2*m2*r^2*diff(theta(t), t)*cos(phi(t))*sin(phi(t))*diff(phi(t), t)-1/3...ecua := [1/3*m2*h^2*diff(theta(t), `$`(t, 2))+1/2*diff(theta(t), `$`(t, 2))*m1*r^2-m2*r^2*diff(theta(t), `$`(t, 2))*cos(phi(t))^2+2*m2*r^2*diff(theta(t), t)*cos(phi(t))*sin(phi(t))*diff(phi(t), t)-1/3...ecua := [1/3*m2*h^2*diff(theta(t), `$`(t, 2))+1/2*diff(theta(t), `$`(t, 2))*m1*r^2-m2*r^2*diff(theta(t), `$`(t, 2))*cos(phi(t))^2+2*m2*r^2*diff(theta(t), t)*cos(phi(t))*sin(phi(t))*diff(phi(t), t)-1/3...ecua := [1/3*m2*h^2*diff(theta(t), `$`(t, 2))+1/2*diff(theta(t), `$`(t, 2))*m1*r^2-m2*r^2*diff(theta(t), `$`(t, 2))*cos(phi(t))^2+2*m2*r^2*diff(theta(t), t)*cos(phi(t))*sin(phi(t))*diff(phi(t), t)-1/3...ecua := [1/3*m2*h^2*diff(theta(t), `$`(t, 2))+1/2*diff(theta(t), `$`(t, 2))*m1*r^2-m2*r^2*diff(theta(t), `$`(t, 2))*cos(phi(t))^2+2*m2*r^2*diff(theta(t), t)*cos(phi(t))*sin(phi(t))*diff(phi(t), t)-1/3...ecua := [1/3*m2*h^2*diff(theta(t), `$`(t, 2))+1/2*diff(theta(t), `$`(t, 2))*m1*r^2-m2*r^2*diff(theta(t), `$`(t, 2))*cos(phi(t))^2+2*m2*r^2*diff(theta(t), t)*cos(phi(t))*sin(phi(t))*diff(phi(t), t)-1/3...ecua := [1/3*m2*h^2*diff(theta(t), `$`(t, 2))+1/2*diff(theta(t), `$`(t, 2))*m1*r^2-m2*r^2*diff(theta(t), `$`(t, 2))*cos(phi(t))^2+2*m2*r^2*diff(theta(t), t)*cos(phi(t))*sin(phi(t))*diff(phi(t), t)-1/3...

Damos valores a los parámetros que quedaban por indicar para obtener una representación del sistema en un instante determinado y para poder integrar numéricamente.

> m1:=10: m2:=1: g:=9.81: r:=2: h:=3:

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

[Plot]

Integración numérica dando los valores de las coordenadas y velocidades iniciales.

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

Gráficas de las coordenadas en el tiempo y enfrentadas.

> odeplot(res,[t,phi(t)],0..2.64);

[Plot]

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

[Plot]

> odeplot(res,[t,beta(t)],0..2.64);

[Plot]

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

[Plot]

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

[Plot]

> odeplot(res,[beta(t),phi(t)],0..2.64);

[Plot]

Gráficos  de las velocidades generalizadas en función del tiempo y enfrentadas.

> odeplot(res,[t,diff(phi(t),t)],0..2.64);

[Plot]

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

[Plot]

> odeplot(res,[t,diff(beta(t),t)],0..2.64);

[Plot]

> odeplot(res,[theta(t),diff(phi(t),t)],0..2.64);

[Plot]

> odeplot(res,[beta(t),diff(phi(t),t)],0..2.64);

[Plot]

> odeplot(res,[phi(t),diff(theta(t),t)],0..2.64);

[Plot]

> odeplot(res,[diff(beta(t),t),diff(theta(t),t)],0..2.64);

[Plot]

> odeplot(res,[diff(phi(t),t),diff(theta(t),t)],0..2.64);

[Plot]

> odeplot(res,[diff(beta(t),t),diff(phi(t),t)],0..2.64);

[Plot]

Animación del movimiento.

> dibu3(2.64,20);

[Plot]

>