Cómo hacer simulaciones de partículas en MatLab. IV: Péndulo (1 part. – 2D)

< III. Órbita gravitacional (1 part. – 2D)

 

En los dos tutoriales anteriores habíamos visto simulaciones de una sola partícula, a las que se les aplicaba una fuerza que era constante (tiro parabólico) o que dependía de la posición de la partícula (órbita). En esta ocasión vamos a ver un caso similar a los anteriores (también una sola partícula en 2D), pero ahora la fuerza sobre la partícula depende tanto de su posición como de su velocidad.

 

  F = const. F (x) F (v) F (x,v)

1D

2D

II (1)

III (1)

IV (1)

3D

 

*(1): Una sola partícula. (N): Un sistema de N partículas independientes. (C): Colisiones.

 

Mecánica del péndulo

El sistema que he escogido para este tipo de simulación es un péndulo simple contenido en un plano, y que puede presentar rozamiento con el aire. Un péndulo simple tiene dos partes: una masa que ocupa un punto (nuestra partícula a simular), y una cuerda que se considera siempre tensa (forma una línea recta), no elástica (siempre tiene la misma longitud) y sin masa.

En el siguiente dibujo pueden verse las fuerzas que actúan sobre la partícula, que son tres: el peso (P), la fuerza que hace la cuerda (T) y la fuerza de fricción (Fr). La posición viene determinada por el ángulo θ. La velocidad y aceleración angular son ω y α respectivamente. La longitud de la cuerda es L

 

Pendulo_dinamica

 

También aparecen las aceleraciones en azul. La partícula sigue un movimiento circular de radio L. Eso significa que su aceleración tendrá dos componentes: una tangencial (perpendicular a la cuerda), que se puede calcular a partir de la aceleración angular como α·L; y una normal (en la dirección de la cuerda) que se calcula en función de la velocidad angular como ω2 ·L. Así, el vector aceleración de la partícula es:

image

 

El peso (P) es una fuerza constante, porque así consideramos la masa (m) y el valor del campo gravitatorio (g), y siempre actúa en la misma dirección y sentido. Se calcula como:

P = m·g

 

La fuerza de fricción (Fr) es la responsable de la novedad de esta simulación con respecto a las anteriores. Es una fuerza que depende de la velocidad de la partícula. La mecánica de fluidos es una ciencia complicada, y la realidad no es tan simple, pero supondremos que la fuerza de fricción viscosa es proporcional a la velocidad, y por supuesto siempre tiene el sentido contrario de esta. He llamado C a la constante de rozamiento viscoso. La velocidad lineal (v) se puede escribir en función de la velocidad angular (ω) conociendo que v = ω·L:

Fr = –C·v = –C·ω·L

 

La fuerza que hace la cuerda (T) tiene siempre la dirección de la cuerda y el sentido que se muestra en el dibujo. Su valor depende tanto del ángulo como de la velocidad angular, y lo calcularemos en cuanto apliquemos las leyes de Newton a continuación. No obstante el cálculo de esta fuerza no es necesario para la simulación.

Si sumamos todas las fuerzas, la resultante (F) es:

image

 

Aplicando la segunda ley de Newton, F = m·a, tenemos la siguiente ecuación:

image

 

De donde obtenemos dos ecuaciones. Una nos servirá para saber la aceleración angular α conociendo la posición y la velocidad. La otra nos permitiría calcular la fuerza T, si nos interesara.

image

 

Nótese que en el caso de no tener rozamiento (C=0), la dinámica de la partícula es independiente de la masa.

 

Simulaciones

Ahora que sabemos cuál es la aceleración angular en cada posición y para cada velocidad, podemos hacer las simulaciones de forma similar a los tutoriales anteriores. El siguiente es un ejemplo de código que simula el movimiento de un péndulo simple dadas las condiciones iniciales del ángulo θ y de velocidad angular ω. También guarda (en theta_graf) y representa al final el recorrido del ángulo θ a lo largo de la simulación.

 

clear
figure(1)

theta_graf=[];

%variables
theta=[];
w=[];

%parametros
m = 1;
g = 9.81;
L = 1; %longitud del péndulo
C = 0.75; %constante de fricción viscosa
h=0.01;

%condiciones iniciales
theta = 30*(pi/180); %entre –pi/2 y pi/2 en radianes 
x=[L*sin(theta);-L*cos(theta)];
w = 0/L;

%entrada
alpha = -(L*w*C + m*g*sin(theta))/(L*m);

