Identificación Experimental del Comportamiento Dinámico de Estructuras Aeronáuticas mediante Sensores de Bajo Coste

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Identificación Experimental del Comportamiento Dinámico de Estructuras Aeronáuticas mediante Sensores de Bajo Coste Autor: Carlos Alonso Vélez Tutor: Javier Fernando Jiménez Alonso Dpto. de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Identificación Experimental del Comportamiento Dinámico de Estructuras Aeronáuticas mediante Sensores de Bajo Coste Autor: Carlos Alonso Vélez Tutor: Javier Fernando Jiménez Alonso Profesor Titular de Universidad Dpto. de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 Trabajo Fin de Grado: Identificación Experimental del Comportamiento Dinámico de Estructuras Aeronáuticas mediante Sensores de Bajo Coste Autor: Carlos Alonso Vélez Tutor: Javier Fernando Jiménez Alonso El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos A mi familia, amigos y pareja, que con su apoyo y cariño han hecho los años de carrera mucho más amenos A mi tutor Javier Jiménez por estar siempre dispuesto a ayudar y haberme enseñado tanto durante estos meses I Índice de Figuras 1.1 GVT realizado en un vehículo aéreo no tripulado 2 1.2 Diagrama esquemático de los pasos que componen la Identificación Estructural 4 1.3 Esquema de problemas básicos de la dinámica estructural 6 1.4 Esquema de la organización por capítulos del documento 8 2.1 Diagrama de técnicas de análisis de prototipos 9 2.2 Imagen de un modelo a escala reducida de un UAV en túnel de viento 10 2.3 Diagrama conceptual de la predicción del comportamiento de un prototipo a través del de un modelo escalado 11 3.1 Esquema de los componentes de un sistema de adquisición de datos 17 3.2 Esquema del modelo simplificado de un acelerómetro 18 3.3 Esquema del funcionamiento y las partes de un acelerómetro MEMS 20 3.4 Placa Arduino Uno R3 22 3.5 Esquema de la arquitectura de un sistema de adquisición de datos para identificación modal 23 3.6 Diagrama de flujo de la utilización del timer para la lectura de las medidas en Arduino 24 3.7 Ejemplo de efecto aliasing en la reconstrucción de una señal tras digitalización 27 3.8 Efecto del leakage en una señal periódica de 50 Hz ( fs=1 kHz) y corrección mediante una ventana 28 3.9 Distorsión generada en la función PSD por efecto del leakage y del aliasing 29 3.10 Ejemplos de arquitecturas para sincronizar canales de datos del ensayo. Se muestran con línea fina los cables digitales y con línea gruesa los analógicos 33 4.1 Representación esquemática de la identificación modal basada en respuesta 37 4.2 Principales tipos de métodos OMA en el dominio de la frecuencia, y algunos ejemplos significativos 45 4.3 Principales tipos de métodos OMA en el dominio del tiempo, y algunos ejemplos significativos 46 4.4 Esquema del proceso de estimación de la función PSD a través del método de Welch 48 4.5 Valores singulares de la matriz PSD para un sistema con dos modos de vibración 50 4.6 Valores singulares de la matriz PSD de respuesta para tres mediciones con 1 y 2 modos 50 4.7 Esquema general del algoritmo FDD para el caso de dos canales de medidas 52 4.8 Ejemplo de diagrama de estabilización en la aplicación del método SSI-COV 56 4.9 Esquema general del algoritmo SSI-COV para el caso de dos canales de medidas 58 4.10 Esquema general del algoritmo p-LSCF para el caso de dos canales de medidas 65 5.1 Diagrama de las fases de la implementación de un filtro de Kalman 70 5.2 Diagrama de la propagación de media y covarianza, en caso de muestreo real y en caso de UT 72 III 2Capítulo 1. Introducción Los GVT tienen como principal objetivo actualizar los modelos matemáticos de la aeronave para mejorar la calidad de las predicciones de flameo de cara a los ensayos en vuelo. Sin embargo, dado que estos test no tienen en cuenta la interacción entre la estructura y el fluido, son necesarios ensayos en túnel de viento y en vuelo para identificar los modos aeroelásticos. Una descripción del estado del arte y el futuro de los GVT en aeronaves de gran tamaño se puede encontrar en [1], mientras que una aplicación a aeronaves flexibles se estudia en [2]. Por su parte, en [3] se lleva a cabo una identificación modal en túnel de viento. Por último, ejemplos de investigación en test en vuelo son [4–6]. Figura 1.1 GVT realizado en un vehículo aéreo no tripulado. Tomada de [2]. Además de la caracterización del flameo, existen otras aplicaciones de la identificación modal en estructuras aeronáuticas. Una de ellas es calibrar y actualizar los modelos (por ejemplo de elementos finitos) de las mismas, como se explicará en el siguiente apartado. A su vez, la obtención de las características modales en servicio permite que se cumplen los requisitos de seguridad e integridad en cuanto a estados límite de vibraciones, o detectar y evaluar daños (Structural Health Monitoring). Otra aplicación es el control activo y semi-activo de estructuras, mediante el empleo de dispositivos que disipen energía, como los amortiguadores de masa activa. Existen casos en los que, antes de la producción de un diseño, resulta complicado generar un prototipo para experimentación, o simplemente no resulta rentable en cuanto dedicación de recursos. Sin embargo, la creación de modelos y su verificación mediante experimentación es un proceso esencial en cualquier fase de diseño. En los casos mencionados, resulta de gran utilidad ensayar prototipos a escala reducida, en los que se realizan simplificaciones. Dentro de las aplicaciones de la monitorización para predecir la salud estructural o Structural Health Monitoring, una de las principales limitaciones suele ser el alto coste de los sensores empleados para medir la respuesta estructural [7]. Sin embargo, en los últimos años se ha empezado a introducir en estas aplicaciones el uso de sensores basados en electrónica de bajo coste, como, por ejemplo, los acelerómetros micro electro-mecánicos o sus siglas en inglés MEMS (Micro Electro-Mechanical Systems), los cuales son mucho más asequibles y están sufriendo continuos desarrollos y evoluciones tecnológicas. 1.3 Fundamentos de la identificación estructural 3 1.3 Fundamentos de la identificación estructural La base fundamental del presente trabajo es la identificación experimental de parámetros estructurales, o Structural Identification (St-Id). Este concepto fue introducido en 1977 por G. Hart y J. Yao [8], y es una división de la identificación de sistemas aplicada a estructuras tanto civiles como mecánicas. Consiste en la correlación paramétrica entre la respuesta estructural esperada mediante predicciones realizadas con un modelo, y la respuesta estructural medida experimentalmente. El objetivo fundamental de dicha comparación reside en la aproximar el modelo al sistema real para mejorar las estimaciones sobre rendimiento y puntos débiles del mismo. Los modelos se pueden clasificar en los siguientes tipos [9]: 1. Modelos físicos o physics-based (PB). Se pueden considerar predictivos y útiles a priori, en el sentido de que, dado que están basados en leyes generalizadas de la estática, mecánica, dinámica, etc., no requieren datos de la respuesta del sistema real. Pueden ser continuos o discretos. • Modelos continuos: teoría de la elasticidad, ecuaciones de campo y de ondas, modelos de la mecánica de materiales (Bernoulli, Vlasov), etc. • Modelos discretos: conjuntos resorte-masa-amortiguador, coeficientes K, M y C, modelos de elementos finitos (FEM), modelos modales (propiedades modales, vectores Ritz), etc. 2. Modelos no-físicos o non-physics-based (NPB). Son de naturaleza descriptiva. En lugar de estar basados en leyes físicas, se derivan a partir de datos de la respuesta del sistema, a través de métodos de modelado, interpretación y reducción. Por ello, no se pueden emplear a priori, aunque, en algunos casos, el entrenamiento de estos modelos puede llevar a considerarlos predictivos al ser capaces de estimar la respuesta futura identificando patrones en el sistema. Se pueden dividir en semánticos, probabilísiticos y numéricos. Debido a las características mencionadas anteriormente acerca de la empleabilidad de los modelos a priori, en este trabajo se va a tratar esencialmente con modelos físicos. La identificación estructural permite crear uno de ellos e ir actualizándolo periódicamente con la adquisición de datos. El comité de St-Id de la ASCE (American Society of Civil Engineers) llegó a un consenso sobre los pasos básicos a integrar para lograr satisfactoriamente la identificación de una estructura en operación [9 – 11]. Estos pasos presentan cierta flexibilidad, dependiendo del sistema estudiado, y vienen representados esquemáticamente en la Figura 1.2. A continuación, se proporciona una explicación detallada de cada una de dichas fases: 1. Justificación del problema. Consiste en establecer las bases y los objetivos de la identificación, y justificar la necesidad de proceder con el proyecto. Es importante familiarizarse con la estructura a analizar e investigar acerca del conocimiento previo disponible sobre la misma. Otro aspecto a destacar en este paso es la identificación de posibles fuentes de incertidumbre en su comportamiento, de cara a evitar un modelado impreciso que afectará al resto de pasos del proceso. 2. Creación de modelos a priori. En esta fase, se estudia la información previa disponible; se capturan tantas dimensiones, propiedades mecánicas y características estructurales como sea posible, y se observa su abanico de condiciones de carga, tanto en operación como ambientales, con el objetivo de generar un 4Capítulo 1. Introducción Identificación Estructural Justificación del problema, objetivos 1 Dimensionamiento de la estructura, modelos a priori 2 Experimentación controlada, monitorización 3 Procesamiento e interpretación de los datos 4 Selección, calibración y validación de modelos 5 Toma de decisiones a través del modelo calibrado 6 Figura 1.2 Diagrama esquemático de los pasos que componen la Identificación Estructural. Editada a partir de [11]. modelo físico a priori con el máximo criterio. Cuanto mayor es la cantidad de información recogida e interpretada preliminarmente, mayor será la reducción de incertidumbre y errores en el modelado y los pasos siguientes. Las distintas aproximaciones que se pueden tomar a la creación del modelo dependen principalmente de los objetivos del proceso de identificación y de la complejidad de la estructura. En algunos casos, el modelo a priori se mantiene durante el proceso de actualización con la identificación de los parámetros (paso 5). 3. Experimentación controlada Los métodos experimentales aplicados para capturar el comportamiento del sistema real sirven como enlace entre el modelo y dicho sistema. Pueden ser de distintos tipos: de excitación natural, de excitación forzada, de carga controlada, y monitorización de eventos operacionales y ambientales. La experimentación comprende la elección del caso de carga; la instrumentación; la adquisición, preprocesado, trasmisión y almacenamiento de datos, y la documentación general del proceso. Una buena práctica durante este paso puede reducir en gran medida la incertidumbre en el procesado e interpretación de los datos posteriores. 4. Procesamiento e interpretación de los datos. La primera de las actividades, el procesamiento, consiste en conseguir que el conjunto de datos extraídos de los experimentos sea adecuado para su interpretación. Para ello, se minimizan errores debidos a picos en las señales, fallos en los sensores, etc., y se mejora la calidad de los datos, a través de filtros, ventanas (windowing). En algunos casos, también se precisa de una transformación de las señales para facilitar la interpretación, por ejemplo, a través de la extracción de propiedades modales o coeficientes de influencia. Además, dicho procesamiento puede ir acompañado de una interpretación directa de los datos, con el objetivo de ajustar modelos no-físicos, tales como redes neuronales, modelos auto-regresivos, modelos de espacio de estados, etc. Este tipo de modelos detectan cambios en 1.3 Fundamentos de la identificación estructural 5 el sistema estableciendo umbrales de lo que considerar normalidad en los datos, sin capturar la explicación física detrás de dicho cambio. Por este motivo, requieren poca interacción con el usuario, siendo ventajosos para monitorización continua. En cambio, presentan dificultad para detectar variaciones lentas y prolongadas en el tiempo, provocadas, por ejemplo, por el ambiente, ya que la no es posible definir umbrales que excluyan dichos cambios tan leves a corto plazo. 5. Selección, calibración y validación de modelos físicos. En esta fase, se realiza una búsqueda del origen físico del cambio en el sistema, calibrando o actualizando (model updating) el modelo. Generalmente, se selecciona un conjunto de parámetros del modelo para posteriormente iniciar un proceso de optimización en el que minimizar la diferencia entre las predicciones del modelo y los resultados experimentales. Los enfoques de calibración se pueden clasificar, entre otros, según cómo se realice la selección de parámetros, la formulación de las funciones objetivo a minimizar, o el método de optimización que emplean (basado en gradiente o no). En el caso de haber realizado una interpretación directa de los datos mediante el ajuste de modelos no-físicos, este paso se puede ver como opcional. Sin embargo, es siempre recomendable realizar una labor de investigación sobre la explicación física detrás del cambio y la influencia que puede tener en la operación de la estructura, así como el empleo de varios tipos de modelos para comparar y elegir el de mejores resultados. 6. Toma de decisiones a través del modelo calibrado. Para justificar la aplicación desde el punto de vista de los recursos, tanto económicos como humanos y temporales, es necesario emplear los modelos, calibrados en el caso de físicos y entrenados en el caso de los no-físicos, para tomar decisiones. Estas decisiones pueden estar enfocadas a mejorar la operación del diseño en ciertas condiciones, introducir mejoras y actualizaciones al sistema, evaluar la seguridad y la integridad de la estructura, o predecir el comportamiento ante condiciones futuras o hipotéticas. Para ello, se hace uso del modelo validado para realizar simulaciones o estudios paramétricos. 1.3.1 Identificación de propiedades modales Como se nombró anteriormente, el presente trabajo se centra en la identificación de las propiedades modales de una estructura, como son las frecuencias naturales, los modos de vibración y los ratios de amortiguamiento modales de la misma. Para ello, se cubren los pasos 1 a 4 de la secuencia explicada. Se tratarán con especial detalle los métodos empleados para el procesamiento de los datos. La identificación modal emplea datos de mediciones de las vibraciones del sistema (por ejemplo, aceleraciones). En el problema dinámico directo, se parte de unas propiedades estructurales especificadas, a partir de las cuales se calculan resolviendo un problema de autovalores las propiedades modales. Estas últimas permiten, a su vez, obtener una predicción de la respuesta estructural mediante la superposición de la respuesta de los distintos modos de vibración. Sin embargo, tal y como se puede apreciar en la Figura 1.3, el problema también permite ser recorrido en sentido inverso, a través de la identificación modal y estructural. De esta manera, la identificación modal presenta significativas ventajas, ya que constituye una manera tanto de calibrar las propiedades de la estructura, como de predecir de manera aceptable la respuesta ante situaciones futuras. Dicha predicción se puede realizar sin conocer entonces las 6Capítulo 1. Introducción Propiedades estructurales Masa, amortiguamiento, rigidez, etc. Propiedades modales Frecuencias naturales, modos de vibración, ratios de amortiguamiento, masas modales Respuesta estructural Desplazamientos, velocidades, aceleraciones, tensiones, etc. Prob. de autovalores Identif. estructural Superposición modal Identif. modal + modelo de amortiguamiento Fuerzas nodales Fuerzas modales Modo 1 Modo 2 Freq. natural Ratio de amort. Masa modal Freq. natural Ratio de amort. Masa modal + Figura 1.3 Esquema de problemas básicos de la dinámica estructural. Editada a partir de [12]. propiedades estructurales del sistema. Además, es posible identificar únicamente las propiedades modales que influyen en mayor medida en la dinámica de la estructura simplificando en gran medida el problema. Esta selección de los parámetros más determinantes en el problema no es generalmente posible en el caso de las propiedades estructurales. Se van a estudiar dos tipos de identificación de propiedades modales: • Identificación offline. Esta consiste en la obtención de los propiedades modales de la estructura a partir del conjunto de datos correspondientes a la respuesta dinámica de la estructura ante una determinada excitación en un periodo prolongado de tiempo. Dependiendo de la naturaleza de la carga aplicada, se utilizan métodos como el análisis modal experimental o EMA (Experimental Modal Analysis) en el caso de excitación forzada, o el análisis modal operacional u OMA (Operational Modal Analysis) en el caso de excitación ambiental. La particular de estos métodos respecto al otro tipo de identificación reside en que esta requiere el procesado de los datos experimentales a posteriori, una vez se han almacenado registros de un ensayo completo. Debido a esto, la obtención de las propiedades modales se realiza de manera periódica. También se llaman métodos no-recursivos [13]. • Identificación online. Consiste en la implementación de un observador adaptativo que permite la identificación conjunta del estado del sistema y de las propiedades modales de la estructura, mediante el filtro de Kalman en la versión unscented o sus siglas en inglés UKF (Unscented Kalman Filter). Esta técnica supone una identificación en tiempo real, sin necesidad de esperar a que finalice el ensayo. Esto permite su uso con sistemas de control de vibraciones semi-activos. También se llaman métodos recursivos. 1.4 Objetivos y contribuciones del trabajo Los objetivos perseguidos con la realización de este trabajo son los siguientes: 1.5 Organización del texto 7 • Aproximación a la identificación de las propiedades modales de estructuras aeronáuticas mediante métodos basados en la respuesta del sistema. • Empleo de un modelo a escala reducida para simular el comportamiento dinámico de un prototipo real. • Definición de un sistema de adquisición de datos de bajo coste adecuado para la identificación modal. •Configuración de un ensayo dinámico de una estructura. • Comparación de diferentes técnicas de identificación de propiedades modales, tanto offline como en tiempo real. •Análisis crítico de los resultados. Por otro lado, las principales contribuciones del mismo son: • Uso de modelos a escala reducida para identificar parámetros de prototipos cuyo ensayo es de gran dificultad o coste. • Utilización de un sistema de adquisición de bajo coste para obtener una estimación aceptable, mediante distintos métodos, de las propiedades modales de una estructura. • Empleo de diferentes técnicas offline de Análisis Modal Operacional (OMA) para la identificación de estructuras aeronáuticas con sensores de bajo coste. • Estimación de propiedades modales de forma continua y en tiempo real mediante filtros de Kalman unscented, abriendo paso a la creación de modelos digitales o digital twins que pueden simular y predecir el comportamiento del sistema y actuar de acuerdo a ello. 1.5 Organización del texto El presente documento se articula en siete capítulos (ver Figura 1.4). El primero de ellos constituye la introducción del tema, los cuatro siguientes son capítulos teóricos donde se tratan las bases necesarias para la identificación modal en modelos a escala reducida. El Capítulo 6 está dedicado a una aplicación práctica a un caso real, mientras que las conclusiones del trabajo se exponen en el último capítulo. Para finalizar, se incluye una bibliografía. Particularmente, el contenido de cada capítulo es el siguiente: • Capítulo 1: se introduce el tema de la identificación estructural y, más concretamente, la de propiedades modales, haciendo especial énfasis en las aplicaciones en la industria aeronáutica. • Capítulo 2: se exponen los fundamentos de la teoría de la similitud y el análisis dimensional para la creación de modelos a escala reducida. • Capítulo 3: se describen los componentes de un sistema de adquisición de datos y la configuración de un ensayo dinámico, además del preprocesado de datos para adecuarlos al análisis modal posterior. • Capítulo 4: se centra en el Análisis Modal Operacional, sus fundamentos y tipos, y se profundiza en tres técnicas, la FDD, el SSI-COV y el p-LSCF. • Capítulo 5: se debate la identificación online de propiedades de un sistema mediante una versión modificada del filtro de Kalman unscented para estimación de parámetros en tiempo real. 8Capítulo 1. Introducción • Capítulo 6: se aplica la identificación de propiedades modales a un prototipo de ala real, exponiendo y analizando los resultados. •Capítulo 7: se enuncian las conclusiones del trabajo y posibles líneas futuras de desarrollo. Marco teórico Capítulo 1 Introducción Capítulo 2 Modelos a escala reducida Capítulo 3 Adquisición de datos Capítulo 4 Identificación offline Capítulo 5 Identificación online Capítulo 6 Aplicación práctica Capítulo 7 Conclusiones, líneas futuras Figura 1.4 Esquema de la organización por capítulos del documento. 2 Modelos a escala reducida 2.1 Introducción En el Capítulo 1 se ha discutido acerca de la importancia del modelado de estructuras. En particular, durante el proceso de diseño y antes de la producción de un nuevo sistema, son necesarias predicciones sobre su comportamiento mediante análisis teóricos, así como la verificación de dichos comportamientos mediante ensayos experimentales [14]. Esto último se hace especialmente indispensable si la complejidad de la estructura impide generar modelos matemáticos adecuados para simulación. Sin embargo, debido a restricciones como el tamaño del diseño, en algunas ocasiones no es posible o no es rentable destinar los recursos económicos, los recursos humanos o el tiempo necesario a la fabricación de un prototipo real para emplear en experimentación. Observación directa Métodos experimentales Prototipo real Suposiciones Prototipo idealizado basados en: Métodos de simulación Experimentación con modelo a escala reducida Modelos matemáticos Figura 2.1 Diagrama de técnicas de análisis de prototipos. Tomada de [15]. Una solución para este tipo de casos es la creación de un modelo a escala reducida 1 , también llamado modelo similar, el cual se usa para simular el comportamiento del prototipo. Esta técnica presenta una considerable simplicidad en comparación con la magnitud del problema real, y frecuentemente conlleva menos suposiciones que el empleo de otros métodos como modelos matemáticos, dado que estos requieren conocimiento preciso de las ecuaciones que gobiernan el comportamiento 1 En adelante, a menos que se indique expresamente lo contrario, la palabra modelo se va a emplear para referirse a artefactos físicos, y no se debe confundir con otros tipos de modelos, como los matemáticos. 9 10 Capítulo 2. Modelos a escala reducida del sistema. En el esquema de la Figura 2.1 se muestran los distintos enfoques que existen para el análisis de un diseño a través de técnicas experimentales o de simulación. La Teoría de similitud se encarga de establecer las condiciones de semejanza que se han de dar entre el modelo real y el modelo a escala reducida. Los objetivos que se persiguen con el uso de estos modelos escalados son, entre otros [16]: • Obtener datos experimentales para realizar una evaluación cuantitativa de un análisis teórico. •Observar el comportamiento producido por un fenómeno poco conocido. • Obtener datos cuantitativos para usarlos en problemas de diseño de prototipos, especialmente cuando el modelado matemático de su comportamiento es complicado o imposible. •Generar una relación funcional empíricamente para resolver un problema general. •Evaluar las limitaciones de un prototipo de alto coste ya existente. Algunas de las principales aplicaciones de estos artilugios se dan en la ingeniería civil, automovilística, naval y aeroespacial [15]. Una de ellas es su empleo en túneles de viento, donde es posible generar condiciones similares a las de operación de manera económica, tanto en el caso de coches [17] como de aeronaves [18] (ver Figura 2.2) y edificios, por ejemplo, en el análisis del rascacielos más alto del mundo, la torre Burj Khalifa [19]. En la industria aeroespacial, otra muestra de uso de modelos a escala reducida es el proyecto de investigación de la NASA del UAV X-48, donde se emplean dichos modelos en vuelos de prueba 2 . Son considerados por ingenieros aeronáuticos de dicha organización como una de las herramientas de diseño más viables y valiosas desde la llegada de los vuelos [20]. Figura 2.2 Imagen de un modelo a escala reducida de un UAV en túnel de viento. Fuente: J. Holsten et al. [18]. 2https:// www.nasa.gov/ centers/ armstrong/ news/ FactSheets/ FS-090-DFRC.html 2.2 Teoría de la similitud 11 2.2 Teoría de la similitud 2.2.1 Fundamentos. Tipos de similitud La teoría de la similitud se encarga de establecer las condiciones necesarias y suficientes de semejanza entre dos fenómenos. Estos criterios de similitud son de gran ayuda en el estudio de nuevos sistemas, ya que permiten: •Predecir el comportamiento de un sistema investigando otros que ya han sido estudiados. •Estimar su comportamiento a través del de otro sistema cuyo estudio es más sencillo. La similitud entre sistemas implica, a su vez, semejanza en el comportamiento de ambos. Es por ello que conociendo como responde un determinado sistema ante una entrada, es posible predecir la respuesta de sistemas similares ante entradas semejantes [21]. Esta relación se muestra conceptualmente en la Figura 2.3. Prototipo Comportamiento estructural Modelo escalado Comportamiento estructural Geometría, Material, Cargas Condiciones de similitud Leyes de escalado Geometría, Material, Cargas Figura 2.3 Diagrama conceptual de la predicción del comportamiento de un prototipo a través del de un modelo escalado. Tomado de [14]. La naturaleza de un sistema puede expresarse a través de un modelo matemático dependiente de sus variables y parámetros como son la geometría, las propiedades del material, la respuesta dinámica o la energía característica del sistema. Un prototipo y su modelo escalado tienen distintos parámetros, y las condiciones necesarias y suficientes de similitud requieren que el modelo matemático del modelo escalado pueda ser relacionado con el del prototipo a través de una transformación biunívoca y viceversa [22]. Esto supone que, siendo Xp el vector de características del prototipo y Xm el vector de características del modelo, existe una matriz de transformación Λtal que [14,15,21]: Xp=ΛXmóXm=Λ−1Xp(2.1) El vector de características es aquel cuyas componentes son las variables y parámetros del sistema. El caso más simple, aunque generalmente no sucede, en el que la matriz Λ es una matriz diagonal donde las componentes, λxi son los factores de escala de las correspondientes componentes del vector X: 18 Capítulo 3. Adquisición de datos 3.2 Sensores Como se ha comentado, el papel de los transductores de movimiento es la conversión de desplazamiento, velocidad o aceleración en una magnitud eléctrica, generalmente voltaje, que sea proporcional a la magnitud de origen. En cuanto a la medición de la respuesta dinámica de estructuras, este trabajo se va a centrar en los acelerómetros. Estos sensores, los más comunes en la aplicación tratada, presentan la ventaja de tener unos principios de funcionamiento básicos y menos complejos que los sensores de desplazamiento o de velocidad, por lo que conllevan, generalmente, un menor coste. Sus dos principales tipos son los piezoeléctricos y los acelerómetros Force-Balance o de equilibrio de fuerzas [24]. 3.2.1 Principios generales de funcionamiento A continuación, se va a describir el principio básico de funcionamiento de un acelerómetro, comentando las particularidades propias tanto del transductor piezoeléctrico como del de equilibrio de fuerzas. Un sensor se puede generalmente modelar como un sistema de una masa con un resorte, como se puede observar en la Figura 3.2. Punto de reposo Superficie estructural Figura 3.2 Esquema del modelo simplificado de un acelerómetro. Tomada de [25]. Se asume que la carcasa del acelerómetro (en gris en la figura) está unida de forma rígida a la superficie de la estructura. Se ha denominado x(t) al movimiento base, el de dicha unión respecto a cierto sistema de referencia. Por otro lado, y(t) se refiere al movimiento relativo de la masa interna del sensor, llamada comúnmente masa sísmica, respecto al punto de reposo en la carcasa. El modelo planteado constituye un sistema de un grado de libertad, en el que la masa sísmica M está unida a un resorte de rigidez k que ejerce una fuerza F=ky sobre ella cuando su desplazamiento relativo es y. Dadas las descripciones anteriores, la frecuencia natural de dicho sistema es: fn=1 2πrk M(3.1) Cuando el conjunto estructura-carcasa se mueve en la dirección x , la aceleración a=¨x(t) induce una fuerza de inercia de valor F=M a sobre la masa sísmica. Asumiendo que el amortiguamiento del sistema es despreciable, por equilibrio de fuerzas se obtiene que el desplazamiento de la masa, y por tanto la deformación del resorte, es proporcional a la aceleración de carcasa. 3.2 Sensores 19 Finercia =M a Fresorte =k y )Finercia =Fresorte ⇒M a =k y ⇒y=M ka(3.2) Como la carcasa del sensor está unida rígidamente a la superficie de la estructura en un punto, su aceleración coincide con la de la estructura en dicho punto [25]. Siendo el descrito el funcionamiento de un acelerómetro genérico, se van a exponer ahora brevemente algunas características propias del acelerómetro piezoeléctrico y el Force-Balance. • Acelerómetro piezoeléctrico. Se basa en el efecto del mismo nombre de algunos materiales como el cuarzo para operar. Cuando se aplica una fuerza sobre el cristal, iones positivos y negativos se acumulan en superficies opuestas del mismo, generando una carga eléctrica que es proporcional a la deformación que sufre. En estos acelerómetros, la masa sísmica está unida al material piezoeléctrico, por lo que la deformación que induce esta sobre él es proporcional a la aceleración de la estructura. • Acelerómetro Force-Balance o de equilibrio de fuerzas. Están diseñados con el objetivo de mantener la masa sísmica en la posición de reposo, es decir, el resorte indeformado. Para ello, la electrónica del dispositivo ajusta continuamente la fuerza del resorte para compensar la fuerza de inercia de la masa. Esta es la fuerza que mide directamente el sensor, y dado que la masa sísmica es conocida, se puede deducir el valor de la aceleración. Además de los dos tipos comentados, otros sensores también utilizados en identificación modal son los electromagnéticos, como el sismógrafo y el geófono. En la actualidad, se está extendiendo el empleo de acelerómetros de bajo coste MEMS, en los que este trabajo se va a centrar y que se desarrollan en la Subsección 3.2.3. 3.2.2 Parámetros básicos de un acelerómetro Las principales características del sensor a tener en cuenta a la hora de llevar a cabo un ensayo de identificación modal son la sensibilidad, la resolución y el rango dinámico [24,25]. La sensibilidad del sensor suele venir dada como la ganancia del mismo, por ejemplo, S=5 V/g. Se define a través del mínimo cambio en aceleración que es capaz de provocar un variación detectable en el voltaje de salida. Dado que ha de ser detectable, esa mínima señal está también limitada por el ruido inducido por la electrónica, y se denomina resolución del sensor, δV0 . Cabe mencionar que el valor de sensibilidad como una constante es una aproximación, dado que en realidad un acelerómetro es un sistema dinámico cuya relación lineal entrada-salida se puede expresar en el dominio de la frecuencia como: ˜ V(f) = ˜ S(f)˜a(f)(3.3) donde la sensibilidad se puede interpretar como la función de respuesta en frecuencia o FRF (Frequency Response Function)H(f) = ˜ S(f) = ˜ V(f)/˜a(f)del sensor. El rango dinámico o sus siglas en inglés DR (Dynamic Range) del acelerómetro es el ratio entre el valor más pequeño δV0y más alto Vmque se puede medir: 20 Capítulo 3. Adquisición de datos DR =20 logVm δV0(en dB) (3.4) 3.2.3 Acelerómetros MEMS (Micro Electro-Mechanical Systems) Los sistemas micro electro-mecánicos o sus siglas en inglés MEMS (Micro Electro-Mechanical Systems) son dispositivos de pequeño tamaño, entre 1 µm y 1 mm, basados en silicio que integran generalmente tanto un sensor acelerómetro y giróscopo como un circuito para acondicionamiento de señal [26]. Los MEMS pueden utilizar distintos principios de operación, como el efecto piezoeléctrico o piezorresistivo [27]. Sin embargo, los más comunes son los capacitivos, sensores de equilibrio de fuerzas que emplean un condensador con placas móviles y sistema de control. Dicho sistema produce una fuerza compensadora cuando las placas se desplazan debido a la aceleración de la estructura. En reposo, las placas móviles, que constituyen la masa sísmica del sensor, son simétricas respecto a la placa fija central. Cuando se induce un desplazamiento relativo entre ellas debido a una aceleración, se genera una señal desequilibrada que se compensa a través de una fuerza magnética por realimentación con el objetivo de que las placas capacitivas vuelvan a la posición original. La corriente que circula por la bobina para generar la fuerza es proporcional a la aceleración experimentada por el punto de la estructura [24]. Las distintas partes básicas de este tipo de sensor, así como su funcionamiento, se pueden observar en la Figura 3.3. Masa sísmica REPOSO Electrodo flexible Electrodo fijo Resorte Sustrato ACELERACIÓN Figura 3.3 Esquema del funcionamiento y las partes de un acelerómetro MEMS. Los sensores inerciales MEMS juegan un papel fundamental en la expansión de los dispositivos electrónicos personales, estando presentes, entre muchos otros, en los smartphones. Además, su continua evolución tecnológica ha propiciado su empleo en una gran cantidad de aplicaciones industriales [26], debido a características como las siguientes: • Compactibilidad y ligereza. Como se ha comentado anteriormente, la tecnología de fabricación a escala microscópica permite a los MEMS tener un tamaño extremadamente pequeño y un bajo peso, lo que supone grandes beneficios para aplicaciones donde estos parámetros son críticos, como UAVs. 3.2 Sensores 21 Tabla 3.1 Tabla comparativa de acelerómetros comúnmente usados en la literatura. Tomada de [28]. Nombre Precio (€) Rango de aceleraciones (g) Rango de frecuencias (Hz) Ruido espectral (µg/√Hz) Temperatura de operación (ºC) Tipo de sensor 3713B112G 2070,0 ±2,0 [0,00, 250] 22,90 [-54, +121] 3 ejes, MEMS 356B08 1610,0 ±50,0 [0,50, 5000] 40,00 [-54, +77] 3 ejes, Piezoel. 356A45 1410,0 ±50,0 [0,70, 7000] 125,00 [-54, +85] 3 ejes, Piezoel. KS48C 750,0 ±6,0 [0,25, 130] 0,60 [-20, +120] 1 eje, Piezoel. 393B12 820,0 ±0,5 [0,15, 1000] 1,30 [-54, +82] 1 eje, Piezoel. 352A24 540,0 ±50,0 [1,00, 8000] 80,00 [-54, +121] 1 eje, Piezoel. ADXL335 10,7 ±3,6 [0,50, 550] 300,00 [-40, +85] 3 ejes, MEMS LIS344ALH 12,0 ±2,0 [1,00, 500] 50,00 [-40, +85] 3 ejes, MEMS MPU9250 5,8 ±16,0 [0,24, 500] 300,00 [-40, +85] 3 ejes, MEMS MPU6050 5,4 ±16,0 [0,24, 500] 400,00 [-40, +85] 3 ejes, MEMS •Consumo de energía bajo. Esta es una característica ideal para dispositivos alimentados por batería que requieren una duración prolongada de la misma. • Bajo costo. La producción en masa de estos dispositivos y los avances en la tecnología de fabricación han permitido que sean significativamente más económicos que otros de los acelerómetros más empleados en la industria, como el piezoeléctrico. •Capacidad de integración con facilidad en circuitos electrónicos de mayor tamaño. • Versatilidad. Permiten la medición de distintos tipos de magnitudes físicas: aceleraciones, velocidades angulares, temperatura, etc. Por ello, tienen utilidad en una amplia gama de sectores, como el aeroespacial, automotriz, electrónica de consumo, medicina, civil y robótica, entre otros. • Rendimiento. Presentan una precisión más que razonable en la mayoría de aplicaciones, en relación al precio del dispositivo. Como se comentó en el Capítulo 1, uno de los principales factores limitantes en las aplicaciones de monitorización de estructuras (SHM) es el precio de los sensores. Adicionalmente, en los últimos años ha habido en dicho campo un cambio progresivo de interés de estudios temporales o puntuales (offline) a monitorización continua y permanente en tiempo real (online). Para ello, se requiere proponer soluciones de bajo coste en cuanto a la instrumentación de los procesos. Por esa razón, el uso de MEMS en la monitorización de estructuras está aumentando significativamente [7]. La Tabla 3.1 presenta una comparativa de algunos acelerómetros empleados en SHM en la industria, ordenados por precio (sin IVA), con algunas de sus principales características. Se muestran sombreados aquellos de tipo MEMS. Como se puede observar, los de bajo coste son significativamente más asequibles que los piezoeléctricos, sin conllevar un aumento excesivo del ruido de 22 Capítulo 3. Adquisición de datos fondo y presentando un rango de frecuencias razonable para la gran mayoría de aplicaciones de SHM. Según el número de direcciones en los que son capaces de medir la aceleración, los sensores pueden ser de 1 eje o de 3 ejes. 3.3 El microcontrolador Arduino Arduino es una plataforma electrónica de código abierto basada en hardware y software de uso sencillo y bajo coste. Las placas de Arduino emplean un microcontrolador reprogramable que puede ser conectado a la mayoría de sensores, tanto analógicos como digitales [28]. Para realizar dicha conexión, disponen de una serie de pines de entrada/salida para controlar dispositivos externos de manera independiente. Existe una gran variedad de placas, siendo algunas de las más comunes Arduino Uno (ver Figura 3.4), Arduino Leonardo, Arduino Galileo, Arduino Mega 2560 y Arduino Due. La mayoría de ellas poseen un microcontrolador AVR Atmel-8 bits, como por ejemplo, en el caso de la placa Arduino Uno R3, el ATmega16U2. Figura 3.4 Placa Arduino Uno R3. Fuente: Arduino Store1. Además, las placas Arduino permiten la comunicación serial entre varios dispositivos mediante un puerto TX/RX (generalmente USB) y un bus I 2 C. El estándar I 2 C (Inter-Integrated Circuit) es un protocolo de comunicación serial cuya operación requiere de dos cables: uno para señal de reloj (SCL) y otro para señal de datos (SDA). En el bus, cada dispositivo tiene una dirección, la cual sirve como identificación para acceder a ellos de forma individual. La arquitectura del I 2 C es de tipo maestro-esclavo, lo que significa que es el dispositivo “maestro” el que inicia la comunicación, envía y recibe datos de los dispositivos “esclavos”. Estos últimos no pueden, ni iniciar comunicación, ni interactuar entre ellos de manera directa. Otra característica del protocolo es la sincronía, ya que el maestro proporciona una señal de reloj con la que se mantiene sincronizados a todos los dispositivos del bus. El soporte I 2 C en Arduino está vinculado físicamente a ciertos pines, los cuales varían de un modelo a otro. En el caso del Arduino Uno, estos son el pin A4 para la señal SDA y el pin A5 para la señal SCL [29]. Una de las grandes ventajas del uso de Arduino reside en la gran comunidad de desarrolladores y usuarios existente, la cual es activa en foros y en comunicación constante enfocada a la resolución de problemas. La plataforma es fácil de utilizar, tiene una interfaz atractiva y su diseño es flexible. Además, la mayoría de sensores MEMS pueden interactuar con Arduino y recibir alimentación del 1Arduino Store: https:// store.arduino.cc/ products/ arduino-uno-rev3 3.4 Arquitectura del sistema de adquisición de datos 23 mismo [28]. 3.4 Arquitectura del sistema de adquisición de datos El objetivo al diseñar un sistema de adquisición de datos para llevar a cabo una identificación modal es tener la capacidad de recopilar las señales de las aceleraciones de la estructura, medidas con sensores, y trasmitir dichos datos a un ordenador para ser almacenados en memoria, visualizados y/o procesados. El esquema del sistema estudiado en el presente trabajo es el mostrado en la Figura 3.5. ADC Bus I2C Esclavos I2C Sensor MEMS 1 Placa Arduino Ordenador ADC Sensor MEMS 2 ADC SDA SCL Sensor MEMS N Maestro I2C Timer / Función ISR Módulo UART Puerto USB Puerto COM Software MATLAB - Identificación - Visualización - Almacenamiento Disparador Serial Comunicación serial por USB Señal digital Figura 3.5 Esquema de la arquitectura de un sistema de adquisición de datos para identificación modal. Como se puede apreciar en dicho diagrama, el inicio del flujo de datos comienza en los sensores, 24 Capítulo 3. Adquisición de datos que transforman la aceleración de la estructura en una señal de voltaje interpretable por hardware. Esta se digitaliza mediante un conversor analógico digital (ADC), el cual viene incorporado en el acelerómetro. Este proceso se desarrolla con más amplitud en la Subsección 3.5.1. A continuación, la señal de datos procedente de cada uno de los sensores (esclavos) se transmite mediante el bus de protocolo I 2 C hacia la placa Arduino (maestro). Para controlar el flujo de datos, se emplea una interrupción que pausa la actividad normal del microprocesador para atender una rutina dada. En este caso, se propone la utilización de eventos programados mediante el Timer o temporizador incorporado en Arduino. El temporizador es un dispositivo que contabiliza el tiempo a través del reloj interno del sistema (16 MHz) o una versión preescalada del mismo [30]. Puede ser configurado en modo overflow o en modo comparación. Este último permite generar una interrupción cuando ocurre el evento dado por la igualdad entre el registro del temporizador que almacena la cuenta del tiempo y el registro que contiene un valor fijado por el usuario. La interrupción conlleva la ejecución de una rutina denominada ISR (Interrupt Service Routine) que se puede configurar para provocar la petición de medidas a los sensores. Este proceso viene descrito por el siguiente diagrama de flujo: Inicio Configuración de sensores y timer Lectura de las medidas de los sensores ISR Envío de datos por el puerto USB No Sí Interrupción de timer Figura 3.6 Diagrama de flujo de la utilización del timer para la lectura de las medidas en Arduino. De esta manera, se dispondrá de medidas procedentes de todos los sensores a una frecuencia de muestreo fijada en la configuración del temporizador. Cuando el microcontrolador de Arduino recibe un conjunto de datos correspondiente a las medidas de las aceleraciones de cada uno de los sensores en el momento en el que se produce la interrupción, este se envía de forma serial por el puerto USB presente en la placa. Seguidamente, la información es recibida por un ordenador mediante uno de sus puertos COM. A través de MATLAB, 3.5 Preprocesado de datos 25 es posible detectar el recibimiento de datos en dicho puerto y tratarlos como resulte conveniente: almacenando los datos en memoria, visualizándolos o transformándolos. Así, es posible disponer de secuencias de medidas de la respuesta dinámica de la estructura a estudiar, listas para ser analizadas e interpretadas por un algoritmo de identificación modal. 3.5 Preprocesado de datos Desde que la señal analógica de voltaje es generada por el acelerómetro, hasta que los datos son datos son apropiados para ser analizados mediante técnicas de identificación, la información ha de sufrir una serie de transformaciones. 3.5.1 Conversión analógico-digital En primer lugar, la señal analógica producida por los sensores ha de ser convertida en una secuencia de números digitales correspondientes a su amplitud. Esta transformación es llevada a cabo por un conversor analógico digital o sus siglas en inglés ADC (Analog-Digital Converter), y se basa en la discretización del tiempo, a través del muestreo, y de la amplitud de la señal, a través de la cuantización. La cuantización conlleva siempre un error asociado a la aproximación de la señal de amplitud continua a un número finito de valores discretos. Sin embargo, el error de cuantización suele ser mucho menor que el procedente de otras fuentes como los sensores y, por tanto, generalmente se desprecia salvo en casos de resolución escasa [24,25]. La resolución de un ADC se refiere al menor paso de voltaje que puede ser detectado, lo que está directamente relacionado con el bit menos significativo o sus siglas en inglés LSB (Least Significant Bit) del conversor. Siendo Nbit el número de bits disponibles en el procesador interno del conversor, entonces se dispone de Nval =2Nbit números para digitalizar la señal. Si el voltaje de entrada es simétrico, es decir, ±Vm,ADC (por ejemplo, ±5 V), la resolución digital ∆VADC del conversor se puede expresar como [24,25]: ∆VADC =2Vm,ADC Nval −1(3.5) En ocasiones, también se denomina resolución al número de bits del conversor, Nbit . Cabe remarcar que, dado que la mayoría de conversores tienen asociado un ruido interno mayor que un paso, la resolución efectiva se corresponde con el número de bits libres de ruido. Se llama nivel de ruido al número de bits ocupados por ruido cuando la entrada al ADC es nula. Típicamente, este suele ser de dos bits en conversores de calidad de 24-bits [24]. Este número de bits, LSB, viene dado por la siguiente expresión: 2LSB =δV0,ADC ∆VADC (3.6) donde δV0,ADC es el ruido del conversor para una banda de frecuencia concreta. El número efectivo de bits será: Nbit,e f f =Nbit −LSB (3.7) 26 Capítulo 3. Adquisición de datos El rango dinámico del ADC se define como el ratio entre el máximo voltaje de entrada permitido Vm,ADC y el nivel de ruido δV0,ADC, expresado normalmente en decibelios: DRADC =20 logVm,ADC δV0,ADC =20 log∆VADC (Nval −1)/2 ∆VADC ·2LSB ≃ ≃20 log2Nbit 2LSB+1=20 log(2Nbit,e f f −1) = 20(Nbit,e f f −1)log 2 ≈6(Nbit,e f f −1) (3.8) Por ejemplo, si el ADC tiene 21 bits efectivos, el rango dinámico será aproximadamente de 120 dB [24,25]. Es importante considerar también que el rango dinámico del sistema completo de toma de datos depende también del valor del rango dinámico del sensor. 3.5.2 Frecuencia de muestreo. Aliasing La definición del valor de la frecuencia de muestreo representa uno de los pasos más importantes en la configuración de un ensayo dinámico. El paso temporal ∆t entre dos medidas determina la frecuencia de muestreo (sampling frequency) fs , así como la frecuencia de Nyquist fN , que se corresponde con el límite superior del rango de frecuencias observables [25]: fs=2fN=1 ∆t(3.9) El teorema de Shannon afirma que la frecuencia de muestreo ha de ser al menos dos veces mayor que la frecuencia más alta de interés contenida en la señal. Una selección inadecuada de dicha frecuencia provoca que, en la reconstrucción de las ondas tras la conversión analógico-digital, señales de alta frecuencia aparezcan como señales de baja frecuencia [24]. A este efecto se le denomina aliasing. Por el teorema de Parseval [25] (siendo Gx(f) la función de densidad espectral de la señal x(t)): E[x2] = σ2 x=2ZfN 0 Gx(f)d f (3.10) se observa que, dado que la energía/potencia de la señal se mantiene tras la digitalización, si hay componentes fuera de la banda de Nyquist, la energía de estas componentes se transfiere artificialmente al intervalo [−fN,+fN]. 2Z∞ 0 Gx,1(f)d f =2ZfN 0 Gx,2(f)d f (3.11) Por ejemplo, se considera el caso de una señal analógica consistente en una senoide de freal =4 Hz, con la que se selecciona una frecuencia de muestreo de fs=5 Hz. Dado que la frecuencia de dicha senoide se encuentra fuera de la banda de Nyquist ( freal =4Hz >2,5Hz =fN=fs/2 ) su energía/potencia se transfiere a la frecuencia de 1 Hz. Este fenómeno se puede observar en la Figura 3.7. Para evitar posibles interpretaciones erróneas de los datos debido a este efecto, la solución es el uso de filtros antialiasing que eliminan las componentes frecuenciales de la señal que se encuentran por encima de la frecuencia de Nyquist. Estos son esencialmente filtros paso bajo o sus siglas en inglés LPF (Low-Pass Filter) que atenúan las frecuencias que superan una determinada frecuencia 3.5 Preprocesado de datos 27 Señal real Señal reconstruida Figura 3.7 Ejemplo de efecto aliasing en la reconstrucción de una señal tras digitalización. de corte. En general, los filtros antialiasing están incluidos en los sensores, en caso contrario, es importante incorporar uno antes de la conversión analógico-digital. Este último punto es crítico, puesto que, una vez ha tenido lugar la conversión, las componentes frecuenciales generadas por aliasing no pueden ser reconocidas y, por tanto, no se puede corregir el fenómeno con filtros digitales [24]. 3.5.3 Eliminación de tendencias y manchado espectral La eliminación de tendencias consiste en eliminar las componentes en régimen permanente de las señales (desfases y derivas), de tal manera que tengan media nula. La razón para realizar esto es que normalmente dichas componentes se deben a ruido de baja frecuencia proveniente de los sensores, además de que componentes frecuenciales muy por debajo de la más baja de la estructura no suelen ser de interés en un análisis dinámico. Para llevar a cabo esta operación, existen dos posibilidades, principalmente [25]: • Empleo de filtros digitales de paso alto o sus siglas en inglés HPF (High-Pass Filter). Se emplean generalmente cuando el error en régimen permanente es considerable. Sin embargo, pueden perturbar la señal en la banda de Nyquist hasta cierto punto que pueda no resultar aceptable. • Uso de segmentación. Consiste en dividir la señal en tramos consecutivos de duración T y restar la el valor medio de cada uno de ellos. Otra fuente de error en el espectro frecuencial de una señal es el denominado manchado espectral oleakage. Este fenómeno ocurre al realizar la transformada discreta de Fourier, debido a que el algoritmo FFT (Fast Fourier Transform) es solo adecuado para: 1. Funciones periódicas. 2. Ventanas de señal cuyo primer y último valor es el mismo, es decir, tienen un número entero de periodos. La segunda de las condiciones se cumple en escasas ocasiones. Cuando la duración de la ventana, Nt∆t , no es un múltiplo del periodo una determinada componente frecuencial, T , esa componente no podrá ser asociada a ninguna frecuencia discreta de la FFT (quedará entre dos de ellas), por lo 34 Capítulo 3. Adquisición de datos capítulo, la función de autocorrelación de dicha coordenada modal tiene la siguiente forma, siendo Auna constante: R(τ) = Ae−ξωminτcos(ωmin τ)(3.17) De esta manera, se aprecia como, para que la amplitud de dicha función se reduzca a un 4%, se requiere el siguiente desfase Tmem: |R(Tmem)| |R(0)|=e−ξωminTmem =e−π≈4% ⇒ξ ωmin Tmem =π⇒Tmem =1 2ξfmin (3.18) El tiempo Tmem es el tiempo de memoria del sistema, y asegura la posibilidad de estimar el amortiguamiento del mismo de manera fiable. Sin embargo, para estimar la función de autocorrelación de la respuesta, se requieren unos 20 segmentos con el objetivo de realizar las medias y el promediado adecuados (véase el método de Welch en la Subsección 4.5.1). Por ello, el tiempo total del ensayo Ttot debe ser: Ttot >20 2ξfmin =10 ξfmin (3.19) La anterior ecuación proporciona el tiempo mínimo requerido para un ensayo de identificación modal, conociendo una predicción de la mínima frecuencia natural de interés y del ratio de amortiguamiento del sistema. 4 Identificación experimental offline 4.1 Introducción En este capítulo se va a tratar la identificación estructural de propiedades modales de forma offline, es decir, evaluando las características de forma puntual. Los métodos de identificación offline suponen una mayor precisión en los resultados; sin embargo, requieren un análisis con mayor coste computacional posterior a la realización del ensayo y el almacenamiento de los datos obtenidos, por lo que no permiten estimación y monitorización continua en tiempo real. La posibilidad y el método a utilizar para realizar la extracción de propiedades modales de una estructura a partir de ensayos de vibración dependen en gran medida de la información que se conoce acerca de la excitación que sufre la misma durante el test. De manera general, se puede considerar que la respuesta medida de una estructura proviene de cuatro fuentes distintas [12]: Respuesta medida =(1) Vibración libre a partir de condiciones iniciales desconocidas +(2) Vibración forzada por excitación conocida +(3) Vibración forzada por excitación desconocida +(4) Errores de medida (desconocidos) (4.1) Los términos (1), (2) y (3) dependen de las propiedades modales de la estructura y, por tanto, son útiles para poder identificarlas. Sin embargo, no siempre es posible desarrollar un método para realizar el análisis apropiadamente. El segundo término viene completamente determinado por las propiedades modales y la excitación conocida. El primer término depende, además de de las propiedades modales, de las condiciones iniciales, que generalmente son desconocidas. En cuanto al tercer término, si la excitación es completamente arbitraria, la respuesta incluirá demasiadas incógnitas y, en consecuencia, la identificación es imposible. Sin embargo, para el caso en el que la excitación se puede modelar estadísticamente, se puede proceder con una identificación modal, aunque esta será tan precisa como las suposiciones acerca de dicha excitación [12]. Existen dos principales enfoques ante el análisis modal de estructuras: el análisis modal experimental o EMA (Experimental Modal Analysis) y el análisis modal operacional u OMA (Operational Modal Analysis). El primero de ellos, también conocido como análisis modal clásico por ser el pionero, se basa en la medida de la respuesta estructural ante unas fuerzas dinámicas también conocidas y medidas. A través de ensayos de excitación controlada, relacionan la carga aplicada con la respuesta observada para obtener las propiedades modales. Los análisis pueden ser llevados a cabo en el dominio de la frecuencia o en el dominio del tiempo. En el primer caso, la identificación 35 36 Capítulo 4. Identificación experimental offline se realiza ajustando un modelo a la función de respuesta en frecuencia o FRF de la estructura, mientras en el dominio del tiempo el ajuste se realiza a la función de respuesta a impulso o IRF (Impulse Response Function), que se corresponde con la transformada inversa de Fourier de la FRF [6,32]. El análisis modal operacional, por otro lado, emplea únicamente la respuesta estructural para llevar a cabo la identificación del sistema. Esto permite realizar los ensayos bajo condiciones reales de operación (por ello el nombre del método) y hallar las propiedades modales bajo las cargas de servicio. Puesto que, generalmente, las cargas aerodinámicas a las que se ven sometidas las estructuras aeronáuticas son difíciles de medir con precisión, los métodos operacionales son muy apropiados para aplicaciones en dicha industria. Por ello, este capítulo se va a centrar en la descripción de los análisis de este tipo [6]. 4.2 Fundamentos del OMA El análisis modal operacional u OMA responde a la necesidad de identificación de propiedades modales de estructuras sometidas a excitaciones que no pueden ser evaluadas [33]. En estos casos, solo la respuesta del sistema ante las fuerzas ambientales está disponible para llevar a cabo el proceso de análisis. Las propiedades estimadas son típicamente las frecuencias naturales, los modos de vibración y los ratios amortiguamiento modales de la estructura. Como se comentó en el Capítulo 1, este método se comenzó a desarrollar durante la década de 1990 para estudios ambientales en el marco de la ingeniería civil (puentes, edificios, puentes de peatones, estructuras históricas, turbinas eólicas, presas) y, debido a sus numerosas ventajas, se ha extendido a áreas como la ingeniería aeroespacial (identificación modal en vuelo, estudios de flameo) o la mecánica (barcos, camiones, chasis de vehículos, motores) [24,34]. El OMA se basa en las siguientes suposiciones [24]: • Linealidad. La respuesta total del sistema a una combinación de entradas o excitaciones es igual a la suma de las correspondientes respuestas. • Estacionariedad. Las propiedades dinámicas de la estructura se mantienen constantes, es decir, los coeficientes de las ecuaciones diferenciales que gobiernan la respuesta dinámica del sistema son independientes del tiempo. • Observabilidad. La colocación de sensores ha sido diseñada para ser capaz de observar adecuadamente todos los modos de interés. Adicionalmente, dado que los ensayos se realizan sin conocer las excitaciones a las que se somete la estructura, son necesarias algunas hipótesis sobre dichas entradas. Si estas están definidas por ruido blanco, es decir, su densidad espectral de potencia (PSD) es constante, todos los modos se encuentran idénticamente excitados, conteniendo la respuesta información completa sobre la estructura. Sin embargo, esto raramente sucede, ya que, en general, la excitación de entrada tiene su propia PSD, por lo que los modos están ponderados por las propiedades tanto de la estructura como de la entrada. Por ello, pueden aparecer en la respuesta armónicos procedentes de la excitación y, tal y como se detalló en el Capítulo 3, ruido o incluso armónicos espurios procedentes de equipamiento externo. La respuesta medida se puede ver de esta manera como la salida de un sistema combinado compuesto por el sistema de excitación (cargado con ruido blanco) y la estructura ensayada [24] (ver Figura 4.1). 4.3 Modelos de dinámica estructural 37 Salida Sistema combinado indentificado Sistema de excitación Sistema estructural Respuesta dinámica Entrada Ruido blanco Figura 4.1 Representación esquemática de la identificación modal basada en respuesta. 4.3 Modelos de dinámica estructural La mayoría de métodos de identificación modal se basan en modelos experimentales del sistema ensayado que se ajustan a los datos obtenidos del test y de los cuales es posible extraer propiedades modales. En esta sección se van a describir aquellos modelos y conceptos necesarios para entender los métodos de identificación que se tratarán con mayor detalle en el resto del capítulo. Los modelos son idealizaciones de sistemas dinámicos reales, por lo que conllevan necesariamente simplificaciones y suposiciones, como las enunciadas anteriormente: linealidad y estacionariedad. Dichos modelos pueden estar establecidos en el dominio del tiempo o en el dominio de la frecuencia, con ecuaciones en forma continua o en forma discreta [35]. 4.3.1 Respuesta frecuencial y respuesta impulsiva El comportamiento dinámico de una estructura se puede representar en el dominio del tiempo mediante un conjunto de ecuaciones diferenciales, o en el dominio de la frecuencia por un conjunto de ecuaciones algebraicas. Las ecuaciones del movimiento que definen el comportamiento dinámico de un sistema de múltiples grados de libertad o MDOF (Multiple Degrees Of Freedom) son generalmente expresadas en el dominio del tiempo mediante el siguiente conjunto de ecuaciones diferenciales lineales de segundo orden: M¨ x(t)+C˙ x(t)+Kx(t) = f(t)(4.2) donde ¨ x(t) , ˙ x(t) y x(t) son los vectores de aceleración, velocidad y desplazamiento, respectivamente; M , C y K definen las matrices de masa, amortiguamiento y rigidez, y f(t) es el vector de fuerza. Estas ecuaciones se pueden desacoplar resolviendo un problema de autovalores, suponiendo amortiguamiento proporcional. De esta manera, la respuesta completa del sistema se puede expresar como superposición de las respuestas modales. Aplicando la transformada de Fourier a la ecuación diferencial matricial anterior, se obtiene un conjunto de ecuaciones lineales algebraicas: (−ω2K+iωC+K)˜ x(ω) = ˜ f(ω)(4.3) 38 Capítulo 4. Identificación experimental offline siendo ˜ x(ω) y ˜ f(ω) las transformadas de Fourier de x(t) y f(t) . Como sistema lineal estacionario, puede ser representado por su FRF, que se define como el ratio entre las tranformadas de la entrada y la salida del sistema. Sobre la ecuación anterior: ˜ Z(ω) = −ω2K+iωC+K⇒˜ Z(ω)˜ x(ω) = ˜ f(ω)(4.4) donde ˜ Z(ω)es la matriz de rigidez dinámica. La FRF se define entonces como: ˜ H(ω) = ˜ Z−1(ω) = ˜ x(ω) ˜ f(ω)(4.5) Esta función representa las propiedades inerciales, disipativas y elásticas de la estructura. Como se demuestra en [36], la FRF se puede expresar en términos de las propiedades modales como una expansión en fracciones parciales (denotando el superíndice [−]∗el conjugado complejo): ˜ H(ω) = N ∑ m=1 Rm iω−λm +R∗ m iω−λ∗ m (4.6) donde N es el número de modos, Rm es la matriz de residuos del modo m , y λm es polo del modo m, el cual se define como: λm=−αm+iωd,m=−ξmωn,m+iωn,mq1−ξ2 m(4.7) siendo ωn,m y ωd,m las frecuencias natural y amortiguada del modo m , y ξm el ratio de amortiguamiento de dicho modo. La matriz de residuos Rm viene dada por el producto tensorial del modo de vibración φmconsigo mismo (denotando el superíndice [−]Tla traspuesta): Rm=QmφmφT m(4.8) incluyendo el escalar Qm información sobre el factor de escala modal. Sustituyendo la definición 4.8 en la expresión 4.6, se llega a ( [−]H se refiere al operador hermitiano, que es equivalente a [−]∗T ): ˜ H(ω) = N ∑ m=1 Qmφmφ∗ m iω−λm +Qmφ∗ mφH m iω−λ∗ m (4.9) La expresión anterior de la FRF tiene importantes implicaciones en el análisis modal. En primer lugar, indica que cada uno de los modos contribuye a la respuesta del sistema, independientemente de la frecuencia, por lo que no se puede excitar solamente un modo mediante una onda de una sola frecuencia. Sin embargo, dado que los amortiguamientos son pequeños, cerca de una frecuencia natural (resonancia) el sumatorio se puede aproximar por el término del correspondiente modo [24]. Teniendo en cuenta la propiedad de convolución de la correlación Rxy(τ) , y que esta última forma un par de transformadas de Fourier con la función de densidad espectral Gxy(ω)[25]: Rxy(τ) = x(−t)∗y(t) Gxy(ω) = F{Rxy(τ)})⇒Gxy(ω) = X∗(ω)Y(ω)(4.10) Generalizando la expresión anterior para las matrices de densidad espectral o matrices PSD (Power Spectral Density)Gx(ω)de la respuesta x(t)del sistema y Gf(ω)de la fuerza f(t): 4.3 Modelos de dinámica estructural 39 Gx(ω) = ˜ x∗(ω)˜ xT(ω) =      Gx1x1(ω)Gx1x2(ω)··· Gx1xL(ω) Gx2x1(ω)Gx2x2(ω)··· Gx2xL(ω) . . .. . ..... . . GxLx1(ω)GxLx2(ω)··· GxLxL(ω)      (4.11) Gf(ω) = ˜ f∗(ω)˜ fT(ω)(4.12) Dada la definición de FRF, ˜ x(ω) = ˜ H(ω)˜ f(ω) , se llega a la siguiente relación entre esta y la matriz PSD: Gx(ω) = ˜ H∗(ω)Gf(ω)˜ HT(ω)(4.13) Tal y como se comentó en la Sección 4.2, una de las hipótesis de partida del OMA es que la entrada al sistema es ruido blanco Gaussiano. Por tanto, su matriz PSD es una matriz constante, permitiendo que la matriz PSD de la salida se puede expresar en forma de polos y residuos, de manera similar a la matriz FRF: Gx(ω) = N ∑ m=1 φmγT m −iω−λm +φ∗ mγH m −iω−λ∗ m +γmφT m iω−λm +γ∗ mφH m iω−λ∗ m (4.14) donde γm es el vector operacional de referencia asociado al modo m . A diferencia del vector de participación modal Qmφm que aparece en la expresión polo-residuo de la matriz FRF, este vector depende de todas las propiedades modales del sistema, la posición de las entradas, y su matriz de correlación Rf(τ) [24]. El vector γm no es proporcional al modo de vibración φm ; sin embargo, sí que es aproximadamente proporcional para casos en los que el amortiguamiento es pequeño [25]. En tal situación: γm∼ =cmφm(4.15) La Ecuación 4.14 tiene el inconveniente de que las filas o columnas de Gx(ω) no son proporcionales a los modos, lo que dificulta la extracción de parámetros a partir de dicha matriz. Por ello, se ha popularizado el uso de una forma en la que se emplea únicamente la mitad positiva del espectro frecuencial, Gxx(ω), reduciendo el número de polos por modo de cuatro a dos: Gxx(ω) = N ∑ m=1 γmφT m iω−λm +γ∗ mφH m iω−λ∗ m (4.16) Añadiendo, a su vez, la aproximación indicada en la Ecuación 4.15, se llega a la siguiente forma de la matriz PSD: Gxx(ω) = N ∑ m=1 cmφmφT m iω−λm +c∗ mφ∗ mφH m iω−λ∗ m (4.17) En la matriz anterior, que constituye la descomposición modal de la matriz de densidad espectral, tanto las columnas como las filas son proporcionales a los modos de vibración, teniendo la misma forma que una oscilación libre amortiguada en el dominio de la frecuencia [25]. Se comprueba que la expresión contiene tanto información sobre las frecuencias naturales y ratios de amortiguamiento 40 Capítulo 4. Identificación experimental offline del sistema en los polos, como sobre los modos de vibración en los residuos. Dado que, tanto la FRF con la función de respuesta impulsiva (IRF), como la función de densidad espectral con la función de correlación, forman pares de transformadas de Fourier, las analogías expresadas entre FRF y densidad espectral se pueden realizar de manera similar entre las IRF y las funciones de correlación, dentro del contexto del OMA. La IRF determina en el dominio del tiempo la respuesta de un sistema dinámico ante una entrada genérica. Su expresión matemática es la siguiente [24]: h(t) = N ∑ m=1 Rmeλmt+R∗ meλ∗ mt(4.18) Como se puede apreciar, la expresión anterior coincide con la transformada inversa de Fourier de la FRF, ˜ H(ω) , en la Ecuación 4.6. Al igual que en los casos anteriores, los polos λm contienen información sobre las frecuencias naturales y los ratios de amortiguamiento, mientras que los modos de vibración están indicados por los residuos Rm(véase la Ecuación 4.8). En el contexto del OMA, en el que solo se conoce la salida del sistema y, al igual que antes, se asume como entrada un ruido blanco gaussiano, la matriz de correlación Rx(τ) se puede expresar también como suma de exponenciales complejas: Rx(τ) =            N ∑ m=1 γmφT meλmτ+γ∗ mφH meλ∗ mττ≥0 N ∑ m=1 φmγT me−λmτ+φ∗ mγH me−λ∗ mττ<0 (4.19) Una vez más, se puede observar la correspondencia entre la expresión anterior y la de la matriz de densidad espectral Gx(ω) (Ecuación 4.14), siendo la transformada inversa de Fourier. Es interesante mencionar como la parte causal de la función ( τ>0 ) contiene los polos estables, es decir, α>0 (ver Ecuación 4.7), mientras que la parte no causal contiene los inestables. Los métodos de identificación temporales estiman generalmente la parte causal, reduciendo el número de polos por modo de cuatro a dos, de manera análoga al uso de la matriz Gxx(ω) (reducida) en lugar de la matriz completa, Gx(ω) , en los métodos frecuenciales. De nuevo, las filas de la parte positiva (causal) de la matriz Rx(τ) tienen la forma de una oscilación libre amortiguada en el dominio del tiempo, al igual que ocurría con las filas y columnas de Gxx(ω) en el dominio de la frecuencia. Esta propiedad es clave en el empleo de algunos métodos OMA del dominio temporal. 4.3.2 Modelos de espacio de estados Los modelos de espacio de estados (state-space models) se emplean para convertir la ecuación diferencial de segundo orden de 4.2 en dos ecuaciones de primer orden, que toman la siguiente forma: ˙ x(t) = Acx(t)+Bcu(t) ˙ y(t) = Ccx(t)+Dcu(t)(4.20) La primera de las anteriores ecuaciones es la ecuación de estado ode transición, donde la matriz Ac es la matriz de estado o de transición y la matriz Bc es la matriz de entrada. La segunda de las 4.3 Modelos de dinámica estructural 41 ecuaciones se corresponde con la ecuación de observación. La matriz Cc es la matriz de salida y Dc es la matriz de transmisión directa. El subíndice [−]c denota tiempo continuo. El vector de estado, x(t) , es una concatenación del vector de desplazamientos y del vector de velocidades (en ese orden) de cada grado de libertad. El vector de respuesta y(t) , en cambio, indica la respuesta medida del sistema. La definición de las matrices del modelo es: Ac=0 I −M−1K−M−1C(4.21) Bc=0 M−1S(4.22) Cc=Cd−CaM−1K Cv−CaM−1C(4.23) Dc=CaM−1S(4.24) La matriz S forma parte de la factorización del vector de fuerzas, tal que f(t) = Su(t) . En cuanto a las matrices Ca , Cv y Cd , indican la localización de las salidas observadas para aceleraciones, velocidades y desplazamientos, respectivamente. La deducción de las anteriores expresiones se puede encontrar en [32]. Cabe destacar que, aunque una estructura real viene caracterizada por infinitos grados de libertad, en la práctica de ensayos este número se reduce a varias docenas o menos [24]. Las ecuaciones expresadas en 4.20 constituyen el modelo de espacio de estados en tiempo continuo, que permite calcular la respuesta x(t) de la estructura dada una entrada o excitación u(t) . Una propiedad importante del modelo es que un único sistema puede ser representado de infinitas formas distintas, cada una de ellas conociéndose como realización. En los ensayos experimentales el objetivo es obtener una de estas realizaciones. Cada una de ellas se obtiene a partir de una transformación de semejanza dada por una matriz T cuadrada no-singular, tal que x(t) = Tz(t) . Como las matrices de las ecuaciones de distintas realizaciones se relacionan por estas transformaciones de semejanza, los autovalores, y por tanto las propiedades modales del sistema, se preservan. Como se detalló en el Capítulo 3, la señal analógica medida por los sensores durante los ensayos experimentales es digitalizada por un convertidor analógico-digital para que los datos puedan ser almacenados y procesados por un ordenador. Por ello, la información disponible sobre la respuesta del sistema bajo estudio es discreta en tiempo. Consecuentemente, es necesario transformar las ecuaciones de estado en tiempo continuo a una versión en tiempo discreto para poder ajustarse a los datos experimentales. Dado un tiempo de muestreo ∆t , las ecuaciones se pueden discretizar y resolver para todo instante de tiempo tk=k·∆t , con k∈N . Es necesario, primeramente, hacer suposiciones sobre el comportamiento de las variables entre dos muestras [24,35]. Aplicando el método de discretización conocido como retenedor de orden cero o, en inglés, Zero-Order Hold (ZOH), se supone que las funciones dependientes del tiempo son constantes entre muestras discretas consecutivas [32], tal que, por ejemplo: x(t) = xk,∀t∈[k∆t,(k+1)∆t](4.25) Bajo esta suposición, el modelo de espacio de estados en tiempo continuo se puede convertir en modelo de espacio de estados en tiempo discreto: 42 Capítulo 4. Identificación experimental offline xk+1=Axk+Buk yk=Cxk+Duk (4.26) Las variables x(t) , y(t) y u(t) se han sustituido por sus valores en los instantes discretos k∆t . Las matrices A , B , C y D son, respectivamente, la matriz de estado discreta, la matriz de entrada discreta, la matriz de salida discreta y la matriz de transmisión directa. Bajo la suposición de ZOH, las relaciones entre las matrices del modelo continuo y del modelo discreto son: A=eAc∆t(4.27) B= (A−I)A−1 cBc(4.28) C=Cc(4.29) D=Dc(4.30) Como se puede observar, el muestreo con ZOH no afecta a las últimas dos matrices. Al igual que esta, existen otras suposiciones para el comportamiento de las variables entre dos muestras consecutivas, como es el retenedor de orden uno o First-Order Hold (FOH), el cual resulta en relaciones más complejas entre las matrices en forma continua y discreta. El modelo expresado en la Ecuación 4.26 se corresponde con un modelo determinista, ya que no se asume ningún tipo de incertidumbre en las entradas que dirigen el sistema. Por ello, es necesario incluir componentes estocásticos para explicar las medidas tomadas realmente. Esto da lugar al modelo determinista-estocástico de espacio de estados en tiempo discreto: xk+1=Axk+Buk+wk yk=Cxk+Duk+vk (4.31) donde wk es el denominado ruido de proceso, debido a perturbaciones e imprecisiones en el modelo, y vk es el ruido de medición debido a errores del los sensores. La ecuación de estado o de transición modela el comportamiento dinámico del sistema, la ecuación de observación expresa que parte de la respuesta dinámica puede ser medida como salida del modelo. En el contexto del OMA, la entrada al sistema es desconocida, por lo que no es posible establecer valores para uk . La respuesta del sistema viene generada entonces únicamente por los procesos estocásticos wk y vk , que incluyen los efectos que tendría el vector uk en la misma. De esta manera se obtiene el modelo estocástico de espacio de estados: xk+1=Axk+wk yk=Cxk+vk (4.32) Es importante remarcar que el rol que tiene el ruido de proceso es de modelar la entrada que dirige el comportamiento dinámico del sistema, mientras que el ruido de medición modela las perturbaciones sobre la respuesta del mismo. Cuando se adopta un modelo de espacio de estados estocástico, el objetivo es, además de determinar el orden del sistema, obtener una estimación de 4.3 Modelos de dinámica estructural 43 una realización de las matrices A y C a partir de un conjunto extenso de medidas yk de la respuesta del sistema. La matriz de estado, A transforma un estado xk al estado siguiente xk+1 junto con la contribución del ruido de proceso, a la vez que la matriz C da, a partir del estado del sistema, la parte del mismo que es observable. Concretamente, el vector de respuesta yk viene dado por la parte observable del estado, más el ruido de medición [24]. Ambos ruidos son imposibles de medir. Se asume que son procesos estacionarios de ruido blanco con media nula, con matrices de covarianza dadas por: Ewp vpwT qvT q=     Q S STRp=q 0p=q (4.33) donde los índices p y q representan instantes de tiempo genéricos [35]. La estimación de las matrices Q , R y S forma también parte del proceso de identificación. En los métodos de identificación de espacios de estados estocásticos o sus siglas en inglés SSI (Subspace System Identification) la suposición de ruido blanco para wk y vk es de vital importancia. Como se indicó en la Sección 4.2, si la entrada al sistema contiene componentes frecuenciales dominantes, en lugar de tener una densidad espectral constante, estas aparecen como polos de la matriz A junto con los autovalores reales de la estructura. Esto es equivalente a realizar la identificación del sistema combinado de la estructura más el sistema de excitación [24]. Además de las propiedades descritas sobre los términos de ruido, los modelos de espacio de estados estocásticos presentan otras propiedades que resultan de gran utilidad para los algoritmos de identificación. Estas se basan en que el vector de estado xk se puede caracterizar como un proceso gaussiano estacionario de media nula. Debido a su estacionariedad, la matriz de covarianza del estado Σes independiente del instante k: Σ=E[xk·xT k](4.34) Dado que wkyvkson ruidos blancos de media nula independientes de xk: E[xk·wk] = 0(4.35) E[xk·vk] = 0(4.36) La respuesta del sistema también se puede representar como un proceso gaussiano de media nula, cuya matriz de correlación para un retraso de tiempo de τ=i·∆tviene dada por (para simplificar la notación, se va a expresar la matriz de correlación Rx(τ=i·∆t)como Ri): Ri=E[xk+i·xT k](4.37) La matriz de covarianza nuevo estado-salida se define como: G=E[xk+1·yT k](4.38) Teniendo en cuenta las ecuaciones 4.32, 4.34, 4.35 y 4.36, se puede realizar la siguiente deducción: 50 Capítulo 4. Identificación experimental offline para definir múltiples sistemas equivalentes de un grado de libertad (SDOF) caracterizados por las mismas propiedades modales que el sistema completo de varios grados de libertad (MDOF) bajo estudio. Dado que la SVD proporciona los valores singulares ordenados de mayor a menor, cerca de un pico de resonancia el primer valor singular contiene información sobre el modo dominante a dicha frecuencia (ver Figura 4.5) . Esta propiedad puede ser empleada para identificar modos cercanos, lo que no era posible con el método Peak-Picking. El número de valores singulares dominantes a una determinada frecuencia se corresponde con el número de modos que contribuyen significativamente a la respuesta del sistema en esa frecuencia. La Figura 4.6 incluye un ejemplo de valores singulares para el caso de tres mediciones y un solo modo en 1 Hz y otro para dos modos en 1 y 1,05 Hz. Modo 2 Valores singulares Modo 1 ( ) ( ) ( ) ( ) Figura 4.5 Valores singulares de la matriz PSD para un sistema con dos modos de vibración. Tomada de [12]. Frecuencia [Hz] 0 0.5 1 1.5 2 10 -1 10 0 10 1 10 2 10 3 10 4 (a) (b) Figura 4.6 Valores singulares de la matriz PSD de respuesta para tres mediciones con: (a) Un modo en 1 Hz. Las dos líneas inferiores ( s2 2 y s2 3 ) son prácticamente constantes e iguales al ruido de fondo, 1 (µg) 2 /Hz. (b) Dos modos en 1 Hz y en 1,05 Hz. s2 1 viene dominado por el modo 1 por debajo de 1 Hz y por el modo 2 por encima de 1,05 Hz, mientras que para s2 2ocurre lo contrario. s2 3representa el ruido de fondo. Tomada de [12]. Si solo hay un modo dominante en torno a la frecuencia ωm (asociada a un pico de resonancia), es decir, solo es significativo un término en la Ecuación 4.17, la matriz PSD de la salida se aproxima a una matriz de rango unidad: Gx(ω) = s2 1u1uH 1,ω−→ωm(4.55) En tal caso, el primer vector singular u1es una estimación del modo de vibración m: 4.5 Frequency Domain Decomposition (FDD) 51 ˆ φm=u1(ωm)(4.56) En cuanto al valor singular s2 1 , representa la función de densidad espectral del sistema de un grado de libertad equivalente correspondiente al modo de interés. La función PSD del sistema SDOF equivalente se identifica como el conjunto de valores singulares que están caracterizados por un vector singular parecido en torno a un determinado pico. Este proceso se realiza comparando la estimación del modo de vibración ˆ φm en el pico de resonancia con los vectores singulares uj asociados a los valores singulares significativos en torno a dicho pico. La métrica para evaluar la correlación entre los vectores modales se denomina MAC (Modal Assurance Criterion), y se calcula de la siguiente forma: MAC(uj,ˆ φm) = |uH jˆ φm|2 (uH juj)( ˆ φH mˆ φm)(4.57) El MAC toma valores entre 0 y 1, siendo igual a 0 cuando los dos vectores comparados son ortogonales, y 1 cuando son proporcionales. Un valor singular se encuentra asociado al modo m en una frecuencia particular si su vector singular asociado uj posee una correlación con la estimación del modo ˆ φmsupera un umbral de MAC, típicamente sobre 0,8. La función PSD del sistema SDOF equivalente, identificada mediante el procedimiento anterior, se puede emplear para evaluar el ratio de amortiguamiento modal a través de la versión avanzada del método, denominada EFDD (Enhanced Frequency Domain Decomposition). Esto se lleva a cabo aplicando la transformada inversa de Fourier a dicha función, obteniendo la función de correlación aproximada del sistema SDOF equivalente. El ratio de amortiguamiento se puede obtener mediante regresión lineal en el decremento logarítmico de la función (Ecuación 4.58) [34], o mediante el ajuste exponencial de la envolvente. δ=log x1 x2 =2π ξ p1−ξ2⇒ξ=δ √δ2+4π2(4.58) Adicionalmente, otra manera de obtener la frecuencia natural del modo es, a su vez, mediante regresión lineal de los tiempos en los que la función cruza el cero, teniendo en cuenta la siguiente relación entre frecuencia natural y frecuencia natural amortiguada [24,41]: fd,m=fn,mq1−ξ2 m(4.59) En la Figura 4.7 se puede observar esquemáticamente la implementación del algoritmo descrito en esta sección. 52 Capítulo 4. Identificación experimental offline Descomposición en valores singulares Figura 4.7 Esquema general del algoritmo FDD para el caso de dos canales de medidas. 4.6 Covariance-driven Stochastic Subspace Identification (SSI-COV) 53 4.6 Covariance-driven Stochastic Subspace Identification (SSI-COV) La identificación de espacio de estados estocástico basada en covarianzas, SSI-COV, es un método paramétrico en el dominio del tiempo que tiene como objetivo la identificación de un modelo estocástico de espacio de estados del sistema a partir de datos de su respuesta, concretamente mediante la matriz de correlaciones de la salida. 4.6.1 Base del método La técnica parte con el cálculo de las matrices de correlación de salida, Ri , de acuerdo con la Ecuación 4.37. Dado que en casos prácticos solo se dispone de una secuencia finita de muestras de la respuesta, se emplea la siguiente expresión: ˆ Ri=1 Nt−1 Nt−i−1 ∑ k=0 xk+ixT k(4.60) donde ˆ Ri constituye una estimación insesgada de la matriz de correlación de las salidas con un desfase de tiempo positivo i·∆t , y Nt es el número de muestras de la secuencia. Cada matriz tiene una dimensión de Nx×Nx , siendo Nx el número de salidas. A continuación, se agrupan las correlaciones con desfases entre 1·∆ty(2jb−1)·∆ten una matriz Toeplitz por bloques, tal que: T1|jb=     ˆ Rjbˆ Rjb−1··· ˆ R1 ˆ Rjb+1ˆ Rjb··· ˆ R2 . . .. . ..... . . ˆ R2jb−1ˆ R2jb−2··· ˆ Rjb      (4.61) Dicha matriz tiene un tamaño de Nxjb×Nxjb , siendo jb el número de matrices de correlación por fila y por columna de la matriz por bloques, parámetro que se puede ajustar libremente para obtener resultados óptimos. Aplicando la propiedad de factorización de Ri descrita en la Ecuación 4.43 (Ri=CAi−1G) a la matriz de Toeplitz: T1|jb=     C CA . . . CAjb−1     ·Ajb−1G··· AG G=OjbΓjb(4.62) Cabe recordar que A , C y G son, respectivamente, las matrices de estado, de salida y de covarianza nuevo estado-salida del sistema. La matriz de observabilidad se define como: Ojb=     C CA . . . CAjb−1      (4.63) mientras que la matriz de controlabilidad inversa: Γjb=Ajb−1G··· AG G(4.64) 54 Capítulo 4. Identificación experimental offline Como se puede observar, la primera de ellas es una columna de jb bloques de dimensiones Nx×n , siendo n el orden del sistema. La matriz Γjb , en cambio, es una fila de jb bloques de dimensiones n×Nx . En caso de que Nxjb≥n , y dado que la matriz de Toeplitz está formada por el producto de una matriz de ncolumnas por una de nfilas, su rango es igual a n. Aplicando la descomposición en valores singulares (SVD) a la matriz de Toeplitz: T1|jb=USVT=U1U2·S10 0 0·VT 1 VT 2(4.65) se obtiene el rango de la misma, que coincide con el número de valores singulares no nulos. Omitiendo los valores nulos: T1|jb=OjbΓjb=U1S1VT 1(4.66) Las matrices U1 y VT 1 tienen dimensiones Nxjb×n y n×Nxjb , respectivamente, mientras que la matriz S1 contiene los n valores singulares positivos en orden descendiente. Dividiendo la SVD en dos partes, se pueden calcular las matrices OjbyΓjbcomo: Ojb=U1S1/2 1(4.67) Γjb=S1/2 1VT 1(4.68) Teniendo en cuenta la Ecuación 4.63 y la Ecuación 4.64, la matriz C puede ser estimada a partir de las Nx primeras filas de la matriz de observaciones, mientras que la matriz G se obtiene con las últimas Nx columnas de la matriz de controlabilidad inversa. El procedimiento para obtener la matriz de estado Ase basa en la resolución de un problema de mínimos cuadrados:      C CA . . . CAjb−2     ·A=     CA CA2 . . . CAjb−1     ⇒A=     C CA . . . CAjb−2      + ·     CA CA2 . . . CAjb−1      = (O jb)+O jb(4.69) El superíndice [−]+ representa la pseudoinversa de Moore-Penrose de una matriz. Las matrices O jb y O jb se obtienen eliminando las últimas y las primeras Nx filas de Ojb , respectivamente. Una vez estimadas las matrices A y C , se pueden extraer las propiedades modales. Los autovalores de la matriz A dan los pares conjugados complejos de polos discretos µm del sistema. Para convertir dichos polos de tiempo discreto a continuo, se emplea la siguiente relación: λm=log µm ∆t(4.70) A continuación, los autovalores con parte imaginaria positiva se emplean para obtener las frecuencias naturales fn,m y los ratios de amortiguamiento modales ξm del sistema, mediante las expresiones: fn,m=|λm| 2π(4.71) 4.6 Covariance-driven Stochastic Subspace Identification (SSI-COV) 55 ξm=−Re(λm) |λm|(4.72) Por último, los modos de vibración del sistema, Φ , se pueden calcular a través de los autovectores de A,Ψ, que son observables, es decir: Φ=CΨ(4.73) Las matrices de correlación de las salidas y las matrices del modelo de espacio de estados son estimaciones de los valores reales de dichas cantidades, al estar basadas en un número finito de muestras. La presencia de ruido puede acarrear errores en las estimaciones, siendo típicamente sus principales fuentes: • Imprecisiones en el modelado. Se refiere, por ejemplo, al caso en el que el sistema real no puede ser representado con un modelo de espacio de estados estocástico. •Errores de medición, debidos a los sensores y al sistema de adquisición de datos en general. •Ruido computacional, debido a la precisión del hardware de análisis. •El número finito de muestras de salidas. Estos errores provocan que, en lugar de tener la SVD de la matriz de Toeplitz tantos valores singulares no nulos como el orden n del sistema, todos los valores singulares sean distintos de cero. Por tanto, son necesarias medidas para identificar el orden real del sistema. 4.6.2 Diagramas de estabilización. Automatización de OMA Como se ha mencionado, el método SSI y, en general, todos los métodos OMA paramétricos presentan la desventaja de que requieren un conocimiento previo del orden aproximado del sistema, por ejemplo, observando la gráfica de valores singulares de la matriz PSD de la salida, aplicando la técnica FDD y obteniendo el número de picos. En la mayoría de aplicaciones prácticas se emplean órdenes del modelo muy superiores al orden físico del sistema real, con el objetivo de asegurar la identificación de todos los modos de interés. Sin embargo, esta sobreestimación del orden genera polos espurios, que han de ser separados tras la identificación de los polos físicos, complicando el análisis. Los polos espurios pueden deberse a: 1. Modos generados por ruido debidos a fuentes físicas, como por ejemplo, los modos del sistema de excitación. 2. Modos generados por el algoritmo matemático para ajustarse al orden superior que se le ha impuesto, debido a ruido e imperfecciones en los datos. La separación entre polos físicos y espurios es comúnmente realizada mediante diagramas de estabilización, que representan las polos detectados para cada orden impuesto en función de la frecuencia. Los criterios para considerar un polo estable se basan generalmente en comparar los polos asociados a un determinado orden con los asociados a un orden inferior. Estas comprobaciones se realizan respecto variaciones relativas de frecuencia natural, amortiguamiento y modo de vibración (MAC), siendo unas condiciones típicas de estabilidad de un polo las siguientes: 56 Capítulo 4. Identificación experimental offline |f(n)−f(n+1)| f(n)<0,01 (4.74) |ξ(n)−ξ(n+1)| ξ(n)<0,05 (4.75) 1−MAC(φ(n),φ(n+1))<0,02 (4.76) Un polo sería, por tanto, considerado estable si cumpliera las tres condiciones anteriores. Los polos físicos se diferencian de los matemáticos en que aparecen como columnas de polos estables (ya comiencen en órdenes bajos o altos), ya que los polos espurios matemáticos no suelen ser estables o, en caso de serlos, tienden a dispersarse. Por otro lado, los polos debidos a modos de ruido se pueden descartar según propiedades como los modos de vibración o el amortiguamiento esperados [24]. La Figura 4.8 muestra un ejemplo de diagrama de estabilización en el método SSI-COV. Figura 4.8 Ejemplo de diagrama de estabilización en la aplicación del método SSI-COV. Uno de los inconvenientes de los diagramas de estabilización, al igual que la selección manual de picos en el FDD, es la necesidad de interactuar con el usuario para realizar la identificación. Por ello, para poder monitorizar los parámetros de una estructura identificándolos en cuanto un conjunto de medidas está disponible, se han desarrollado métodos para procesar todos los datos de manera automática, lo que se denomina OMA automático (Automated Operational Modal Analysis). Uno de ellos es la selección automática de picos en el FDD. Sin embargo, el mayor esfuerzo de investigación está focalizado en la interpretación automática de la información presentada en los diagramas de estabilización, tanto para la distinción de modos espurios de modos físicos, como para la agrupación de los modos físicos de diferentes órdenes en un solo modo. El objetivo es desarrollar algoritmos que imiten las decisiones que toma un analista 4.6 Covariance-driven Stochastic Subspace Identification (SSI-COV) 57 experimentado en el proceso de examinación. Una de las aproximaciones más empleadas para dicha automatización son los clusters, cuyo uso consiste en crear grupos de modos que tienen propiedades como frecuencias o amortiguamientos que son similares y que, por tanto, se pueden asociar a un mismo modo físico. Los métodos de clusters pueden clasificarse en jerárquicos y no jerárquicos. Un desarrollo más amplio con aplicaciones prácticas se puede encontrar en [35]. La Figura 4.9 muestra un esquema de la aplicación del algoritmo SSI-COV descrito en la presente sección. 58 Capítulo 4. Identificación experimental offline Medidas de la respuesta Cálculo de la matriz de Toeplitz Matriz de correlación: Matriz de Toeplitz por bloques: Descomposición en valores singulares Diagrama de estabilización Modos de vibración M1 M2 M3 Figura 4.9 Esquema general del algoritmo SSI-COV para el caso de dos canales de medidas. 4.7 Poly-reference Least Squares Complex Frequency Domain (p-LSCF) o PolyMAX 59 4.7 Poly-reference Least Squares Complex Frequency Domain (p-LSCF) o PolyMAX El método poly-reference Least Squares Complex Frequency Domain (p-LSCF), también conocido por su nombre comercial PolyMAX, es una técnica paramétrica en el dominio de la frecuencia creada hace dos décadas para identificar propiedades modales a partir de matrices FRF. Dada la similitud entre la descomposición modal de la matriz FRF (Ecuación 4.9) y de la matriz de densidad espectral positiva de la respuesta (Ecuación 4.16) en el caso de un sistema excitado por ruido blanco, el método se pudo adaptar de manera significativamente directa a análisis en los que solo se dispone de la información sobre la respuesta (OMA), introduciéndose dicho uso en 2005 [42]. 4.7.1 Ajuste del modelo Los modelos modales introducidos en la Subsección 4.3.1 son útiles para entender el problema y el significado de las matrices desde un punto de vista físico. Sin embargo, no son adecuados para realizar un ajuste con los datos experimentales obtenidos de los ensayos. Por ello, este método emplea para describir la matriz PSD un modelo de fracción matricial en el dominio de la frecuencia en tiempo discreto (dominio z ), el cual se basa en la división de dos matrices polinómicas A y B (no confundir con las matrices AyBdel modelo de espacio de estados), tal que: Gxx(zj) = B(zj)A−1(zj)(4.77) Teniendo en cuenta que la variable z se puede expresar como z=es∆t=eiω∆t (si se restringe en dominio de s al eje imaginario, s=iω , como se ha hecho hasta ahora), y tomando un mismo orden pde polinomios para denominador y numerador, AyB, ambas matrices se definen como: A(zj) = p ∑ r=0 Arzr j= p ∑ r=0 Areiωj∆t r =A(ωj)(4.78) B(zj) = p ∑ r=0 Brzr j= p ∑ r=0 Breiωj∆t r =B(ωj)(4.79) donde Ar y Br son las matrices reales de coeficientes de los polinomios que es necesario estimar, puesto que contienen la información de los parámetros del sistema. Las dimensiones de la matriz de medio espectro Gxx son Nx×Nre f , donde Nx es el número total de salidas medidas del sistema y Nre f es el número de salidas de referencia ( Nre f ≤Nx ). Por tanto, las matrices A y B tendrán dimensiones Nre f ×Nre f yNx×Nre f , respectivamente. Para llevar a cabo la estimación comentada de Ar y Br , se trata de minimizar el error entre la matriz PSD que se obtiene por el modelo teórico anterior, y la matriz PSD estimada a través de las medidas experimentales de la respuesta del sistema, ˆ Gxx(ωj): ENLS(ωj) = p ∑ r=0 Breiωj∆t r p ∑ r=0 Areiωj∆t r−1 −ˆ Gx(ωj)(4.80) Dado que la anterior ecuación para los errores a ser minimizados implica la resolución de un problema de mínimos cuadrados no lineal (NLS), se propone la definición de los errores de una forma alternativa para simplificar el problema a uno lineal. Multiplicando por la derecha por A(ωj) : 5 Identificación experimental online 5.1 Introducción Este último capítulo de la parte teórica del trabajo va a tratar una segunda tipología de identificación estructural, la identificación online o en tiempo real. Como se ha comentado a lo largo del trabajo, durante las últimas décadas se ha incrementado constantemente la importancia de la identificación estructural en áreas como la ingeniería aeroespacial, la mecánica y la civil, con el objetivo de hallar parámetros críticos de la misma como son las frecuencias naturales o los ratios de amortiguamiento. Sin embargo, no ha sido hasta recientemente que la identificación de estructuras no lineales, incluyendo daños del sistema, han ganado atención. Efectos, a simple vista pequeños, que desvían el sistema de la linealidad, como tales daños, pueden ocasionar un deterioro grave en el comportamiento dinámico de la estructura. Por ello, han surgido dos tipologías de métodos de identificación, offline yonline. La primera de ellas, tratada en el capítulo anterior, requiere un conjunto completo de datos de la respuesta del sistema para poder realizar la identificación, la cual tiene lugar a posteriori. Los métodos online, en cambio, son necesarios para identificar cambios no lineales en la estructura en tiempo real. Por ello, se denominan también técnicas recursivas, ya que la estimación se realiza paso a paso basándose en estados e información del sistema a priori [45]. A diferencia de lo que ocurre para los métodos de identificación offline, los online únicamente emplean análisis en el dominio del tiempo, debido a que las técnicas frecuenciales son, generalmente, no recursivas. Uno de los ejemplos más destacados es la estimación por mínimos cuadrados o LSE (Least Squares Estimation) [46]. La identificación de parámetros se realiza en este caso minimizando el error cuadrático entre la medida tomada experimentalmente de la respuesta del sistema, y el valor esperado de la misma. La principal desventaja de este método es, sin embargo, la necesidad de disponer de medidas de los desplazamientos y las velocidades de la estructura, lo cual, en la mayoría de los casos prácticos, presenta una gran dificultad. Como se ha hecho mención en anteriores capítulos, en general las medidas de la respuesta dinámica se toman mediante acelerómetros, y aunque las señales de aceleración pueden ser integradas para la obtención de velocidades y desplazamientos, dicha integración conlleva errores y ruido adicional. Otra técnica de identificación recursiva u online es el método secuencial de Monte Carlo o SMC (Sequential Monte Carlo, también denominado filtro de partículas o PF (Particle Filter) [47]. Dado que se trata de un método bayesiano, emplea una serie de muestras temporales (partículas) para representar la probabilidad de que se de un determinado estado o parámetro en el futuro. En el caso de sistemas complejos, el número de muestras requeridas aumenta considerablemente, provocando 67 68 Capítulo 5. Identificación experimental online a su vez que el coste computacional crezca, lo que supone un inconveniente considerable en su aplicación [13]. Además de los anteriores métodos, uno de los más populares para identificación de sistemas es el filtro de Kalman o KF (Kalman Filter). Este permite una estimación de tanto el estado del sistema como de parámetros del mismo. Sin embargo, sólo puede tener en cuenta comportamientos lineales, propiciando el empleo de su evolución, el filtro de Kalman extendido o EKF (Extended Kalman Filter). Su base es la misma que la del KF, solo que implementa en cada paso de tiempo una linealización de las ecuaciones no lineales de estado y observación, a través del cálculo del Jacobiano. Esto último supone uno de sus inconvenientes, ya que el cálculo de dicha matriz en cada paso conlleva un elevado coste computacional. Por otro lado, para sistemas cuyo comportamiento se desvía en gran medida de la linealidad, las aproximaciones realizadas son poco precisas. Por ello, Julier et al. [48] propuso en 1995 un tercer desarrollo, el filtro de Kalman unscented o UKF (Unscented Kalman Filter). Este método corrige los inconvenientes del EKF, ya que, en lugar de la linealización de las ecuaciones, emplea una serie de puntos de muestra que son propagados en el tiempo mediante una transformación unscented (UT) con las ecuaciones no lineales, reduciendo su coste computacional [13,45]. A pesar de la capacidad del UKF para identificar sistemas no lineales, las variaciones abruptas de parámetros no pueden ser cubiertas por algoritmos tradicionales en tiempo real. Por ello, es necesario un enfoque mediante un observador adaptativo para permitir el funcionamiento del método tanto antes como tras el cambio [13]. Una de las aplicaciones de este tipo de identificación es durante el empleo de amortiguadores de masa sincronizada o TMD (Tuned Mass Dampers), dispositivos de control de estructuras que generan amortiguamiento y fuerzas de recuperación adicionales en las mismas. Los sistemas con TMD presentan retos adicionales debido a la precisión que requieren a la hora de identificar cambios en las propiedades de las estructuras, ya que estos pueden deteriorar gravemente su actuación [45]. A continuación, se van a pasar a describir los algoritmos del filtro de Kalman y el filtro de Kalman unscented, así como la incorporación del observador adaptativo para la estimación conjunta del estado y los parámetros de un sistema. 5.2 Filtro de Kalman El filtro de Kalman es una técnica recursiva de aproximación gaussiana en el dominio del tiempo que permite estimar el estado, por ejemplo, desplazamientos y velocidades, de un sistema no lineal ante cualquier tipo de excitación. Para ello, pondera la diferencia entre la respuesta medida y la respuesta estimada con la denominada ganancia de Kalman, que se calcula a través de las covarianzas del error de estimación, del error de mediciones y del modelo [45]. La base del filtro de Kalman parte del modelo determinista-estocástico de espacio de estados en tiempo discreto de un sistema lineal, Ecuación 4.31. Los ruidos de proceso y de medición, w y v , se suponen, al igual que en el capítulo anterior, procesos estacionarios de ruido blanco gaussiano con media nula. Adicionalmente, se consideran independientes entre sí, por lo que sus matrices de covarianza presentan las siguientes propiedades [13]: Ewp vpwT qvT q=     Q 0 0 Rp=q 0p=q (5.1) 5.2 Filtro de Kalman 69 donde Q y R , generalmente constantes en el tiempo, son las matrices de covarianza del ruido de proceso y el ruido de medición, respectivamente, y quedan definidas como: Q=E[wkwT k](5.2) R=E[vkvT k](5.3) El KF consta de una fase de innovación y una fase de predicción del estado, que se ejecutan secuencialmente. Durante la primera de ellas, la salida del sistema viene dada por (denotando ˆ· estimación a priori o predicción y ˜·estimación a posteriori o corrección): ˆ yk=Cˆ xk+Dˆ uk(5.4) siendo la estimación a priori del estado, ˆ xk , y ˆ uk conocidos. A continuación, se actualiza la estimación del estado a una estimación a posteriori ˜ xk , a partir de ˆ xk y de la diferencia o error entre la estimación de la salida ˆ yky la salida medida yk, a través de una ponderación con la ganancia de Kalman Kk: ˜ xk=ˆ xk+Kk(yk−ˆ yk)(5.5) Tal diferencia entre la salida medida y la estimada se conoce comúnmente como error de innovación, que se define por tanto como: ek=yk−ˆ yk(5.6) La expresión para el cálculo de la ganancia de Kalman se obtiene a través de la minimización del error cuadrático medio de estimación del estado. Para ello, se emplean la matriz de covarianza del estado ˆ Pk , la matriz de covarianza del ruido de medición R y la matriz de salida C , dando lugar a la siguiente fórmula: Kk=ˆ PkCT(Cˆ PkCT+R)−1(5.7) Durante el proceso de innovación, también tiene lugar la actualización de la matriz de covarianzas del estado, tal que: ˜ Pk= (I−KkC)ˆ Pk(5.8) Por otro lado, en la fase de predicción se calculan tanto la estimación a priori del estado para el siguiente instante (k+1): ˆ xk+1=A˜ xk+Buk(5.9) como la estimación a priori de la matriz de covarianzas del estado: ˆ Pk+1=A˜ PkAT+Q(5.10) 70 Capítulo 5. Identificación experimental online En la Figura 5.1 se puede observar esquemáticamente la arquitectura de un filtro de Kalman. Cabe remarcar que la ganancia de Kalman toma la vital función de distribuir el peso en la corrección del estado entre la predicción y la medición de la salida. Si la incertidumbre de la predicción a través de la ecuación de estado es alta, la ganancia de Kalman será mayor, dándole más importancia al término de la medición (error de innovación) en la Ecuación 5.5. Por lo contrario, en caso de que la incertidumbre de la medición sea mayor, el valor de Kk será bajo y, por tanto, tendrá más peso la predicción a la hora de corregir el estado. Salida medida Estado corregido Fase de predicción 1) Predicción del estado 2) Predicción de la matriz P: Fase de innovación 1) Cálculo de la estimación de la salida: 2) Corrección del estado: 3) Corrección de la matriz P: Figura 5.1 Diagrama de las fases de la implementación de un filtro de Kalman. Adaptado de [13]. 5.3 Filtro de Kalman unscented Tal y como se comentó durante la sección introductoria del capítulo, el filtro de Kalman no es suficiente para poder estimar el estado de un sistema con comportamiento no lineal. Para poder incorporar no linealidades sin perder la eficiencia computacional del KF, se emplea la evolución del filtro de Kalman unscented. Este método crea una serie de puntos de muestra de forma determinista, que representan estadísticamente el estado y los parámetros del sistema. Posteriormente, dichos puntos son propagados mediante las ecuaciones no lineales del sistema, a través de un proceso denominado transformación unscented o UT, de tal manera que los puntos resultantes también representan las propiedades estadísticas del sistema correctamente. El fundamento del UKF parte de las ecuaciones del modelo determinista-estocástico de espacio de estados de un sistema no lineal en tiempo discreto: xk+1=f(xk,uk)+ wk yk+1=h(xk,uk)+ vk (5.11) donde fyhson, respectivamente, las funciones discretas de transición y de salida no lineales. 5.3 Filtro de Kalman unscented 71 Al igual que el filtro de Kalman lineal, el UKF se divide en una fase de innovación y otra de predicción. La primera comienza con la creación de la serie de puntosσ que se aproximan a la distribución gaussiana, es decir, representan fielmente las propiedades estadísticas, del estado en el instante k: ˆ σk=ˆ X0,k··· ˆ Xi,k··· ˆ X2n,kW(m) 0W(c) 0···W(m) iW(c) i···W(m) 2nW(c) 2n(5.12) donde n es la dimensión del vector de estado, ˆ Xi,k es el i -ésimo puntoσ , W(m) i es el factor de peso correspondiente para el cálculo de la media y W(c) i para el cálculo de la covarianza. Como se puede observar, hay un total de 2n+1 puntos en la distribución, estando situado ˆ X0,k en el centro de la misma y los restantes 2n puntosσ alrededor, según el valor de la matriz ˆ Pk , de manera equidistante. Los puntos-σy sus correspondientes pesos se calculan de la siguiente forma: ˆ X0,k=ˆ xk,ˆ Xi,k=         ˆ xk+αq(n+λ)ˆ Pki i=1, ... ,n ˆ xk−αq(n+λ)ˆ Pki−n i=n+1, ... ,2n (5.13)        W(m) 0=λ n+λ W(c) 0=λ n+λ+1−α2+β W(m) i=W(c) i=1 2(n+λ)i=1, ... ,2n(5.14) La expresión [−]i se refiere a la i -ésima columna de la matriz interior a los corchetes. La estimación a priori o predicción del estado ˆ xk se supone conocida para el primer instante de tiempo. El parámetro λviene definido por: λ=α2(n+κ)−n(5.15) siendo α un parámetro que, como se puede apreciar en la Ecuación 5.13, escala la distancia de los puntosσ al centro de la distribución, al que generalmente se le da un valor comprendido entre 10−4 y 1. Por su parte, κ se define óptimamente como κ=n−3 , mientras que β=2 para el caso de suposición de distribución gaussiana. Ambos son parámetros de ajuste de la transformación. A continuación, se pueden calcular los puntosσ de la salida del sistema a través de la UT, propagando los puntos-σdel estado con la función de salida no lineal: ˆ Yi,k=hˆ Xi,k,uk(5.16) Dado que se dispone de los pesos de cada uno de los puntosσ en las medias y covarianzas, se puede calcular, en primer lugar, la predicción de la salida del sistema como la media de los puntosσ de la salida: ˆ yk= 2n ∑ i=0 W(m) iˆ Yi,k(5.17) 72 Capítulo 5. Identificación experimental online Por otro lado, se calculan la matriz de covarianza de las medidas o salidas ˆ Pyy,k y la matriz de covarianza cruzada ˆ Pxy,k: ˆ Pyy,k= 2n ∑ i=0 W(c) iˆ Yi,k−ˆ ykˆ Yi,k−ˆ ykT+R(5.18) ˆ Pxy,k= 2n ∑ i=0 W(c) iˆ Xi,k−ˆ xkˆ Yi,k−ˆ ykT(5.19) La Figura 5.2 muestra, para una dimensión del estado de n=2 , el proceso de propagación de puntosσ para la obtención de la estimación de la salida ˆ yk a partir de la predicción del estado ˆ xk y de la matriz de covarianzas de dicho estado ˆ Pk . Además, se pueden comparar los resultados de la UT con métodos de muestreo generales como el filtro de partículas, en el que se obtienen medias y covarianzas muy precisas, pero se requiere de una gran nube de puntos. Se aprecia que los resultados de covarianza y media para la UT son muy similares a los reales, siendo necesarios únicamente cinco puntos, lo que conlleva un importante ahorro de coste computacional [49]. Muestreo real ESTADO SALIDA Transformación unscented (UT) ESTADO SALIDA Media Covarianza Covarianza del estado Puntos Sigma Covarianza de la salida Media real Covarianza real Predicción del estado Predicción de la salida Figura 5.2 Diagrama de la propagación de media y covarianza, en caso de muestreo real y en caso de UT. A través de las matrices de covarianza de las salidas y cruzada, se puede obtener la ganancia de Kalman Kka partir de la expresión: 5.4 Identificación conjunta de estado y parámetros del sistema 73 Kk=ˆ Pxy,kˆ P−1 yy,k(5.20) A partir de ella, es posible el cálculo de la corrección del estado ˜ xk una vez se dispone de las medidas de la respuesta yk: ˜ xk=ˆ xk+Kk(yk−ˆ yk)(5.21) A su vez, la matriz de covarianzas del estado se corrige de acuerdo con la siguiente expresión, dando lugar al final de la fase de innovación del filtro: ˜ Pk=ˆ Pk−Kkˆ Pyy,kKT k(5.22) En cuanto a la fase de predicción del UKF, esta comienza con la actualización de los puntosσ , empleando las Ecuaciones 5.13 y 5.14 y sustituyendo las estimaciones a priori del estado y la matriz de covarianzas ˆ xk y ˆ Pk por sus valores a posteriori ˜ xk y ˜ Pk . Una vez se dispone del conjunto de puntos corregidos ˜ Xi,k , se propagan a través de la función de transición no lineal f para obtener los puntos-σde la predicción del siguiente estado: ˆ Xi,k+1=f˜ Xi,k,uk(5.23) La predicción del estado ˆ xk+1 se obtiene a través de la siguiente media ponderada de los puntosσ : ˆ xk+1= 2n ∑ i=0 W(m) iˆ Xi,k+1(5.24) Para finalizar, la matriz de covarianzas del estado se predice como: ˆ Pk+1= 2n ∑ i=0 W(c) iˆ Xi,k+1−ˆ xk+1ˆ Xi,k+1−ˆ xk+1T+Q(5.25) En la Figura 5.3 se muestra un esquema de la arquitectura del filtro UKF. 5.4 Identificación conjunta de estado y parámetros del sistema En esta sección, se va a introducir el empleo de un observador adaptativo para la estimación tanto del estado como de parámetros del sistema, teniendo el mismo la capacidad de detectar cambios bruscos en las propiedades del mismo, como las frecuencias naturales. Como se ha comentado en el apartado anterior, el filtro de Kalman unscented es capaz de estimar el estado de sistemas no lineales. Sin embargo, para la estimación de parámetros, el vector de estado xk , que contiene desplazamientos, velocidades o aceleraciones, ha de ser ampliado con un vector de parámetros θk que contenga los parámetros a identificar. De esta manera, las ecuaciones del modelo de espacio de estados han de ser modificadas para operar con un vector de estado aumentado xa k definido como: 74 Capítulo 5. Identificación experimental online Salida medida Estado corregido Fase de innovación 1) UT: 2) Cálculo de las predicciones de salida, correlación de salidas y correlación cruzada: 3) Corrección del estado: 4) Corrección de la matriz P: Fase de predicción 1) UT: 2) Cálculo de las predicciones de estado y correlación de estado: Figura 5.3 Diagrama de las fases de la implementación de un filtro de Kalman unscented. xa k=xk θk(5.26) Dado que la dimensión del vector de estado es 2n , la dimensión de xa k será de 2n+s , donde s es el número de parámetros del sistema a identificar. La base de la arquitectura adaptativa a emplear es la actualización de la matriz corregida de covarianzas del estado, ˜ Pk . Dicha matriz es indicativa del nivel de confianza que se tiene en la estimación a posteriori del estado, ˜ xk , y se puede utilizar para influir en las estimaciones futuras de parámetros, ya que un alto valor de covarianza es señal de un sistema muy sensible, haciendo posible la identificación de variaciones bruscas por no linealidades. Con el anterior objetivo, se define un parámetro disparador γ basado en el error de innovación ek+1normalizado con la matriz de covarianza del ruido de medición R: γ=eT k+1R−1ek+1(5.27) En el caso en el que se cuenta con Nx sensores con covarianzas del ruido de medición idénticas y constantes Ri=R, la anterior expresión se puede desarrollar como: γ=   e1 . . . eNx    T   R1··· 0 . . ..... . . 0··· RNx    −1   e1 . . . eNx   =e2 1 R1 +...+e2 Nx RNx =e2 1+...+e2 Nx R(5.28) 5.4 Identificación conjunta de estado y parámetros del sistema 75 Cuando el parámetro disparador γ supera un umbral γ0 , se procede a realizar la adaptación para detectar cambios en el sistema: γ≥γ0⇒adaptación (5.29) El error de innovación econtiene tanto errores del sistema como errores de medidas: Error de innovación =Error del sistema+Error de medidas Cambios en los parámetros del sistema se corresponden con errores del sistema. Por ello, el umbral γ0 establecido debe de cubrir la porción del error de innovación debida a mediciones, para que cada vez que se supere dicho umbral y se proceda con la adaptación, sea porque ha ocurrido un cambio en el sistema. Para la definición del umbral, se considera en primer lugar un sistema lineal sin cambios, es decir, el error de innovación se debe completamente a errores de medición: ei=yi−ˆyi=yi−(Cˆ x)i−(Du)i(5.30) Tal y como se discutió en la Subsección 4.3.2, el ruido de medición incluye tanto errores de medición de la salida yi , como errores de medición de la entrada u . En el caso de empleo de acelerómetros, provocando que la entrada al sistema tenga efecto directo sobre la salida ( D=I ), el error de innovación ei se puede aproximar como la diferencia entre el error de medición de la salida nyy el error de medición de la entrada nu: ei≈ny−nu(5.31) Considerando ambos ruidos como gaussianos y de media nula, con varianza Ri=R , el error de innovación será también gaussiano con media nula y una varianza acumulada de 2Ri=2R . Definiendo la variable aleatoria ze,icomo: ze,i=ei √2R(5.32) y sustituyendo en la Ecuación 5.28, el umbral γviene dado por la siguiente expresión: γ0≈2Nxz2 0(5.33) Puesto que Nx depende de la configuración del sistema, el parámetro z0 determina la probabilidad de que se supere el umbral debido exclusivamente a errores de medición, asumiendo que el sistema no ha sufrido cambios y, por tanto, no es una fuente de error. Tras la detección de una variación en el sistema, se ha de continuar con un algoritmo de localización del parámetro que causa el cambio. Para ello, cabe describir la matriz corregida de correlación del estado como: ˜ Pk=˜ Pxx,k˜ Pxθ,k ˜ Pθx,k˜ Pθθ,k(5.34) 82 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala Por otro lado, la cuerda del modelo se ha seleccionado de 25 cm, con el objetivo de disponer de margen respecto a las máximas dimensiones en planta de la impresora ( 30 ×30 cm). Finalmente, el grosor tm del PLA se ha elegido de 2 mm para evitar roturas durante la fabricación, a la vez que el porcentaje de infill o relleno en la impresión ha sido del 70%. En cuanto a las dimensiones del prototipo, se conoce tanto la envergadura b como el alargamiento Λ=b2/Sdel mismo: Tabla 6.1 Datos conocidos del avión estudiado. Fuente: Jane’s All the World’s Aircraft 2004-2005 [51]. b(m) 11,53 m Λ(-) 7,6 A partir de la envergadura y el alargamiento, se pueden obtener tanto la longitud del semiala, como la cuerda media de la misma (suponiendo que el estrechamiento del ala se puede considerar unidad): Lp=b 2=5,7671 m cp=S b=b Λ=1,5177 m (6.6) De esta manera, todos los factores de escala de la Ecuación 6.5 quedan determinados, a excepción de λL, conociendo las propiedades de los materiales (ver Tabla 6.2). Tabla 6.2 Propiedades de los materiales en prototipo y modelo. Fuentes: MatWeb, Xakata 3D2. Propiedad Aluminio Acero PLA (70% relleno) E(GPa) 70,3 200 1,85 ρ(kg/m3) 2810 7800 1620 Calculando los distintos parámetros de las dos secciones definidas, se obtienen factores de escala mostrados en la Tabla 6.3. Tabla 6.3 Factores de escala obtenidos mediante el dimensionamiento de las secciones. λEI 6263,2 λρA6,5223 Si, además, se impone igualdad de frecuencias naturales como se ha comentado anteriormente: λL=λEI λρAλf1/4 =5,5667 ⇒Lm=Lp λL =1,0356 m (6.7) 2Matweb: https:// www.matweb.com/ index.aspx Xakata 3D PLA 850: https:// www.sakata3d.com/ descargables/ TDS-PLA-850.pdf 6.2 Creación del modelo a escala reducida 83 Así, las dimensiones finales del modelo a fabricar vienen resumidas en la Tabla 6.4. Tabla 6.4 Resumen de parámetros del modelo a escala reducida. L(m) 1,036 cm(m) 0,25 tm(mm) 2 b(mm) 30 h(mm) 10 em(m) 0,124 Aplicando la Ecuación 6.2 a ambos diseños, se comprueba que se obtienen valores iguales para la primera frecuencia natural, f1=6,50 Hz. Como aparece en referencias como [52,53], este valor es común en aeronaves de tales alargamientos. 6.2.2 Modelado 3D y fabricación La impresora 3D con la que se ha llevado a cabo la fabricación del modelo es la Creality CR 10V3 3 (ver Figura 6.3), con extrusora directa. A su vez, el material empleado ha sido el PLA 850 4 de Xakata 3D, cuya impresión se realiza a una temperatura de boquilla de unos 200 ºC. Los filamentos de PLA se caracterizan por unas buenas propiedades de resistencia y rigidez, a la vez que mantienen un bajo coste, de unos 26 €/kg, cuando se comparan con otros materiales de impresión 3D del mercado. Figura 6.3 Impresora 3D empleada para la fabricación. La primera consideración a tener en cuenta es la restricción que presenta la impresora en cuanto a dimensiones la fabricación, ya que el tamaño de impresión es de 30 ×30 ×40 mm, siendo esta última medida la dirección de la extrusión (vertical). Por tanto, a pesar de que la cuerda del modelo no supone un problema a la hora de fabricar, no es posible imprimir todo el modelo de una sola vez a causa de su longitud. Debido a esto, se han diseñado distintas partes o módulos de 20 cm que 3 Creality CR 10V3: https:// www.creality3dofficial.com/ es/ products/ creality-cr-10-v3-3d-printer-with-genuine-e3d- direct-drive-extruder-2020-latest-version 4PLA 850 Xakata 3D: https:// www.sakata3d.com/ es/ content/ 14-pla-850 84 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala al unirse a la barra de acero mediante tornillos, y entre ellas mediante epoxi, conformarán el ala completa. El ensamblaje de 5 módulos permite alcanzar la longitud de 1 m, mientras que los 3,6 cm restantes para igualar Lmse dejarán de margen para ajustar el empotramiento de la barra. El software empleado para el diseño 3D ha sido CATIA V5, donde cada uno de los módulos del ala se ha dividido en una cubierta superior y una inferior para poder ser unidas de manera sencilla mediante tornillos a la barra de acero. Las Figuras 6.4, 6.5 y 6.6 muestran distintas fotos de este modelado. Una vez creados los modelos CAD, se procede a exportarlos en formato .stl para ser procesados por el programa UltiMaker Cura. Esta herramienta se encarga de dividir el modelo 3D en una serie de capas en dos dimensiones sobre las que la impresora puede mover su cabezal y extrusor, y llevar a cabo la fabricación capa por capa. De esta manera, tiene el papel de generar el G-code que máquina es capaz de interpretar. Además, permite configurar los distintos parámetros de la impresión, como el porcentaje de relleno, el grosor de las capas, la temperatura del extrusor, la velocidad de impresión, etc. Figura 6.4 Vista superior de ambas cubiertas de un módulo del ala en CATIA V5. Figura 6.5 Vista inferior de ambas cubiertas de un módulo del ala en CATIA V5. Una vez unida y ensamblada el ala, ha sido instalada en el laboratorio del Departamento de Mecánica de Medios Continuos y Teoría de Estructuras de la ETSI. El empotramiento se ha realizado fijando mediante un sargento la parte sobrante de la barra de acero a una viga que se ha considerado de rigidez mucho mayor a la del ala bajo estudio. En las Figuras 6.7, 6.8 y 6.9 se observa el resultado de dicha instalación. 6.2 Creación del modelo a escala reducida 85 Figura 6.6 Vista frontal del módulo en CATIA V5. Figura 6.7 Foto general de la instalación del modelo del ala en el laboratorio. Figura 6.8 Empotramiento de la barra mediante sargento. 86 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala Figura 6.9 Vista frontal del modelo instalado en el laboratorio. 6.3 Sistema de adquisición de datos En esta sección se va describir la arquitectura empleada para recopilar datos de la respuesta dinámica del modelo fabricado. Como se explicó en el Capítulo 3, un sistema DAQ para identificación modal requiere de una estructura física sobre la que realizar los experimentos, unos sensores de movimiento, un sistema de adquisición de datos que convierta la señal analógica en digital, y un sistema de procesamiento de datos para recopilar las medidas y adecuarlas al análisis. Anteriormente, se analizaron las ventajas que presentan los acelerómetros MEMS frente a los acelerómetros piezoeléctricos tradicionalmente empleados en aplicaciones de SHM. En particular, se observó como el precio de los mismos es significativamente más bajo, sin aumentar excesivamente el ruido de fondo. Por ello, se ha optado por el empleo de sensores MPU6050, un modelo de acelerómetro y giróscopo MEMS muy popular en aplicaciones de medida de vibraciones de bajo coste. Este dispositivo tiene incorporado un conversor ADC que digitaliza la señal analógica generada por el movimiento la estructura, lo que elimina la necesidad de uso de un sistema de conversión externo. El sistema de procesamiento de datos elegido consta tanto de una placa con microcontrolador Arduino como de un ordenador con software MATLAB. El modelo de placa empleada ha sido Arduino UNO, que contiene un microcontrolador ATmega328P 5 . Tal y como se comentó en la Sección 3.3, las placas Arduino son de bajo coste, y pueden ser conectadas a la mayoría de sensores, como el MPU6050, a través de pines de entrada y salida. Concretamente, la comunicación entre los elementos tiene lugar mediante un bus I 2 C, siendo la placa la que actúa como maestro en la transmisión de datos. Dado que el bus está vinculado físicamente a los pines A4, para la señal de datos SDA, y A5, para la señal de reloj SCL, es necesario conectar ambos con los pines SDA y SCL de los acelerómetros empleados. El Arduino también se encarga de proporcionar alimentación de 5V a cada sensor. 5Arduino UNO: https:// store.arduino.cc/ products/ arduino-uno-rev3 6.3 Sistema de adquisición de datos 87 La primera restricción que impone el sistema de adquisición en la toma de medidas viene dada por el número de direcciones de I 2 C que puede llevar un MPU6050, que está limitado a dos, 0x68 y 0x69. Dado que una placa Arduino UNO dispone de un solo bus I 2 C, únicamente es posible conectar dos sensores a la misma, ya que cada uno ha de tener asociada una dirección distinta para ser identificado correctamente por el microcontrolador. El esquema de conexiones realizadas entre placa y dos sensores viene mostrado en la Figura 6.12, y la foto real en la Figura 6.13. Como se puede observar, uno de los acelerómetros recibe la alimentación en el pin AD0, en lugar de en el pin VCC. Esto sirve para asignar a dicho MPU6050 la dirección 0x69, en lugar de la 0x68 que posee por defecto. Además, para llevar a cabo todas las conexiones necesarias se ha requerido el uso de una protoboard o placa de pruebas. Por último, también se ha empleado un ordenador portátil con MATLAB para permitir un almacenamiento, visualización y procesamiento más avanzado de los datos adquiridos durante los ensayos. De esta manera, la arquitectura del sistema ha estado basada en la propuesta de la Figura 3.5, y se puede observar de forma más detallada y completa en la Figura 6.15 (una foto del sistema real se puede ver en la Figura 6.14). El flujo seguido por los datos desde que son generados por el sensor hasta que se interpretan en MATLAB es el descrito en la Sección 3.4. Figura 6.10 Acelerómetro y giróscopo MEMS MPU6050. Figura 6.11 Placa de Arduino UNO empleada. 88 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala Figura 6.12 Esquema de conexiones entre Arduino UNO y dos MPU6050. 6.3 Sistema de adquisición de datos 89 Figura 6.13 Foto de las conexiones entre Arduino UNO y dos MPU6050. Figura 6.14 Foto del sistema de adquisición de datos durante los ensayos. 90 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala Medida por: ADC MPU6050 Mov. ADC MPU6050 Ref. Respuesta dinámica Modelo a escala reducida Excitación distribuida Ruido blanco Bus I2C dir. 0x68 dir. 0x69 Placa Arduino UNO Puerto USB Disparador (ISR) Sincronización con TIMER1 Ordenador Puerto COM - Lectura - Almacenamiento Visualización ID offline ID online - Preprocesado - Aplicación método OMA UKF con estimación de parámetros Comunicación serial por USB Figura 6.15 Esquema general del proceso de adquisición de medidas de la respuesta dinámica del sistema. 6.4 Análisis modal de un modelo numérico 91 6.4 Análisis modal de un modelo numérico Con el objetivo de contar con unos resultados de referencia en las primeras fases del desarrollo de la identificación modal, se ha creado un modelo de elementos finitos de la estructura del ala ensayada, para posteriormente ser analizada con ANSYS APDL para la obtención de una estimación numérica de las frecuencias naturales y modos de vibración de la misma. 6.4.1 Creación del modelo de elementos finitos La geometría en la que se ha basado el modelo es la descrita en la Tabla 6.4. Se han empleado dos tipos de elementos, por un lado un elemento tipo chapa, el SHELL181 , y un elemento de volumen, el SOLID185 , ambos comunes en análisis estructurales modales como el que se ha realizado. Para el primero de los elementos, se han definido dos tipos de sección: por un lado, una de 2 mm de espesor para las áreas del revestimiento y, por otro, una de 4 mm de espesor para las áreas coincidentes con la parte superior de la cubierta inferior y la parte inferior de la cubierta superior. Los volúmenes se corresponden con la barra central de acero y con las zonas sólidas de PLA que se sitúan por encima y por debajo de ella. En cuanto a los materiales, se han empleado los módulos elásticos y las densidades de acero y PLA al 70% de relleno definidos en la Tabla 6.2. El tamaño de los elementos se eligió como 5 mm de lado debido a limitaciones de licencia del software. Por último, se ha impuesto una condición de empotramiento en la raíz del ala, a través de la restricción de desplazamientos y giros en los nodos situados en la base de la barra de acero. Figura 6.16 Áreas del modelo geométrico del ala generado para el análisis numérico. 6.4.2 Resultados del análisis modal numérico Una vez creado el modelo de elementos finitos del ala a ensayar, se ha llevado a cabo un análisis modal estimando numéricamente los diez primeros modos de vibración de la misma. Los resultados para las frecuencias naturales se pueden encontrar en la Tabla 6.5. Como se puede apreciar, los modos con una mayor participación en el desplazamiento en la dirección vertical son el 1, 3 y 5, por lo que se van a tomar dichos modos como modos de referencia para compararlos con los métodos experimentales. 98 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala En algunos de los ensayos realizados se ha presenciado la aparición de un modo espurio alrededor de 60 Hz (Figura 6.25), con gran probabilidad producido por máquinas o equipamiento eléctrico cercanos en el laboratorio. Por ello, en dichos casos era necesario aplicar un filtro rechazo banda en torno a esas frecuencias (véase la Figura 6.24) para evitar la detección de dicho modo como modo físico. Figura 6.25 Presencia de un modo espurio en torno a 60 Hz. 6.7 Aplicación de los métodos de identificación offline En esta sección, se van a mostrar los resultados de aplicar cada uno de los métodos OMA desarrollados en el Capítulo 4 para la identificación de los tres primeros modos de flexión del ala del modelo a escala reducida, comparándolos con los obtenidos en el análisis numérico presentado en la sección anterior. En cuanto a la técnica FDD, el algoritmo empleado, bajo la misma idea que el mostrado en la Figura 4.7, ha sido desarrollado por E. Cheynet [54] y su implementación es posible a través de un código de MATLAB 6 compartido en la web de File Exchange de MathWorks. Este algoritmo incorpora una selección automática de los picos del primer valor singular de la matriz PSD de la respuesta, lo que evita la necesidad de interacción con el usuario para la identificación. En lo que respecta al método SSI-COV, el algoritmo implementado ha sido también creado por E.Cheynet [55] y ejecutado en MATLAB mediante el código compartido en la web de MathWorks 7 . Su base es la esquematizada en la Figura 4.9 y, adicionalmente, incorpora una automatización de la agrupación de modos físicos de diferentes órdenes del modelo en modos estructurales de interés, mediante un algoritmo de clusters (inspirado en el descrito en [56]), como se describió en la Subsección 4.6.2. El método p-LSCF ha sido implementado con un algoritmo como el esquematizado en la Figura 4.10 gracias al código desarrollado por R. Trumpp [44] y subido a la plataforma GitHub 8 . Dicho algoritmo permite seleccionar entre emplear un método EMA aportando la función de respuesta en 6 Automated Frequency Domain Decomposition (AFDD) por E. Cheynet: https:// es.mathworks.com/ matlabcentral/ fileexchange/ 57153-automated-frequency-domain-decomposition-afdd 7 Operational modal analysis with automated SSI-COV algorithm por E. Cheynet: https:// es.mathworks.com/ matlabcentral/ fileexchange/ 69030-operational-modal-analysis-with-automated-ssi-cov-algorithm?s_tid=prof_ contriblnk 8Modal Analysis OMA/EMA por R. Trumpp: https:// github.com/ KIT-FAST/ modalAnalysis-OMA-EMA 6.7 Aplicación de los métodos de identificación offline 99 frecuencia del sistema estudiado, o la técnica p-LSCF mediante los datos de la respuesta del sistema al ruido blanco como excitación. Además, una gran ventaja de este código respecto a los anteriores es que permite de manera directa integrar los cuatro ensayos realizados por identificación, reduciendo las cuatro matrices PSD de dimensión 2 estimadas en total, a una sola matriz de dimensión 5. Este algoritmo también incluye un método con clusters para identificar los modos físicos y estructurales entre los estables. 6.7.1 FDD. Valores singulares de las matrices PSD A continuación, se muestran las gráficas de los valores singulares de las matrices PSD de la respuesta estructural para cada uno de los cuatro ensayos en función de la frecuencia. 20 40 60 80 100 120 10-10 10-5 100 (a) Medida 1. 20 40 60 80 100 120 10-10 10-5 100 (b) Medida 2. 20 40 60 80 100 120 10-10 10-5 100 (c) Medida 3. 20 40 60 80 100 120 10-10 10-5 100 (d) Medida 4. Figura 6.26 Gráficas de valores singulares de las matrices PSD en función de la frecuencia. Como se puede observar, los modos están claramente diferenciados en el primer valor singular, permaneciendo el segundo prácticamente constante en toda la banda de frecuencias. Esto indica que este segundo valor singular es aproximadamente igual al ruido de fondo generado durante las mediciones, es decir, unos 3·10−5 g 2 /Hz =10 (mg) 2 /Hz. La gran diferencia en decibelios entre el valor de s2 1 en los picos y el valor de s2 2 significa una buena relación señal ruido o SNR (Signal to Noise Ratio): SNR =20 log10−1 3·10−5=70,45 dB (6.9) Lo cual supone, de acuerdo con C. Ventura et al. [25], un buen ratio para ensayos OMA, ya que sitúa el mínimo aceptable en 30-40 dB. 100 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala 6.7.2 SSI-COV. Diagramas de estabilización En cuanto a la aplicación del algoritmo de identificación modal SSI-COV, se van a representar los diagramas de estabilización generados para cada una de los cuatro ensayos. 0 20 40 60 80 100 0 2 4 6 8 10 12 14 16 18 20 (a) Medida 1. 0 20 40 60 80 100 120 0 2 4 6 8 10 12 14 16 18 20 (b) Medida 2. Figura 6.27 Diagramas de estabilización del método SSI-COV de ensayos 1 y 2 realizados para la identificación. 6.7 Aplicación de los métodos de identificación offline 101 0 20 40 60 80 100 120 0 2 4 6 8 10 12 14 16 18 20 (a) Medida 3. 0 20 40 60 80 100 120 0 2 4 6 8 10 12 14 16 18 20 (b) Medida 4. Figura 6.28 Diagramas de estabilización del método SSI-COV de ensayos 3 y 4 realizados para la identificación. Como se aprecia en las figuras anteriores, la selección de modos físicos resulta significativamente clara a partir de los polos estables alineados a lo largo de 5-10 órdenes consecutivos obtenidos en los diagramas, ya que, a pesar de que aparecen algunos fuera del conjunto de frecuencias naturales buscado, estos son aislados y se dispersan o dejan de ser estables cuando aumenta el orden. 102 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala 6.7.3 p-LSCF. Diagrama de estabilización Por último, se va a graficar el diagrama de estabilización resultante de la aplicación del método p-LSCF, en el que, gracias al empleo de clusters, se han eliminado los modos computacionales e inestables. Además, se incluye una gráfica con la curva auto-PSD experimental y la generada mediante el ajuste del modelo matemático con los polos seleccionados. 0 20 40 60 80 100 120 10-10 10-8 10-6 10-4 10-2 100 6 8 10 12 14 16 18 20 Diagrama de estabilización 0 20 40 60 80 100 120 10-10 10-8 10-6 10-4 10-2 100 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 Curva de PSD Figura 6.29 Diagrama de estabilización y curvas PSD generadas con tras la aplicación del método p-LSCF. En este caso, se identifica de manera más clara la presencia de los modos estructurales que con los diagramas obtenidos del SSI-COV. 6.7.4 Comparación de los resultados En primer lugar, se muestra una comparación a través de una gráfica de los distintos valores de frecuencia natural (Figura 6.30) y ratio de amortiguamiento modal (Figura 6.31) obtenidos con cada uno de los métodos experimentales. A continuación, otra gráfica con los modos de vibración, en la que se añade el resultado del análisis numérico (Figura 6.32). Segundo, se presentan tres tablas con las frecuencias naturales y ratios de amortiguamiento obtenidos mediante cada método (en FDD y SSI-COV se ha promediado entre todas las medidas), junto el porcentaje de variación relativa respecto al resultado numérico (Tablas 6.6, 6.7 y 6.8). Además, se han comparado los modos de vibración obtenidos experimentalmente con la estimación mediante el análisis numérico a través del MAC (Modal Assurance Criterion), índice cuya expresión se definió en la Ecuación 4.57. Por último, se incluye una visualización 3D de los modos de vibración estimados con el método FDD, en las Figuras 6.33, 6.34 y 6.35. Como se aprecia en dichos resultados, las frecuencias experimentales obtenidas son bastante similares, sin embargo, en el caso de los modos de vibración, el método SSI-COV se aleja lige- 6.7 Aplicación de los métodos de identificación offline 103 ramente del resto de resultados. En cuanto a los ratios de amortiguamiento, la variación de una técnica a otra es mayor, siendo el p-LSCF el más alejado de la media en todos los modos. Figura 6.30 Gráficas de comparación entre las frecuencias naturales obtenidas con cada método experimental. Figura 6.31 Gráficas de comparación entre los ratios de amortiguamiento modales obtenidos con cada método experimental. 104 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala Figura 6.32 Gráficas de comparación entre los modos de vibración obtenidos con cada método de análisis modal empleado. Tabla 6.6 Resultados del método FDD para fnyξ, y comparación con la estimación numérica. Modo fn,num (Hz) fn(Hz) ξn,num (%) Error fn(%) 1 7,1086 6,0182 0,343 -15,339 2 43,6208 38,1832 0,2582 -12,466 3 90,5576 99,5157 0,6731 9,892 Tabla 6.7 Resultados del método SSI-COV para fn y ξ , y comparación con la estimación numérica. Modo fn,num (Hz) fn(Hz) ξn,num (%) Error fn(%) 1 7,1086 6,0237 0,331 -15,262 2 43,6208 38,2116 0,4215 -12,401 3 90,5576 99,5696 0,6929 9,952 Tabla 6.8 Resultados del método p-LSCF para fn y ξ , y comparación con la estimación numérica. Modo fn,num (Hz) fn(Hz) ξn,num (%) Error fn(%) 1 7,1086 6,0279 0,0055 -15,203 2 43,6208 38,2664 0,0991 -12,275 3 90,5576 99,365 0,4466 9,726 6.7 Aplicación de los métodos de identificación offline 105 Tabla 6.9 MAC comparando los modos de vibración obtenidos experimentalmente con el resultado numérico. Modo FDD SSI-COV p-LSCF 1 0,9993 0,9812 0,9967 2 0,9964 0,9893 0,9960 3 0,9786 0,9121 0,9840 Figura 6.33 Visualización 3D del modo 1. Figura 6.34 Visualización 3D del modo 2. 106 Capítulo 6. Aplicación práctica: identificación modal del modelo a escala reducida de un ala Figura 6.35 Visualización 3D del modo 3. 6.7.5 Estudio para corregir el error de sincronización Por último, se va a comentar una técnica que ha empleado con el objetivo de tratar de corregir el desfase comentado anteriormente entre muestras de un acelerómetro y de otro (error de sincronización). Esta consistió en emplear la función synchronize de MATLAB, que interpola una secuencia de medidas con una determinada secuencia de tiempos asignada, a otra secuencia de tiempos distinta. Se ha empleado dicha técnica de tres maneras diferentes: 1. Interpolando las aceleraciones del primer sensor para que sus marcas de tiempo coincidan con las del segundo. 2. Interpolando las aceleraciones del segundo sensor para que sus marcas de tiempo coincidan con las del primero. 3. Interpolando ambas aceleraciones para que sus marcas de tiempo coincidan con una secuencia temporal intermedia entre las dos. Comparando las tres opciones con la suposición sin sincronización de que los pares de aceleraciones están medidos en el mismo instante, se ha obtenido que las frecuencias naturales y los ratios de amortiguamiento detectado no sufren variación, mientras que los modos de vibración difieren como se muestra en la Figura 6.36. Como se puede apreciar, la mayor variación aparece en el tercer modo (ya que la sincronización afecta más a las frecuencias más altas), encontrándose la opción de no sincronización y de sincronización de a1 y a2 en la solución intermedia. Esto se puede explicar debido a que, generalmente, cuando se interpola se está disminuyendo la amplitud de los picos en la respuesta, por lo que, por ejemplo, si se interpola a1 (sensor de la punta) cobrará más importancia en el modo detectado la coordenada de posiciones más cercanas a la raíz (posición 6). 6.8 Aplicación del método de identificación online 107 Modos de vibración Figura 6.36 Comparación de los modos de vibración obtenidos con la técnica FDD cuando se trata de interpolar las secuencias de medidas. Se muestra en azul el modo 1, en rojo el modo 2 y en verde el modo 3. 6.8 Aplicación del método de identificación online Para finalizar la parte práctica del trabajo, se ha procedido a tratar de identificar las dos primeras frecuencias naturales del ala mediante un filtro de Kalman unscented con estimación conjunta de estado y parámetros del sistema, tal y como se explicó en el Capítulo 5. Dicho filtro se ha incorporado en un código de MATLAB, haciendo uso del objeto unscentedKalmanFilter que posee dicho software. El primer paso en la implementación es la construcción de las funciones de transición y de salida discretas, como se ha descrito en la Subsección 5.4.1. A continuación, es necesario seleccionar los valores de las matrices de covarianzas del ruido de proceso, Q , y del ruido de medición, R . La segunda de ellas se puede calcular experimentalmente a través de las medidas obtenidas de los ensayos, ya que este error se debe al fallo de los sensores a la hora de medir correctamente tanto la respuesta del sistema como la carga, que en este apartado se supone nula (cuando no lo es). De esta manera, se obtiene un valor de: R=0,0691 0 0 0,0167g2(6.10) El valor de Qse ha elegido mediante prueba y error, resultando finalmente: Q=         10−8m20 0 0 0 0 0 10−8m20 0 0 0 0 0 10−8(m/s)20 0 0 0 0 0 10−8(m/s)20 0 0 0 0 0 10−10 (rad/s)20 0 0 0 0 0 10−14 (rad/s)2         (6.11) Se han realizado dos cambios respecto al algoritmo desarrollado anteriormente. En primer lugar, se han seleccionado valores diferentes de P0 para los casos en los que se cambia la primera frecuencia 114 Bibliografía [12] S.-K. Au, Operational modal analysis: Modeling, Bayesian Inference, uncertainty laws. Springer Singapore, 2017. [13] O. Altay, Vibration Mitigation Systems in Structural Engineering. CRC Press, 2021. [14] C. P. Coutinho, A. J. Baptista, and J. Dias Rodrigues, “Reduced scale models based on similitude theory: A review up to 2015,” Engineering Structures, vol. 119, pp. 81–94, 2016. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0141029616301274 [15] C. P. Coutinho, “Structural reduced scale models based on similitude theory,” Ph.D. dissertation, Universidade do Porto, Porto, Portugal, 2017. [16] W. E. Baker, P. S. Westine, and F. T. Dodge, Similarity Methods in Engineering Dynamics: Theory and Practice of Scale Modeling, rev sub ed., ser. Fundamental Studies in Engineering 12. Elsevier Science Ltd, 1991. [17] I. Horváth, J. van Beeck, and B. Merci, “Full-scale and reduced-scale tests on smoke movement in case of car park fire,” Fire safety journal, vol. 57, pp. 35–43, 2013. [18] J. Holsten, T. Ostermann, and D. Moormann, “Design and wind tunnel tests of a tiltwing uav,” CEAS aeronautical journal, vol. 2, pp. 69–79, 2011. [19] J. Feblowitz, “Confusing the wind: The burj khalifa, mother nature, and the modern skyscraper,” Inquiries Journal, vol. 2, 2010. [Online]. Available: https://api.semanticscholar.org/CorpusID: 109109729 [20] J. R. Chambers, “Modeling flight: The role of dynamically scaled free-flight models in support of nasa’s aerospace programs,” 2010. [Online]. Available: https://api.semanticscholar. org/CorpusID:106935152 [21] G. J. Simitses and J. Rezaeepazhand, “Structural similitude and scaling laws for laminated beam-plates,” Tech. Rep., 1992. [22] E. Szücs, Similitude and modelling, 1st ed. Elsevier, 1980. [23] A. Sonin, The Physical Basis of Dimensional Analysis, 2nd ed. Department of Mechanical Engineering, MIT, Cambridge, 2001. [24] C. Rainieri and G. Fabbrocino, Operational Modal Analysis of Civil Engineering Structures: An Introduction and Guide for Applications. Springer New York, 2014. [25] R. Brincker and C. Ventura, Introduction to Operational Modal Analysis, 1st ed. John Wiley & Sons, 2015. [26] M. Looney, “An introduction to mems vibration monitoring,” Analog Dialogue, vol. 48, no. 2, pp. 3–5, 2014. [Online]. Available: https://www.mouser.com.gt/pdfdocs/intro-to-mems- vibration-monitoring.pdf [27] S. Komarizadehasl, F. Lozano, J. A. Lozano-Galant, G. Ramos, and J. Turmo, “Low-cost wireless structural health monitoring of bridges,” Sensors, vol. 22, no. 15, 2022. [Online]. Available: https://www.mdpi.com/1424-8220/22/15/5725 [28] S. Komarizadehasl, B. Mobaraki, H. Ma, J.-A. Lozano-Galant, and J. Turmo, “Development of a low-cost system for the accurate measurement of structural vibrations,” Sensors, vol. 21, no. 18, 2021. [Online]. Available: https://www.mdpi.com/1424-8220/21/18/6191 Bibliografía 115 [29] L. Llamas, “El bus i2c en arduino,” https://www.luisllamas.es/arduino-i2c/, accessed August 17, 2023. [30] N. Abdukarimov, “Easy daq system development with arduino ®and tri-axis accelerometer - model adxl335,” Master’s thesis, Turin Polytechnic University, Tashkent, Uzbekistan, 2019. [31] N. Audio, “Transformación rápida de fourier fft - conceptos básicos,” https://www.nti-audio. com/es/servicio/conocimientos/transformacion-rapida-de-fourier-fft, accessed August 17, 2023. [32] G. Goitom Tekeste, “Dynamics of footbridges through operational modal analysis and vibration control using tuned mass dampers,” Master’s thesis, Instituto Superior Técnico de Lisboa, Lisboa, Portugal, 2015. [33] J. Debille, S. Pauwels, and B. Peeters, “The benefits of operational modal analysis of aircraft and spacecraft structures,” in European Test & Telemetry Conference, June 2005, pp. 1–6. [34] E. Cárdenas and L. Medina U., “Non-parametric operational modal analysis methods in frequency domain: A systematic review,” International Journal of Engineering and Technology Innovation, vol. 11, pp. 34–44, Jan. 2021. [35] F. Magalhães, “Operational modal analysis for testing and monitoring of bridges and special structures,” Ph.D. dissertation, Universidade do Porto, Porto, Portugal, 2017. [36] D. J. Ewins, Modal testing: theory, practice and application, 2nd ed. John Wiley & Sons, 2000. [37] G. James, T. Carne, J. Lauffer, and A. Nord, “Modal testing using natural excitation,” in Proceedings of the 10th International Modal Analysis Conference. Sem Society for Experimental Mechanics Inc, Feb. 1992. [38] P. Van Overschee and B. De Moor, Subspace identification for linear systems: Theory—Implementation—Applications. Springer Science & Business Media, 2012. [39] R. Brincker, L. Zhang, and P. Andersen, “Modal identification from ambient responses using frequency domain decomposition,” in IMAC 18: Proceedings of the International Modal Analysis Conference (IMAC), San Antonio, Texas, USA, Feb. 2000, pp. 625–630. [40] P. Welch, “The use of fast fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms,” IEEE Transactions on Audio and Electroacoustics, vol. 15, no. 2, pp. 70–73, 1967. [41] R. Brincker, L. Zhang, and P. Andersen, “Modal identification of output-only systems using frequency domain decomposition,” Smart Materials and Structures, vol. 10, no. 3, p. 441, June 2001. [Online]. Available: https://dx.doi.org/10.1088/0964-1726/10/3/303 [42] B. Peeters and H. Van der Auweraer, “Polymax: a revolution in operational modal analysis,” in Proceedings of the 1st International Operational Modal Analysis Conference, vol. 820, Copenhagen, Denmark, Apr. 2005, pp. 1–12. [43] E. Reynders, G. Degrande, and G. De Roeck, “System identification and modal analysis in structural mechanics,” Ph.D. dissertation, Katholieke Universiteit, Leuven, Belgium, 2009. [44] R. F. Trumpp, “Implementierung des poly-reference least square complex frequency (p-lscf) algorithmus zur operational modal analysis,” Bachelor’s Thesis, Karlsruher Institut für Technologie (KIT), 2017. 116 Bibliografía [45] S. Schleiter and O. Altay, “Identification of abrupt stiffness changes of structures with tuned mass dampers under sudden events,” Structural Control and Health Monitoring, vol. 27, no. 6, p. e2530, 2020. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/stc.2530 [46] A. W. Smyth, S. F. Masri, A. G. Chassiakos, and T. K. Caughey, “On-line parametric identification of mdof nonlinear hysteretic systems,” Journal of Engineering Mechanics, vol. 125, no. 2, pp. 133–142, 1999. [Online]. Available: https://ascelibrary.org/doi/abs/10.1061/ %28ASCE%290733-9399%281999%29125%3A2%28133%29 [47] I. Yoshida and T. Sato, “Health monitoring algorithm by the monte carlo filter based on non-gaussian noise,” Journal of Natural Disaster Science, vol. 24, no. 2, pp. 101–107, 2002. [48] S. J. Julier, J. K. Uhlmann, and H. F. Durrant-Whyte, “A new approach for filtering nonlinear systems,” in Proceedings of 1995 American Control Conference-ACC’95, vol. 3, Seattle, WA, USA, 1995, pp. 1628–1632. [49] E. Wan and R. Van Der Merwe, “The unscented kalman filter for nonlinear estimation,” in Proceedings of the IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium (Cat. No.00EX373), Lake Louise, AB, Canada, 2000, pp. 153–158. [50] T. Irvine, “Bending frequencies of beams, rods and pipes,” 2022. [51] P. Jackson, Jane’s All the World’s Aircraft, ser. IHS Jane’s All the World’s Aircraft. Jane’s Information Group, 2004. [Online]. Available: https://books.google.es/books?id= eS2yMQEACAAJ [52] H. A. Abdulameer and H. R. Wasmi, “Vibration control analysis of aircraft wing by using smart material,” Innovative Systems Design and Engineering, vol. 6, no. 8, pp. 7–42, 2015. [53] S. Fichera, S. Jiffri, and J. Mottershead, “Design and wind tunnel test of a modular aeroelastic flexible wing (modflex),” in Proceedings of ISMA2016 International Conference on Noise and Vibration Engineering and USD2016 International Conference on Uncertainty in Structural Dynamics, 2016, pp. 445–456. [54] E. Cheynet, J. B. Jakobsen, and J. Snæbjörnsson, “Damping estimation of large wind-sensitive structures,” Procedia Engineering, vol. 199, pp. 2047–2053, 2017. [Online]. Available: https://doi.org/10.1016%2Fj.proeng.2017.09.471 [55] E. Cheynet, “Operational modal analysis with automated ssi-cov algorithm,” 2020. [Online]. Available: https://zenodo.org/record/3774061 [56] F. Magalhães, A. Cunha, and E. Caetano, “Online automatic identification of the modal parameters of a long span arch bridge,” Mechanical Systems and Signal Processing, vol. 23, no. 2, pp. 316–329, 2009. [Online]. Available: https://www.sciencedirect.com/science/article/ pii/S0888327008001301