ENUNCIADO EJEMPLO 7

     Un sólido rígido está formado por un disco pesado con centro en el punto A, de radio R y masa M y una varilla OA de masa m de longitud 2r que se encuentra articulada a un punto fijo O.

    El punto A se puede mover libremente en todas las direcciones del espacio. Además, el disco puede girar libremente alrededor de OA.

   Paso 0. Reiniciación de las variables del sistema y llamada a los paquetes linalg, plots y plottools.

> 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

> libname:="C:/",libname:

> with(mecapac3d):

   Paso 1. Definimos las coordenadas generalizadas del sistema en una lista que se denominará cg.. Al moverse A en todo el espacio y O estar fijo tendremos 3 grados de libertad que definiremos a continuación:

> cg := [alpha,theta,beta] :

   Paso 2. Definición mediante variables de los elementos que forman el sistema mecánico, la varilla y el disco.

Para la varilla:

> xgv := [r*sin(alpha)*sin(theta),-r*sin(alpha)*cos(theta),r*cos(alpha)];

xgv := [r*sin(alpha)*sin(theta), -r*sin(alpha)*cos(theta), r*cos(alpha)]

> rotv1 := rota(theta,3);

rotv1 := matrix([[cos(theta), -sin(theta), 0], [sin(theta), cos(theta), 0], [0, 0, 1]])

> rotv2 := rota(alpha,1);

rotv2 := matrix([[1, 0, 0], [0, cos(alpha), -sin(alpha)], [0, sin(alpha), cos(alpha)]])

> rottotv := evalm(rotv1 &* rotv2);

rottotv := matrix([[cos(theta), -sin(theta)*cos(alpha), sin(theta)*sin(alpha)], [sin(theta), cos(theta)*cos(alpha), -cos(theta)*sin(alpha)], [0, sin(alpha), cos(alpha)]])

> v1 := [varilla,xgv,rottotv,mvar,2*r] :

A continuación definimos igualmente el disco y sus correspondientes giros:

> xgd := [2*r*sin(alpha)*sin(theta),-2*r*sin(alpha)*cos(theta),2*r*cos(alpha)];

xgd := [2*r*sin(alpha)*sin(theta), -2*r*sin(alpha)*cos(theta), 2*r*cos(alpha)]

> rotd1 := rota(theta,3);

rotd1 := matrix([[cos(theta), -sin(theta), 0], [sin(theta), cos(theta), 0], [0, 0, 1]])

> rotd2 := rota(alpha,1);

rotd2 := matrix([[1, 0, 0], [0, cos(alpha), -sin(alpha)], [0, sin(alpha), cos(alpha)]])

> rotd3 := rota(beta,3);

rotd3 := matrix([[cos(beta), -sin(beta), 0], [sin(beta), cos(beta), 0], [0, 0, 1]])

> rottotd := evalm(rotd1 &* rotd2 &* rotd3);