for step=1:1000
   
    hold off
    plot(x(1),x(2),‘o’, ‘MarkerFaceColor’,‘b’,‘MarkerSize’,10);
    hold on
    plot([0;x(1)],[0;x(2)]);
   
    title([‘step : ‘ num2str(step)]);
    axis([-L L –L 0]);
    set(gca,‘dataAspectRatio’,[1 1 1])
    pause(0.001);
   
    theta_a = theta;
    wa = w;
   
    %paso integración
    wpm = wa + (h/2)*alpha;
    theta_pm = theta_a + (h/2)*wa;
    alpha_pm = -(L*wpm*C + m*g*sin(theta_pm))/(L*m);
   
    w = wa + h*alpha_pm;
    theta = theta_a + h*wpm;
    x=[L*sin(theta);-L*cos(theta)];
    alpha = -(L*w*C + m*g*sin(theta))/(L*m);
   
    theta_graf = [theta_graf theta];
       
end

figure(2)
plot(theta_graf,‘b’);

 

Paso a explicar algunas novedades o cosas que quizás no quedan claras solo con ver el código:

Variables: Esta vez no uso la posición y velocidad lineales como en los anteriores tutoriales. En un péndulo, basta con saber el ángulo y la velocidad angular para determinar la posición y velocidad. Las ecuaciones de esta forma son más sencillas.

Parámetros: La única novedad es el coeficiente de rozamiento viscoso C. Su valor siempre debe ser positivo o cero. Si es cero no hay rozamiento. Si es menor que un cierto valor crítico la oscilación del péndulo será sub-amortiguada: oscilará varias veces reduciendo la amplitud. Si el coeficiente es superior al valor crítico la oscilación será sobre-amortiguada: no llegará a hacer una oscilación completa y tenderá a la posición de ángulo cero. 

Puede calcularse el valor crítico de C suponiendo que el sistema es lineal (haciendo la aproximación sin(θ) ≈ θ). No quiero desviarme mucho justificándolo, pero se puede comprobar con la simulación que el valor crítico para el coeficiente C es:

image

 

donde ωn es la frecuencia natural (del sistema lineal aproximado). Para los parámetros de la simulación, el valor crítico es de 6.26 aproximadamente. Recomiendo al lector comprobar diferentes valores por encima y por debajo del crítico para ver las distintas clases de amortiguamiento.

Condiciones iniciales: El ángulo theta debe estar en radianes (rad), pero he incluido el cálculo para pasar los grados a radianes. En el código de ejemplo el ángulo inicial es de 30º. Solo se tiene que cambiar el 30 por cualquier otro ángulo en grados, pero debe estar comprendido entre -180º y 180º.

La velocidad inicial ω es angular (rad/s), pero como con el ángulo theta, puede darse en m/s, substituyendo el valor de la fórmula (el numerador) por el que se desee. En el ejemplo el valor es de 0 m/s. Basta con cambiar el 0 por cualquier otro valor de velocidad lineal.   

Cálculo de alpha: Cada vez que conozcamos el ángulo y la velocidad angular, podremos calcular la aceleración angular, que nos servirá para predecir el siguiente ángulo y velocidad. La fórmula para el cálculo es la que hemos deducido en la parte teórica.

Dibujo: Lo primero que se hace en el bucle es dibujar la partícula y la cuerda, según el ángulo theta actual. hold off sirve para eliminar el dibujo anterior y hold on para dibujar tanto la cuerda como la partícula. El primer plot dibuja la partícula, indicando la forma, el tamaño y el color de relleno.

plot(x(1),x(2),‘o’, ‘MarkerFaceColor’,‘b’,‘MarkerSize’,10);

 

Utilizo plot en lugar de scatter porque he descubierto que lo hace más rápido, así que lo recomiendo. El segundo plot es el de la cuerda, simplemente una línea entre dos puntos. Podría cambiarse el color, grosor, o lo que sea si se desea. Se dibuja el número de paso y se establecen los ejes. Una novedad es que utilizo la siguiente función para mantener la escala en los ejes:

set(gca,‘dataAspectRatio’,[1 1 1])

 

Paso de integración: Utilizo el método del punto medio como en los tutoriales anteriores. La única novedad es que ahora necesito también la velocidad del punto medio para calcular la aceleración.

 

El resultado del programa es una simulación como la que se ve aquí abajo. Es un ejemplo de oscilación sub-amortiguada, porque la constante C está por debajo de la crítica. 

 

Tutorial03_roz

 

Y este es el resultado del plot al final del programa, el que representa el ángulo theta a lo largo de la simulación.

 

Amort_1

 

Pueden cambiarse los parámetros de masa, ángulo inicial, velocidad inicial, longitud de la cuerda, o constante de fricción para comprobar cómo afectan a la dinámica de un péndulo. Por ejemplo puede comprobarse que para ángulos pequeños el periodo de oscilación es independiente de la masa y del ángulo inicial, pero sí depende de la longitud de la cuerda. Esto es algo muy poco intuitivo, pero ya lo demostró Galileo. La frecuencia natural ωn que mostré antes solo depende de la gravedad y de la longitud de la cuerda.

