ENUNCIADO DEL EJEMPLO 23

Un disco de radio R y masa M se mueve en todo momento en un plano vertical y rueda sin deslizar sobre una recta horizontal.

En el disco existe una ranura diametral lisa en la que se mueve una partícula de masa m, que está unida a un punto de la periferia mediante un resorte (muelle) de constante elástica k y longitud natural R.

Se pide:

            1.Definir los parámetros que definen el movimiento del sistema.

            2.Expresar las energías cinética y potencial.

            3.Expresar la ecuación de Lagrange.

            4.Ecuaciones del movimiento.

            5.Gráficas con la variación de los parámetros en funcion del tiempo.

            6.Animación del sistema.

Empezaremos cargando los paquetes necesarios para realizar el ejercicio, entre ellos el mecapac3d, para el cual tendremos que indicar la ubicación de la librería en el disco duro.

> 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):

Seguiremos ahora definiendo lo primero que se nos ha pedido en el enunciado,es decir, los parámetros que definen el sistema,que serán dos, ya que el sistema cuenta con dos grados de libertad. Un primer parámetro angular (theta) para definir la posicion del disco, y otro parámetro (s) que defina la posicion de la partícula según la elongación del muelle. Definimos por tanto esas coordenadas generalizadas.

> cg:=[theta,s];

cg := [theta, s]

A continuacion vamos a ir definiendo las componentes del sistema:

Primero empezaremos por el disco, indicando las rotaciones a las que está sometido, que serán dos, que al considerarlas absolutas (referidas a los ejes fijos), las multiplicaremos por la izquierda. Después definimos el disco, dejando como parámetros sus masa y su radio.

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

> rot2:=rota(theta,1):

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

> d1:=[disco,[0,R*theta,R],rottot,M,R]:

Habiendo definido ya el disco pasaremos a definir lo que nos queda del sistema,es decir, el muelle  y la partícula.

> m2:=[muelle,[0,R*theta+R*sin(theta),R+R*cos(theta)],[0,R*theta-s*sin(theta),R-s*cos(theta)],k,l0]:

> p3:=[punto,0,R*theta-s*sin(theta),R-s*cos(theta),m]:

Ahora para completar el sistema incluiremos una representación de los ejes fijos mediante componentes de tipo segmento y texto para nombrarlos.

> ejey:=[segmento,[0,-20,0],[0,20,0],green]:

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

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

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

Creamos ahora el sistema definitivo que englobará a todos los anteriores es decir tanto al disco, muelle y partícula como a los ejes fijos que acabamos de definir.

> sistema:=[d1,m2,p3,ejey,ejez,TY,TZ]:

A partir de este momento al tener todo definido empezaremos a ir respondiendo a las demás preguntas que formulaba el ejercicio,la segunda pregunta formulada era la expresión de las energías cinética y potencial, que serán las siguientes.

Primero el potencial:

> V:=fV(sistema);

V := M*g*R+1/2*k*(((-s*sin(theta)-R*sin(theta))^2+(-s*cos(theta)-R*cos(theta))^2)^(1/2)-l0)^2+m*g*(R-s*cos(theta))

Ahora la energía cinética

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

T := -m*R*theta1*s1*sin(theta)-m*R*theta1^2*s*cos(theta)+3/4*theta1^2*M*R^2+1/2*m*R^2*theta1^2+1/2*m*s1^2+1/2*m*s^2*theta1^2

Habiendo ya definido las energías cinética y potencial formaremos la Lagrangiana, y contestaremos a la tercera pregunta.

> L:=simplify(T-V);

L := -m*R*theta1*s1*sin(theta)-m*R*theta1^2*s*cos(theta)+3/4*theta1^2*M*R^2+1/2*m*R^2*theta1^2+1/2*m*s1^2+1/2*m*s^2*theta1^2-M*g*R-1/2*k*R^2-1/2*k*s^2-k*s*R+k*csgn(R+s)*l0*R+k*csgn(R+s)*l0*s-1/2*k*l0^...

En el cuarto apartado nos piden las ecuaciones del movimiento,que serán las siguientes.

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

ecua := [-m*R*diff(s(t), `$`(t, 2))*sin(theta(t))-2*m*R*diff(s(t), t)*cos(theta(t))*diff(theta(t), t)-2*m*R*diff(theta(t), `$`(t, 2))*s(t)*cos(theta(t))+m*R*diff(theta(t), t)^2*s(t)*sin(theta(t))+3/2*...ecua := [-m*R*diff(s(t), `$`(t, 2))*sin(theta(t))-2*m*R*diff(s(t), t)*cos(theta(t))*diff(theta(t), t)-2*m*R*diff(theta(t), `$`(t, 2))*s(t)*cos(theta(t))+m*R*diff(theta(t), t)^2*s(t)*sin(theta(t))+3/2*...

Empezaremos ahora a asignar valores numéricos a todos los parámetros definidos anteriormente,para poder obtener una representación gráfica del sistema en un instante determinado y realizar la integración numérica.

> g:=9.81:M:=20:m:=10:R:=5:l0:=10:k:=25:

> fG([evalf(Pi/3),-1]);

[Plot]

> res:=fint([(evalf(Pi/3)),.5,-1,0.1]):

En el apartado quinto nos piden la evolución de los parámetros, tanto el angular como el lineal, con respecto al tiempo.

Primero representaremos el angular.

> odeplot(res,[t,theta(t)],0..50,numpoints=50);

[Plot]

Ahora haremos la evolución del parámetro s.

> odeplot(res,[t,s(t)],0..20);

[Plot]

Y ahora la animación del movimiento, durante 10 segundos con 50 imágenes.

> dibu3(10,50);

[Plot]

>