ENUNCIADO DEL EJEMPLO 22

Sea un sistema formado por una varilla de masa m1 y longitud a que tiene uno de sus extremos fijo en el centro de gravedad y por otra varilla de masa m2 y longitud a que tiene un extremo enlazado con el otro extremo de la anterior y que está siempre en el mismo plano vertical que la anterior.  En el extremo que comparten ambas varillas hay una partícula de masa mp1 y en el extremo libre de la segunda partícula hay otra partícula de masa m2. El punto de enlace entre varillas se encuentra unido además a una deslizadera que se mueve por el eje Z por medio de un muelle de constante k y longitud natural l. Obtener las ecuaciones del movimiento del sistema.

> restart;

Cargamos los paquetes de Maple que vamos a emplear, entre ellos el Mecapac3d, para lo cuál es necesario indicar previamente donde se encuentra la librería dentro del disco duro.

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

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

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

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

> with(mecapac3d):

El sistema formado por las dos varillas tiene tres grados de libertad, ya que con dos parámetros (por ejemplo angulares) podemos determinar la posición de la primera varilla, mientras que con otro podemos determinar la posición de la segunda, ya que se cumple que están en el mismo plano vertical. Por ello emplearemos tres coordenadas generalizadas, theta,phi y psi. Theta representa el ángulo que forma la primera varilla con la vertical. Phi representa el ángulo que forma la segunda varilla con la vertical. Psi representa el ángulo que ha girado el conjunto alrededor del eje Z partiendo del plano (X,Z). Por tanto, definimos dichas coordenadas generalizadas.

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

cg := [theta, phi, psi]

Definimos a continuación los puntos o partículas que indica el enunciado.

> m1 := [punto,a*sin(theta)*cos(psi),a*sin(theta)*sin(psi),a*cos(theta),masa1]:

> m2:=[punto,a*(sin(theta)+sin(phi))*cos(psi),a*(sin(theta)+sin(phi))*sin(psi),a*cos(theta)-a*cos(phi),masa2]:

Definimos ahora las coordenadas de los centros de gravedad de las varillas.

> xgvar1:=[a*(sin(theta)/2)*cos(psi),a*(sin(theta)/2)*sin(psi),a*cos(theta)/2]:

> xgvar2:=[a*(sin(theta)+sin(phi)/2)*cos(psi),a*(sin(theta)+sin(phi)/2)*sin(psi),a*cos(theta)-a*cos(phi)/2]:

Definimos la matriz de rotación de la primera varilla considerando el giro alrededor de su eje propio x como giro relativo y multiplicando por lo tanto por la derecha. A continuación se defina la varilla.

> r1:=evalm(rota(Pi/2+psi,3)&*rota(theta,1)):

> var1:=[varilla,xgvar1,r1,mvar1,a]:

Definimos la matriz de rotación de la segunda varilla y la varilla

> r2:=evalm(rota(-Pi/2+psi,3)&*rota(phi,1)):

> var2:=[varilla,xgvar2,r2,mvar2,a]:

Definimos a continuación los extremos del muelle y el muelle.

> ext1:=[0,0,a*cos(theta)]:ext2:=[a*sin(theta)*cos(psi),a*sin(theta)*sin(psi),a*cos(theta)]:

> mue:=[muelle,ext1,ext2,k,lonnat]:

Definimos a continuación elementos gráficos para representar los ejes.

> s1 := [segmento,[0,0,0],[10,0,0],red] :

> s2 := [segmento,[0,0,0],[0,10,0],green] :

> s3 := [segmento,[0,0,0],[0,0,10],blue] :

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

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

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

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

Definimnos el sistema global

> sistema:=[m1,m2,mue,var2,var1,s1,s2,s3,TX,TY,TZ,TO];

sistema := [[punto, a*sin(theta)*cos(psi), a*sin(theta)*sin(psi), a*cos(theta), masa1], [punto, a*(sin(theta)+sin(phi))*cos(psi), a*(sin(theta)+sin(phi))*sin(psi), a*cos(theta)-a*cos(phi), masa2], [mu...sistema := [[punto, a*sin(theta)*cos(psi), a*sin(theta)*sin(psi), a*cos(theta), masa1], [punto, a*(sin(theta)+sin(phi))*cos(psi), a*(sin(theta)+sin(phi))*sin(psi), a*cos(theta)-a*cos(phi), masa2], [mu...sistema := [[punto, a*sin(theta)*cos(psi), a*sin(theta)*sin(psi), a*cos(theta), masa1], [punto, a*(sin(theta)+sin(phi))*cos(psi), a*(sin(theta)+sin(phi))*sin(psi), a*cos(theta)-a*cos(phi), masa2], [mu...sistema := [[punto, a*sin(theta)*cos(psi), a*sin(theta)*sin(psi), a*cos(theta), masa1], [punto, a*(sin(theta)+sin(phi))*cos(psi), a*(sin(theta)+sin(phi))*sin(psi), a*cos(theta)-a*cos(phi), masa2], [mu...

Calculamos las energías cinética y potencial para obtener después la Lagrangiana.

> T:= fT(sistema):

> V:=fV(sistema):

Lagrangiana del sistema

> L:=simplify(T-V);

L := masa2*a^2*sin(theta)*psi1^2*sin(phi)+masa2*a^2*cos(theta)*theta1*cos(phi)*phi1-masa2*g*a*cos(theta)+masa2*g*a*cos(phi)+1/6*phi1^2*mvar2*a^2+1/6*theta1^2*mvar1*a^2-1/2*mvar2*a^2*sin(theta)*theta1*...L := masa2*a^2*sin(theta)*psi1^2*sin(phi)+masa2*a^2*cos(theta)*theta1*cos(phi)*phi1-masa2*g*a*cos(theta)+masa2*g*a*cos(phi)+1/6*phi1^2*mvar2*a^2+1/6*theta1^2*mvar1*a^2-1/2*mvar2*a^2*sin(theta)*theta1*...L := masa2*a^2*sin(theta)*psi1^2*sin(phi)+masa2*a^2*cos(theta)*theta1*cos(phi)*phi1-masa2*g*a*cos(theta)+masa2*g*a*cos(phi)+1/6*phi1^2*mvar2*a^2+1/6*theta1^2*mvar1*a^2-1/2*mvar2*a^2*sin(theta)*theta1*...

Ecuaciones del movimiento.

> ecua:=simplify(ec_lag());

ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...ecua := [-1/6*a*(2*diff(psi(t), t)^2*mvar1*a*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*masa2*a*diff(psi(t), t)^2*cos(theta(t))*sin(theta(t))*(-a^2*(-1+cos(theta(t))^2))^(1/2)+6*m...

Damos valores a los parámetros para poder representar el sistema en una posición genérica dando valor a las coordenadas generalizadas y para poder realizar después la integración numérica.

> masa1:=2:masa2:=1:mvar1:=1:mvar2:=1:a:=10:lonnat:=0:k:=100:g:=9.8:

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

[Plot]

Realizamos la integración numérica indicando los valores iniciales de las coordenadas generalizadas y sus velocidades.

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

> odeplot(res,[t,theta(t)],0..2,color=blue);

[Plot]

> odeplot(res,[t,phi(t)],0..2,color=blue);

[Plot]

> odeplot(res,[t,psi(t)],0..2,color=blue);

[Plot]

> dibu3(2,25);

[Plot]

>