Para ángulos pequeños (es cuando nuestro sistema se aproxima a uno lineal), el periodo de oscilación también es independiente de la constante de fricción. He alterado un poco el programa para simular diferentes valores de C, desde cero (sin amortiguación) hasta 1. Como el valor crítico está por encima de 1, todas las oscilaciones son sub-amortiguadas (menos la de C=0, claro). Este es el resultado:

 

Amort

 

Para el que quiera un reto

Mientras hacía este tutorial me vino a la mente un video en el que se hacían oscilar varios péndulos de diferentes longitudes pero con el mismo ángulo inicial. El resultado era un fantástico baile sincronizado que formaba figuras geométricas hasta que volvían a coincidir todos en la misma fase inicial. Este es el video:

 

 

Me pregunté si podría hacer una simulación que imitara esa sincronía con varios péndulos. Debo admitir que no fue una tarea fácil, no tanto por escoger los periodos de cada péndulo como por hacer que se viera bien. Intentaré explicar un poco las matemáticas detrás de ese efecto tan curioso.

Lo primero que supuse fue que el sistema era lineal (ángulos pequeños), y por tanto, como ya he ido diciendo, el periodo solo depende de la longitud de la cuerda. Así que toda la dificultad reside en encontrar las longitudes adecuadas. Usé una amplitud de 10 grados para las oscilaciones. También consideré que no había fricción (C=0), y que la velocidad inicial es cero. 

Para empezar podemos pensar qué tienen de especial los periodos de los péndulos para lograr el efecto del vídeo. Lo que debe pasar es que, comenzando todos con el mismo ángulo (misma amplitud), al cabo de un cierto número de oscilaciones (diferente para cada péndulo), deben volver a coincidir todos en el mismo ángulo inicial. Pongamos el ejemplo de dos péndulos solamente. Llamemos T1 al periodo del péndulo 1 y T2 al periodo del péndulo 2. El péndulo 1 tarda T1 segundos en dar una oscilación completa, 2* T1 en dar 2 oscilaciones, 3* T1 en dar 3, etc. Y lo mismo para el péndulo 2 con su periodo. Si vuelven a coincidir en la posición inicial será porque los dos péndulos han tardado lo mismo en hacer un cierto número (diferente para cada uno) de oscilaciones completas. Por tanto, siendo n1 y n2 números naturales que representan las oscilaciones completas de cada péndulo, entonces debe pasar:

n1·T1 = n2·T2

 

Si los dos periodos cumplen lo de arriba (la relación entre ellos es una fracción de naturales), entonces seguro que volverán a coincidir donde empezaron. Nótese también que el momento en que coincidan será el mínimo común múltiplo de sus periodos. Si por ejemplo el péndulo 1 tiene periodo T1 = 3 segundos, y el péndulo 2 tiene periodo T2 = 2 segundos, se volverán a encontrar al cabo de 6 segundos (el mínimo común múltiplo de 2 y 3). Los periodos no tienen por qué ser naturales. Pueden ser incluso irracionales, siempre y cuando se cumpla la condición de arriba. 

Nosotros no elegimos directamente el periodo del péndulo, si no que lo hacemos a través de la longitud de la cuerda. Como ya he ido diciendo, si el sistema fuese lineal la frecuencia de su oscilación (no confundir esta frecuencia con velocidad angular del péndulo, aunque use la misma letra) solo depende de L y de g. Podemos expresar la frecuencia en función del periodo:

image

 

La longitud depende de una constante y del cuadrado del periodo. Ahora veamos la relación que debe haber entre dos longitudes para que sus periodos cumplan la condición expresada antes. Al dividir dos longitudes usando el resultado anterior, las constantes se cancelan y el resultado es:

image

 

Esta es la condición que deben cumplir las longitudes de la cuerdas. Si las longitudes cumplen esta condición, los periodos cumplirán la de n1·T1 = n2·T2. Esto sirve para un número cualquiera de péndulos. Si la pareja de péndulos 1 y 2 cumplen la condición, y por otro lado los péndulos 1 y 3 también la cumplen, entonces la pareja 2 y 3 también la cumple. Por eso podemos preocuparnos solo de que esa condición se cumpla entre el péndulo 1 y los demás (1-2, 1-3, 1-4,…), y automáticamente los demás la cumplirán entre si.

