El Oscilador Amortiguado en Matlab y Simulink

Maestría en Ciencias Físicas Fundamentos Teóricos de la Computación

[31/01/2018]

EL OSCILADOR AMORTIGUADO EN MATLAB Y SIMULINK

MATOS FRANCO, Richard 1; PEREYRA DIAZ, Jose2 . Universidad Nacional de Trujillo, Universidad Cesar Vallejo 1; Universidad Nacional de Trujillo-Perú 2

Resumen: El presente trabajo, busca resolver el problema del oscilador armónico de tres maneras diferentes en el MATLAB. La primera forma es conociendo la solución de la ecuación diferencial para el oscilador amortiguado y tabulándola. La segunda forma fue desarrollada usando el SIMULINK que que es una herramienta incorporada dentro del propio MATLAB, y en la tercera forma, empleamos el método de Runge-Kutta4, conocido como RK4 . Del análisis se puede decir que los tres métodos son soluciones a nuestro problema de estudio de manera aproximada.

I.

INTRODUCCION A. Péndulo Simple Consideremos un móvil sobre una pista sin s in rozamiento enganchado a un muelle fijo tal como muestra la figura 01. La fuerza aplicada sobre el móvil viene dada por la ley de Hooke F x x(x) ( x) = -kx . Entonces la ecuación de movimiento es o

⃛=  =  =  

⃛=  =

0101

Donde hemos introducido la contante

La cual como veremos es la frecuencia angular con la que el móvil oscilará. Aunque la ecuación (01 (01)) está en el contexto de un móvil que está sujeto a un muelle y que se mueve a lo largo del eje x, es aplicable a muchas situaciones oscilantes diferentes en muchos sistemas de coordenadas diferentes. En el caso de un monopatín que oscila en la parte interna de un cascaron esférico inferior, el ángulo con relación al centro del casquete, está gobernado por la misma ecuación, . La ecuación (01 (01)) es llamada la ecuación del oscilador armónico simple.

  ̈ = 

Figura 01: Un 01: Un móvil de masa m que oscila en el extremo de un muelle . Las soluciones exponenciales La ecuación (01 (01)) es una ecuación diferencial de segundo orden, lineal y homogénea, y por ello tiene dos soluciones independientes. Estas soluciones se pueden escoger de varias maneras, pero quizás la más conveniente sea esta:

 =    = − 1


[31/01/2018]

Como se puede comprobar fácilmente, estas dos funciones satisfacen ((01 01). ). Mas aun cualquier múltiplo constante de cualquiera de las soluciones es también una solución y también lo es cualquier suma de tales múltiplos. Por tanto, la función

 =   −

0202

Es también una solución para cualesquiera dos constantes C 1 y C2. (Que cualquier combinación lineal de soluciones como estas sea solución se conoce como principio de superposición y y juega un papel crucial en muchas ramas de la física). Cualquier solución se puede expresar de la forma (02 (02)) mediante una elección adecuada de los coeficientes C1 y C2. Las soluciones seno y coseno Las funciones exponenciales (02 (02)) son tan sencillas de manejar que (02 (02)) muchas veces es la mejor forma de la solución. Sin embargo, esta forma si tiene una desventaja. Sabemos, por supuesto, que x(t) es real, mientras que las dos exponenciales en (02 ( 02)) son complejas. Esto significa que los coeficientes C1 y C2 deben escogerse con cuidado para asegurar que el propio x(t) sea x(t) sea real. De la fórmula de Euler , sabemos que las dos exponenciales se pueden escribir como:

 =  =  =    cos cos      =  cos cos    =      =  

Sustituyendo en (02 (02)) y reagrupando encontramos que:

0303 0404

donde B1 y B2 son simplemente nuevos nombres para los coeficientes de la ecuación anterior,

La forma (03 (03)) puede tomarse como la definición del movimiento armónico simple (o simple (o MAS): MAS): cualquie8tr movimiento que sea una combinación de un seno y un coseno de esta forma se llama armónico simple. Como las funciones cos(ωt) y sen(ωt) son reales, el requerimiento de que x(t) sea x(t) sea real significa simplemente que los coeficientes B1 y B2 deben ser reales. Podemos identificar fácilmente los coeficientes B1 y B2 en términos de las condiciones iniciales del problema. Claramente t=0 Claramente t=0 , en (03 (03)) implica que x(0) = B 1. Esto es, B1 es justo la posición inicial x(0) = x o. Análogamente, diferenciando (03 ( 03), ), podemos identificar que ωB2 como la velocidad inicial v o. Si comenzamos la oscilación tirando del móvil hacia x = xo y dejándolo libre desde el reposo (v (v o = 0 ), ), entonces B2 = 0 en (03 (03)) y solo el termino del coseno perdura, así que:

 =  cos cos  =  sen sen

0505

Si lanzo el móvil desde el origen (x (x o = 0 ) dándole un toque en t = 0 , solo perdura el término del coseno, y

0505

Estos dos casos simples están ilustrados en la figura 02. Advertimos que ambas soluciones, como la solución general (03 (03), ), son periódicas porque tanto el seno como el coseno lo son. 2


[31/01/2018]

Como el argumento de ambos, seno y coseno, es ωt , la función x(t) se repite después del tiempo τ para el cual ωτ = 2 π. Esto es, el periodo es

= 2 =2 

06

Figura 02: (a) Cuando un móvil se libera desde x o en t = 0 las oscilaciones siguen una curva coseno. (b) Si el móvil es empujado desde el origen en t = 0, las oscilaciones siguen una curva con pendiente inicial v o. En cualquiera de los dos casos el periodo de las oscilaciones es τ . B. Péndulo Amortiguado Consideremos ahora que en nuestro oscilador hay una fuerza resistiva que amortiguara las oscilaciones. Existen varias posibilidades para la fuerza resistiva. El rozamiento de deslizamiento ordinario es aproximadamente constante en magnitud, pero siempre dirigido en el sentido opuesto a la velocidad. La resistencia que opone un fluido, tal como aire o agua, depende de la velocidad de una manera complicada. Sin embargo, es mejor suponer que la fuerza resistiva es proporcional a la velocidad. Aquí supusimos que la fuerza resistiva es proporcional a v; específicamente, . Consideremos entonces un objeto en una dimensión, tal como un móvil enganchado a un muelle, que esta sometido a la fuerza de la ley de Hooke, , y a una fuerza resistiva . La fuerza neta sobre el objeto es , y la segunda ley de Newton dice

 = ̇ ̈  ̇  =0



̇ 07

Para resolver esta ecuación es conveniente dividir entre m y después introducir otras dos contantes que renombramos la constante b/m como 2 γ

2= 

08

 =  

09

Este parámetro γ, se que puede denominar constante de amortiguamiento, es simplemente una forma conveniente de caracterizar la resistencia de la fuerza de amortiguamiento. Renombramos la constante k/m como ωo2, esto es,

Desde ahora en adelante usaremos la notación ωo para denotar la frecuencia natural del sistema, la frecuencia a la cual oscilaría si no hubiera una fuerza resistiva, como la que da (09). Con estas notaciones, la ecuación (07) para el oscilador amortiguado se convierte en

̈  2̇   =0 3

10


[31/01/2018]

La ecuación (08) es otra ecuación de segundo orden, lineal, homogénea. Por lo tanto, si por cualquier motivo podemos encontrar dos soluciones independientes, x 1(t) y x 2 (t) por ejemplo, entonces cualquier solución debe tener la forma C 1x 1(t) + C 2 x2 ( t). Se puede encontrar una solución de la forma

Para la cual

= =  =  2 =0  =±√   =√ =√      =∙ ∙ =−∙  −∙ − −∙

Sustituyendo en (10) vemos que nuestra suposición (11) satisface (10) si y solo si

11

12

[una ecuación auxiliar a veces llamada ecuación auxiliar para la ecuación diferencial ( 10)]. Las soluciones de esta ecuación son, por supuesto . Por consiguiente, definimos las dos contantes

Entonces las dos funciones general es:

y

13

son dos soluciones independientes de (10) y la solución

14 15

Oscilación Subamortiguada ( γ = 0) Si no hay amortiguamiento, entonces la constante de amortiguamiento γ es cero, la raíz cuadrada en los exponentes de (15) es justo i ωo, y nuestra solución se reduce a

Si reemplazamos

−∙  − −∙   =    √  = ∙ −∙ =  , se obtiene

La solución familiar para el oscilador armónico.

Amortiguamiento débil ( γ < 0) Supongamos ahora que la constante de amortiguamiento supongamos que

<

γ

16

es pequeña. Específicamente,

17

Una condición a veces denominada subamortiguamiento. En este caso, la raíz cuadrada en os exponentes de (15) es de nuevo imaginaria, y podemos escribir

   =   = Donde

4


[31/01/2018]

 =  

18

El parámetro ω1 es una frecuencia, la cual es menor que la frecuencia natural ωo. En el caso importante de amortiguamiento muy débil , ω1 está muy cerca de γ de ωo. Con esta notación, la solución (15) se convierte en

≪ =−∙ −∙ − 

19

Esta solución es el producto de dos factores. El primero, es una exponencial decreciente, que va disminuyendo gradualmente hacia cero. El segundo factor tiene exactamente la forma (16) de oscilaciones subamortiguadas, excepto que la frecuencia natural ωo, se sustituye por la frecuencia ω1, de alguna manera mas baja. Podemos reescribir el segundo factor como:

=−cos −

20

Esta solución describe claramente el movimiento armónico simple de la frecuencia ω1 con una amplitud que decrece exponencialmente , como se muestra en la Figura 03.

Figura 03 . Las oscilaciones subamortiguadas se pueden entender como oscilaciones armónicas simples con una amplitud que decrece exponencialmente . Las curvas a trazos son las envolventes, .

−

±−

Amortiguamiento Fuerte Supongamos que la constante de amortiguamiento es grande. Específicamente supongamos que

>

21

Una condición que a veces denominamos sobreamortiguamiento. En este caso, la raíz cuadrada en los exponentes de (15) es real y nuestra solución es

−∙ −+ −∙ −−  = 

22

Aquí tenemos dos funciones exponenciales reales, las cuales disminuyen a medida que pasa el tiempo (porque los coeficientes de t en ambos exponentes son negativas). En este caso, el movimiento esta tan amortiguado que no completa oscilaciones auténticas. La figura 04 muestra un caso típico en el que el oscilador es empujado desde O en t = 0, se mueve hasta un desplazamiento máximo y luego se mueve otra vez hacia atrás aún más despacio; solo vuelve al origen en el limite . El primer término del segundo miembro de (22) disminuye tan lentamente que el segundo, porque el coeficiente en su ex ponente es el menor

→∞

5


[31/01/2018]

de los dos. Por consiguiente, el movimiento a largo plazo está dominado por este primer término.

Figura 04: El movimiento sobreamortiguado en el que el oscilador es empujado desde el origen en el instante t = 0. Se mueve hasta un desplazamiento máximo y después vuelve atrás hacia O asintóticamente cuando .

→∞

Amortiguamiento Critico ( γ = ω ) o La frontera entre subamortiguamiento y superamortiguamiento se denomina amortiguamiento crítico y ocurre cuando la constante de amortiguamiento es igual a la frecuencia natural, γ = ωo. Cuando γ = ωo las dos soluciones que encontramos en (15) son la misma solución, a saber

=−

23

=−

24

=− −

25

[esto sucedió porque las dos soluciones de la ecuación auxiliar (12) coinciden cuando γ = ωo]. Este es un caso en el que nuestra intuición, al buscar una solución de la forma x(t) = e rt , falla a la hora de encontrar dos soluciones de la ecuación de movimiento, y tenemos que usar algunos otros métodos para encontrar una segunda solución. Afortunadamente, en este caso, nos es difícil hallar una segunda solución: cómo se puede comprobar fácilmente, la solución Es también una solución de la ecuación de movimiento (10) en el caso especial en que γ = ωo. Por lo tanto, la solución general para el caso de amortiguamiento critico es

Observamos que ambos términos contienen el mismo factor exponencial

−

.

C. Simulink Simulink proporciona un entorno grafico al usuario que facilita enormemente el análisis, diseño y simulación de sistemas (de control, electrónicos, dinámicos, etc.), al incluir una serie de rutinas que resuelven los cálculos matemáticos de fondo, junto con una sencilla interfaz para su uso. Proporciona un entorno de usuario grafico que permite dibujar los sistemas como diagramas de bloques tal y como se haría sobre un papel. Es un entorno de programación de más alto nivel de abstracción que el lenguaje interpretado Matlab (archivos con extensión .m). Simulink genera archivos con extensión .mdl (de model"). Simulink viene a ser una herramienta de simulación de modelos o sistemas, con cierto grado de abstracción de los fenómenos físicos involucrados en los mismos. Se hace hincapié en el análisis de sucesos, a través de la concepción de sistemas (cajas negras que realizan alguna operación). El conjunto de componentes incluidos junto al programa Simulink, incluye bibli otecas de fuentes de señal, dispositivos de presentación de datos, sistemas lineales y no lineales, conectores y

6


[31/01/2018]

funciones matemáticas. En caso de que sea necesario, se pueden crear nuevos bloques a medida por el usuario.

Figura 05: Ventana principal de trabajo de Simulink y la ventana de la Librería Browser, donde se encuentran los bloques para armar dibujar los circuitos que permiten simular los sistemas de estudio.

Figura 06: Circuito elaborado en Simulink formado por un generador de Señal (Signal Generator), un sumador (Add), integrador (Integrator), una ganancia (Gain con valor de -2 y un Osciloscopio (Scope). A la derecha se observa la grafica obtenida en el Oscilador al aplicarle una señal al circuito.

7


[31/01/2018]

D. Runge-Kutta Probablemente uno de los procedimientos numéricos mas populares, así como mas preciso, usado para obtener soluciones aproximadas para un problema con valores iniciales y’ = f(x, y), y(x o) = y o es el método de Runge-Kutta de cuarto orden. Como el nombre lo indica existen métodos de Runge Kutta de diferentes órdenes. Método de Runge Kutta En esencia, los métodos de Runge-Kutta son generalizaciones de la forma básica de Euler, en que la función pendiente se reemplaza por un promedio ponderado de pendientes en el intervalo . Es decir,

 ≤≤+ + = ℎ⏟    ⋯

26

Aquí los pesos w i , i = 1, 2, …, m , son constantes que generalmente satisfacen w 1+w 2 +…+w m = 1, y cada k i , i = 1, 2, …, m es la función f evaluada en un punto seleccionado (x,y) para el que . Veremos que las k i se definen recursivamente. El numero m se llama el orden del método. Observe que al tomar m=1, w 1 = 1 y k 1 = f(x n , y n ), se obtiene la conocida formula de Euler y n+1 = y n + hf(x n, y n) . Por esta razón se dice que el método de Euler es un método de Runge-Kutta de primer orden.

 ≤≤+

Método de Runge-Kutta de Cuarto Orden Un procedimiento de Runge-Kutta de cuarto orden consiste en determinar parámetros de modo que la formula Donde

+ = ℎ     =,  = ℎ, ℎ  = ℎ, ℎ ℎ  = ℎ, ℎ ℎ ℎ

27

Concuerda con el polinomio de Taylor de cuatro. Esto da como resultado un sistema de 11 ecuaciones con 13 incógnitas. El conjunto de valores usado con más frecuencias para los parámetros produce el siguiente resultado.

+ =  ℎ6  2 2   =,  =(  12 ℎ,  12 ℎ)  =(  12 ℎ,  12 ℎ)  =( ℎ,   12 ℎ)

28

Mientras que las otras formular de cuarto orden se deducen con facilidad, el algoritmo resumido en (28) que es muy usado y reconocido como una invaluable herramienta de 8


[31/01/2018]

cálculo, se denomina el método de Runge-Kutta de cuarto orden o método clásico de RungeKutta. De aquí en adelante, se debe considerar a (28), cuando se use la abreviatura método RK4. Método de Runge-Kutta de Cuarto Orden en Ecuaciones de Orden Superior Un problema con valores iniciales de segundo orden

"=,,′; =; ′=

29

Se puede expresar como un problema con valores iniciales para un sistema de ecuaciones diferenciales de primer orden. Si y’ = u, la ecuación dif erencial en (29) se convierte en el sistema

 =   =,,

30

Puesto que y’(x o) = u(x o), las condiciones iniciales correspondientes para (30) son El sistema (30) se puede resolver de forma numérica mediante la simple aplicación de un método numérico a cada ecuación diferencial de primer orden en el sistema. Así el método de RungeKutta de cuarto orden o método RK4, seria:

Donde

+ =  ℎ6 ℎ 2 2  31 + =  6  2 2   =;  =,  =ℎ  12 ℎ;  =(  12 ℎ,  12 ℎ,  12 ℎ )  =ℎ  12 ℎ;  =(  12 ℎ,  12 ℎ,  12 ℎ)  =ℎ ℎ;  = ℎ,  ℎ, ℎ

En general, se puede expresar cada ecuación diferencial de n-ésimo orden y (n) = f(x, y, y’, …, y (n-1)) como un sistema de n ecuaciones diferenciales de primer orden usando sustituciones y = u 1, y’ = u 2, y’’=u 3, … , y (n-1) = u n.

II. MATERIALES Y METODOS A. Materiales. Para el presente trabajo, se necesitó los siguientes materiales. ✓ Una laptop Marca : LENOVO Modelo : Z51 Sistema Operativo : Windows 8.1 Pro, 64 bits. Procesador : Intel ® Core™ i3-5005U CPU @ 2.00 GHz Memoria Ram : 4096 MB Disco Duro : 512 GB • • • • • •

9

Maestría en Ciencias Físicas Fundamentos Teóricos de la Computación ✓

✓ ✓

El programa MATLAB R2014a ▪ Sistemas Operativos ▪ Procesador ▪ Espacio en disco ▪ RAM Libros. Hojas y lapiceros.

[31/01/2018]

: : : :

De Windows XP (x64) en Adelante. Cualquier procesador Intel 1 GB solo para MATLAB, 3-4 GB 1024 MB (se recomiendan al menos 2048 MB)

B. Métodos Para desarrollar este trabajo, seguimos los siguientes pasos para cada modo de solución: Conociendo la Solución: El caso más común para oscilaciones con amortiguamiento débil, corresponde a la solución de mediante la ecuación

=−cos

20

Para ello, como ya conocemos la solución, simplemente tabularemos, dándole los valores para cada valor de t.

1. Inicio 2. Var real: tf, h, Yo, m, b, k. 3. Escribir Datos: (tiempo final, pasos, Amplitud Inicial, masa, coeficiente disipativo, constante de resorte) Leer (tf, h, Yo, m, b, k). 4. Calcular Gamma←b/(2*m). 5. Calcular matriz t = [0: h: tf]. 6. Calcular Matriz Exponencial. funexp←exp(-gamma*t). 7. Calcular frecuencia de amortiguamiento w. w←sqrt((k/m) - ((b*b)/(4*m*m))). 8. Calcular matriz solución Y. Y←Yo.*funexp.*cos(w*t). 9. Calcular matriz Velocidad V V←-Yo*(gamma*cos(w*t))-w*sen(w*t)).*funexp. 10. Calcular matriz Momento P←m*V. 11. Graficar t vs. Y 12. Graficar t vs. V 13. Graficar Y vs. P 14. Fin

10


[31/01/2018]

function OscilacionAmortiguada() %tf tiempo final %h tamaño de pasos %Y0 amplitud inicial %m masa del cuerpo oscilante %b constante disipativa %k constante del resorte o longitud del péndulo %phi ángulo de fase tf=input('Ingrese tiempo final en segundos: ' ); h=input('Ingrese intervalos de tiempo en segundos: ' ); Y0=input('Ingrese amplitud inicial en cm: ' ); m=input('Ingrese masa del cuerpo oscilante kg: ' ); b=input('Ingrese constante disipativa: ' ); k=input('Ingrese constante del resorte o longitud del péndulo N/m: ' ); gamma=b/(2*m); t=[0:h:tf]'; funexp=exp(-gamma*t); w=sqrt((k/m)-(b^2)/(4*m^2)); Y=Y0*funexp.*cos(w*t); %Matriz velocidad V=-Y0*(gamma*cos(w*t)+w*sin(w*t)).*funexp; %Matriz Momentos P=m*V; subplot(1,3,1) plot(t,Y) title('Grafico 01: Posición vs. Tiempo' ) xlabel('Tiempo (s)') ylabel('Posición (cm)') subplot(1,3,2) plot(t,V) title('Grafico 02: Velocidad vs. Tiempo' ) xlabel('Tiempo (s)') ylabel('Velocidad (cm/s)') subplot(1,3,3) plot(Y,P) title('Grafico 03: Posición vs. Momento' ) xlabel('Posición (cm)') ylabel('Momento (kg cm/s)' )

Código Fuente del Programa conociendo la solución para el caso de amortiguamiento débil

Figura 07: Valores solicitados por el programa. Podemos verificar que, por los valores ingresados, el comportamiento es de tipo de amortiguamiento débil.

11


[31/01/2018]

Grafico A. Las gráficas muestran el comportamiento de un oscilador amortiguado en el caso de tener un amortiguamiento débil, su Posición vs. Tiempo; Velocidad vs. Tiempo y Posición vs. Momento (diagrama de fase).

Grafico B . Para las mismas condiciones anteriores, con la diferencia en que la constante disipativa b es igual a cero (Sin efecto disipativo) lo cual seria un movimiento oscilatorio subamortiguado ( γ=0) o armónico simple .

12


[31/01/2018]

Desarrollo por SIMULINK Para desarrollar la ecuación diferencial (07) empleando SIMULINK, diseñamos un circuito considerando cada uno de los términos. 1. En la barra de Menú, seleccionamos la opción Tools, pestaña Library Browser.

Figura 08: Ventana principal de Simulink, la barra superior es a barra de Menu. 2. En la ventana de la Library Browser, verificamos que, en la parte izquierda, se encuentre en la opción Simulink ► Commonly Used Block. Una vez allí, seleccionamos la opción CONSTANT (constante) y trasladamos tres veces el bloque (son tres constantes en nuestra ecuación diferencial; masa, coeficiente de amortiguamiento y la constante de rigidez del resorte)

Figura 09: La figura muestra la ventana de la Librería Browser y su lista de bloques . Seleccionamos los bloques y los ubicamos en la ventana principal.

̈ =   ̇   

3. La ecuación (07) se puede escribir despejando

̈

07.

Por lo tanto, debemos considerar los bloques que nos ayuden a hacer la división de cada par de constantes, entonces consideramos dos bloques de multiplicación para las constantes y unimos las constantes con cada entrada de los bloques. 13


[31/01/2018]

Figura 10: La unión de los bloques simplemente se efectuó al colocar el puntero sobre la entrada de cada bloque multiplicador hasta la salida de cada bloque constante. 4. Por defecto, cada bloque toma a cada numero y los multiplica, pero podemos antes invertir un numero (en este caso el valor de la masa divide a los otros coeficientes) y luego multiplicar con el sobrante, para ello hacemos doble clic en cada bloque y luego le ponemos el signo de división “/” segundo el orden del valor donde se ingresó dicho número.

Figura 11: La ventana de configuración del bloque producto.

14


[31/01/2018]

̇ 

5. Como ambas constantes se multiplican por sus funciones ( y ) consideramos dos bloques mas para multiplicar a cada nueva constante obtenida en los dos productos anteriores

Figura 12: La unión con los nuevos bloques productos

̈

6. Una vez obtenido los nuevos productos debemos sumar y obtener la función ,

Figura 13: ubicación del bloque sumador, al igual que en el bloque productos, también se puede configurar el bloque sumador haciendo doble clic , el cambio de un signo “+” por un “ -” implic a en el recuadro de Number Imputs cambia el signo en la parte entrante.

̇

̈ , 

entonces a la salida del sumador debemos hacerlo pasar por un bloque Integrador; y será la integral de por lo tanto

7. La función es la primera integral de



debemos hacerlo pasar por otro bloque integral para obtener .

15

̇


[31/01/2018]

Figura 14: Los bloques Integrator se encuentran en la librería Simulink ►Continous. 8. Para poder rotar un bloque, le dimos click izquierdo y seleccionamos la opción Rotate & Flip y luego se puede aplicar una rotación horaria o anti horaria. Otra forma es seleccionar el bloque y presionar control R para una rotación horaria.

Figura 15: Colocamos los bloques Integradores y luego los rotamos para que se nos hiciera más fácil poder conectarlos entre sí y con los demás bloques. Cada rotación se hace a 90° de su posición original. 9. Unimos los integradores y multiplicamos a cada uno de ellos con sus contantes uniéndolos en las entradas de los bloques de multiplicación

16


[31/01/2018]

Figura 16: Unimos las partes faltantes de nuestro circuito. 10. Para poder graficar la posición vs. tiempo, debemos incluir un bloque que nos ayude a graficarlo, en este caso es SCOPE. Un SCOPE para la posición y otro para la velocidad. Un bloque SCOPE grafica la dependencia de una variable en el tiempo.

Figura 17: Los bloques Scope se encuentran en la librería de Simulink ► Sinks a estos bloques los conectamos a la salida del integrador de la velocidad y de la posición. 11. En la implementación del diagrama de fase, necesitamos graficar la posición vs. el momento, para calcular el momento, multiplicaremos la masa por la velocidad añadiendo el bloque de masa y el bloque de multiplicación

17


[31/01/2018]

Figura 18: Ubicación de los bloques constante y multiplicador para obtener el momento (masa por velocidad). 12. Luego añadimos el bloque XY GRAPH que grafica dos variables diferentes al tiempo.

Figura 19: Ubicación del bloque XYGraph, el cual unimos con la salida del momento y la posición. 13. Nuestro circuito que describe al sistema quedo del siguiente modo.

18


[31/01/2018]

Figura 20: Acabado final del sistema. 14. Para darle las condiciones iniciales hacemos doble clic en cada bloque integrador, recordando que la primera integral es la velocidad y la segunda es la posición 15. Ahora presionamos el botón verde en la parte superior para poder correr programa, y luego hicimos doble clic en cada uno de los Scope.

Figura 21: Graficas correspondientes a los valores de masa 0.5, coeficiente de disipación -0.1, constante de rigidez -0.5; y condiciones iniciales de velocidad Inicial 0 y posición inicial 3. Observamos que las graficas son las mismas que las obtenidas con la función solución para el caso de amortiguamiento débil. 19


[31/01/2018]

16. Luego consideramos el caso cuando el coeficiente de disipación es cero.

Figura 22: Graficas correspondientes a los valores de masa 0.5, coeficiente de disipación 0, constante de rigidez -0.5; y condiciones iniciales de velocidad Inicial 0 y posición inicial 3. Observamos que las gráficas son las mismas que las obtenidas con la función solución para el caso de subamortiguamiento. Solución por Método de Runge-Kutta de Cuarto Orden Para desarrollar el problema del oscilador amortiguado por el método de Runge-Kutta de cuarto orden, empleamos el siguiente algoritmo. 1. Inicio 2. Var real: tf, h, Yo, Vo, m , b, k. 3. Escribir Datos: (tiempo final, pasos, Amplitud Inicial, velocidad inicial, masa, coeficiente disipativo, constante de resorte) Leer (tf, h, Yo, Vo, m, b, k). 4. Calcular Gamma←b/(2*m). 5. Calcular matriz t = [0: h: tf]. 6. Calcular numero de elementos de la matriz t n←cantidad de elementos de t 7. Definir matriz Posicion Y←zeros(n,1) 8. Definir matriz V V←zeros(n,1) 9. Asignar primero elemento a Y Y(1,1) ←Yo 10. Asignar primer elemento a V V(1,1) ←Vo 11. Desde i←1 hasta n 20


[31/01/2018]

k1←h*Y(i,1); m1←V(i,1); k2←-h*(2*g*(V(i,1)+0.5*m1)+wo*(Y(i,1)+0.5*k1)); m2←h*(V(i,1)+0.5*k1); k3←-h*(2*g*(V(i,1)+0.5*m2)+wo*(Y(i,1)+0.5*k2)); m3←h*(V(i,1)+0.5*k2); k4←-h*(2*g*(V(i,1)+m3)+wo*(Y(i,1)+k3)); m4←h*(V(i,1)+k3); k←(k1+2*k2+2*k3+k4)/6; m← (m1+2*m2+2*m3+m4)/6; Y(i+1,1) ←Y(i,1)+m; V(i+1,1) ←V(i,1)+k;

12. Fin desde 13. Calcular Matriz Momento P←masa*V 14. Graficar Posición vs. Tiempo 15. Graficar Velocidad vs. Tiempo 16. Graficar Posición vs. Momento 17. Fin function RungeKuttaOscilAmortiguadoDZ() %tf tiempo final %h tamaño de pasos %Y0 amplitud inicial %V0 velocidad inicial %m masa del cuerpo oscilante %b constante disipativa %k constante del resorte o longitud del péndulo tf=input('Ingrese tiempo final: ' ); h=input('Ingrese intervalos de tiempo: ' ); Y0=input('Ingrese amplitud inicial: ' ); V0=input('Ingrese velocidad inicial: ' ); masa=input( 'Ingrese masa del cuerpo oscilante: ' ); b=input('Ingrese constante disipativa: ' ); ko=input('Ingrese constante del resorte o longitud del péndulo: ' ); g=b/(2*masa); %factor gamma wo=ko/masa; t=[0:h:tf]'; n=length(t); Y=zeros(n,1); Y(1,1)=Y0; V=zeros(n,1); V(1,1)=V0; for i=1:n-1 k1=h*Y(i,1); m1=V(i,1); k2=-h*(2*g*(V(i,1)+0.5*m1)+wo*(Y(i,1)+0.5*k1)); m2=h*(V(i,1)+0.5*k1); k3=-h*(2*g*(V(i,1)+0.5*m2)+wo*(Y(i,1)+0.5*k2)); m3=h*(V(i,1)+0.5*k2); k4=-h*(2*g*(V(i,1)+m3)+wo*(Y(i,1)+k3)); m4=h*(V(i,1)+k3); k=(k1+2*k2+2*k3+k4)/6; m=(m1+2*m2+2*m3+m4)/6;

21


[31/01/2018]

Y(i+1,1)=Y(i,1)+m; V(i+1,1)=V(i,1)+k; end P=masa*V;%Matriz de Momentos masa por velocidad subplot(1,3,1) plot(t,V(:,1)) title('Tiempo vs Velocidad' ) xlabel('Tiempo (s)') ylabel('Velocidad (cm/s)') subplot(1,3,2) plot(t,Y(:,1)) title('Tiempo vs Posición' ) xlabel('Tiempo (s)') ylabel('Posición (cm)') subplot(1,3,3) plot(Y(:,1), P) title('Tiempo vs Posición' ) xlabel('Posición (cm)') ylabel('Momento (cm/s)')

Figura 23: Valores solicitados por el programa. Podemos verificar que, por los valores ingresados, el comportamiento es de tipo de amortiguamiento débil.

Grafico C. Las gráficas muestran el comportamiento de un oscilador amortiguado en el caso de tener un amortiguamiento débil, su Posición vs. Tiempo; Velocidad vs. Tiempo y Posición vs. Momento (diagrama de fase) empleando el método de Runge-Kutta. 22


[31/01/2018]

Grafico D . Para las mismas condiciones anteriores, con la diferencia en que la constante disipativa b es igual a cero (Sin efecto disipativo) lo cual sería un movimiento oscilatorio subamortiguado ( γ=0) o armónico simple .

III. DISCUSIONES Consideramos los casos amortiguamiento débil, porque fijamos la solución en nuestro programa de referencia, el cual tiene como caso limite cuando b=0 (γ = 0) la oscilación subamortiguada. Comparamos el Gráficos C y las figura 21 con la figura A, para el caso de amortiguamiento débil (esto es para, γ < 0), nos dimos cuenta por inspección rápida que el grafico C y la figura 21 difieren del grafico A, porque las crestas son de mayor tamaño que las crestas de la solución esto es debido a que los otros dos gráficos son obtenidos por métodos numéricos, los cuales llevan son soluciones aproximadas y van sumando los errores debido a los procesos recursivos como es en el bucle for para el caso del Método Ringe-Kutta de cuarto orden. Sin embargo, cuando comparamos el Gráficos D y la figura 22 con la figura B, en el caso de subamortiguamiento (γ = 0) observamos que la figura 22 (empleando el Simulink) es prácticamente idéntica a la grafica A, pero el grafico D, muestra un comportamiento inusual donde la amplitud en vez de mantenerse constante, se incrementa en cada periodo de oscilación, este comportamiento es el típico caso cuando existe una fuerza adicional llamada fuerza de restitución. Estos resultados se deben a que a medida que se va haciendo mas interacciones, los errores van incrementando porque ya no tiene un factor adicional que le ayude a disminuir como era el coeficiente de restitución, es decir, que ahora su incremento sería mucho más alto. En el simulink, empleamos el método Runge-Kutta45 (por defecto) a fin de poder comparar con la gráfica de la función solución.

23


[31/01/2018]

IV. CONCLUSIONES A pesar que el Simulink es una herramienta incluida en el Matlab, nos da una mayor aproximación a la función solución que el Método de RK4 elaborado para este caso, donde se aprecia la falla de dicho método en el caso de la oscilación subamortiguada. Mientras mas iteraciones se haga (caso del bucle for, if , while) mayor será el incremento de los errores y por lo tanto nuestros resultados difieren del valor esperado. Se sugiere para otros estudios, considerar los otros casos en cada algoritmo que nos permita comprar dichos métodos a fin de poder verificar, así como también diferentes valores de parámetros incluidos los tamaños de pasos.

V. REFERENCIAS [01] Burden, R. Douglas, J.; (2002), Análisis Numérico. Editorial Thomson; Séptima edición; Ciudad de México, México. [02] Gekeler, A.; (2008). Mathematical Methods for Mechanics A Handbook with MATLAB Experiments. Editorial Springer; Primera edición; Heidelberg, Alemania. [03] Gilat, A.; (2005); Matlab: Una introducción con ejemplos prácticos; Editorial Reverte; Segunda edición, Barcelona, España. [04] Joyanes, L.; (2008); Fundamentos de Programación: Algoritmos, estructura de datos y objetos; Editorial Mc Graw Hill, Cuarta edición; Madrid, España. [05] Raffo, E.; (2011); Métodos Numéricos para ciencias e ingeniería con Matlab. Editorial; Editorial Raffo; Segunda edición; Lima, Perú. [06] Ramírez J.; Fecha de consulta: 29/01/2018. Recuperado de: http://www.ugr.es/~javierrp/master_files/Seminario%20de%20Matlab.pdf . [07] Universidad de Granada; Fecha de consulta: 28/01/2018. Recuperado de: http://www.esi2.us.es/~jaar/Datos/RegAuto/Practica3.pdf . [08] Taylor, J.; (2013); Mecánica Clásica; Editorial Reverte; Primera edición; Barcelona, España. [09] Vásquez, L. Jiménez S. Aguirre, C. Pascual, J.; (2009); Métodos Numéricos para la Física y la Ingeniería, Editorial Mc Graw Hill; Primera edición; Doosan, Corea del Sur. [10] Zill D; (2009); Ecuaciones Diferenciales con aplicaciones de modelado; Editorial Cengage Learning; Novena edición; Ciudad de México, México.

24

El Oscilador Amortiguado en Matlab y Simulink

Recommend Documents