ENUNCIADO DEL EJEMPLO 28

Se considera un sistema formado por dos discos iguales de masa m y radio r, cuyos centros se encuentran unidos perpendicularmente mediante una varilla de longitud l y masa despreciable. Uno de los discos (A) puede girar libremente alrededor del eje de la varilla, mientras que el otro (B), se encuentra soldado a la varilla por su centro. Se considera que en todo instante el giro del disco B alrededor de la varilla es nulo, y que a su vez la varilla no gira alrededor de su propio eje.

El conjunto formado por la varilla y los discos se apoya sin rozamiento sobre un soporte fijo articulado que permite tanto el balanceo de la varilla (inclinación respecto del plano horizontal) como su movimiento de precesión alrededor del eje de apoyo. El punto de apoyo O en la varilla se encuentra en todo momento a una distancia l/3 del disco A y 2l/3 del disco B.

Se pide:

1. Expresar las ecuaciones de segundo orden que gobiernan el movimiento del sistema.

2. Representar las ecuaciones primeras del movimiento.

3.  Representar el caso en el que el movimiento resultante sea tal que la precesión del conjunto sea constante con la varilla siempre horizontal.

> restart;

Cargamos los paquetes del Maple que vamos a necesitar, entre ellos, mecapac3d, para el cuál es necesario indicar la situación de la librería en el 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 tiene tres grados de libertad, por lo que introduzco como coordenadas generalizadas los ángulos de Euler; precesión(y), nutación (Q) y rotación propia(F).

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

Defino las coordenadas de los centros de gravedad (hay que multiplicar por constantes para cada elemento):

> xgv:=[l/3*cos(theta)*cos(psi),l/3*cos(theta)*sin(psi),l/3*sin(theta)];

xgv := [1/3*l*cos(theta)*cos(psi), 1/3*l*cos(theta)*sin(psi), 1/3*l*sin(theta)]

Defino todas las rotaciones que van a experimentar los dos discos y la varilla, empleando matrices abslutas y locales para el disco A en su giro de precesión alrededor de la varilla.

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

> r2:=rota(-psi,1):

> r3:=rota(-theta,2):

> r4:=rota(phi,3):

> rtot:=evalm(r1 &* r2 &* r3):

> rotda:=evalm(r1 &* r2 &* r3 &* r4):

Defino los elementos que van a conformar el sistema:

> v1:=[varilla,(1/2)*xgv,rtot,0,l]:

> da:=[disco,-xgv,rotda,m,r]:

> db:=[disco,2*xgv,rtot,m,r]:

En el sistema incluimos unos elementos gráficos para representar los ejes del sistema de referencia fijo.

> ejex:=[segmento,[0,0,0],[2,0,0],green]:

> ejey:=[segmento,[0,0,0],[0,2,0],red]:

> ejez:=[segmento,[0,0,0],[0,0,2],blue]:

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

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

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

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

> sistema:=[v1,da,db,ejex,ejey,ejez, TO, TX,TY,TZ];

sistema := [[varilla, [1/6*l*cos(theta)*cos(psi), 1/6*l*cos(theta)*sin(psi), 1/6*l*sin(theta)], rtot, 0, l], [disco, [-1/3*l*cos(theta)*cos(psi), -1/3*l*cos(theta)*sin(psi), -1/3*l*sin(theta)], rotda,...sistema := [[varilla, [1/6*l*cos(theta)*cos(psi), 1/6*l*cos(theta)*sin(psi), 1/6*l*sin(theta)], rtot, 0, l], [disco, [-1/3*l*cos(theta)*cos(psi), -1/3*l*cos(theta)*sin(psi), -1/3*l*sin(theta)], rotda,...sistema := [[varilla, [1/6*l*cos(theta)*cos(psi), 1/6*l*cos(theta)*sin(psi), 1/6*l*sin(theta)], rtot, 0, l], [disco, [-1/3*l*cos(theta)*cos(psi), -1/3*l*cos(theta)*sin(psi), -1/3*l*sin(theta)], rotda,...

Cálculo de la lagrangiana del sistema para sacar las ecuaciones del movimiento:

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

> V:= simplify(fV(sistema)):

> L:=simplify(T-V);

L := 1/36*m*(10*l^2*cos(theta)^2*psi1^2+10*l^2*theta1^2-9*cos(theta)^2*psi1^2*r^2+18*r^2*sin(theta)*psi1*phi1+18*psi1^2*r^2+9*theta1^2*r^2+9*r^2*phi1^2-12*g*l*sin(theta))

Ecuaciones del movimiento

> ecua:= ec_lag();

ecua := [1/36*m*(-40*l^2*cos(theta(t))*diff(psi(t), t)*sin(theta(t))*diff(theta(t), t)+20*l^2*cos(theta(t))^2*diff(psi(t), `$`(t, 2))+36*cos(theta(t))*diff(psi(t), t)*r^2*sin(theta(t))*diff(theta(t), ...ecua := [1/36*m*(-40*l^2*cos(theta(t))*diff(psi(t), t)*sin(theta(t))*diff(theta(t), t)+20*l^2*cos(theta(t))^2*diff(psi(t), `$`(t, 2))+36*cos(theta(t))*diff(psi(t), t)*r^2*sin(theta(t))*diff(theta(t), ...ecua := [1/36*m*(-40*l^2*cos(theta(t))*diff(psi(t), t)*sin(theta(t))*diff(theta(t), t)+20*l^2*cos(theta(t))^2*diff(psi(t), `$`(t, 2))+36*cos(theta(t))*diff(psi(t), t)*r^2*sin(theta(t))*diff(theta(t), ...

Para poder integrar numéricamente, asigno valores a los distintos parámetros:

> m:=1: r:=1: l:=3: g:=9.8:

Realizamos la integración indicando los valores de las coordenadas generalizadas en el instante inicial y las velocidades generalizadas. Obtenemos después la representación de las coordenadas en el tiempo.

> res:=fint([0,0.5,0,0,0,5]):

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

[Plot]

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

[Plot]

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

[Plot]

Animación del movimiento.

> dibu3(10,100);

[Plot]

Establezco una ligadura anholónoma, de forma que Q sea 0 en todo el movimiento. (En realidad, este paso indica que theta es constante)

> Phi:=[theta1]:

Nuevas ecuaciones del movimiento.

> ecua:= simplify(ec_lag());

ecua := [-9*cos(theta(t))*diff(psi(t), t)*sin(theta(t))*diff(theta(t), t)+9/2*cos(theta(t))^2*diff(psi(t), `$`(t, 2))+1/2*cos(theta(t))*diff(theta(t), t)*diff(phi(t), t)+1/2*sin(theta(t))*diff(phi(t),...ecua := [-9*cos(theta(t))*diff(psi(t), t)*sin(theta(t))*diff(theta(t), t)+9/2*cos(theta(t))^2*diff(psi(t), `$`(t, 2))+1/2*cos(theta(t))*diff(theta(t), t)*diff(phi(t), t)+1/2*sin(theta(t))*diff(phi(t),...ecua := [-9*cos(theta(t))*diff(psi(t), t)*sin(theta(t))*diff(theta(t), t)+9/2*cos(theta(t))^2*diff(psi(t), `$`(t, 2))+1/2*cos(theta(t))*diff(theta(t), t)*diff(phi(t), t)+1/2*sin(theta(t))*diff(phi(t),...

Realizamos la integración numérica de nuevo.

> res:=fintr([0,0.5,0,0,0,5]);

res := proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 14; if ...

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

[Plot]

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

[Plot]

Animación del movimiento.

> dibu3(12,100);

[Plot]

>