Yo lo que hice fue establecer una “longitud de referencia” de 1, y hacer que el resto de longitudes cumplieran la condición de arriba con respecto a la de referencia. Si en la ecuación de arriba se substituye L2 = 1, tenemos una forma de calcular el resto de longitudes simplemente inventando valores de n1 y n2. Pero, ¿qué valores son mejores?

Responder a esa pregunta ya no es cuestión de que se cumpla la condición o no, si no de que se vea bien la simulación. Si por ejemplo unas longitudes son muy largas y otras muy cortas, los péndulos cortos no se verán bien. O si el mínimo común múltiplo de todos los periodos es demasiado grande, tardaríamos mucho en ver cómo vuelven a coincidir los péndulos. De hecho lo más difícil para mí no fue descubrir la condición que debían cumplir las longitudes, sino encontrar el conjunto de longitudes que hicieran que la simulación fuese visualmente buena.

Al final conseguí un conjunto bastante aceptable de valores para n2 / n1 . Si N es el número de péndulos, lo que hice fue asignar a la fracción n2 / n1  una fracción entera de N, con el denominador sucesivamente más grande, empezando por N. Por ejemplo si hubiera 5 péndulos (N=5), el conjunto de valores n2 / n1  sería:

n2 / n1  =  { 5/5 , 5/6 , 5/7 , 5/8, 5/9 }

 

Recordemos que n2 y n1 nos decían el número de oscilaciones completas que debe dar cada péndulo para encontrarse de nuevo. Esta forma de escoger los valores (siempre el mismo n2) impone que el número de oscilaciones del péndulo de referencia (el de longitud 1) sea N cuando se vuelven a encontrar todos. Cuando el primer péndulo hace 5 oscilaciones completas, el segundo habrá hecho 6, el tercero 7, el cuarto 8, y el quinto 9.

Para encontrar el conjunto de longitudes solo hay que usar la formula que relaciona longitud con la fracción n2 / n1 . En nuestro ejemplo serían:

Li = { 1 , 25/36 , 25/49 , 25/64 , 25/81 }

 

Esta es la simulación con mi forma particular de repartir las longitudes que cumplen con la condición de sincronía.

 

Tutorial03

 

Seguro que existen otras formas de encontrar conjuntos de fracciones n2 / n1, no hay nada que nos obligue a usar este en particular. Me encantaría si algún lector encontrara otras formas de calcular estos conjuntos de longitudes. Nótese que no están distribuidas nunca de forma lineal. Esto se puede ver también en el vídeo, al inicio los péndulos en reposo no forman una línea recta, sino una ligera curva (me aventuraría a decir que del tipo 1/x2). Puede que otro tipo de curvas en la distribución de las longitudes den como resultado efectos distintos. Sería interesante comprobarlo.

 


Esta entrada pertenece a la colección Animación con MatLab

6 comentarios

  1. Información Bitacoras.com

    Valora en Bitacoras.com: No hay resumen disponible para esta anotación

  2. Me parece muy buena la explicación y el reto que te pusiste, pero a mi me pusieron un pendulo que es un prisma rectangular y me piden que el centro de masa del prisma influya en el movimiento del mismo, sinceramente no se como hacerlo. Si pudieras ayudarme te lo agradeceria.

    1. Siempre que la gravedad sea uniforme el centro de masas coincide con el de gravedad. Y supongo que la densidad será uniforme así que el centro de masas coincidirá también con el centroide (centro geométrico del prisma).
      Si pasa eso (el centro de gravedad coincide con el centroide), entonces el movimiento del péndulo es el mismo que sigue una partícula que este situada en el centro de gravedad, con toda la masa del prisma en ese punto.

      Si el centro de masa cambia de lugar para el mismo prisma (mismo tamaño me refiero), y la gravedad es uniforme, entonces es que el prisma no tiene densidad uniforme. Pero igualmente se puede argumentar que el centro de masas es el de gravedad, e influirá en el movimiento en el sentido de que será como tener un péndulo puntual con diferente longitud.

      En resumen: el movimiento es el del centro de gravedad con la masa total concentrada ahí, y ese es un problema de una partícula. Y el centro de gravedad es el masas si la gravedad es uniforme (que supongo que sí).

    2. El CG de ese prisma describirá el mismo movimiento que una partícula suspendida, pero el prisma oscilará con la posición en la que está en equilibrio cuando está suspendido de la cuerda. En esa posición de equilibrio, cada punto del sólido describirá el mismo movimiento, tratando de conservar la misma distancia relativa respecto a los demás puntos del sólido.

      Saludos

  3. […] Cómo hacer simulaciones de partículas en MatLab. IV: Péndulo (1 part. – 2D): Nueva entrega de la colección sobre simulaciones en Matlab. En esta reproduzco un video sobre […]

¿Qué te parece?

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

A %d blogueros les gusta esto: