Implementación en Matlab de métodos de 1 g.d.l. para realizar análisis modal experimental de estructuras

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental de Estructuras Autor: Sergio Alducin Vázquez Tutora: Rosario Chamorro Moreno Dep. de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 iii Trabajo Fin de Grado Grado en Tecnología de las Tecnologías Industriales Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental de Estructuras Autor: Sergio Alducin Vázquez Tutora: Rosario Chamorro Moreno Profesor contratado doctor Dep. de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 v Trabajo Fin de Grado: Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental de Estructuras Autor: Sergio Alducin Vázquez Tutora: Rosario Chamorro Moreno El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2017 El Secretario del Tribunal vii Resumen Desde mediados del siglo XX se han llevado a cabo importantes avances relativos a la comprensión del comportamiento dinámico de los sistemas estructurales. La aplicación de estos conocimientos ha sido fundamental a la hora de poder garantizar la integridad estructural y el correcto funcionamiento en condiciones de operación de cualquier estructura. Los análisis dinámicos analíticos y numéricos no suelen ser suficientes para llevar a cabo el diseño de una estructura. En el primer caso porque solo son aplicables a sistemas muy simples, y en el segundo caso porque lo que se analiza no deja de ser un modelo, que puede diferir de la realidad. Es por esto por lo que el análisis modal experimental juega un papel fundamental en el estudio de vibraciones. El objetivo del presente Trabajo Fin de Grado es realizar unas funciones en Matlab que permitan llevar a cabo el análisis modal de cualquier sistema a partir de las medidas experimentales de sus Funciones de Respuesta en Frecuencia (FRFs). De esta manera se puede realizar una descripción del comportamiento dinámico de una estructura, conociendo los parámetros modales (frecuencia natural, amortiguamiento y constante modal) de todos los modos dentro del rango de frecuencias de interés. El presente trabajo de centra en los métodos de un grado de libertad, que, a pesar de lo que su nombre parece indicar, son aplicables a sistemas de varios grados de libertad siempre y cuando se verifiquen ciertas hipótesis. En el capítulo 2 se resume la teoría necesaria para definir todos los conceptos fundamentales relativos a las vibraciones y para obtener las ecuaciones elementales que se usarán más adelante. En el capítulo 3 se realiza una clasificación de los distintos tipos de análisis modal experimental, y se desarrollan en detalle tres métodos de un grado de libertad. Con esta base teórica ya es posible realizar los archivos de Matlab para obtener las propiedades dinámicas de un sistema estructural mediante los tres métodos previamente descritos. Para ilustrar el proceso de obtención de resultados se van a aplicar estos métodos a las FRFs de un voladizo empotrado en su base. Esto permitirá comprobar las ventajas e inconvenientes de cada uno de los métodos, ya que, como se puede presuponer, ninguno es mejor que los otros para todos los posibles casos de estudio. La condición para poder aplicar los métodos de un grado de libertad se puede resumir de forma simplificada como que un modo no debe verse influido significativamente por el resto de modos, y, de forma aún más simplificada, como que los modos deben estar suficientemente alejados entre sí. Sin embargo, no existe una definición clara de cuando dos modos se pueden considerar lo suficientemente alejados y cuando no. Los ficheros de Matlab pueden aplicarse también a curvas de FRFs generadas de forma teórica, imponiendo los valores que las definen. De esta manera se puede estudiar cómo afecta la variación de los parámetros que definen los modos de la FRF de un sistema cualquiera a la diferencia entre los resultados que se obtienen respecto a los que se esperaba obtener. ix Índice Resumen vii Índice ix Índice de Tablas xi Índice de Figuras xiii Notación xv 1 Introducción 1 2 Fundamento Teórico de las Vibraciones 5 2.1. Sistemas de un grado de libertad 5 2.1.1. Sistemas no amortiguados 5 2.1.1.1. Vibración libre 5 2.1.1.2. Vibración forzada 6 2.1.2. Sistemas con amortiguamiento viscoso 6 2.1.2.1. Vibración libre 6 2.1.2.2. Vibración forzada 8 2.1.3. Amortiguamiento estructural o histerético 8 2.1.4. Amortiguamiento seco o de Coulomb 9 2.2. Funciones de respuesta en frecuencia (FRF) en sistemas de 1 g.d.l. (SDOF) 10 2.3. Sistemas de varios grados de libertad (MDOF) no amortiguados 10 2.3.1. Vibración libre 11 2.3.2. Propiedades de ortogonalidad 11 2.3.3. Vibración forzada 11 2.4. Sistemas de varios grados de libertad (MDOF) con amortiguamiento proporcional 12 2.4.1. Amortiguamiento proporcional a la rigidez 13 2.4.2. Forma general del amortiguamiento proporcional 13 2.5. Sistemas de varios grados de libertad (MDOF) con amortiguamiento histerético 14 2.5.1. Vibración libre 14 2.5.2. Vibración forzada 14 2.6. Sistemas de varios grados de libertad (MDOF) con amortiguamiento viscoso 15 2.6.1. Vibración libre 15 2.6.2. Vibración forzada 16 3 Análisis Modal Experimental 19 3.1. Clasificación de los tipos de análisis 19 3.2. Dificultades en el análisis modal experimental 20 3.3. Métodos de un grado de libertad (SDOF) 20 3.3.1. Introducción 20 3.3.2. Peak-Picking 20 3.3.3. Circle-Fit 22 3.3.3.1. Propiedades del círculo modal 22 3.3.3.2. Descripción del método 24 [𝑀] Matriz de masa [𝐾] Matriz de rigidez 𝑑𝑒𝑡| | Determinante de una matriz [ ]𝑇 Transpuesta de una matriz [ ]−1 Inversa de una matriz 𝜔𝑟 Frecuencia natural r-ésima {𝜓}𝑟 Modo de vibración r-ésimo [⋱𝜔𝑟2⋱] Matriz de frecuencias naturales (o de autovalores) [𝛹] Matriz de modos de vibración (o de autovectores) 𝑚𝑟 Masa modal del modo r-ésimo 𝑘𝑟 Rigidez modal del modo r-ésimo [𝑚𝑟] Matriz de masas modales [𝑘𝑟] Matriz de rigideces modales 𝜙𝑟 Autovector del modo r-ésimo normalizado respecto a la masa [𝛷] Matriz de modos (o de autovectores) normalizada respecto a la masa [𝛼(𝜔)] Matriz de receptancia 𝛼𝑗𝑘(𝜔) Elemento de la matriz de receptancia 𝐴𝑗𝑘𝑟 Constante modal [𝐶] Matriz de amortiguamiento viscoso 𝑐𝑟 Amortiguamiento modal del modo r-ésimo [𝑐𝑟] Matriz de amortiguamientos modales ( )∗ Conjugado de una variable compleja 𝜂𝑟 Factor de pérdida del modo r-ésimo 𝜁𝑟 Factor de amortiguamiento del modo r-ésimo { }𝐻 Transpuesta hermitiana {𝑢} Matriz de coordenadas (contiene los desplazamientos y las velocidades) 1 1 INTRODUCCIÓN n movimiento vibratorio es la variación de la configuración de un sistema respecto al tiempo alrededor de una posición de referencia. En general, cualquier sistema que tenga masa y capacidad para la deformación elástica es susceptible de vibrar libremente [1]. Las vibraciones son un fenómeno presente en prácticamente todos los ámbitos de la ingeniería. Salvo algunos casos concretos de máquinas que basan su funcionamiento en la vibración de algunos de sus componentes, las vibraciones son por lo general un fenómeno indeseado, que aparece inevitablemente durante la operación de un sistema mecánico, debido a la aplicación de una excitación en forma de colisión o fuerza variable con el tiempo. Los problemas asociados a la existencia de vibraciones son varios. Por ejemplo, las vibraciones pueden afectar negativamente al confort de las personas próximas a la máquina, ya sea por la propia vibración, o por el ruido que ésta produce. Además, las vibraciones pueden provocar el deterioro progresivo de sus componentes, ya que la fatiga es un fenómeno ampliamente estudiado y conocido, por el que un sistema estructural puede llegar a fallar debido a cargas que no superan la carga máxima admisible, pero que se repiten continuamente durante largos periodos de tiempo. En otros casos, la propia vibración puede ocasionar el colapso de manera inmediata si la amplitud es lo suficientemente grande. Un fenómeno de especial interés asociado a las vibraciones es el de resonancia, que ocurre cuando la frecuencia de la excitación coincide con alguna de las frecuencias naturales del sistema. El comportamiento de una estructura en resonancia puede ser algo complejo de entender a priori, ya que una fuerza variable, pero de amplitud máxima constante, produce un movimiento con una amplitud que aumenta con el tiempo. El resultado de la resonancia es que fuerzas excitadoras de pequeña amplitud, que el sistema podría soportar si se aplicaran de manera estática, pueden provocar el colapso, pero no solo a largo plazo como ocurre con la fatiga, sino también durante la propia resonancia. Uno de los sucesos claves para el desarrollo de los análisis dinámicos, por su espectacularidad y gravedad, fue el colapso del puente de Tacoma en 1940. Durante los pocos meses que estuvo en servicio antes de su derrumbe, era habitual que el puente oscilase debido a que entraba en resonancia. El fenómeno de resonancia no era totalmente desconocido en aquella época, y se intentaron llevar a cabo ciertas modificaciones en el puente que asegurasen su integridad, pero algunas de ellas no fueron efectivas y otras no llegaron a tiempo, y el puente acabó cayendo debido al fenómeno de autoexcitación aerodinámica tras cerca de una hora oscilando (el día del derrumbe las oscilaciones fueron en forma de torsión, mientras que las oscilaciones que se habían observado previamente en el puente eran verticales). [2][3] U Introducción 2 Figura 1-1: Derrumbe del puente de Tacoma [4] Todo esto pone de manifiesto la importancia de conocer el comportamiento dinámico de cualquier máquina o estructura. El estudio dinámico de un sistema mecánico se lleva a cabo tanto en la etapa de diseño, sobre todo para conocer las frecuencias naturales y ajustarlas a los requisitos de diseño, como durante la vida útil del sistema, por ejemplo para verificar la validez de las hipótesis y simplificaciones llevadas a cabo en la etapa de diseño. Para conocer las frecuencias naturales y los modos de vibración de un sistema mecánico antes de su construcción se pueden llevar a cabo análisis analíticos o numéricos. Los primeros están muy limitados, ya que solo son aplicables a sistemas muy simples, con poca aplicación práctica, por lo que es común recurrir a modelos de elementos finitos. Figura 1-2: Primer modo de vibración de una pasarela peatonal obtenido mediante el MEF. Desde el punto de vista experimental, se puede proceder de dos formas. En primer lugar, se pueden tomar medidas de la respuesta de un sistema ante una excitación desconocida, generalmente en condiciones de servicio. Por otra parte, se puede medir la respuesta del sistema ante una excitación conocida, por lo general en condiciones controladas que permiten mayor precisión. Estos dos tipos de análisis se conocen como signal analysis y system analysis, respectivamente, según [5]. El Análisis Modal se puede definir formalmente como el conjunto de procesos involucrados en testar componentes o estructuras con el objetivo de obtener una descripción matemática de su comportamiento dinámico [6]. El Análisis Modal incluye la toma de medidas, así como el análisis, propiamente dicho, de las mismas. Los mayores avances en dinámica de sistemas mecánicos se llevaron a cabo a partir de mediados del siglo XX. Uno de los documentos más importantes sobre la materia fue el artículo de Kennedy y Pancu [9], publicado en 1947, y que es la base del método de Circle Fit para extraer propiedades modales. De aquí en adelante se escribieron multitud de artículos e incluso libros dedicados exclusivamente a la dinámica. Uno de los libros más destacados es Modal Testing: Theory, Practice and Application de D. J. Ewins [6], citado en multitud de ocasiones en artículos publicados con posterioridad. Además, en paralelo a estos avances teóricos, se desarrolló el Método de los Elementos Finitos, que también tiene una clara aplicación a problemas dinámicos, 3 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental y mejoraron los componentes electrónicos y digitales, lo que favoreció el desarrollo del análisis modal experimental. Las aplicaciones del Análisis Modal son varias. La primera de ellas es obtener las propiedades dinámicas del sistema real, y utilizarlas para comparar el sistema con el modelo utilizado, generalmente de elementos finitos. La comparación del sistema real con el sistema modelizado permite comprobar la validez del modelo, así como modificarlo en caso de ser necesario hasta que se corresponda con la realidad. Además, sabemos que las propiedades dinámicas de una estructura son inherentes a la misma. Utilizando esta importante propiedad, podemos usar el Análisis Modal para estudiar variaciones en las propiedades dinámicas, que estarán asociadas a algún cambio en la configuración de la máquina o estructura. Una de las aplicaciones puede ser encontrar fallos en una estructura, de una manera más fiable y sistemática que si se buscasen mediante simple inspección visual [7]. Esto permite detectar a tiempo imperfecciones en la estructura, lo suficientemente significativas como para modificar las propiedades dinámicas, pero no para provocar el fallo, y anticiparse al colapso. La variación de las propiedades dinámicas no siempre está debida a fallos no deseados, también se puede mejorar el comportamiento vibratorio variando o añadiendo elementos estructurales o no estructurales, lo que permite mejorar a posteriori el comportamiento dinámico de una estructura que presentaba alguna irregularidad desde su etapa de diseño, como por ejemplo ocurre en [8]. El análisis experimental permite cuantificar esta mejora, y modificar el modelo utilizado en un principio. Otras posibles aplicaciones pueden ser hallar el valor de la fuerza aplicada en casos en los que esta sea desconocida y sea de interés conocer su valor, o estudiar componentes por separado antes de ser ensamblados. Introducción 4 5 2 FUNDAMENTO TEÓRICO DE LAS VIBRACIONES L objetivo de este capítulo es definir todos los conceptos necesarios para comprender el comportamiento dinámico de cualquier sistema y obtener todas las ecuaciones que definen dicho comportamiento. Se empieza por los sistemas de un único grado de libertad, por su sencillez, y luego se estudian los de varios grados de libertad. En ambos casos se analizan los diferentes tipos de amortiguamiento que pueden darse, así como las situaciones de vibración libre y de vibración forzada. Además, se definen los distintos tipos de respuesta en frecuencia, así como las relaciones entre las mismas, haciendo especial hincapié en la receptancia, sobretodo en las expresiones que permiten expresarla como un sumatorio con tantos términos como modos de vibración, ya que son las ecuaciones básicas para el siguiente capítulo. 2.1. Sistemas de Un Grado de Libertad (SDOF) Aunque pocas estructuras reales pueden modelarse como un sistema de un único grado de libertad, el estudio de este tipo de sistemas resulta interesante porque sistemas más complejos de varios grados de libertad (MDOF) pueden representarse como la superposición lineal de sistemas SDOF. El modelo básico de un sistema de un grado de libertad es el siguiente: Figura 2-1: Modelo de un sistema de un grado de libertad [1] donde f(t) representa la fuerza, x(t) el desplazamiento, m la masa, k la rigidez y c (o d) el amortiguamiento. 2.1.1 Sistemas no amortiguados 2.1.1.1 Vibración libre En este caso se tiene f(t)=0, por lo que la ecuación de movimiento es la siguiente: 𝑚∙ 𝑥󰇘+ 𝑘∙𝑥=0 (2.1) La solución de la ecuación será del tipo 𝑥(𝑡)=𝑋∙𝑒𝑖∙𝜔∙𝑡. Sustituyendo en la ecuación de movimiento se obtiene: (𝑘−𝜔2∙𝑚)=0 (2.2) E Fundamento Teórico de las Vibraciones 6 de donde se obtiene la frecuencia natural del sistema:: 𝜔0     =√𝑘 𝑚 (2.3) 2.1.1.2 Vibración forzada Se considera ahora una excitación de la forma 𝑓(𝑡)=𝐹∙𝑒𝑖∙𝜔∙𝑡, y una ecuación de movimiento: 𝑚∙𝑥󰇘+𝑘∙𝑥=𝑓(𝑡) (2.4) La solución será del mismo tipo que para el caso de vibración libre: 𝑥(𝑡)=𝑋∙𝑒𝑖∙𝜔∙𝑡. Sustituyendo en la ecuación de movimiento: (𝑘−𝜔2∙𝑚)∙𝑋∙𝑒𝑖∙𝜔∙𝑡=𝐹∙𝑒𝑖∙𝜔∙𝑡 (2.5) A partir de esta ecuación se puede obtener la función de respuesta en frecuencia (FRF): 𝛼(𝜔)=𝑋 𝐹=1 (𝑘−𝜔2∙𝑚) (2.6) siendo α(ω) una de las posibles formas de la FRF, llamada Receptancia. Más adelante se presentarán otras posibles formas de expresar la FRF, así como la relación entre las mismas. 2.1.2 Sistemas con amortiguamiento viscoso Se define el amortiguamiento como el proceso de disipación de energía de un sistema durante la vibración. En sistemas sin fuerzas externas aplicadas (vibración libre) el amortiguamiento produce una reducción progresiva de la amplitud de la vibración hasta alcanzar la parada final. El amortiguamiento puede estar producido por diferentes causas, y es común que en un mismo sistema vibratorio coexistan simultáneamente diferentes mecanismos de amortiguamiento. Por lo general, todos estos tipos de amortiguamiento se pueden agrupar en tres bloques: amortiguamiento fluido, amortiguamiento seco (o de Coulomb) y amortiguamiento estructural (o histerético). El amortiguamiento fluido se origina por el paso de un fluido a través de unos orificios o por la interacción del sistema vibrante con el medio fluido en el que se encuentre inmerso. Dentro de este grupo se distinguen dos casos: amortiguamiento viscoso (cuando el movimiento del fluido es de régimen laminar) y amortiguamiento turbulento (cuando el movimiento del fluido es de régimen turbulento). Por lo general son más frecuentes los sistemas en los que el amortiguamiento viscoso predomina sobre el turbulento, y además el viscoso es más sencillo de analizar analíticamente ya que la fuerza de amortiguamiento varía linealmente con la velocidad, algo que no ocurre en el amortiguamiento turbulento, donde la fuerza de amortiguamiento es proporcional a la velocidad al cuadrado. Por ello, este apartado se centra en el estudio del amortiguamiento viscoso. 2.1.2.1 Vibración libre Considerando f(t)=0 y añadiendo el amortiguamiento viscoso a la ecuación de movimiento, ésta queda: 𝑚∙ 𝑥󰇘+ 𝑐∙ 𝑥󰇗+ 𝑘∙ 𝑥=0 (2.7) En este caso se debe tomar una solución de la forma 𝑥(𝑡)=𝑋∙𝑒𝑠∙𝑡, más general que para la estructura sin amortiguamiento, ya que ahora el exponente es complejo con parte real e imaginaria, mientras que antes solo tenía parte imaginaria. Sustituyendo se obtiene: (𝑚∙𝑠2+𝑐∙𝑠+𝑘)=0 (2.8) Resolviendo esta ecuación de segundo grado se llega a la solución: 𝑠1,2=− 𝑐 2∙𝑚±√𝑐2−4∙𝑘∙𝑚 2∙𝑚 =−𝜔0     ∙𝜁±𝑖∙𝜔0     ∙√1−𝜁2 (2.9) 7 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Siendo 𝜔0     =√𝑘 𝑚 y 𝜁= 𝑐 2∙√𝑘∙𝑚 , donde 𝜁representa el factor de amortiguamiento. La solución modal consta de un único modo de vibración con una frecuencia natural compleja: 𝑥(𝑡)=𝑋∙𝑒−𝜔0     ∙𝜁∙𝑡∙𝑒𝑖∙(𝜔0     ∙√1−𝜁2)∙𝑡=𝑋∙𝑒−𝑎∙𝑡∙𝑒𝑖∙𝜔0 ,∙𝑡 (2.10) Siendo 𝑎=𝜔0     ∙𝜁 y 𝜔0,=𝜔0     ∙√1−𝜁2 Analizando esta expresión se observa que la forma de la función de respuesta depende del valor del factor de amortiguamiento. Se pueden presentar los siguientes casos: Figura 2-2: Tipos de respuesta en función del coeficiente de amortiguamiento [1] Figura 2-3: Tipos de respuesta de un sistema amortiguando [10]  Si 𝜁=0 la respuesta es un movimiento armónico simple.  Si 𝜁>1 la respuesta está formada por dos exponenciales decrecientes y tenderá asintóticamente a cero sin que se produzca vibración. Se dice que el sistema está sobreamortiguado. Fundamento Teórico de las Vibraciones 8  Si 𝜁=1 la respuesta también tenderá asintóticamente a cero sin vibración, pero de manera más rápida que en el caso anterior. Se dice que el sistema está críticamente amortiguado. Al valor de c que hace 𝜁=1 se le denomina coeficiente de amortiguamiento crítico. 1= 𝑐𝑐𝑟 2∙𝑚∙𝜔𝑛=𝑐𝑐𝑟 2∙√𝑘∙𝑚 (2.11)  Si 0<𝜁<1 se dice que el sistema está subamortiguado. La respuesta consta de un movimiento armónico, de frecuencia 𝜔𝑑, modulado por una exponencial decreciente. 𝜔𝑑 se denomina frecuencia natural amortiguada y tiene la siguiente expresión: 𝜔𝑑=𝜔0     ∙√1−𝜁2 (2.12) Figura 2-4: Respuesta en el tiempo de un sistema subamortiguado [6] 2.1.2.2 Vibración forzada Tomando una fuerza de excitación 𝑓(𝑡)=𝐹∙𝑒𝑖∙𝜔∙𝑡, la ecuación de movimiento queda de la siguiente manera: 𝑚∙𝑥󰇘+𝑐∙𝑥󰇗+𝑘∙𝑥=𝑓(𝑡) (2.13) En este caso, la solución de la ecuación diferencial es de la forma 𝑥(𝑡)=𝑋∙𝑒𝑖∙𝜔∙𝑡, y sustituyendo se llega a: (−𝜔2∙𝑚+𝑖∙𝜔∙𝑐+𝑘)∙𝑋∙𝑒𝑖∙𝜔∙𝑡=𝐹∙𝑒𝑖∙𝜔∙𝑡 (2.14) de donde se puede obtener fácilmente la expresión de la receptancia: 𝛼(𝜔)=𝑋 𝐹=1 (𝑘−𝜔2∙𝑚)+𝑖∙(𝜔∙𝑐) (2.15) A diferencia de lo que ocurría en el caso en el que el sistema no presentaba amortiguamiento, la receptancia es ahora compleja, y contiene información de la magnitud y de la fase: |𝛼(𝜔)|=|𝑋| |𝐹|=1 √(𝑘−𝜔2∙𝑚)2+(𝜔∙𝑐)2 (2.16) arg(𝛼(𝜔))=arg(𝑋)−arg(𝐹)=tan−1 −𝜔∙𝑐 𝑘−𝜔2∙𝑚 (2.17) 2.1.3 Amortiguamiento estructural o histerético El amortiguamiento estructural se produce debido al rozamiento interno del material al sufrir deformación. Este fenómeno puede ocurrir en cualquier material, incluso en los que se considera que tienen comportamiento elástico. El rozamiento interno provoca que durante la deformación aparezca un desfase entre la fuerza aplicada y el desplazamiento producido, lo que se conoce como ciclo de histéresis. 9 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Figura 2-5: Ciclo de histéresis [1] El análisis de este tipo de sistemas se suele realizar únicamente en situaciones de vibración forzada, y su resolución se lleva acabo definiendo previamente el coeficiente de amortiguamiento estructural como: 𝑑=𝑐𝑒∙ 𝜔 (2.18) Entonces se tiene: (−𝜔2∙𝑚+𝑘+𝑖∙𝑑)∙𝑋∙𝑒𝑖∙𝜔∙𝑡=𝐹∙𝑒𝑖∙𝜔∙𝑡 (2.19) 𝑋 𝐹=𝛼(𝜔)=1 (𝑘−𝜔2∙𝑚)+𝑖∙(𝑑) (2.20) 𝛼(𝜔)=1𝑘 ⁄ (1−(𝜔𝜔𝑜 ⁄)2+𝑖∙𝜂) (2.21) donde η se define como el factor de pérdida 𝜂=2∙𝜁𝑒 (2.22) 2.1.4 Amortiguamiento seco o de Coulomb El amortiguamiento seco es el debido a la disipación de energía producida por las fuerzas de rozamiento que aparecen entre dos superficies que deslizan durante la vibración. El comportamiento de este tipo de sistemas está modelado por la ley de Coulomb, de manera que: 𝐹𝑟=−𝜇∙𝑁∙𝑠𝑖𝑔𝑛𝑜(𝑥󰇗)=−𝜇∙𝑁∙ 𝑥󰇗 |𝑥󰇗| (2.23) siendo N la fuerza normal entre las superficies y μ el coeficiente de rozamiento. Se observa que la fuerza de rozamiento actúa paralela a las superficies de rozamiento y en sentido contrario a la velocidad relativa de las superficies. Planteando equilibrio de fuerzas se obtiene la ecuación de movimiento: 𝑚∙𝑥󰇘+𝑘∙𝑥=−𝜇∙𝑁∙ 𝑥󰇗 |𝑥󰇗| (2.24) Esta ecuación es no lineal, pero se puede transformar en dos ecuaciones lineales, de manera que el estudio de estos sistemas se hace por tramos. Cuando 𝑥󰇗>0: 𝑚 ∙ 𝑥󰇘+ 𝑘∙ 𝑥=−𝜇∙ 𝑁 (2.25) Fundamento Teórico de las Vibraciones 16 Desarrollando matemáticamente el caso en el que los modos r y q sean un par conjugado se obtiene: 2∙𝜔𝑟∙𝜁𝑟={𝜓}𝑟 𝐻∙[𝐶]∙{𝜓}𝑟 {𝜓}𝑟 𝐻∙[𝑀]∙{𝜓}𝑟=𝑐𝑟 𝑚𝑟 (2.92) 𝜔𝑟2={𝜓}𝑟 𝐻∙[𝐾]∙{𝜓}𝑟 {𝜓}𝑟 𝐻∙[𝑀]∙{𝜓}𝑟=𝑘𝑟 𝑚𝑟 (2.93) siendo 𝜔𝑟 la frecuencia natural, { }𝐻 la transpuesta hermitiana y 𝑚𝑟 , 𝑘𝑟 , y 𝑐𝑟 la masa modal, la rigidez modal y el amortiguamiento modal respectivamente. 2.6.2. Vibración forzada La ecuación de movimiento de este tipo de sistemas es la siguiente: [𝑀]∙{𝑥󰇘(𝑡)}+[𝐶]∙{𝑥󰇗(𝑡)}+[𝐾]∙{𝑥(𝑡)}={𝑓(𝑡)} (2.94) Asumiendo que la solución del sistema matricial de ecuaciones es de la forma: {𝑥(𝑡)}={𝑋}∙𝑒𝑖∙𝜔∙𝑡 (2.95) Tras operar y sustituir en la ecuación de movimiento se llega a: {𝑋}=([𝐾]−𝜔2∙[𝑀]+𝑖∙𝜔∙[𝐶])−1∙{𝐹} (2.96) Esta ecuación presenta la relación entre las fuerzas y los desplazamientos, por lo que la obtención de la expresión de la matriz de receptancia es inmediata. Sin embargo, tal y como se discutió en el estudio de sistemas MDOF no amortiguados, esta forma de expresar la matriz de receptancia no es conveniente para su posterior resolución numérica. Es por ello por lo que se vuelve a buscar expresar la receptancia como un sumatorio de términos, solo que en este caso la forma de los autovalores dificulta el proceso. Para ello es necesario definir un vector de coordenadas, {u}, de dimensiones 2Nx1, que contiene la velocidad y el desplazamiento, y posteriormente reescribir la ecuación diferencial (2.81) para el caso de vibración libre del apartado anterior. {𝑢}={𝑥 𝑥󰇗} (2.97) [𝐶:𝑀]∙{𝑢󰇗}+[𝐾:0]∙{𝑢}={0} (2.98) En este punto se tienen N ecuaciones y 2N incógnitas, por lo que es necesario añadir otra igualdad matricial para cerrar el problema: [𝑀:0]∙{𝑢󰇗}+[0:−𝑀]∙{𝑢}={0} (2.99) Combinando ambas ecuaciones: [𝐶 𝑀 𝑀 0]∙{𝑢󰇗}+[𝐾 0 0 −𝑀]∙{𝑢}={0} (2.100) [𝐴]∙{𝑢󰇗}+[𝐵]∙{𝑢}={0} (2.101) Con esta ecuación ya se puede proceder como se había hecho hasta ahora, es decir, suponiendo la forma de la solución (en este caso {𝑢(𝑡)}={𝑈}∙𝑒𝑠∙𝑡), sustituyendo, y resolviendo el problema de autovalores y autovectores. En este caso se obtienen 2N autovalores con 2N autovectores asociados, designados por 𝑠𝑟 y {𝜃}𝑟 respectivamente. Estos atovalores y autovectores serán, por lo general, complejos en pares conjugados. (𝑠𝑟∙[𝐴]+[𝐵])∙{𝜃}𝑟=0 ; 𝑟=1,…,2𝑁 (2.102) Esta solución también posee propiedades de ortogonalidad: [𝛩]𝑇∙[𝐴]∙[𝛩]=[𝑎𝑟] (2.103) [𝛩]𝑇∙[𝐵]∙[𝛩]=[𝑏𝑟] (2.104) 𝑠𝑟=𝑏𝑟 𝑎𝑟 ; 𝑟=1,…,2𝑁 (2.105) 17 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Añadiendo a la ecuación (2.98) el vector {P} de fuerzas aplicadas es posible analizar la vibración forzada. Además, realizando un desarrollo análogo al realizado en el apartado 2.4.3., es posible llegar a la expresión de un término genérico de la matriz de receptancias como un sumatorio de términos. {𝑃}2𝑁𝑥1={𝐹 0} (2.106) {𝑋 … 𝑖∙𝜔∙𝑋}2𝑁𝑥1=∑{𝜃}𝑟 𝑇∙{𝑃}∙{𝜃}𝑟 𝑎𝑟∙(𝑖∙𝜔−𝑠𝑟) 2𝑁 𝑟=1 (2.107) Teniendo en cuenta que tanto los autovalores como los autovectores se obtienen en forma de pares conjugados el sumatorio se puede reescribir como: {𝑋 … 𝑖∙𝜔∙𝑋}=∑({𝜃}𝑟 𝑇∙{𝑃}∙{𝜃}𝑟 𝑎𝑟∙(𝑖∙𝜔−𝑠𝑟)+{𝜃}𝑟 𝐻∙{𝑃}∙{𝜃∗}𝑟 𝑎𝑟∗∙(𝑖∙𝜔−𝑠𝑟)) 𝑁 𝑟=1 (2.108) A partir de esta expresión se puede obtener el valor de un elemento genérico 𝛼(𝜔)𝑗𝑘 de la matriz de receptancia, que relaciona la respuesta 𝑋𝑗 con la fuerza 𝐹𝑘: 𝛼𝑗𝑘(𝜔)=∑(𝑅𝑗𝑘𝑟)+𝑖∙(𝜔 𝜔𝑟)∙(𝑆𝑗𝑘𝑟) 𝜔𝑟2−𝜔2+2∙𝑖∙𝜔∙𝜔𝑟∙𝜁𝑟 𝑁 𝑟=1 (2.109) donde: {𝑅𝑘𝑟 }=2∙(𝜁𝑟∙𝑅𝑒{𝐺𝑘𝑟 }−𝐼𝑚{𝐺𝑘𝑟 }∙√1−𝜁𝑟2) (2.110) {{ 𝑆𝑘𝑟 }}=2∙𝑅𝑒{𝐺𝑘𝑟 } (2.111) {𝐺𝑘𝑟 }=(𝜃𝑘𝑟 𝑎𝑟)∙{𝜃}𝑟 (2.112) Para llegar a esta forma reducida de la expresión se ha tenido en cuenta que: 𝑠𝑟=𝜔𝑟∙(−𝜁𝑟+𝑖∙√1−𝜁𝑟2 ) (2.113) Fundamento Teórico de las Vibraciones 18 19 3 ANÁLISIS MODAL EXPERIMENTAL STE capítulo comienza con una clasificación de todos los tipos de análisis modal experimental, con el que se pretende tener una visión general de los mismos. Seguidamente, nos centraremos en los métodos de un grado de libertad, comentando sus características, dificultades y limitaciones. Dentro de este tipo de análisis se desarrollarán en detalle tres de ellos: el Peak-Picking, el Cicle-Fit y el Line-Fit, describiendo cómo llevarlos a cabo. Siguiendo los procedimientos desarrollados en este capítulo se crearán los archivos de Matlab capaces de extraer las propiedades modales de cualquier sistema a partir de medidas experimentales de sus respuestas en frecuencia, que se pondrán en práctica en los siguientes capítulos con un sistema en concreto. 3.1 Clasificación de los Tipos de Análisis La primera posible clasificación hace referencia al dominio en el que se lleva a cabo el análisis. Aunque lo más común es realizar métodos de curve-fitting de las Funciones de Respuesta en Frecuencia (FRF), también es posible aplicar estos métodos en el dominio del tiempo. Para ello se tiene en cuenta que la transformada de Fourier inversa de la FRF da como resultado la Función de Respuesta a Impulso (IRF), que también es una función característica del sistema. Otra posible clasificación atiende al rango de frecuencias sobre el que se realizará cada análisis, y depende de si se extraen uno o varios modos al mismo tiempo. De esta manera se definen los métodos de un único grado de libertad (SDOF) y los métodos de varios grados de libertad (MDOF). Por último, se distingue entre análisis en los que se considera una única FRF, y análisis en los que se consideran varias FRF, lo cual está directamente relacionado con cómo se lleven a cabo las medidas. La manera más sencilla de obtener medidas experimentales de FRFs es la que se conoce como SISO (una entrada, una salida), de esta manera se obtiene un conjunto de FRFs, cada una de ellas obtenida individualmente (generalmente de forma secuencial). En segundo lugar, se pueden obtener varias FRF simultáneamente en distintos puntos a partir de una excitación aplicada sobre un único punto. Es lo que se conoce como SIMO (una entrada, varias salidas). Por último, se puede medir de forma simultánea la respuesta en varios puntos cuando la estructura es excitada en varios puntos al mismo tiempo. Esta manera de proceder se conoce como MIMO (varias entradas, varias salidas). En resumen, las posibles clasificaciones son:  Dependiendo del dominio considerado: o Dominio de la frecuencia (FRF) o Dominio del tiempo (IRF)  Dependiendo del número de modos que se extraen al mismo tiempo: o Un grado de libertad (SDOF) o Varios grados de libertad (MDOF) E Análisis Modal Experimental 20  Dependiendo del número de FRFs: o Una única FRF o Varias FRFs  Métodos Globales (SIMO)  Polirreferencia (MIMO) 3.2 Dificultades en el Análisis Modal Experimental Al igual que ocurría en el análisis teórico, una de las dificultades que aparecen al realizar un análisis modal experimental está en decidir cómo modelar el efecto del amortiguamiento. Para ello, es necesario hacer algunas suposiciones, que en caso de no ser correctas pueden provocar errores significativos en la estimación de los parámetros modales (y no sólo en los parámetros del amortiguamiento). Otro problema es debido al desconocimiento de cuántos modos hay en la FRF (o en el rango de frecuencias de la misma) que se va a estudiar. Mientras que en curvas dibujadas a partir de las expresiones teóricas se pueden reconocer fácilmente los picos, en curvas experimentales puede no verse tan claramente debido a los problemas típicos que aparecen en la toma de medidas, especialmente al ruido. El objetivo de los algoritmos de curve-fitting es minimizar el error entre la curva original y la resultante, y esto puede provocar que aparezcan más modos de los que hay realmente. La dificultad aparece a la hora de decidir cuáles de estos modos son reales y cuáles son ficticios. 3.3 Métodos de Un Grado de Libertad (SDOF) 3.3.1 Introducción Los métodos SDOF de análisis modal no implican que el sistema bajo estudio deba tener un único grado de libertad, sino que analizan los modos por separado, de manera secuencial. La principal desventaja es que es difícil separar el comportamiento de dos modos de vibración si éstos no están lo suficientemente alejados entre sí. Estos métodos basan su funcionamiento en la suposición de que cerca de una zona de resonancia el comportamiento está dominado por un único modo, siendo está suposición más válida cuanto más alejadas estén las resonancias una de otra. Partiendo de la expresión de la receptancia: 𝛼𝑗𝑘(𝜔)=∑ 𝐴𝑗𝑘𝑠 𝜔𝑠2−𝜔2+𝑖∙𝜂𝑠∙𝜔𝑠2 𝑁 𝑠=1 (3.1) Esta ecuación se puede reescribir como: 𝛼𝑗𝑘(𝜔)=𝐴𝑗𝑘𝑟 𝜔𝑟2−𝜔2+𝑖∙𝜂𝑟∙𝜔𝑟2+ ∑ 𝐴𝑗𝑘𝑠 𝜔𝑠2−𝜔2+𝑖∙𝜂𝑠∙𝜔𝑠2 𝑁 𝑠=1 ; ≠𝑟 (3.2) La suposición que se lleva a cabo para realizar estos métodos es que el segundo término es independiente de la frecuencia en un pequeño rango de frecuencias próximas a la frecuencia natural del modo r, por lo que: 𝛼𝑗𝑘(𝜔)𝜔≈𝜔𝑟≅𝐴𝑗𝑘𝑟 𝜔𝑟2−𝜔2+𝑖∙𝜂𝑟∙𝜔𝑟2+ 𝐵𝑗𝑘𝑟 (3.3) 3.3.2 Peak-Picking El método de Peak-Picking es el más sencillo y directo de todos los disponibles para estimar frecuencias naturales, modos de vibración y parámetros modales a partir de medidas experimentales. Por el contrario, también es el que tiene el rango de aplicación más limitado. En primer lugar, se supone que toda la respuesta 21 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental es debida a un único modo de vibración, por lo que los demás son ignorados. Además, también resulta crítico el valor del amortiguamiento. Si éste es demasiado pequeño, las medidas experimentales no serán lo suficientemente precisas, y si es demasiado grande la respuesta en una zona de resonancia se verá influenciada por varios modos de vibración. Aunque como se puede ver la aplicabilidad de este método es limitada, sigue siendo el método más común, ya que debido a su sencillez se aplica incluso a sistemas que no cumplen todos los requisitos como para que los resultados puedan considerarse válidos. Estos resultados, aún sin ser completamente correctos, se pueden utilizar como estimación inicial antes de iniciar otros métodos más sofisticados. El método de Peak-Picking consta de los siguientes pasos: 1) En primer lugar, partiendo de la representación en forma de Bode de los resultados experimentales, se identifican los picos de la curva en la gráfica módulo frente a frecuencia de la FRF. La frecuencia a la que se produce uno de los picos se toma como la frecuencia natural 𝜔𝑟 de ese modo. 2) La amplitud que alcanza el pico de la función se denota como |𝛼|. A partir de este valor, se caracterizan los half-power points como los puntos, a derecha e izquierda del máximo, cuyo módulo vale |𝛼|√2 ⁄. Figura 3-1: Método de Peak-Picking [11] 3) Una vez conocidas las frecuencias de los dos half-power points se puede obtener el amortiguamiento modal con la siguiente expresión: 𝜂𝑟=𝜔𝑎 2−𝜔𝑏 2 2∙𝜔𝑟2≅∆𝜔 𝜔𝑟 (3.4) 2∙𝜁𝑟=𝜂𝑟 (3.5) 4) Por último, suponiendo que la respuesta en la zona de resonancia es debida a un único modo, se puede obtener el valor de la constante modal de dicho modo. |𝛼|=𝐴𝑟 𝜔𝑟2∙𝜂𝑟 (3.6) 𝐴𝑟=|𝛼|∙𝜔𝑟2∙𝜂𝑟 (3.7) Como se puede observar, la precisión de los resultados obtenidos mediante este método depende de la precisión de los resultados experimentales cerca de la zona de resonancia. Esto supone un problema, ya que precisamente éstas zonas son las más susceptibles de contener errores en la medición, especialmente para sistemas ligeramente amortiguados. Por otra parte, solo se pueden obtener constantes modales reales, lo que implica que se puede suponer el modelo de amortiguamiento proporcional. También existe una alternativa a este método, conocida como método de la parte real de la receptancia, que mejora los resultados para casos en los que las zonas de resonancia no estén claramente alejadas una de otra. Análisis Modal Experimental 22 Figura 3-2: Método de la parte real de la receptancia [6] 3.3.3 Circle-Fit 3.3.3.1 Propiedades del círculo modal La principal característica en la que se basa este método es que al representar el diagrama de Nyquist de la FRF de un sistema de un grado de libertad se obtiene una curva que se asemeja claramente a un círculo. De hecho, dependiendo del tipo de amortiguamiento que se considere, se puede seleccionar una de las posibles formas de la FRF de forma que ésta genere un círculo de forma exacta. Para sistemas de varios grados de libertad, se observa que la representación en forma de Nyquist de la FRF produce una curva con distintas secciones que se asemejan a arcos de circunferencia en las proximidades de la zona de resonancia. Para sistemas con amortiguamiento estructural se representa la receptancia, que es la que genera círculos exactos al ser representada en el plano complejo. Para sistemas con amortiguamiento viscoso la FRF que produce círculos de manera exacta es la movilidad. En ambos casos el proceso a seguir es análogo a pesar de que la apariencia de las curvas sea significativamente diferente al estar una girada 90º respecto a la otra. Se va a realizar el desarrollo para el caso de amortiguamiento estructural, por lo que se parte de la expresión de la receptancia. 𝛼(𝜔)=1 𝜔𝑟2∙(1−(𝜔𝜔𝑟 ⁄)2+𝑖∙𝜂𝑟) (3.8) No se incluye en el numerador la constante modal, que lo que haría sería aumentar el tamaño del círculo (por el módulo de la constante modal) y girarlo (por el argumento de la constante modal). Figura 3-3: El círculo modal [11] 23 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Para cualquier valor de la frecuencia, la curva cumple: tan(𝛾)=𝜂𝑟 1−(𝜔𝜔𝑟 ⁄)2 (3.9) tan(90°−𝛾)=tan(𝜃 2)=1−(𝜔𝜔𝑟 ⁄)2 𝜂𝑟 (3.10) 𝜔2=𝜔𝑟2∙(1−𝜂𝑟∙tan(𝜃 2)) (3.11) Derivando esta ecuación con respecto a 𝜃 y derivando posteriormente respecto a 𝜔 se llega a una propiedad importante a la hora de localizar la frecuencia natural, y es que la máxima separación entre 2 puntos experimentales consecutivos a lo largo del arco de la circunferencia se produce cuando la frecuencia es la frecuencia natural. 𝑑𝜔2 𝑑𝜃 =−𝜔𝑟2∙𝜂𝑟∙(1−(𝜔𝜔𝑟 ⁄)2)2 2∙𝜂𝑟2 (3.12) 𝑑 𝑑𝜔(𝑑𝜔2 𝑑𝜃)=0 𝑐𝑢𝑎𝑛𝑑𝑜 (𝜔𝑟2−𝜔2)=0 (3.13) Tomando dos puntos del círculo, uno con frecuencia 𝜔𝑎 mayor que la frecuencia natural y otro con frecuencia 𝜔𝑏 menor, se puede deducir la expresión de 𝜂𝑟 que se introdujo en el apartado 3.3.2. 𝑡𝑎𝑛(𝜃𝑏 2)=1−(𝜔𝑏𝜔𝑟 ⁄)2 𝜂𝑟 (3.14) 𝑡𝑎𝑛(𝜃𝑎 2)=(𝜔𝑎𝜔𝑟 ⁄)2−1 𝜂𝑟 (3.15) Por lo que: 𝜂𝑟=𝜔𝑎 2−𝜔𝑏 2 𝜔𝑟2∙(𝑡𝑎𝑛(𝜃𝑎 2)+𝑡𝑎𝑛(𝜃𝑏 2)) (3.16) Para el caso particular en el que estos dos puntos considerados sean los half-power points, es decir, 𝜃𝑎=𝜃𝑏= 90°, la ecuación se simplifica. 𝜂𝑟=𝜔𝑎 2−𝜔𝑏 2 2∙𝜔𝑟2 (3.17) Para amortiguamientos pequeños esta ecuación se puede aproximar como: 𝜂𝑟≅𝜔𝑎−𝜔𝑏 𝜔𝑟 (3.18) Una vez conocidos el factor de pérdida del modo y la frecuencia natural, ya sólo queda conocer la constante modal para que el modo quede completamente definido. Como se indicó al inicio del capítulo, la constante modal afecta a la representación del círculo modal escalándolo (por el módulo de la constante modal) y girándolo (por el argumento de la constante modal). De esta manera se puede obtener el módulo de la constante modal a partir del diámetro del círculo como: 𝐷𝑗𝑘=|𝐴𝑗𝑘𝑟| 𝜔𝑟2∙𝜂𝑟 𝑟 (3.19) Y el argumento se obtiene teniendo en cuenta que el círculo completo esta girado de forma que el diámetro principal (el que pasa por el punto de frecuencia natural) está orientado un ángulo arg ( 𝐴𝑗𝑘) 𝑟respecto al eje imaginario negativo. Análisis Modal Experimental 24 Para sistemas con amortiguamiento viscoso, se parte de la expresión de la movilidad: 𝑌(𝜔)=𝑖∙𝜔 (𝑘−𝜔2∙𝑚)+𝑖∙(𝜔∙𝑐) (3.20) Y siguiendo un proceso análogo al del caso con amortiguamiento estructural se llega a: 𝜁𝑟≅𝜔𝑎−𝜔𝑏 2∙𝜔𝑟 (3.21) 3.3.3.2 Descripción del método Los pasos a seguir son los siguientes: 1) Seleccionar los puntos experimentales que se van a usar. Observando el diagrama de Bode de la FRF se puede localizar la zona de resonancia de cada uno de los modos. Para uno de ellos, se toman varios puntos por encima y varios por debajo de la frecuencia a la que se produce el máximo. Se pueden encontrar distintas recomendaciones en la bibliografía, algunas más restrictivas que otras, sobre cuantos puntos tomar. Hay que tener presente que si se toman demasiados puntos es posible que estos estén influidos por otros modos de vibración. Una de las posibles recomendaciones a seguir es tomar al menos 6 puntos, aunque lo óptimo sería poder dibujar tres cuartas partes del círculo. 2) Aplicar un algoritmo de circle-fitting y comprobar la calidad de dicho ajuste. Dentro de los múltiples algoritmos de circle-fitting disponibles se selecciona el que mejor funcione para los datos experimentales seleccionados, ya que por ejemplo algunos algoritmos están optimizados para situaciones en las que los puntos experimentales cubran un pequeño arco de la circunferencia completa, mientras que otros funcionan mejor en otras circunstancias. En este paso se obtiene el centro y el radio del círculo, así como un factor que de una idea de la validez de la aproximación. 3) Localizar la frecuencia natural. Ya se ha comentado que una de las características del círculo modal es que, suponiendo que los puntos experimentales están medidos tomando un incremento constante de frecuencia, cada uno de estos puntos, representados en el plano complejo, estará más alejado del punto que le precede y del que le sucede cuanto más próxima sea la frecuencia de ese punto a la frecuencia natural. Partiendo de la representación gráfica de estos puntos, y del centro del círculo de ajuste, puede hallarse el ángulo existente entre todos los pares de puntos consecutivos, de manera que la frecuencia natural se encontrará entre la frecuencia de los dos puntos que estén separados por el mayor ángulo. Figura 3-4: Localización de la frecuencia natural mediante circle-fit [6] Además de ésta, existen otras formas de hallar la frecuencia natural a partir de la representación del círculo modal, de manera que la frecuencia natural se puede definir como:  La frecuencia de máxima respuesta.  La frecuencia con la mayor receptancia imaginaria (en valor absoluto).  La frecuencia para la que la parte real de la receptancia valga cero. 25 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental En la siguiente figura se muestra la frecuencia natural de un sistema dado obtenida a partir de cada una de estas alternativas: Figura 3-5: Distintas formas de localizar la frecuencia natural mediante circle-fit [6] Se observa que entre cada una de ellas puede haber una diferencia significativa. Además hay que considerar que para calcular tanto el factor de pérdida como el módulo de la constante modal, en la ecuación aparece la frecuencia natural al cuadrado, por lo que el error en el cálculo de estos parámetros puede llegar a ser realmente significativo. 4) Calcular los coeficientes de pérdida. Utilizando la ecuación (3.22) se puede hallar el coeficiente de pérdida empleando cualquier par de puntos siempre que estos cumplan que uno tenga frecuencia mayor que la frecuencia natural y otro menor. Por tanto es posible obtener multitud de estimaciones de 𝜂𝑟 para un único círculo modal. Idealmente todos los valores de 𝜂𝑟 obtenidos deberían ser idénticos, por lo que este hecho se puede emplear como criterio para verificar la validez de las medidas tomadas. Desviaciones menores a un 5% se pueden considerar válidas, y desviaciones mayores del 20% son un indicativo de que las medidas no tienen la validez deseada. 𝜂𝑟=𝜔𝑎 2−𝜔𝑏 2 𝜔𝑟2∙(𝑡𝑎𝑛(𝜃𝑎 2)+𝑡𝑎𝑛(𝜃𝑏 2)) (3.22) 5) Por último, se calcula la constante modal (tanto el módulo como el argumento), a partir del diámetro del círculo y de los parámetros modales ya calculados 𝐷𝑗𝑘=|𝐴𝑗𝑘𝑟| 𝜔𝑟2∙𝜂𝑟 𝑟 (3.23) Además, si se quiere reconstruir la curva de la FRF usando la expresión teórica de la misma (3.3) sustituyendo en ella los valores de frecuencia natural, constante modal y coeficiente de pérdida recién calculados, es necesario hallar también la contribución del resto de modos al modo que se está estudiando. Para ello se calcula el valor de 𝐵𝑗𝑘𝑟 como la distancia entre el origen y el extremo del diámetro principal opuesto al de frecuencia natural. 3.3.4 Line-Fit (o método inverso) Este método es una alternativa al método de circle-fitting, y utiliza alguno de los conceptos presentados en dicho método. En un principio este método se basaba en que si una función representada en el plano complejo genera una curva circular, la inversa de dicha función produce una línea en el plano complejo. Esto se puede verificar partiendo de la ecuación de la receptancia para un sistema con amortiguamiento estructural y observando la forma que toma su inversa: Peak-Picking 32 4.2 Resultados Experimentales La estructura a estudiar es un voladizo de acero empotrado en una base de hormigón, de la que se pretenden obtener los parámetros modales de sus 5 primeros modos. La viga tiene una longitud de 0.9 metros y es de sección rectangular de lados 15 mm y 20 mm. El acero tiene una densidad de 7850 kg/m3 y un módulo elástico de 210∙109N/m2. Para tomar las medidas se utiliza un martillo instrumentado, con el que se golpea el lado más ancho de la sección, en 8 puntos de la viga y se mide la aceleración en el punto número 8. Figura 4-3: Viga en voladizo objeto de estudio [15] Punto de excitación 1 2 3 4 5 6 7 8 Coordenada z (m) 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Tabla 4-3: Ubicación de los puntos de excitación [15] La toma de medidas, y el posterior tratamiento de las mismas, se hará en dos fases, una para las frecuencias bajas (que en este caso comprenden los tres primeros modos de vibración) y otra para las frecuencias altas (modos cuarto y quinto). El proceso es análogo en ambos casos, con la diferencia que para las bajas frecuencias se utiliza una punta de martillo más blanda que para las altas frecuencias. De esta forma se obtienen mejores medidas que si se utilizase el mismo martillo para todo el rango de frecuencias. Los ficheros de Matlab PP_Baja_1FRF.m y PP_Alta_1FRF.m reciben como datos los ficheros de extensión .txt que contienen las frecuencias, las amplitudes y las fases (en radianes) de las acelerancias de cada una de las FRFs. Para aplicar el método de Peak-Picking sólo se necesitan los valores de las amplitudes, cuyas gráficas son las siguientes: 33 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Figura 4-4: Amplitud de la acelerancia del sistema para todos los puntos de excitación hasta 400 Hz Figura 4-5: Amplitud de la acelerancia del sistema para todos los puntos de excitación entre 400 y 1200 Hz A partir de estos datos de las acelerancias es posible pasar a las receptancias, que es la forma de FRF empleada para el método de Peak-Picking, utilizando para ello la ecuación: 𝐴(𝜔)=−𝜔2∙𝛼(𝜔) (4.4) Peak-Picking 34 Figura 4-6: Amplitud de la receptancia del sistema para todos los puntos de excitación hasta 400 Hz Figura 4-7: Amplitud de la receptancia del sistema para todos los puntos de excitación entre 400 y 1200 Hz Para cada una de las 8 FRFs se aplica el método de Peak-Picking, de forma que el fichero de Matlab devuelve como resultado las propiedades modales y la gráfica de la receptancia con los puntos de interés señalados. Por ejemplo, para la FRF del punto 8 los resultados son: Modo 1 Modo 2 Modo 3 Modo 4 Modo 5 Frecuencia Natural 14'75 89'5 250'25 494'65 814'65 Amortiguamiento 0'021829 0'0073501 0'0027604 0'0037078 0'0023114 Constante Modal 2'5922 3'2028 3'8277 2,3064 1,9925 Tabla 4-4: Propiedades modales obtenidas mediante Peak-Picking para la octava FRF 35 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Figura 4-8: Aplicación del Peak-Picking para la octava FRF (hasta 400 Hz) Figura 4-9: Aplicación del Peak-Picking para la octava FRF (entre 400 y 1200 Hz) Haciendo zoom en cada uno de los máximos se pueden ver con más detalle los puntos que se utilizan para calcular las propiedades modales: Peak-Picking 36 Figura 4-10: Detalle de una de las zonas de resonancia donde la amplitud de los half-power points es la amplitud del máximo entre √2 y las frecuencias de cada uno de ellos ha sido calculada por interpolación lineal entre dos puntos experimentales. Procediendo análogamente para las 8 FRFs se obtienen las propiedades modales (frecuencia natural, amortiguamiento y constante modal) de los 5 primeros modos de vibración: FRF 1 FRF 2 FRF 3 FRF 4 FRF 5 FRF 6 FRF 7 FRF 8 Modo 1 14'75 14'75 14'75 14'75 14'75 14'75 14'75 14'75 Modo 2 89'75 89'5 89'5 89'5 89'5 89'5 89'5 89'5 Modo 3 250'5 250'5 250'5 250'5 250'5 250'5 250'25 250'25 Modo 4 494'65 494'65 494'65 494'65 494'65 494'65 494'65 494'65 Modo 5 814'65 814'65 814'65 814'65 814'65 814'65 814'65 814'65 Tabla 4-5: Frecuencias naturales obtenidas mediante el método de Peak-Picking FRF 1 FRF 2 FRF 3 FRF 4 FRF 5 FRF 6 FRF 7 FRF 8 Modo 1 0'40966 0'73154 1'0758 1'3827 1'7049 2'2724 2'0696 2'5922 Modo 2 1'5783 2'3751 2'7413 2'4016 1'4646 0'058844 1'6997 3'2028 Modo 3 2'7259 2'5431 0'90016 1'2777 2'4531 1'6933 0'66826 3'8277 Modo 4 2'0116 0'35797 1'6892 1'4478 0'56027 1'6487 0'34605 2,3064 Modo 5 0'79461 1,2598 0'9107 1,239 1,0329 1'1878 1,1403 1,9925 Tabla 4-6: Constantes modales obtenidas mediante el método de Peak-Picking FRF 1 FRF 2 FRF 3 FRF 4 FRF 5 FRF 6 FRF 7 FRF 8 Modo 1 0'027864 0'025545 0'025668 0'023638 0'022796 0'022573 0'019594 0'021829 Modo 2 0'0076853 0'0077075 0'0077842 0'0075396 0'0076271 0'008212 0'0071034 0'0073501 Modo 3 0'0024685 0'0023769 0'0023327 0'0024862 0'0026775 0'0028654 0'0024763 0'0027604 Modo 4 0'0041281 0'0044683 0'0037917 0'0036442 0'0037 0'0034514 0'003384 0'0037078 Modo 5 0'0021698 0'0021788 0'002231 0'0022626 0'0022493 0'0024639 0'0031015 0'0023114 Tabla 4-7: Coeficientes de amortiguamiento obtenidos mediante el método de Peak-Picking 37 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Se observa que todas las constantes modales calculadas tienen valores reales, ya que el método de Peak- Picking no permite calcular constantes modales con parte imaginaria. En este caso, la suposición puede considerarse válida, ya que una viga en voladizo empotrada es uno de los más claros ejemplos de sistema con modos reales [16]. Para otros sistemas, especialmente aquellos en los que haya amortiguamientos localizados en vez de amortiguamientos proporcionalmente distribuidos, esta suposición no será válida y será necesario usar otro método. Una forma de verificar la consistencia del método para extraer parámetros modales puede ser comparar los valores de las frecuencias naturales y los coeficientes de amortiguamiento de cada uno de los modos obtenidos a partir de cada FRF, ya que al ser propiedades propias de cada uno de los modos, los valores obtenidos tienen que ser iguales en todas las FRFs. Esto se traduce en que en las tablas 4-1 y 4-3 los valores de cada una de las filas deberían ser iguales entre sí. Para estudiar la validez de los resultados experimentales se pueden comparar las frecuencias naturales obtenidas con la solución analítica, ya que al tratarse de un sistema estructural tan sencillo ésta es conocida [17]. En el caso de estructuras más complejas sin solución analítica conocida se podría recurrir a modelos de elementos finitos. 𝜔1=3.5160 2∙𝜋 ∙√𝐸∙𝐼 𝑚∙𝐿4=15.47 𝐻𝑧 (4.5) 𝜔2=22.034 2∙𝜋 ∙√𝐸∙𝐼 𝑚∙𝐿4=96.96 𝐻𝑧 (4.6) 𝜔3=61.697 2∙𝜋 ∙√𝐸∙𝐼 𝑚∙𝐿4=271.50 𝐻𝑧 (4.7) 𝜔4=120.90 2∙𝜋 ∙√𝐸∙𝐼 𝑚∙𝐿4=532.03 𝐻𝑧 (4.8) 𝜔5=199.86 2∙𝜋 ∙√𝐸∙𝐼 𝑚∙𝐿4=879.50 𝐻𝑧 (4.9) En la siguiente tabla se muestra una comparación entre los resultados obtenidos experimentalmente y los resultados analíticos, así como el error relativo entre los mismos. Analítica Experimental Error relativo (%) Frec. natural 1 (Hz) 15'47 14'75 4'65 Frec. natural 2 (Hz) 96'96 89'5 7'69 Frec. natural 3 (Hz) 271'50 250'5 7'73 Frec. natural 4 (Hz) 532'03 494'65 7'03 Frec. natural 5 (Hz) 879'50 814'65 7'37 Tabla 4-8: Comparación de las frecuencias naturales experimentales con las analíticas Se observan errores significativos, de hasta un 7’73%, por lo que es pertinente analizar las posibles fuentes de estas discrepancias entre resultados. En primer lugar, podría ser posible que las propiedades geométricas y/o mecánicas reales no coincidan con los valores tomados para los cálculos, pero al tratarse de un perfil normalizado de acero lo lógico es pensar que las propiedades reales no difieren significativamente de las dadas por el catálogo. Por otra parte, en la solución numérica no se considera el amortiguamiento, lo que puede provocar un error en los resultados, pero al tratarse de un sistema con amortiguamiento muy pequeño este error debería ser despreciable. Por último se puede considerar que el error provenga de la aplicación de las condiciones de contorno, en este caso que el empotramiento real no sea perfecto como el que se considera en Peak-Picking 38 el modelo analítico. De estas opciones la tercera es la que parece más lógica. Además de las propiedades modales obtenidas mediante el método de Peak-Picking, el archivo de Matlab también devuelve por pantalla como resultado las amplitudes y las fases de la acelerancia en los puntos de resonancia, lo que permite dibujar los modos de vibración: FRF 1 FRF 2 FRF 3 FRF 4 FRF 5 FRF 6 FRF 7 FRF 8 Modo 1 14'702 28'637 41'913 58'493 74'787 100'67 105'63 118'75 Modo 2 205'36 308'15 352'17 318'53 192'03 7'1657 239'28 435'75 Modo 3 1104'3 1069'9 385'89 513'94 916'19 590'95 269'86 1386'6 Modo 4 487'29 80'112 445'49 397'29 151'42 477'71 102'26 622'05 Modo 5 366'21 578'21 408'2 547'61 459'21 482'07 367'68 862'05 Tabla 4-9: Amplitudes de las 8 acelerancias en los 5 primeros puntos de resonancia FRF 1 FRF 2 FRF 3 FRF 4 FRF 5 FRF 6 FRF 7 FRF 8 Modo 1 0'248 0'132 0'301 0'346 0'337 0'256 0'416 0'299 Modo 2 -2’882 -2'78 -2'719 -2'702 -2'718 -0'645 0'433 0'438 Modo 3 0'184 0'215 0'111 -2'961 -2'95 -2'846 0'268 0'195 Modo 4 -3'078 -3'037 0'029 0'035 -3'085 -3'092 -3'141 0'023 Modo 5 0'054 -3'071 -3'047 0'108 0'078 -3'015 -2'913 0'045 Tabla 4-10: Fases de las 8 acelerancias en los 5 primeros puntos de resonancia Al igual que las frecuencias naturales, los modos de vibración de un voladizo tienen solución analítica conocida, que puede usarse para compararla con los modos de vibración obtenidos experimentalmente. Figura 4-11: Cinco primeros modos de vibración analíticos de un voladizo empotrado [18] 39 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Los modos de vibración obtenidos a partir de estos datos son los siguientes: Figura 4-12: Cinco primeros modos de vibración obtenidos a partir de las medidas experimentales Las curvas que representan las deformadas de la barra se dibujan a partir de los datos experimentales usando la función spline [19] de Matlab. Dadas las condiciones de contorno de la estructura (extremo izquierdo empotrado, extremo derecho libre), la pendiente de la curva en el primer punto debe ser nula y la pendiente en el último punto es libre. La función spline permite imponer pendientes, pero sólo si se imponen las pendientes en ambos extremos, por lo que esta opción no es válida para este caso. Esto provoca que la deformada en las proximidades del empotramiento no se corresponda con la realidad, especialmente en el quinto modo. Para poder representar los modos con más exactitud habría sido necesario obtener medidas experimentales de al menos un punto más próximo al origen. Peak-Picking 40 41 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental 5 CIRCLE-FIT N este capítulo se desarrolla un fichero de Matlab que aplica el método de Circle-Fit a cualquier sistema estructural, partiendo de las medidas experimentales de las FRFs, poniendo en práctica lo explicado en el apartado 3.3.3. El desarrollo de los apartados 5.1 y 5.2 es análogo al de los apartados 4.1. y 4.2., con la diferencia de que en este capítulo las medidas experimentales son diferentes que las usadas en el capítulo anterior, aunque el sistema estructural es el mismo. 5.1 Funcionamiento del archivo de Matlab Al igual que en el apartado 4.2. se comienza definiendo una receptancia con un único modo dándole valores a los parámetros modales que lo definen: Frecuencia Natural 15 Amortiguamiento 0.025 Constante Modal 0.5 Tabla 5-1: Parámetros modales tomados para definir la receptancia A partir de la expresión analítica de la receptancia se obtiene el diagrama de Bode de la misma: Figura 5-1: Diagrama de Bode de la receptancia E Circle-Fit 48 Figura 5-9: Diagrama de Nyquist de la acelerancia para la octava FRF Conociendo los valores de las partes real e imaginaria de la acelerancia se pueden calcular las amplitudes y fases de la misma, lo que permite representarla en forma de Bode. Figura 5-10: Diagrama de Bode de la acelerancia para la octava FRF A partir de la acelerancia se obtiene la receptancia, que se puede representar tanto en forma de Bode como de Nyquist. 49 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Figura 5-11: Diagrama de Bode de la receptancia para la octava FRF Figura 5-12: Diagrama de Nyquist de la receptancia para la octava FRF Para aplicar el método de Circle-Fit se utiliza el Nyquist de la receptancia, pero, como ya se ha comentado, no en todo el rango de frecuencias, sino que para cada modo hay que seleccionar un intervalo de frecuencias adecuado alrededor de la resonancia. Es conveniente, aunque no necesario, que la frecuencia natural esté lo más centrada posible en el intervalo que se elige. Circle-Fit 50 Para el primer modo se selecciona el intervalo (4.88281 – 34.17967) Hz: Figura 5-13: Diagramas de Bode y Nyquist de la octava FRF alrededor de la primera resonancia Para el segundo modo se selecciona el intervalo (63.47653 – 112.30463) Hz: Figura 5-14: Diagramas de Bode y Nyquist de la octava FRF alrededor de la segunda resonancia Para el tercer modo se selecciona el intervalo (224.60926 – 273.43736) Hz: Figura 5-15: Diagramas de Bode y Nyquist de la octava FRF alrededor de la tercera resonancia 51 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Para el cuarto modo se selecciona el intervalo (468.74976 – 517.57786) Hz: Figura 5-16: Diagramas de Bode y Nyquist de la octava FRF alrededor de la cuarta resonancia Para el quinto modo se selecciona el intervalo (791.01522 – 839.84332) Hz: Figura 5-17: Diagramas de Bode y Nyquist de la octava FRF alrededor de la quinta resonancia Los resultados que se obtienen tras aplicar el Circle-Fit a los 5 primeros modos de la octava FRF son los siguientes: FRF 8 Frec Natural Amortiguamiento Desv Estándar Módulo Cte Modal Argum Cte Modal ( º ) Modo 1 13.99 0.31532 0.029792 3.1268 -6.0432 Modo 2 90.298 0.049989 0.0079512 2.7195 -14.047 Modo 3 252.61 0.01861 0.0013633 2.2205 -5.4797 Modo 4 495.59 0.010416 0.00088542 2.2415 -7.2213 Modo 5 817.89 0.0073632 0.00042141 1.8914 -10.353 Tabla 5-3: Resultados para la octava FRF obtenidos mediante Cicle-Fit Circle-Fit 52 Procediendo análogamente para las 8 FRFs se obtienen todos los resultados. FRF 1 Frec Natural Amortiguamiento Desv Estándar Módulo Cte Modal Argum Cte Modal ( º ) Modo 1 15.704 0.23667 0.020732 0.44685 -19.735 Modo 2 90.025 0.066188 0.0045723 1.2811 173.8 Modo 3 253 0.031672 0.0027731 2.0565 -1.8341 Modo 4 496.17 0.012823 0.00041937 2.0641 176.77 Modo 5 817.89 0.0076842 0.00062288 0.90018 -3.3105 Tabla 5-4: Resultados para la primera FRF obtenidos mediante Cicle-Fit FRF 2 Frec Natural Amortiguamiento Desv Estándar Módulo Cte Modal Argum Cte Modal ( º ) Modo 1 15.16 0.20219 0.16024 0.33813 -38.865 Modo 2 90.123 0.053476 0.00098821 2.053 173.94 Modo 3 252.47 0.01996 0.00095824 1.9325 1.1619 Modo 4 495.31 0.010071 0.0014553 0.31973 -177.21 Modo 5 817.81 0.0074157 1.9519e-005 1.2912 174.91 Tabla 5-5: Resultados para la segunda FRF obtenidos mediante Cicle-Fit FRF 3 Frec Natural Amortiguamiento Desv Estándar Módulo Cte Modal Argum Cte Modal ( º ) Modo 1 13.989 0.30043 0.047496 0.97286 -14.131 Modo 2 90.355 0.052983 0.0012962 2.4081 176.44 Modo 3 252.58 0.01957 0.0023179 0.72958 7.2723 Modo 4 495.44 0.011761 0.00025494 1.6825 -4.2664 Modo 5 817.71 0.0079045 0.00035452 0.99897 177.82 Tabla 5-6: Resultados para la tercera FRF obtenidos mediante Cicle-Fit FRF 4 Frec Natural Amortiguamiento Desv Estándar Módulo Cte Modal Argum Cte Modal ( º ) Modo 1 13.966 0.31245 0.0025892 2.6771 -28.801 Modo 2 90.101 0.060516 0.0023104 2.2131 177.91 Modo 3 252.91 0.020102 0.0012517 0.76964 176.33 Modo 4 495.65 0.010282 0.00028067 1.4793 -1.0916 Modo 5 818.08 0.0072695 0.00013997 1.2231 -4.2499 Tabla 5-7: Resultados para la cuarta FRF obtenidos mediante Cicle-Fit FRF 5 Frec Natural Amortiguamiento Desv Estándar Módulo Cte Modal Argum Cte Modal ( º ) Modo 1 no no no no no Modo 2 89.74 0.038231 0.024007 0.80262 -161.33 Modo 3 252.83 0.019011 0.00035611 1.4264 178.61 Modo 4 495.58 0.0091617 0.00098014 0.48089 173.49 Modo 5 817.65 0.007956 0.00019305 0.95372 -3.6664 Tabla 5-8: Resultados para la quinta FRF obtenidos mediante Cicle-Fit 53 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental FRF 6 Frec Natural Amortiguamiento Desv Estándar Módulo Cte Modal Argum Cte Modal ( º ) Modo 1 no no no no no Modo 2 no no no no no Modo 3 252.25 0.021342 0.0031104 1.2525 -174.96 Modo 4 495.54 0.010992 0.00019966 1.7288 179.31 Modo 5 817.79 0.007742 0.0001297 1.1717 173.52 Tabla 5-9: Resultados para la sexta FRF obtenidos mediante Cicle-Fit FRF 7 Frec Natural Amortiguamiento Desv Estándar Módulo Cte Modal Argum Cte Modal ( º ) Modo 1 13.777 0.44291 0.095944 3.6349 -26.363 Modo 2 90.392 0.053902 0.017999 1.3356 -18.663 Modo 3 254.59 0.021587 0.011383 0.42683 -22.011 Modo 4 495.12 0.0087451 0.0039267 0.30972 -165.99 Modo 5 817.89 0.0077291 0.0006134 0.83661 -179.36 Tabla 5-10: Resultados para la séptima FRF obtenidos mediante Cicle-Fit FRF 8 Frec Natural Amortiguamiento Desv Estándar Módulo Cte Modal Argum Cte Modal ( º ) Modo 1 13.99 0.31532 0.029792 3.1268 -6.0432 Modo 2 90.298 0.049989 0.0079512 2.7195 -14.047 Modo 3 252.61 0.01861 0.0013633 2.2205 -5.4797 Modo 4 495.59 0.010416 0.00088542 2.2415 -7.2213 Modo 5 817.89 0.0073632 0.00042141 1.8914 -10.353 Tabla 5-11: Resultados para la octava FRF obtenidos mediante Cicle-Fit Se observa que los resultados del primer modo para las distintas FRFs (y también el del segundo modo para la sexta FRF) no son correctos o no se pueden calcular. Esto ya era esperable sólo con ver los diagramas de Bode o de Nyquist para frecuencias bajas. En vista de esto, parece que habría sido recomendable haber tomado las medidas para baja y alta frecuencia por separado como se hizo con las medidas utilizadas en el capítulo 4.2. Por lo demás, para cada modo se obtienen valores de la frecuencia natural y del amortiguamiento muy parecidos para las diferentes FRFs, y además para los amortiguamientos se tienen desviaciones estándar pequeñas, lo que sirve como indicador de que las medidas experimentales son buenas. Circle-Fit 54 55 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental 6 LINE-FIT N este capítulo se desarrolla un fichero de Matlab que aplica el método de Line-Fit a cualquier sistema estructural, partiendo de las medidas experimentales de las FRFs, poniendo en práctica lo explicado en el apartado 3.3.4. El desarrollo de los apartados 6.1 y 6.2 es análogo al de los apartados 5.1. y 5.2., y las medidas experimentales usadas en este capítulo son las mismas que se usaron en el capítulo 5. 6.1 Funcionamiento del archivo de Matlab Como se hizo en los apartados 4.1. y 5.1. se le aplica el método a una curva de receptancia con un solo modo de vibración definida a partir de su expresión analítica, imponiendo arbitrariamente sus parámetros modales. El objetivo es poder explicar cómo funciona el archivo de Matlab, y saber así en qué pasos hay que prestar especial atención. Los parámetros modales que se han tomado son los siguientes: Frecuencia Natural 15 Amortiguamiento 0.025 Constante Modal 0.5+0.01i Tabla 6-1: Parámetros modales tomados para definir la receptancia Con estos parámetros, la receptancia, representada en forma de Bode, es la siguiente: Figura 6-1: Diagrama de Bode de la receptancia E Line-Fit 56 También se podría representar el diagrama de Nyquist como se hizo en el apartado 5.1., pero para este método no tiene ninguna aplicación. Para aplicar el archivo de Matlab a unos datos experimentales habrá que obtener los valores de la receptancia a partir de las medidas que se hayan tomado, generalmente de la acelerancia. Una vez se tenga la receptancia, el proceso es análogo al que se explica a continuación. Lo siguiente es elegir el intervalo de frecuencias al que se le va a aplicar el método. Para ello hay que tener en cuenta las mismas consideraciones que ya se comentaron en el apartado 5.1. Conviene intentar elegir el mayor número de puntos experimentales posible, siempre y cuando el intervalo de frecuencias no sea lo suficientemente grande como para que la zona resonante que se está estudiando esté influida por otros modos. A diferencia del método de Circle-Fit, en el método de Line-Fit no hay que prestarle tanta atención a que el punto de máxima frecuencia esté centrado en el intervalo elegido. Por otra parte, hay que ser más cuidadoso con lo grande que es el intervalo, ya que los puntos alejados de la resonancia tienen una mayor influencia negativa en los resultados que en el método de Circle-Fit. Esto, que puede parecer algo complejo, no lo es realmente. En las figuras que devuelve el programa será fácil comprobar si el intervalo elegido es adecuado o no, y en el caso de no serlo es muy rápido modificarlo y volver a obtener resultados, ya que el tiempo de ejecución es despreciable. Para este ejemplo el intervalo elegido es (12.5 – 17.5) Hz, aunque al haber sólo un modo se obtendrían los mismos resultados para cualquier intervalo. Si hubiera más modos habría que resolverlos uno a uno, eligiendo un intervalo apropiado para cada uno de ellos. Esto es lo único que hay que modificar cada vez que se aplique el método, todo lo demás es automático. El diagrama de Bode de este intervalo es: Figura 6-2: Diagrama de Bode de la receptancia en las inmediaciones de la resonancia Se calculan ahora la nuevas FRFs 𝛼𝑗𝑘 ′ y ∆(𝜔) que se definieron en el capítulo 3.3.4. como: 𝛼𝑗𝑘 ′(𝜔)=𝛼𝑗𝑘(𝜔)−𝛼𝑗𝑘(𝛺) (6.1) ∆(𝜔)=𝜔2−𝛺2 𝛼𝑗𝑘 ′(𝜔) (6.2) siendo 𝛺 una variable conocida como frecuencia fijada, que es una de las frecuencias comprendidas en el intervalo que se está estudiando. En este ejemplo, el intervalo va desde 12.5 hasta 17.5 Hz, con una resolución de 0.5 Hz, por los que las frecuencias que se tienen son 12.5, 13, 13.5, 14, …, 17.5. 57 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental La manera de proceder es la siguiente: Primero se toma como frecuencia fijada 12.5 Hz y se calcula 𝛼𝑗𝑘 ′(𝜔) para el rango 12.5—17.5, que no es más que restarle a todos los valores de la receptancia el valor de la receptancia en 12.5 Hz. En segundo lugar se calcula ∆(𝜔), y se obtiene su parte real e imaginaria. Con esto ya se pueden representar 2 curvas: 𝑅𝑒(∆) frente a 𝜔2, e 𝐼𝑚(∆) frente a 𝜔2. Una vez hecho esto tomando 12.5 Hz como frecuencia fijada, se vuelve a repetir para 13 Hz, y luego para 13.5, y así sucesivamente hasta 17.5 Hz. El resultado son dos figuras, 𝑅𝑒(∆) frente a 𝜔2, e 𝐼𝑚(∆) frente a 𝜔2, con tantas curvas como posibles valores de frecuencias fijadas haya. Figura 6-3: 𝑅𝑒(∆) frente a 𝜔2 para los distintos valores de frecuencia fijada Figura 6-4: 𝐼𝑚(∆) frente a 𝜔2 para los distintos valores de frecuencia fijada Line-Fit 64 FRF 3 intervalo frec Frec Natural Amortiguamiento Módulo Cte Modal Argum Cte Modal (rad) Modo 1 14.648 - 34.18 12.941 0.3487 1.3877 0.069028 Modo 2 83.008 - 107.42 91.145 0.053382 2.6291 2.9691 Modo 3 229.49 - 268.55 254.65 0.019365 0.64375 -0.15401 Modo 4 478.52 - 507.81 494.87 0.013089 1.7371 -0.041437 Modo 5 800.78 - 834.96 817.5 0.0088487 1.0304 3.0835 Tabla 6-6: Resultados para la tercera FRF obtenidos mediante Line-Fit FRF 4 intervalo frec Frec Natural Amortiguamiento Módulo Cte Modal Argum Cte Modal (rad) Modo 1 14.648 - 34.18 13.983 0.28027 2.9302 -0.29255 Modo 2 83.008 - 107.42 90.114 0.06769 2.3294 3.1014 Modo 3 234.38 - 258.79 253.34 0.021128 0.79129 3.1024 Modo 4 478.52 - 507.81 496.08 0.0099525 1.4516 -0.072837 Modo 5 800.78 - 834.96 818.25 0.007347 1.2238 -0.058823 Tabla 6-7: Resultados para la cuarta FRF obtenidos mediante Line-Fit FRF 5 intervalo frec Frec Natural Amortiguamiento Módulo Cte Modal Argum Cte Modal (rad) Modo 1 14.648 - 34.18 13.432 0.58767 3.7035 -1.3316 Modo 2 83.008 - 107.42 90.555 0.042742 1.0749 3.0109 Modo 3 234.38 - 253.91 252.59 0.015665 1.1385 -2.9624 Modo 4 478.52 - 502.93 495.21 0.0096462 0.52778 3.1168 Modo 5 800.78 - 834.96 817.81 0.008702 0.99522 -0.088911 Tabla 6-8: Resultados para la quinta FRF obtenidos mediante Line-Fit FRF 6 intervalo frec Frec Natural Amortiguamiento Módulo Cte Modal Argum Cte Modal (rad) Modo 1 14.648 - 34.18 13.439 0.38667 2.7836 -0.27948 Modo 2 no no no no no Modo 3 234.38 - 258.79 253.01 0.022199 1.2253 3.0362 Modo 4 468.75 - 498.05 495.82 0.010629 1.7331 3.0601 Modo 5 795.9 - 830.08 817.68 0.0085209 1.2071 3.0224 Tabla 6-9: Resultados para la sexta FRF obtenidos mediante Line-Fit FRF 7 intervalo frec Frec Natural Amortiguamiento Módulo Cte Modal Argum Cte Modal (rad) Modo 1 9.7656 - 34.18 12.089 0.43485 3.3815 0.031553 Modo 2 73.242 - 102.54 89.119 0.056519 1.4361 0.012704 Modo 3 234.38 - 268.55 252.3 0.017141 0.37455 0.088299 Modo 4 468.75 - 502.93 496.48 0.010927 0.3795 2.8488 Modo 5 791.02 - 825.2 818.83 0.0079029 0.84794 2.9152 Tabla 6-10: Resultados para la séptima FRF obtenidos mediante Line-Fit FRF 8 intervalo frec Frec Natural Amortiguamiento Módulo Cte Modal Argum Cte Modal (rad) Modo 1 9.7656 - 34.18 14.2 0.34872 3.5545 -0.071361 Modo 2 73.242 - 102.54 89.759 0.046237 2.7882 -0.11103 Modo 3 239.26 - 273.44 252.91 0.019104 2.3125 -0.039674 Modo 4 468.75 - 517.58 495.28 0.01028 2.2538 -0.047498 Modo 5 791.02 - 839.84 817.5 0.0073856 1.9028 -0.10679 Tabla 6-11: Resultados para la octava FRF obtenidos mediante Line-Fit 65 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental Como ocurría en el método de Circle-Fit, para el primer modo de vibración se obtienen valores del coeficiente de amortiguamiento excesivamente grandes. Todo parece indicar que es debido a que la toma de medidas no era óptima para las frecuencias más bajas. Aún así, con el método de Line-Fit al menos se pueden obtener resultados para el primer modo a partir de las 8 FRFs, algo que no ocurrió con el Circle-Fit. Sin embargo, los resultados del segundo modo para la sexta FRF no se han considerado válidos, ya que se obtenía un amortiguamiento negativo. En las tablas de resultados se incluye el intervalo de frecuencias seleccionado para estudiar cada modo, que en la mayoría de los casos ha sido menor que el usado al aplicar Circle-Fit, con el objetivo de minimizar la influencia de otros modos de vibración en los resultados. Por lo demás, se observa que para cada modo de vibración se obtienen valores de frecuencia natural y amortiguamiento parecidos para las distintas FRFs, y que los argumentos de las constantes modales son próximos a 0 y a ±𝜋 (indicativo de modos de vibración reales). Además, los resultados obtenidos mediante Circle-Fit son por lo general parecidos a los de Line-Fit. Line-Fit 66 67 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental 7 CONCLUSIONES Y DESARROLLOS FUTUROS 7.1 Conclusiones El objetivo del presente Trabajo Fin de Grado ha sido realizar en Matlab unos archivos que permitan llevar a cabo el análisis modal de cualquier sistema estructural a partir de las medidas experimentales de sus FRFs. El análisis modal experimental se puede definir como el conjunto de procesos involucrados en testar componentes o estructuras con el objetivo de obtener una descripción matemática de su comportamiento dinámico [6]. El análisis modal experimental incluye también el proceso de toma de medidas, aunque en este TFG se ha partido de dos conjuntos de medidas experimentales tomadas previamente durante el desarrollo de una práctica impartida en la Escuela. En los capítulos más teóricos del Trabajo se han definido varios conceptos fundamentales, siendo uno de ellos el de Función de Respuesta en Frecuencia (FRF), que es una relación entre la fuerza y la respuesta del sistema (posición, velocidad o aceleración). Por lo general las medidas experimentales recogidas serán de la acelerancia del sistema, a partir de la cual es sencillo obtener la receptancia, que es la forma de FRF que se utilizará posteriormente en la aplicación de los métodos de análisis modal. Para manejar la receptancia, en vez de usar la expresión que se puede obtener de manera más inmediata, se pone en forma de sumatorio, con tantos términos como modos de vibración, de forma que es más sencillo comprender su comportamiento y su uso es mucho más eficiente desde un punto de vista numérico. El Trabajo se centra en tres métodos de los llamados SDOF (Single Degree Of Motion): el Peak-Picking, el Circle-Fit y el Line-Fit. La aplicación de estos métodos está limitada a sistemas con frecuencias naturales suficientemente alejadas entre sí, de forma que un modo de vibración no se vea afectado por los demás. A pesar de lo que el nombre puede parecer indicar, su aplicación no se limita a sistemas de un grado de libertad, sino que el sistema puede tener varios modos de vibración, pero cada uno de estos modos se estudia por separado, de manera secuencial, a diferencia de lo que ocurre en los métodos MDOF. El primero de los métodos desarrollados ha sido el de Peak-Picking. Se ha visto que es el más sencillo y rápido, ya que apenas requiere interacción por parte del usuario. Sin embargo, también es el que está más limitado. En primer lugar, los resultados dependen mucho de los puntos experimentales más cercanos a la resonancia, que son los más susceptibles de contener errores. Se ha visto también que sólo permite obtener constantes modales reales, esto no ha supuesto un gran problema para el sistema que se ha estudiado, ya que es un claro ejemplo de sistema con modos de vibración reales, pero se puede perder mucha información si se le aplica este método a un sistema con modos de vibración complejos. Se ha comprobado también la gran incidencia que tiene resolución de las medidas en los resultados obtenidos, ya que disminuir ligeramente la resolución puede provocar cambios en el orden de magnitud de los resultados. Por lo tanto, la validez de los resultados no sólo depende de que los puntos experimentales sean correctos, sino también de lo separados que están entre sí. El siguiente método es el de Circle-Fit, y lo primero que se ha visto es que con las medidas experimentales usadas para el Peak-Picking no era posible obtener resultados válidos. Esto pone de manifiesto que el tipo de medidas tomadas es un factor tan importante como el tipo de sistema a la hora de elegir un método respecto a otros. En este método tiene menos importancia la resolución, ya que es capaz de tener en cuenta lo que ocurre entre dos puntos experimentales adyacentes, pero tiene más importancia la medida de la fase, ya que errores en ésta afectan al diagrama de Nyquist que se utiliza. Por último, se ha desarrollado el método de Line-Fit. Para este método ha sido necesario definir una nueva FRF, y se ha comprobado que, en comparación con el método de Circle-Fit, los puntos del intervalo más Conclusiones y Desarrollos Futuros 68 alejados de la resonancia tienen mayor incidencia en los resultados. Esto quiere decir que los resultados son más sensibles a modos adyacentes al que se esté estudiando, por lo que hay que ser más conservador a la hora de elegir el intervalo de frecuencias. En cuanto a los resultados, se observa que son consistentes para cada método, ya que se cumple que para un mismo método, los valores de la frecuencia natural y del amortiguamiento de un determinado modo de vibración que se obtienen a partir de las distintas FRFs suelen ser bastante parecidos entre sí. Para el método de Circle-Fit tenemos un indicador de la validez de las medidas experimentales, que es la desviación típica de los valores del coeficiente de amortiguamiento que se obtienen para un mismo modo. Se obtienen desviaciones pequeñas, que es lo adecuado. Además, mediante los métodos de Line-Fit y Circle-Fit se obtienen argumentos de la constante modal próximos a 0 o a ±𝜋 (excepto para el primer modo de vibración, ya que las medidas no son óptimas para frecuencias bajas), lo que era de esperar en un sistema con modos de vibración reales. Mediante Circle-Fit y Line-Fit se obtienen resultados parecidos, aunque estos son algo diferentes a los obtenidos por Peak-Picking, sobretodo en el caso del amortiguamiento. El hecho de que no se hayan usado las mismas medidas medidas para el Peak-Picking que para los otros dos métodos dificulta sacar conclusiones al respecto, ya que en algunos casos se puede ver, simplemente observando las curvas experimentales, que estas presentan diferencias entre sí, seguramente porque no se llevaron a cabo en las mismas condiciones. Se pueden comparar las frecuencias naturales obtenidas mediante análisis modal con las obtenidas analíticamente, ya que al tratarse de un sistema estructural sencillo éstas son conocidas. Se observan diferencias de varios Hz. La hipótesis más lógica parece que el empotramiento del voladizo no sea perfecto, como el que se considera en el modelo. 7.2 Desarrollos Futuros Para poder comparar entre sí adecuadamente los resultados obtenidos mediante los distintos métodos, sería recomendable realizar el ensayo de forma que un mismo conjunto de medidas se pueda analizar mediante los tres métodos. De esta forma se sabría que las discrepancias en los resultados son provocadas por los propios métodos, y no por las condiciones en las que se tomaron las medidas. Además, estos resultados podrían compararse también con los obtenidos por algún programa comercial de análisis modal con métodos MDOF. Se ve que el sistema estudiado tiene unos modos de vibración claramente diferenciados, por lo que, en principio, no debería haber una gran diferencia. En los apartados 4.1., 5.1., y 6.1., en vez de aplicarle los métodos de análisis modal experimental a unas medidas tomadas en el laboratorio, se utilizan curvas de receptancias generadas a partir de su expresión teórica. Esto se ha hecho para poder explicar el funcionamiento de cada método, y también para estudiar cómo afecta la resolución de las medidas a los resultados del Peak-Picking. Siguiendo esta línea, se puede estudiar en profundidad como afectan diferentes factores a la validez de los resultados. Por ejemplo, es posible generar una receptancia que, en vez de tener un modo, tenga dos, y estudiar como varían las propiedades modales de uno de los modos, obtenidas mediante alguno de los tres métodos, respecto a variaciones en los parámetros que definen el otro modo. Uno de los parámetros que se puede ir variando es, por ejemplo, la frecuencia natural del segundo modo, lo que permitiría ver como varían los resultados a medida que la hipótesis de que los modos están lo suficientemente separados se va dejando de cumplir. De esta manera se podría estimar cuando se puede considerar que dos modos están lo suficientemente alejados entre sí como para poder considerar válidos los resultados de métodos SDOF, ya que lo cerca o lejos que están dos modos no es algo que esté claramente definido en la bibliografía, a pesar de ser un concepto fundamental. 69 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental REFERENCIAS [1] Jaime Domínguez Abascal, Juana Mayo Núñez «Tecnología de Máquinas: Teoría.» Universidad de Sevilla [2] Arturo Quirantes «La resonancia bien entendida: el puente de Tacoma Narrows.» 2012 [3] «https://en.wikipedia.org/wiki/Tacoma_Narrows_Bridge_(1940)» [4] Patrick Guillaume «Modal Analysis» Department of Mechanical Engineering, Vrije Universiteit Brussel [5] Ole Døssing «Structural Testing. Part 1: Mechanical Mobility Measurements.» Brüel & Kjær, 1988 [6] D. J. Ewins «Modal Testing: Theory, practice and application. Second Edition» Research Studies Press Ltd., 2000 [7] Ramón Rojas Díaz (tutor: Pedro Galvín Barrera) «Proyecto Fin De Carrera: Detección de fallos en estructuras mediante la medida de la variación de sus propiedades dinámicas.» Universidad de Sevilla, 2006 [8] Jorge Pérez Aracil, Alejandro M. Hernández Díaz, Javier F. Jiménez Alonso, Andrés Sáez Pérez «Experimental and numerical assessment of the unexpected damping in a slender footbridge at Murcia (Spain).» International Operational Modal Analysis Conference, 2015 [9] Charles C. Kennedy, C. D. P. Pancu «Use of Vectors in Vibration Measurement and Analysis.» Journal of the Aeronautical Sciences Volume 14, 1947 [10] Matthew T. Clay «Differential Equations and Laplace Transforms. Free Damped Motion.» University of Arkansas, 2013 [11] S. Ziaei-Rad «Modal Analysis and Testing: Modal Parameter Extraction Methods.» [12] Nuno Manuel Mendes Maia «Extraction of valid modal properties from measured data in stuctural vibrations.» Imperial College of Science, Technology and Medicine, 1988 [13] José M. Goicolea «Análisis sísmico de estructuras: Dinámica estructural.» Universidad Politécnica de Madrid, 2004 [14] «https://es.mathworks.com/help/signal/ref/findpeaks.html» [15] Rosario Chamorro Moreno, Jaime Domínguez Abascal, María Esther Reina Romo «Ampliación de Cinemática y Dinámica de Máquinas. Practica nº 4: 2ª Parte.» Universidad de Sevilla [16] Ole Døssing «Structural Testing. Part 2: Modal Analysis and Simulation.» Brüel & Kjær, 1988 [17] Manuel Hidalgo Martínez «Teoría de Vibraciones.» Universidad de Córdoba, 2009 [18] Ivan Mustakerov, Daniela Borissova «Multi-Criteria Model for Optimal Number and Placement of Sensors for Structural Health Monitoring: Lexicographic Method Implementation.» Institute of Information and Communications Technologies at Bulgarian Academy of Sciences, Department of Information Processes and Decision Support Systems, 2014 [19] «https://es.mathworks.com/help/matlab/ref/spline.html» [20] Nikolai Chernov «Fitting ellipses, circles, and lines by least squares: Matlab code for circle fitting algotithms.» 2012 Referencias 70 [21] David K. Anthony «Ventajas del uso de la parte real del espectro para la medida de los parámetros modales frente a errores de fase.» 40ª Edición del Congreso Español de Acústica (TecniAcústica), 2009 [22] «https://es.mathworks.com/help/matlab/ref/polyfit.html» 71 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental ANEXOS 1. PEAK_PICKING_BAJA.m %************************************************************************** %APLICACION DEL METODO DE PEAK-PICKING A LAS MEDIDAS EXPERIMENTALES DE BAJA %FRECUENCIA (HASTA 400 Hz) %************************************************************************** %LINEAS QUE HAY QUE MODIFICAR: En las lineas 19 y 22 se seleccionan los %archivos .txt que contienen las medidas experimentales. clc clear close all format short g %***************************************************************** %DEFINICION DE VARIABLES A PARTIR DE LOS RESULTADOS EXPERIMENTALES %***************************************************************** %Cargo el fichero .txt de amplitudes y frecuencias A1 = load ('H88_baja_frec.txt'); %Cargo el fichero .txt de fases B1 = load ('G8_baja_frec.txt'); %Defino el vector de amplitudes, con tantos terminos como medidas %experimentales amp_acel_comp = A1(:,2); %Defino el vector de frecuencias, con tantos terminos como medidas %experimentales frec_acel_comp = A1(:,1); %Defino el vector de fases, con tantos terminos como medidas %experimentales fase_acel_comp(:,1) = B1(:,1); %Me quedo con las 1601 primeras componentes de los vectores de Amplitudes, Fases %y Frecuencias (los primeros 400 Hz) amplitud_acelerancia = amp_acel_comp(1:1601); frecuencia_acelerancia = frec_acel_comp(1:1601); fase_acelerancia = fase_acel_comp(1:1601); %Represento en una grafica la amplitud de la acelerancia de la FRF %elegida, solo hasta 400 Hz figure(1); loglog(frecuencia_acelerancia,amplitud_acelerancia); title('Amplitud de la Acelerancia hasta 400 Hz') xlabel('Frecuencia [Hz]') ylabel('Amplitud de la Acelerancia') Anexos 72 %Dibujo el Bode de la acelerancia (amplitud y fase vs frecuencia) hasta 400 %Hz %figure(2) %amplitud_acel_grafica = subplot(2,1,1); %loglog(amplitud_acel_grafica,frecuencia_acelerancia,amplitud_acelerancia) %title('Diagrama de Bode de la Acelerancia') %xlabel('Frecuencia [Hz]') %ylabel('Amplitud de la Acelerancia') %fase_acel_grafica = subplot(2,1,2); %semilogx(fase_acel_grafica,frecuencia_acelerancia,fase_acelerancia) %xlabel('Frecuencia [Hz]') %ylabel('Fase de la Acelerancia [rad]') %************************************************************************ %LA ACELERANCIA VIENE DEFINIDA POR SU AMPLITUD Y SU FASE, LA PASO A FORMA %COMPLEJA (PARTE REAL + PARTE IMAGINARIA) %************************************************************************ n = length(frecuencia_acelerancia) for i=1:1:n Parte_Real_Acelerancia(i) = amplitud_acelerancia(i) * cos(fase_acelerancia(i)); Parte_Imag_Acelerancia(i) = amplitud_acelerancia(i) * sin(fase_acelerancia(i)); end for i=1:1:n Acelerancia_Compleja (i) = Parte_Real_Acelerancia(i) + j*Parte_Imag_Acelerancia(i); end %************************************************* %CALCULO LA RECEPTANCIA A PARTIR DE LA ACELERANCIA %************************************************* for i=1:1:n Receptancia_Compleja (i) = Acelerancia_Compleja(i) / - (frecuencia_acelerancia(i))^2; amplitud_receptancia(i) = abs(Receptancia_Compleja(i)); fase_receptancia (i) = angle(Receptancia_Compleja(i)); end %******************************** %DIBUJO EL BODE DE LA RECEPTANCIA %******************************** %Bode %figure(3) %amplitud_recep_grafica = subplot(2,1,1); %loglog(amplitud_recep_grafica,frecuencia_acelerancia,amplitud_receptancia) %title('Diagrama de Bode de la Receptancia') %xlabel('Frecuencia [Hz]') %ylabel('Amplitud de la Receptancia') %fase_recep_grafica = subplot(2,1,2); %semilogx(fase_recep_grafica,frecuencia_acelerancia,fase_receptancia) %xlabel('Frecuencia [Hz]') %ylabel('Fase de la Receptancia [rad]') 73 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental %Solo amplitud vs frecuencia figure(4) loglog(frecuencia_acelerancia,amplitud_receptancia) title('Amplitud de la Receptancia hasta 400 Hz') xlabel('Frecuencia [Hz]') ylabel('Amplitud de la Receptancia') %***************************************************** %APLICACION DEL PEAK PICKING A LA CURVA DE ACELERANCIA %***************************************************** %La función findpeaks encuentra los máximos del vector %"amplitud_acelerancia" y los almacena en el vector %"Amplitud_Acelerancia_Maximos". Además, almacena en el vector %"Componente_de_los_Maximos" las componentes de "amplitud_acelerancia" en %las que se producen los máximos. Está hecho con la versión R2008a de %Matlab, en versiones posteriores creo que la función cambia un poco su %funcionamiento %Uso la funcion findpeaks para encontrar los máximos y las componentes de %los máximos Amplitud_Acelerancia_Maximos = []; Componente_de_los_Maximos = []; [Amplitud_Acelerancia_Maximos,Componente_de_los_Maximos]=findpeaks(amplitud_a celerancia,'MinPeakDistance',50,'Threshold',1); Componente_de_los_Maximos; Amplitud_Acelerancia_Maximos; num_maximos = length(Componente_de_los_Maximos); %Uso el vector "Componente_de_los_Maximos" para encontrar las frecuencias %naturales en el vector "frecuencia_acelerancia" Frecuencias_Naturales=[]; for j=1:1:num_maximos Frecuencias_Naturales(j)=frecuencia_acelerancia(Componente_de_los_Maximos(j)) ; end Frecuencias_Naturales; %Uso el vector "Componente_de_los_Maximos" para encontrar las fases de los %picos en el vector "fase_acelerancia" Fase_Acelerancia_Maximos=[]; for j=1:1:num_maximos Fase_Acelerancia_Maximos(j)=fase_acelerancia(Componente_de_los_Maximos(j)); end Fase_Acelerancia_Maximos; %******************************************************************* %CALCULO DE LAS CONSTANTES MODALES Y LOS FACTORES DE AMORTIGUAMIENTO %******************************************************************* %Uso el vector "Componente_de_los_Maximos" para encontrar las amplitudes de %la receptancia en los maximos Amplitud_Receptancia_Maximos=[]; for j=1:1:num_maximos Anexos 80 epsilon=1e-12; ynew=1e+20; IterMax=20; xnew = 0; % Newton's method starting at x=0 for iter=1:IterMax yold = ynew; ynew = A0 + xnew*(A1 + xnew*(A2 + 4.*xnew*xnew)); if (abs(ynew)>abs(yold)) disp('Newton-Pratt goes wrong direction: |ynew| > |yold|'); xnew = 0; break; end Dy = A1 + xnew*(A22 + 16*xnew*xnew); xold = xnew; xnew = xold - ynew/Dy; if (abs((xnew-xold)/xnew) < epsilon), break, end if (iter >= IterMax) disp('Newton-Pratt will not converge'); xnew = 0; end if (xnew<0.) fprintf(1,'Newton-Pratt negative root: x=%f\n',xnew); xnew = 0; end end % computing the circle parameters DET = xnew*xnew - xnew*Mz + Cov_xy; Center = [Mxz*(Myy-xnew)-Myz*Mxy , Myz*(Mxx-xnew)-Mxz*Mxy]/DET/2; Par = [Center+centroid , sqrt(Center*Center'+Mz+2*xnew)]; end % CircleFitByPratt 81 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental 4. CIRCLE_FIT.m %******************************************************************* %APLICACIÓN DEL MÉTODO DE CIRCLE-FIT PARA LA OBTENCIÓN DE PARÁMETROS %MODALES %******************************************************************* %LINEAS QUE HAY QUE MODIFICAR % LINEA 40: Hay que poner la resolución, en Hz, de las medidas % experimentales % LINEA 42: Hay que poner el nombre del fichero .txt que contiene las % medidas experimentales % LINEAS 135-137: Hay que seleccionar qué intervalo de frecuencias se va % a tomar para estudiar un modo en particular. En las LINEAS 129-133 % vienen las que se han usado para este ejemplo en concreto. Hay que % elegir el intervalo de forma que la frecuencia natural esté incluida % en él, y a ser posible que esté lo más centrada posible en el intervalo % LINEA 266: Hay que darle un valor a la variable num_puntos. Se calculan % num_puntos^2 valores del coeficiente de amortiguamiento, y se da como % resultado la media de todos ellos. También se incluye la desviación % típica, ya que es una medida de la calidad de las medidas % experimentales %POSIBLES ERRORES % -Si no hay al menos 1 punto experimental a cada lado de la frecuencia % natural, no la puede calcular => Ampliar el intervalo o centrarlo % mejor % -Si no hay a cada lado de la frecuencia natural, como mínimo, tantos % puntos experimentales como el valor que se le haya dado a num_puntos, % no puede calcular el coeficiente de amortiguamiento => Ampliar el % intervalo, centrarlo mejor, o disminuir el valor de num_puntos close all clear all clc format short g %***************************** %CARGO LOS DATOS (ACELERANCIA) %***************************** incr_Hz = 4.88281; A = load ('x8_copia.txt'); Frecuencias = A(:,1); Parte_Real_Acelerancia = A(:,2); Parte_Imag_Acelerancia = A(:,3); %*********************************** %DIBUJO EL NYQUIST DE LA ACELERANCIA %*********************************** figure(1) hold on plot(Parte_Real_Acelerancia,Parte_Imag_Acelerancia) scatter(Parte_Real_Acelerancia, Parte_Imag_Acelerancia) title('Diagrama de Nyquist de la Acelerancia') xlabel('Parte Real Acelerancia') ylabel('Parte Imaginaria Acelerancia') Anexos 82 hold off %******************************** %DIBUJO EL BODE DE LA ACELERANCIA %******************************** [n,m] = size(A); for i=1:1:n Acelerancia_Compleja (i) = Parte_Real_Acelerancia(i) + j*Parte_Imag_Acelerancia(i); Amplitud_Acelerancia(i) = abs(Acelerancia_Compleja(i)); Fase_Acelerancia (i) = angle(Acelerancia_Compleja(i)); end figure(2) amplitud_acel_grafica = subplot(2,1,1); loglog(amplitud_acel_grafica,Frecuencias,Amplitud_Acelerancia) title('Diagrama de Bode de la Acelerancia') xlabel('Frecuencia [Hz]') ylabel('Amplitud de la Acelerancia') fase_acel_grafica = subplot(2,1,2); semilogx(fase_acel_grafica,Frecuencias,Fase_Acelerancia) xlabel('Frecuencia [Hz]') ylabel('Fase de la Acelerancia [rad]') %********************** %CALCULO LA RECEPTANCIA %********************** for i=1:1:n Receptancia_Compleja (i) = Acelerancia_Compleja(i) / -(Frecuencias(i))^2; Amplitud_Receptancia(i) = abs(Receptancia_Compleja(i)); Fase_Receptancia (i) = angle(Receptancia_Compleja(i)); end %******************************** %DIBUJO EL BODE DE LA RECEPTANCIA %******************************** figure(3) amplitud_recep_grafica = subplot(2,1,1); loglog(amplitud_recep_grafica,Frecuencias,Amplitud_Receptancia) title('Diagrama de Bode de la Receptancia') xlabel('Frecuencia [Hz]') ylabel('Amplitud de la Receptancia') fase_recep_grafica = subplot(2,1,2); semilogx(fase_recep_grafica,Frecuencias,Fase_Receptancia) xlabel('Frecuencia [Hz]') ylabel('Fase de la Receptancia [rad]') %************************************************************************* %CALCULO LAS PARTE REAL E IMAGINARIA DE LA RECEPTANCIA Y DIBUJO EL NYQUIST %************************************************************************* for i=1:1:n Parte_Real_Receptancia(i) = Amplitud_Receptancia(i) * cos(Fase_Receptancia(i)); Parte_Imag_Receptancia (i) = Amplitud_Receptancia(i) * sin(Fase_Receptancia(i)); end 83 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental figure(4) hold on plot(Parte_Real_Receptancia,Parte_Imag_Receptancia) scatter(Parte_Real_Receptancia, Parte_Imag_Receptancia) title('Diagrama de Nyquist de la Receptancia') xlabel('Parte Real Receptancia') ylabel('Parte Imaginaria Receptancia') hold off %*************************************************************** %SELECCIONO EL INTERVALO DE FRECUENCIAS EN EL QUE VOY A TRABAJAR %*************************************************************** %Primer Modo ==> 3:8 o 2:8 %Segundo Modo ==> 14:24 %Tercer Modo ==> 47:57 %Cuarto Modo ==> 97:107 %Quinto Modo ==> 163:173 Amplitud_Receptancia = Amplitud_Receptancia(163:173); Fase_Receptancia = Fase_Receptancia(163:173); Frec = Frecuencias(163:173); clear i for j=1:1:length(Frec) Real_Receptancia(j) = Amplitud_Receptancia(j) * cos(Fase_Receptancia(j)); Imag_Receptancia(j) = Amplitud_Receptancia(j) * sin(Fase_Receptancia(j)); Receptancia(j) = Real_Receptancia(j) + i*Imag_Receptancia(j); end %GRAFICA EN FORMA DE BODE DEL INTERVALO QUE SE ESTUDIA figure(5) amplitud_recep_grafica_trozo = subplot(2,1,1); loglog(amplitud_recep_grafica_trozo,Frec,Amplitud_Receptancia) title('Diagrama de Bode de la Receptancia en el tramo de interés') xlabel('Frecuencia [Hz]') ylabel('Amplitud de la Receptancia') fase_recep_grafica_trozo = subplot(2,1,2); semilogx(Frec,Fase_Receptancia) xlabel('Frecuencia [Hz]') ylabel('Fase de la Receptancia [rad]') %GRAFICA EN FORMA DE NYQUIST DEL INTERVALO QUE SE ESTUDIA figure(6) scatter(Real_Receptancia,Imag_Receptancia,'r') hold on plot(Real_Receptancia,Imag_Receptancia,'r') title('Diagrama de Nyquist de la Receptancia en el tramo de interés') axis equal xlabel('Parte Real Receptancia') ylabel('Parte Imaginaria Receptancia') hold off %Uso la funcion CircleByPratt para hallar el circulo que se ajusta mejor a %los puntos Coordenadas=[Real_Receptancia',Imag_Receptancia']; CIRCLE_FIT = CircleFitByPratt (Coordenadas) %Calculo las coordenadas de los puntos del circulo de ajuste Anexos 84 th = 0:pi/500:2*pi; xunit = CIRCLE_FIT(3) * cos(th) + CIRCLE_FIT(1); yunit = CIRCLE_FIT(3) * sin(th) + CIRCLE_FIT(2); %*********************************** %OBTENCION DE LOS PARAMETROS MODALES %*********************************** %********************** %1º) FRECUENCIA NATURAL %********************** %En este apartado se calcula la frecuencia natural como la frecuencia a la %que la respuesta tiene la mayor amplitud. Es la segunda de las 4 formas %que propone el libro de Ewins %Recorro los puntos del circulo de ajuste, y hallo el punto que está mas %alejado del origen (es decir, el punto del circulo con una mayor amplitud %de la respuesta) Coord_Punto_Circ=[xunit',yunit']; for i=1:1:length(xunit) dist(i)=norm(Coord_Punto_Circ(i,:)); end dist; [dist_max,coord_max]=max(dist); Coord_Max=[xunit(coord_max),yunit(coord_max)]; %Coord_Max tiene las componentes del punto DE EL CIRCULO con la mayor %respuesta %Ahora tengo que encontrar los dos puntos experimentales entre los que esta %comprendido el punto en el que se produce la máxima respuesta (que es el %de la frecuencia natural) Coord_Centro_Circ=[CIRCLE_FIT(1),CIRCLE_FIT(2)]; Coord_Punto_Exp=[Real_Receptancia',Imag_Receptancia']; for i=1:1:(length(Real_Receptancia)) Punto1=Coord_Max; Punto2=Coord_Punto_Exp(i,:); Linea1=Punto1-Coord_Centro_Circ; Linea2=Punto2-Coord_Centro_Circ; Modulo1=norm(Linea1); Modulo2=norm(Linea2); producto_esc=Linea1*Linea2'; cos_beta=producto_esc/(Modulo1*Modulo2); beta(i)=acosd(cos_beta); end %Beta es un vector que contiene el angulo que forma la linea que pasa por %el centro del circulo y por el punto de mayor respuesta con la linea que %pasa por el centro del circulo y por cada uno de los puntos %experimentales. Los dos puntos que estoy buscando son los dos que tengan %un beta mas pequeño. [beta_min,Coord_beta_min]=min(beta); Primer_Punto=[Frec(Coord_beta_min),Real_Receptancia(Coord_beta_min),Imag_Rece ptancia(Coord_beta_min),beta_min]; 85 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental if beta(Coord_beta_min+1) >= beta(Coord_beta_min-1) Segundo_Punto=[Frec(Coord_beta_min-1),Real_Receptancia(Coord_beta_min- 1),Imag_Receptancia(Coord_beta_min-1),beta(Coord_beta_min-1)]; end if beta(Coord_beta_min+1) < beta(Coord_beta_min-1) Segundo_Punto=[Frec(Coord_beta_min+1),Real_Receptancia(Coord_beta_min+1),Imag _Receptancia(Coord_beta_min+1),beta(Coord_beta_min+1)]; end Primer_Punto; Segundo_Punto; %Primer_Punto y Segundo_Punto son los dos puntos entre los que esta %comprendido el punto de frecuencia natural. En cada uno de los dos %vectores la primera componente es la frecuencia, la segunda la parte real, %la tercera la parte imaginaria y la cuarta el angulo que forman con el %punto de maxima respuesta %Con esto ya puedo calcular la frecuencia natural if Primer_Punto(1) > Segundo_Punto(1) FRECUENCIA_NATURAL = Segundo_Punto(1) + ((Segundo_Punto(4)*incr_Hz)/(Segundo_Punto(4)+Primer_Punto(4))); end if Primer_Punto(1) < Segundo_Punto(1) FRECUENCIA_NATURAL = Primer_Punto(1) + ((Primer_Punto(4)*incr_Hz)/(Segundo_Punto(4)+Primer_Punto(4))); end FRECUENCIA_NATURAL %********************************** %2º) COEFICIENTE DE AMORTIGUAMIENTO %********************************** num_puntos=3; % Se calculan (num_puntos)^2 valores del Coef de Amort. Aux=[Frec,Real_Receptancia',Imag_Receptancia',beta']; for i=1:1:length(Frec) if Frec(i)<=FRECUENCIA_NATURAL Puntero = i; end end Puntero; Puntos_Coef_Amort=Aux(Puntero-num_puntos+1:Puntero+num_puntos,:); [filas,columnas]=size(Puntos_Coef_Amort); %Calculo los num_puntos^2 valores del coeficiente de amortiguamiento for i=1:1:num_puntos for j=num_puntos+1:1:num_puntos+num_puntos i; j; Num(i,j-num_puntos) = Puntos_Coef_Amort(j,1)^2 - Puntos_Coef_Amort(i,1)^2; Den(i,j-num_puntos) = FRECUENCIA_NATURAL^2 * (tand(Puntos_Coef_Amort(j,4)/2)+tand(Puntos_Coef_Amort(i,4)/2)); Anexos 86 COEFICIENTE_AMORTIGUAMIENTO(i,j-num_puntos) = Num(i,j-num_puntos) / Den(i,j-num_puntos); end end Num; Den; COEFICIENTE_AMORTIGUAMIENTO MEDIA_COEFICIENTE_AMORTIGUAMIENTO = mean(mean(COEFICIENTE_AMORTIGUAMIENTO)) %Calculo la desviacion tipica de todos los valores obtenidos, ya que esta %es un indicador de la calidad de los datos [filas_matriz,columnas_matriz] = size(COEFICIENTE_AMORTIGUAMIENTO); vector_COEFICIENTE_AMORTIGUAMIENTO = reshape(COEFICIENTE_AMORTIGUAMIENTO,[1 filas_matriz*columnas_matriz]); DESVIACION_TIPICA_COEF_AMORT = std(vector_COEFICIENTE_AMORTIGUAMIENTO) %******************************** %3º) MODULO DE LA CONSTANTE MODAL %******************************** Radio_Circulo = CIRCLE_FIT(3); MODULO_CTE_MODAL = (2*Radio_Circulo) * FRECUENCIA_NATURAL^2 * MEDIA_COEFICIENTE_AMORTIGUAMIENTO %*********************************** %4º) ARGUMENTO DE LA CONSTANTE MODAL %*********************************** Linea_1 = Coord_Max; Linea_2 = [0,-5]; Modulo_1 = norm(Linea_1); Modulo_2 = norm(Linea_2); Prod_Esc = Linea_1*Linea_2'; cos_arg_cte_modal = Prod_Esc/(Modulo_1*Modulo_2); ARGUMENTO_CTE_MODAL=acosd(cos_arg_cte_modal); if Coord_Max(1) < 0 ARGUMENTO_CTE_MODAL = - ARGUMENTO_CTE_MODAL; end ARGUMENTO_CTE_MODAL %en grados %****************************** %5º) GRAFICA CON LOS RESULTADOS %****************************** th = 0:pi/500:2*pi; xunit = CIRCLE_FIT(3) * cos(th) + CIRCLE_FIT(1); yunit = CIRCLE_FIT(3) * sin(th) + CIRCLE_FIT(2); figure(7) hold on plot(xunit, yunit,'r:') plot(Real_Receptancia,Imag_Receptancia,'b') scatter(Real_Receptancia,Imag_Receptancia,'xb', 'LineWidth',2) scatter(CIRCLE_FIT(1),CIRCLE_FIT(2),'xr') plot([Coord_Centro_Circ(1) Coord_Max(1)] , [Coord_Centro_Circ(2) Coord_Max(2)],'-r') title('CIRCLE FITTING DE LA RECEPTANCIA') xlabel('Real') ylabel('Imag') axis equal hold off 87 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental 5. LINE_FIT.m %************************************************************************ %APLICACIÓN DEL MÉTODO DE LINE-FIT PARA LA OBTENCIÓN DE PARÁMETROS MODALES %************************************************************************ %LINEAS QUE HAY QUE MODIFICAR % LINEA 28: Hay que poner el nombre del fichero .txt que contiene las % medidas experimentales % LINEAS 121-123: Hay que seleccionar qué intervalo de frecuencias se va % a tomar para estudiar un modo en particular. En las LINEAS 115-119 % vienen los intervalos que se han usado para los 5 primeros modos para % este ejemplo. A diferencia del método del Circle-Fit, aquí es muy % posible que para un mismo modo haya que tomar intervalos diferentes % para las distintas FRFs. En el método de Circle-Fit cerrar o abrir un % poco el intervalo apenas afecta a los resultados, en el método de % Line-Fit sí, así que hay que tener más cuidado. La mejor forma de ver % cómo ir modificando el intervalo para conseguir los mejores resultados % posibles es ir viendo las figuras 8 y 9, y ver si algún punto ''sobra'' close all clear all clc format short g %***************************** %CARGO LOS DATOS (ACELERANCIA) %***************************** A = load ('x8_copia.txt'); Frecuencias = A(:,1); Parte_Real_Acelerancia = A(:,2); Parte_Imag_Acelerancia = A(:,3); %*********************************** %DIBUJO EL NYQUIST DE LA ACELERANCIA %*********************************** figure(1) hold on plot(Parte_Real_Acelerancia,Parte_Imag_Acelerancia) scatter(Parte_Real_Acelerancia, Parte_Imag_Acelerancia) title('Diagrama de Nyquist de la Acelerancia') xlabel('Parte Real Acelerancia') ylabel('Parte Imaginaria Acelerancia') hold off %******************************** %DIBUJO EL BODE DE LA ACELERANCIA %******************************** [n,m] = size(A); for i=1:1:n Acelerancia_Compleja (i) = Parte_Real_Acelerancia(i) + j*Parte_Imag_Acelerancia(i); Amplitud_Acelerancia(i) = abs(Acelerancia_Compleja(i)); Fase_Acelerancia (i) = angle(Acelerancia_Compleja(i)); end Anexos 88 figure(2) amplitud_acel_grafica = subplot(2,1,1); loglog(amplitud_acel_grafica,Frecuencias,Amplitud_Acelerancia) title('Diagrama de Bode de la Acelerancia') xlabel('Frecuencia [Hz]') ylabel('Amplitud de la Acelerancia') fase_acel_grafica = subplot(2,1,2); semilogx(fase_acel_grafica,Frecuencias,Fase_Acelerancia) xlabel('Frecuencia [Hz]') ylabel('Fase de la Acelerancia [rad]') %********************** %CALCULO LA RECEPTANCIA %********************** for i=1:1:n Receptancia_Compleja (i) = Acelerancia_Compleja(i) / -(Frecuencias(i))^2; Amplitud_Receptancia(i) = abs(Receptancia_Compleja(i)); Fase_Receptancia (i) = angle(Receptancia_Compleja(i)); end %******************************** %DIBUJO EL BODE DE LA RECEPTANCIA %******************************** figure(3) amplitud_recep_grafica = subplot(2,1,1); loglog(amplitud_recep_grafica,Frecuencias,Amplitud_Receptancia) title('Diagrama de Bode de la Receptancia') xlabel('Frecuencia [Hz]') ylabel('Amplitud de la Receptancia') fase_recep_grafica = subplot(2,1,2); semilogx(fase_recep_grafica,Frecuencias,Fase_Receptancia) xlabel('Frecuencia [Hz]') ylabel('Fase de la Receptancia [rad]') %************************************************************************* %CALCULO LAS PARTE REAL E IMAGINARIA DE LA RECEPTANCIA Y DIBUJO EL NYQUIST %************************************************************************* for i=1:1:n Parte_Real_Receptancia(i) = Amplitud_Receptancia(i) * cos(Fase_Receptancia(i)); Parte_Imag_Receptancia (i) = Amplitud_Receptancia(i) * sin(Fase_Receptancia(i)); end figure(4) hold on plot(Parte_Real_Receptancia,Parte_Imag_Receptancia) scatter(Parte_Real_Receptancia, Parte_Imag_Receptancia) title('Diagrama de Nyquist de la Receptancia') xlabel('Parte Real Receptancia') ylabel('Parte Imaginaria Receptancia') hold off %*************************************************************** %SELECCIONO EL INTERVALO DE FRECUENCIAS EN EL QUE VOY A TRABAJAR %*************************************************************** 89 Implementación en Matlab de Métodos de 1 g.d.l. Para Realizar Análisis Modal Experimental %Primer Modo ==> 3-4 : 8 %Segundo Modo ==> 16-18 : 22-24 %Tercer Modo ==> 48-50 : 53-57 %Cuarto Modo ==> 97-99 : 103-107 %Quinto Modo ==> 163-165 : 170-173 Amplitud_Receptancia = Amplitud_Receptancia(163:173); Fase_Receptancia = Fase_Receptancia(163:173); Frec = Frecuencias(163:173); clear i for j=1:1:length(Frec) Real_Receptancia(j) = Amplitud_Receptancia(j) * cos(Fase_Receptancia(j)); Imag_Receptancia(j) = Amplitud_Receptancia(j) * sin(Fase_Receptancia(j)); Receptancia(j) = Real_Receptancia(j) + i*Imag_Receptancia(j); end figure(5) amplitud_recep_grafica_trozo = subplot(2,1,1); loglog(amplitud_recep_grafica_trozo,Frec,Amplitud_Receptancia) title('Diagrama de Bode de la Receptancia en el tramo de interés') xlabel('Frecuencia [Hz]') ylabel('Amplitud de la Receptancia') fase_recep_grafica_trozo = subplot(2,1,2); semilogx(fase_recep_grafica_trozo,Frec,Fase_Receptancia) xlabel('Frecuencia [Hz]') ylabel('Fase de la Receptancia [rad]') %************************************************************************** %CALCULO DE LA NUEVA FRF (alfa_prima), DE LA INVERSA DE LA FRF (Delta) Y DE %LA PARTE REAL E IMAGINARIA DE DELTA. HAY QUE HACERLO UNA VEZ POR CADA %FRECUENCIA QUE SE FIJE, POR ESO EL BUCLE DE j=1:1:length(Frec) %************************************************************************** for j=1:1:length(Frec) Componente_Frec_Fijada = j; Frec_Fijada(j) = Frec(Componente_Frec_Fijada); %********************************************************************** %IMPORTANTE: CUANDO LA FRECUENCIA VALE LO MISMO QUE LA FRECUENCIA QUE %HE FIJADO, DELTA RESULTA NaN (cero partido por cero), POR ESO DEFINO %UN VECTOR FREC_AUX QUE CONTIENE TODAS LAS FRECUENCIAS MENOS LA %FRECUENCIA QUE HE FIJADO. A ESTAS FRECUENCIAS SON A LAS QUE LE CALCULO %ALFA_PRIMA Y LUEGO DELTA %********************************************************************** for i=1:1:Componente_Frec_Fijada-1 Frec_Aux_1(i) = Frec(i); end for i=Componente_Frec_Fijada+1:1:length(Frec) Frec_Aux_2(i-Componente_Frec_Fijada) = Frec(i); end if Componente_Frec_Fijada == 1 Frec_Aux(j,:) = Frec_Aux_2; end if Componente_Frec_Fijada == length(Frec) Frec_Aux(j,:) = Frec_Aux_1;