ENUNCIADO DEL EJEMPLO 24

Este ejecicio consta de una varilla unida por un extremo al eje ox y por el otro a un disco perpendicular a ella en todo momento. La varilla gira alredor del eje x, el disco gira a su vez con la varilla y puede girar alrededor del eje de la varilla. Obtener las ecuaciones del movimiento.

> restart;

Cargamos los paquetes de Maple que vamos a emplear, entre ellos el Mecapac3d, para lo cuál es necesario indicar antes la dirección en la que se encuentra dentro del disco duro.

> 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 formado por la varilla y el disco tiene tres grados de libertad, que identificaremos con las coordenadas generalizadas [x, theta, phi], que representan, respectivamente, la coordenada x en la que se encuentra el extremo de la varilla (y el resto de ella también), el ángulo girado por la varilla alrededor del eje x (partiendo de la vertical descendente y en sentido antihorario) y el ángulo girado por el disco alrededor de la dirección de la varilla. Definimos por tanto dichos grados de libertad.

> cg:=[x,theta,phi]:

Pasamos a continuación a definir los elementos del sistema, comenzando por la varilla, con la posición de su centro de masas y la matriz de rotación entre sus ejes propios y los ejes fijos, que es únicamente la que representa un giro alrededor del eje X de valor theta.

> v1:=[varilla,[x,(l/2)*sin(theta),-(l/2)*cos(theta)],rota(theta,1),mv,l]:

Definimos ahora el disco, calculando primero su matriz de rotación. Dado que el giro propio se indica de forma relativa a los ejes del disco, se realiza multiplicando por la derecha las matrices de rotación.

> r1:=rota(theta,1):r2:=rota(phi,3):

> rotp:=evalm(r1&*r2):

> d1:=[disco,[x,l*sin(theta),-l*cos(theta)],rotp,md,r]:

Definimos también unos elementos gráficos para representar los ejes.

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

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

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

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

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

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

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

> radio:=[segmento,[x,l*sin(theta),-l*cos(theta)],[x+r*cos(phi),l*sin(theta)+r*sin(phi)*cos(theta),-l*cos(theta)+r*sin(phi)*sin(theta)],red]:

> sistema:=[v1,d1,ejeX,ejeY,ejeZ,TO,TY,TX,TZ,radio]:

Calculamos ahora la Lagrangiana del sistema, calculando previamente la energía potencial y la cinética.

> V:=fV(sistema):

> T:=simplify(fT(sistema)):

> L:=simplify(T-V);

L := 1/2*mv*x1^2+1/6*theta1^2*mv*l^2+1/2*md*x1^2+1/2*md*l^2*theta1^2+1/8*theta1^2*md*r^2+1/4*phi1^2*md*r^2+1/2*mv*g*l*cos(theta)+md*g*l*cos(theta)

Obtenemos ahora las tres ecuaciones del movimiento.

> ecua:=map(simplify,ec_lag());

ecua := [diff(x(t), `$`(t, 2))*(md+mv), 1/3*diff(theta(t), `$`(t, 2))*mv*l^2+md*l^2*diff(theta(t), `$`(t, 2))+1/4*diff(theta(t), `$`(t, 2))*md*r^2+1/2*mv*g*l*sin(theta(t))+md*g*l*sin(theta(t)), 1/2*di...

Damos ahora unos valores a los parámetros que antes dejamos como constantes sin valor para poder obtener una representación del sistema para unos valores determinados de las coordenadas generalizadas, y para poder llevar a cabo la integración numérica.

> mv:=1:md:=2:l:=3:r:=1:g:=9.8:

> fG([1,evalf(5*Pi/3),evalf(Pi/2)]);

[Plot]

Realizamos ahora la integración numérica, para lo cual es necesario indicar los valores de las coordenadas generalizadas y sus velocidades en el instante inicial.

> res:=fint([1,0.2,evalf(5*Pi/3),0.1,evalf(Pi/2),0.1]):

Obtenemos ahora una representación animada del movimiento y las gráficas de las coordenadas generalizadas en función del tiempo.

> dibu3(7,30);

[Plot]

> g1:=odeplot(res,[t,theta(t)],0..5,color=red):

> g2:=odeplot(res,[t,phi(t)],0..5,color=black):

> g3:=odeplot(res,[t,x(t)],0..5,color=blue):

> display([g1]);

[Plot]

> display([g2]);

[Plot]

> display([g3]);

[Plot]

>