rottotd := matrix([[cos(theta)*cos(beta)-sin(theta)*cos(alpha)*sin(beta), -cos(theta)*sin(beta)-sin(theta)*cos(alpha)*cos(beta), sin(theta)*sin(alpha)], [sin(theta)*cos(beta)+cos(theta)*cos(alpha)*sin...

> d1:= [disco,xgd,rottotd,mdis,R]:

   Paso 3. Definición de los elementos gráficos que definiran nuestro sistema de ejes.

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

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

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

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

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

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

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

   Paso 4. Definición de la variable sistema que agrupa en una lista todos los elementos anteriores.

> sistema := [v1,ejex,ejey,ejez,d1,TO,TX,TY,TZ]:

   Paso 5. Obtención de la energía cinética del sistema mediante fT asignándola a la variable T.

> T := fT(sistema);

T := 1/2*mvar*((r*cos(alpha)*alpha1*sin(theta)+r*sin(alpha)*cos(theta)*theta1)^2+(-r*cos(alpha)*alpha1*cos(theta)+r*sin(alpha)*sin(theta)*theta1)^2+r^2*sin(alpha)^2*alpha1^2)+1/6*alpha1^2*mvar*r^2+1/6...T := 1/2*mvar*((r*cos(alpha)*alpha1*sin(theta)+r*sin(alpha)*cos(theta)*theta1)^2+(-r*cos(alpha)*alpha1*cos(theta)+r*sin(alpha)*sin(theta)*theta1)^2+r^2*sin(alpha)^2*alpha1^2)+1/6*alpha1^2*mvar*r^2+1/6...

   Paso 6. Obtención de la energía potencial del sistema mediante fV asignándola a la variable V

> V := fV(sistema);

V := mvar*g*r*cos(alpha)+2*mdis*g*r*cos(alpha)

   Paso 7. Obtención de la lagrangiana como diferencia de energías entre la energía cinética y la potencial.

> L := simplify(T-V);

L := 1/2*mdis*R^2*beta1*cos(alpha)*theta1+2/3*alpha1^2*mvar*r^2+2*mdis*r^2*theta1^2-2/3*mvar*r^2*theta1^2*cos(alpha)^2-2*mdis*r^2*theta1^2*cos(alpha)^2+2*mdis*r^2*alpha1^2+1/8*mdis*R^2*alpha1^2+1/4*md...

   Paso 8. Obtención de las ecuaciones de lagrange para las dos coordenadas generalizadas mediante el operador Ec_lag

> ecua := simplify(ec_lag());

ecua := [4/3*diff(alpha(t), `$`(t, 2))*mvar*r^2+4*mdis*r^2*diff(alpha(t), `$`(t, 2))+1/4*mdis*R^2*diff(alpha(t), `$`(t, 2))+1/2*mdis*R^2*diff(beta(t), t)*sin(alpha(t))*diff(theta(t), t)-4/3*mvar*r^2*d...ecua := [4/3*diff(alpha(t), `$`(t, 2))*mvar*r^2+4*mdis*r^2*diff(alpha(t), `$`(t, 2))+1/4*mdis*R^2*diff(alpha(t), `$`(t, 2))+1/2*mdis*R^2*diff(beta(t), t)*sin(alpha(t))*diff(theta(t), t)-4/3*mvar*r^2*d...ecua := [4/3*diff(alpha(t), `$`(t, 2))*mvar*r^2+4*mdis*r^2*diff(alpha(t), `$`(t, 2))+1/4*mdis*R^2*diff(alpha(t), `$`(t, 2))+1/2*mdis*R^2*diff(beta(t), t)*sin(alpha(t))*diff(theta(t), t)-4/3*mvar*r^2*d...ecua := [4/3*diff(alpha(t), `$`(t, 2))*mvar*r^2+4*mdis*r^2*diff(alpha(t), `$`(t, 2))+1/4*mdis*R^2*diff(alpha(t), `$`(t, 2))+1/2*mdis*R^2*diff(beta(t), t)*sin(alpha(t))*diff(theta(t), t)-4/3*mvar*r^2*d...ecua := [4/3*diff(alpha(t), `$`(t, 2))*mvar*r^2+4*mdis*r^2*diff(alpha(t), `$`(t, 2))+1/4*mdis*R^2*diff(alpha(t), `$`(t, 2))+1/2*mdis*R^2*diff(beta(t), t)*sin(alpha(t))*diff(theta(t), t)-4/3*mvar*r^2*d...ecua := [4/3*diff(alpha(t), `$`(t, 2))*mvar*r^2+4*mdis*r^2*diff(alpha(t), `$`(t, 2))+1/4*mdis*R^2*diff(alpha(t), `$`(t, 2))+1/2*mdis*R^2*diff(beta(t), t)*sin(alpha(t))*diff(theta(t), t)-4/3*mvar*r^2*d...

   Paso 9. Asignación de valores numéricos a los parámetros que queden sun asignar para poder proceder a la integración numérica.

> mvar:=1: g:=9.8: mdis:=10: r:=1: R:=1:

   Paso 10. Visualizamos la configuracion del sistema para distintas posiciones:

Posición en el instante inicial

> fG([evalf(0),evalf(0),evalf(0)]);

[Plot]

En una posición generica

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

[Plot]

   Paso 11. Integración numérica del problema mediante la función fint asignando el resultado a la variable res.

> res := fint([0.1,0,0,0,0,50]):

   Paso 12. Procedemos a realizar una animación del movimiento del conjunto por medio de la función dibu3.

> dibu3(6,70);

[Plot]

>