Actualización Bayesiana de Modelos de Elementos Finitos de Estructuras de Ingeniería Civil

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería Civil Actualización Bayesiana de Modelos de Elementos Finitos de Estructuras de Ingeniería Civil Autor: Arturo Novales Salamero Tutor: Javier Fernando Jimenez Alonso Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 2 Trabajo de Fin de Grado Grado en Ingeniería Civil Actualización Bayesiana de Modelos de Elementos Finitos de Estructuras de Ingeniería Civil Autor: Arturo Novales Salamero Tutor: Javier Fernando Jimenez Alonso Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 4 Trabajo de Fin de Grado: Actualización Bayesiana de Modelos de Elementos Finitos de Estructuras de Ingeniería Civil Autor: Arturo Novales Salamero Tutor: Javier Fernando Jimenez Alonso El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2024 El Secretario del Tribunal 6 A mi familia por confiar en mi a cada paso que he dado A mi pareja por su apoyo incondicional A mis amigos por remar conmigo en esta ardua travesía A mis maestros por su ayuda y comprensión 8 Agradecimientos Tras unos largos años de trabajo y dedicación, toca poner el broche al camino que tanto esfuerzo ha requerido. Nunca antes sin agradecer a las personas que han tenido un papel crucial en toda esta etapa. Gracias a mis amigos, a Iciar y a mi familia por horas y horas de comprensión, empatía, desconexión y consejos que han resultado necesarios. Han sido tantas las horas dedicadas e innumerables momentos de sacrificio que me han condicionado que, sin el apoyo de estas personas, nada hubiera sido lo mismo. Este satisfactorio resultado es tanto mío como de ellos. Así como también agradecer a maestros y compañeros por la ayuda siempre brindada y a todo aquel que me ha ayudado en el camino y ha empujado mi afán de aprender. Arturo Novales Salamero Estudiante del Grado en Ingeniería Civil Sevilla, 2024 16 Tabla de ilustraciones Ilustración 1. Esquema conceptual Digital Twin enunciado en 2002 por Dr. Michael Grieves (Fuente: (Grieves, 2016)) 21 Ilustración 2. Extensómetro de terraplén usado para medir asientos en el terreno (Fuente: SISGEO) 23 Ilustración 3. Sensor de fibra óptica (OFS) - Foto de la fibra e ilustración 3D de su sección (Fuente: (Bado et al., 2022)) 23 Ilustración 4. Interfaz de la plataforma Space Out con información, a tiempo real, de la densidad de población en cada zona urbana (Fuente: https://www.developer.tech.gov.sg/products/categories/digital-solutions-to- address-covid-19/space-out/overview.html) 24 Ilustración 5. Descripción del flujo de trabajo del modelado DTE Hydrology (Fuente: (Brocca et al., 2024)) 26 Ilustración 6. Distribución de publicaciones científicas sobre los gemelos digitales en función del área de conocimiento (Fuente: (Bado et al., 2022) ) 27 Ilustración 7. Esquema conceptual de funcionamiento de una Smart Factory (Fuente: (Friederich et al., 2022) ) 28 Ilustración 8. Flujo de trabajo de un proceso de actualización de un modelo de elementos finitos 29 Ilustración 9. Esquema de funcionamiento - Actualización en el tiempo 30 Ilustración 10. Diagrama de flujo de la actualización de un modelo de elementos finitos bajo el método bayesiano 36 Ilustración 11. Situación gran escala - Municipios - Google Maps 46 Ilustración 12. Situación pequeña escala - Google Maps 46 Ilustración 13. Pasarela - Google Maps 47 Ilustración 14. Tablero compuesto real 48 Ilustración 15. Geometría de un nervio homogeneizada a hormigón 49 Ilustración 16. Modelado - Vigas cuadrangulares longitudinales 52 Ilustración 17. Modelado - Vigas cuadrangulares transversales y verticales 53 Ilustración 18. Modelado - Vigas transversales en dirección diagonal con sección en L – Visión global 53 Ilustración 19. Modelado - Vigas transversales en dirección diagonal con sección en L – Visión local 53 Ilustración 20. Modelado – Tablero – Visión global 53 Ilustración 21. Modelado – Tablero – Visión local 53 Ilustración 22. Modelo de elementos finitos completo 54 Ilustración 23. Modelado – Condiciones de contorno en direcciones transversal y vertical 55 Ilustración 24. Modelado - Condiciones de contorno en dirección longitudinal - Elemento muelle unido a los extremos del tablero y a un empotramiento 55 Ilustración 25. Análisis modal - Primer modo de vibración vertical – Correspondiente a una frecuencia de 3,78904 Hz 56 Ilustración 26. Análisis modal - Segundo modo de vibración vertical – Correspondiente a una frecuencia de 12,6789 Hz 56 Ilustración 27. Análisis modal - Tercer modo de vibración vertical – Correspondiente a una frecuencia de 21,9883 Hz 56 Ilustración 28. Análisis modal - Cuarto modo de vibración vertical – Correspondiente a una frecuencia de 53,3202 Hz 56 Ilustración 29. Análisis modal - Elección de nodos representativos para la determinación de la forma de los modos de vibración 57 Ilustración 30. Visualización energía de deformación - Primer modo de vibración vertical 60 Ilustración 31. Visualización energía de deformación - Segundo modo de vibración vertical 60 Ilustración 32. Visualización energía de deformación - Tercer modo de vibración vertical 60 Ilustración 33. Visualización energía de deformación - Cuarto modo de vibración vertical 60 Ilustración 34. Desplazamientos verticales nodales para los modos verticales de vibración 62 Ilustración 35. Resultado del estudio del coeficiente de Spearman 63 Ilustración 36. Valores numéricos y reales de los parámetros a actualizar 64 Ilustración 37. Evolución de las soluciones mediante el algoritmo genético en el método determinista 66 Ilustración 38. Coeficientes resultantes – Método determinista 66 Ilustración 39. Comparativa coeficientes teóricos y numéricos – Método determinista 67 Ilustración 40. Comparativa frecuencias numéricas y experimentales – Método determinista 67 Ilustración 41. Resultado de los indicadores MAC – Método determinista 67 Ilustración 42. Resultados – Variables actualizadas – Método determinista 68 Ilustración 43. Modelo de edificio a cortante bajo estudio 69 Ilustración 44. Matriz de masa - Edificio a cortante 70 Ilustración 45. Matriz de rigidez - Edificio a cortante 70 Ilustración 46. Edificio a cortante - Valores experimentales de las cuatro primeras frecuencias 72 Ilustración 47. Edificio a cortante - Histogramas de las cuatro primeras frecuencias 72 Ilustración 48. Resultados de la optimización mediante el método de enjambre de partículas - Aceptación obtenida para el 𝜎𝑃 optimizado 74 Ilustración 49. Edificio a cortante - MCMC - 1000 iter; 0 burn in 75 Ilustración 50. Edificio a cortante - MCMC - 1000 iter; 200 burn in 75 Ilustración 51. Edificio a cortante - MCMC - COV y aceptación 76 Ilustración 52. Edificio a cortante - MCMC – Histogramas 76 Ilustración 53. Edificio a cortante - MCMC - Medias de las iteraciones de los anchos de los pilares para cada planta 76 Ilustración 54. Edificio a cortante - MCMC - Tiempo invertido 77 Ilustración 55. Edificio a cortante - TMCMC - 1000 iter; 0 burn in 77 Ilustración 56. Edificio a cortante - TMCMC - 1000 iter; 100 burn in 78 Ilustración 57. Edificio a cortante - TMCMC - COV 78 18 Ilustración 58. Edificio a cortante - TMCMC - Histogramas 78 Ilustración 59. Edificio a cortante - TMCMC - Medias de las iteraciones de los anchos de los pilares para cada planta 79 Ilustración 60. Edificio a cortante - TMCMC - Tiempo invertido 79 Ilustración 61. Edificio a cortante – SMC – 1000 iter 79 Ilustración 62. Edificio a cortante – SMC - COV 80 Ilustración 63. Edificio a cortante – SMC - Histogramas 80 Ilustración 64. Edificio a cortante – SMC – Medias de las iteraciones de los anchos de los pilares para cada planta 80 Ilustración 65. Edificio a cortante – SMC – Tiempo invertido 81 Ilustración 66. Edificio a cortante - Comparativa métodos 81 Ilustración 67. Pasarela - Histogramas - Mediciones experimentales de las propiedades dinámicas 82 Ilustración 68. Pasarela - 1000 iteraciones; 0 burn-ins 83 Ilustración 69. Pasarela – 1000 iteraciones; 100 burn-ins 84 Ilustración 70. Pasarela – COV 84 Ilustración 71. Pasarela – Medias de los parámetros estimadas 84 Ilustración 72. Pasarela – Histogramas de los valores de los parámetros actualizados 85 Tabla de ecuaciones Ecuación 1. Resolución del método de máxima verosimilitud 35 Ecuación 2. Expresión del teorema de Bayes 37 Ecuación 3. Simplificación del teorema de Bayes para su aplicación práctica 37 Ecuación 4. Teorema de Bayes simplificado 37 Ecuación 5. Expresión de la función de probabilidad para Dk independientes e idénticamente distribuidas 38 Ecuación 6. Expresión de la función de probabilidad para una distribución normal 38 Ecuación 7. Función de probabilidad logarítmica 38 Ecuación 8. Criterio de aceptación o rechazo del algoritmo M-H 41 Ecuación 9. Criterio de aceptación o rechazo del algoritmo M-H simplificado 41 Ecuación 10. Distribuciones transicionales del método TMCMC 42 Ecuación 11. Coeficiente de variación del método TMCMC 43 Ecuación 12. Expresión de 𝛽𝑗 del método TMCMC 43 Ecuación 13. Función de peso del método TMCMC 43 Ecuación 14. Tamaño efectivo de muestra del método SMC 44 Ecuación 15. Coeficiente de homogeneización 48 Ecuación 16. Obtención del ancho equivalente en la homogeneización 48 Ecuación 17. Igualdad necesaria para la veracidad del modelado 49 Ecuación 18. Sistema de ecuaciones definitivo a cumplir para la veracidad del modelado 50 Ecuación 19. Primera condición definida para el modelado 50 Ecuación 20. Condición consecuente impuesta por la definición anterior 50 Ecuación 21. Obtención del espesor equivalente del tablero modelado 50 Ecuación 22. Obtención de la última propiedad mecánica del material para su modelado 50 Ecuación 23. Cálculo de volúmenes y masas reales de cada material que compone el tablero compuesto 51 Ecuación 24. Cálculo de la masa del tablero modelado 51 Ecuación 25. Obtención de la densidad equivalente modelada del tablero para el ajuste de masas 51 Ecuación 26. Condiciones de igualdad a cumplir referentes a los coeficientes de Poisson y módulos de deformación modelados y obtención de los primeros 52 Ecuación 27. Obtención del módulo de torsión del material modelado 52 Ecuación 28. Expresión del indicador estadístico MACG 57 Ecuación 29. Expresión de las frecuencias naturales para una viga aislada biapoyada en 2D 58 20 Ecuación 30. Igualdad que seguirá el resultado obtenido por el método determinista 64 Ecuación 31. Residuo de las frecuencias í-esimas numéricas y experimentales ponderadas al 90% 65 Ecuación 32. Residuo de los modos de vibración i-ésimos numéricos y experimentales ponderados al 10% 65 Ecuación 33. Función objetivo que minimizar 66 Ecuación 34. Rigidez de la planta i-ésima 70 Ecuación 35. Momento de inercia de un pilar de sección cuadrada 70 Ecuación 36. Expresión de las frecuencias naturales de un sistema 70 Ecuación 37. Frecuencias naturales de un sistema en hercios 71 Ecuación 38. Ecuación para determinar los autovectores del sistema - Modos de vibración asociados a la frecuencia wi 71 Ecuación 39. Ecuación para la obtención de las medidas experimentales para la propiedad dinámica i – Edificio a cortante 71 Ecuación 40. Obtención del valor superior o inferior de los parámetros a actualizar 73 1 INTRODUCCIÓN unque la terminología ha cambiado con el tiempo, el concepto básico del Digital Twin (DT), o Gemelos Digitales, se ha mantenido constante desde su primera aplicación por parte del Dr. Michael Grieves, que expuso este concepto de manera oficial en la Universidad de Michigan en 2002. Este consiste en el desarrollo de una construcción digital basada en la información que conforma un sistema físico tangible. De esta manera, la información digital será “gemela” de la información proveniente del sistema físico y estará enlazada con él a lo largo de toda la vida útil de éste. Así fue como, en esta ocasión, el doctor expuso una presentación sobre la formación de un centro de gestión del ciclo de vida del producto (PLM por sus siglas en inglés) cuyo esquema conceptual se englobó en la imagen a continuación. En esta se tratan sus agentes principales, como son el espacio físico real, espacio virtual digital y el enlace entre ambos: Ilustración 1. Esquema conceptual Digital Twin enunciado en 2002 por Dr. Michael Grieves (Fuente: (Grieves, 2016)) En esta ilustración se representa la posibilidad de tener espacios virtuales múltiples (VS) para un mismo espacio físico (producto), consiguiendo una aplicación práctica tan funcional como la posibilidad de realizar cualquier tipo de simulación sobre el modelo virtual considerando su verosimilitud con el real. Hecho que implica la A 22 posibilidad de realización de acciones, de ninguna otra manera no-intrusiva posibles, tales como ensayos destructivos múltiples sin coste directo alguno. Con esta rama de estudio se pueden analizar objetivos del producto (predicciones deseadas) y problemas (predicciones no deseadas e incidencias ni si quiera concebidas (no predichas)), así como, mediante la variación de los parámetros de simulación, a lo largo de los rangos que pueden tomar, es posible una investigación de comportamientos no lineales en sistemas complejos. Esta amplitud de aptitudes no se traduce en poder modelar y simular todas las posibilidades disponibles en la imaginación, pues el análisis de todas las permutaciones y combinaciones posibles de un sistema complejo puede no ser viable desde un punto de vista de economía de tiempo. Sin embargo, los avances exponenciales en la capacidad computacional van acercándonos a análisis cada vez más exhaustivos y que será materia de mención en este documento. Como se puede observar, la premisa consiste en que cada sistema lo conforman dos subsistemas, el físico y el virtual, que contiene toda la información del primero. El título de PLM que trató el doctor hacía hincapié en la idea de que esos dos sistemas estuvieran enlazados a lo largo de toda la vida útil del sistema, diferenciando ésta en 4 fases: creación, fabricación, operación y desecho. Finalmente, recogiendo estas ideas, se puede desarrollar una definición aproximada del gemelo digital como una serie de construcciones informáticas virtuales que consiguen describir un producto físico fabricado, potencial o actual, en su totalidad, desde un nivel micro-atómico a un nivel macro-geométrico. Siendo, el ideal del gemelo digital, aquel modelo en el que existe la posibilidad de consulta de cualquier información que pueda ser inspeccionada en el producto físico (Grieves, 2016). 1.1 Gemelos digitales en la ingeniería civil 1.1.1 Ingeniería estructural Vivimos en un entorno con una demanda constantemente creciente de infraestructuras de transporte, pues se relacionan con el funcionamiento y crecimiento de la economía del lugar. Por otro lado, además, también vivimos en un entorno en el que sus estructuras llevan puestas en funcionamiento varios años y, algunas, van alcanzando el final de su vida útil. De esta manera, mientras la sustitución de estas estructuras no siempre es viable debido a factores económicos, logísticos y medioambientales, el único acercamiento restante es el de mantenimiento: rehabilitación, restauración y regeneración. Así, nace el concepto de durabilidad estructural, denominado como Monitorización de la Salud Estructural (conocido universalmente en la lectura científica como Structural Health Monitoring), que puede entenderse como la medición continua y análisis de parámetros estructurales y medioambientales clave bajo condiciones de operación, con el objetivo de alertar de un estado anormal o de un accidente en su etapa inicial (Bado et al., 2022). Como puede empezar a visualizarse entre estas ideas, el concepto se entrelaza con el relativo al BIM (siglas en inglés de Building Information Modeling) donde, de acuerdo a la ISO 19650:2019, es el uso de una representación digital de una entidad para facilitar el diseño y construcción de esta, así como la generación de unas bases fiables para la toma de decisiones. El potencial del BIM reside, a su vez, en la posibilidad de gestión, en un entorno digital, de toda la información y documentación de una infraestructura a lo largo de toda su vida útil. Así, para conseguir e, incluso, exceder los objetivos de la construcción, nace un modelo digital inteligente y paralelo al real, el gemelo digital. No obstante, este último concepto de gemelo digital va un paso más adelante en el momento en el que se integra un sistema de reconocimiento y lectura de datos mediante un sensor que actúa como un sistema nervioso de la estructura. Consiguiendo, de esta manera, el intercambio de información (como pueden ser valores de deformación, temperatura, etc..) constante y automático, es decir, obteniendo datos a tiempo real y enlazándolos entre la construcción tangible física y la estructura digital desarrollada a partir de ella. Denominando, a estos sistemas, Smart Infrastructure o Smart Buildings (Bado et al., 2022). Estos sensores pueden ser cualquier dispositivo con la capacidad de medir una respuesta estructural e intercambiar los datos leídos, algunos más tradicionales como los inclinómetros, acelerómetros, extensómetros o células de carga y otros más novedosos como Sensores de Fibra Óptica (OFS), Sistemas de Posicionamiento Global (GPS), radares, Sistemas Microelectromecánicos (MEMS) y técnicas de procesamiento de imágenes como la correlación de imágenes digitales. Ilustración 2. Extensómetro de terraplén usado para medir asientos en el terreno (Fuente: SISGEO) Destacando, entre los más avanzados tecnológicamente, los OFS, pues una de sus destacables capacidades es la de medir deformaciones derivadas de variaciones mecánicas y térmicas a lo largo de la fibra mediante la dispersión y retrodispersión de la luz al ser, estos dispositivos, usados para confinar y guiar ésta. Ilustración 3. Sensor de fibra óptica (OFS) - Foto de la fibra e ilustración 3D de su sección (Fuente: (Bado et al., 2022)) 1.1.2 Ingeniería de planificación urbana, hidráulica, hidrológica, geotécnica y otras Se ha dedicado un apartado a la ingeniería estructural por ser el foco de este documento. Sin embargo, la aplicación de técnicas de calibración de modelos es una práctica utilizada en todos los campos que abarca la ingeniería civil. De esta manera, en la actualidad, las principales aplicaciones de los DT, en los distintos sectores, abarcan el diseño/planteamiento, optimización, mantenimiento, seguridad, toma de decisiones, acceso remoto y entrenamiento entre otros. Algunos ejemplos generales son: - Elementos de generación de energía: gestión de mantenimiento de aerogeneradores gracias a la monitorización de las distintas piezas que lo conforman. - Estructuras: optimización de grandes conjuntos estructurales. 24 - Operaciones de fabricación: guía de productos desde su diseño hasta el cumplimiento de su vida útil. - Industria automovilística: mejora del rendimiento y aumento de la eficiencia sobre su producción. - Planificación urbana: elaboración de espacios 3D y 4D en tiempo real. Además, la aplicación de esta tecnología, con el paso del tiempo, va aumentando cada vez más su alcance. Entrando en ejemplos particulares: 1. La empresa TESLA, como ejemplo, apunta a desarrollar un gemelo digital para cada coche fabricado y, mediante la sincronización de transmisión de datos entre el coche y la fábrica, realizar una gestión de mantenimiento detallada al nivel de usuario particular, optimizando recursos. 2. Por otro lado, compañías energéticas como General Electric o Chevron, usan esta tecnología para monitorear la operatividad de los aerogeneradores. 3. Anteriormente se ha hablado de lo conocido como Smart Buildings; relativo a ese concepto se encuentra, entre las aplicaciones más complejas de esta tecnología, el concepto de Smart Cities. El DT de un modelo digital de una ciudad permite, a planificadores urbanos, ingenieros y arquitectos, simular soluciones a problemas urbanos sin necesidad de aplicarlos en la realidad, una práctica denominada como Experimentación Virtual. Este campo de trabajo suele analizar una gran variedad de aspectos, desde la mejora de la movilidad hasta uno de los ejemplos más famosos, el conocido como: Virtual Singapore. Este proyecto implicó la creación, por parte del gobierno de Singapur junto con Dassault Systèmes, de un modelo 3D de la ciudad apoyado por datos estáticos y dinámicos recogidos a partir de una extensa red de sensores que actúan como recolectores de información. Con esta acción querían cumplir el objetivo de llevar, a la máxima expresión, el mundo de la planificación urbana, y optimizar la gestión urbana aplicada, desde una evacuación de emergencia, a una forma cómoda de vivir en la ciudad. Uno de sus ejemplos de uso fue en Abril de 2020 cuando, en un intento de gestión y contención de la propagación del virus COVID-19, el gobierno lanzó la plataforma Space Out, que mostraba el gemelo digital de la ciudad representando los niveles de población en ella, con el fin de que los propios ciudadanos pudieran tomar decisiones informadas con conocimiento a tiempo real en base a las recomendaciones de distanciamiento social (Bado et al., 2022). Ilustración 4. Interfaz de la plataforma Space Out con información, a tiempo real, de la densidad de población en cada zona urbana (Fuente: https://www.developer.tech.gov.sg/products/categories/digital-solutions-to- address-covid-19/space-out/overview.html) Como se puede observar, la monitorización constante tiene sus ventajas a la hora de toma de decisiones por parte del usuario. Sin embargo, por otro lado, también se utilizó este gemelo digital para aplicaciones como: - Calcular la cantidad de energía que podría generar la instalación de paneles solares en determinadas azoteas - Desarrollo y optimización de modelos de evacuación para gestión de desastres - Visualización del paisaje existente ante presentes y futuros proyectos de construcción con el fin de armonizar y promover una planificación urbana holística - Identificación de rutas sin obstáculos para aquellos usuarios que tengan una movilidad más reducida Una vez vistas todas estas aplicaciones, podemos sacar la conclusión de que los DT resultan una herramienta muy útil para cualquier ámbito y, por recorrer más usos relativos a la ingeniería civil, el campo de la hidráulica e hidrología no iba a ser menos en la implementación de este concepto. 4. Se ha usado este tecnología para la diagnosis de fuga en una red de distribución de agua. El gemelo digital permitió la operación remota midiendo, online, la presión en cabeza y el caudal de flujo. De esta manera, mediante un modelo bien calibrado del software EPANET, bien conocido en la escuela, llevaron a cabo simulaciones para estimar la presión en cabeza sin sensores y, mediante un algoritmo genético, fueron capaces de detectar y estimar el tamaño de esta discontinuidad (Gómez-Coronel et al., 2023). 5. Por otro lado, en el mundo de la hidrología, se ha avanzado la herramienta hasta un nuevo concepto: modelos Digital Twin Earth (DTE) que consiste en una innovadora solución que ofrece réplicas digitales para monitorear y simular procesos de la Tierra con una resolución espacio-temporal sin precedentes. La nube de datos del 4D DTE Hydrology fusiona los datos de alta resolución provenientes de la observación terrestre vía tecnología satelital y modelos avanzados de precipitación, humedad del suelo, evaporación y descarga fluvial. 32 2.2.1 Métodos directos Son métodos no iterativos que se usan para actualizar la estructura del modelo de elementos finitos mediante la modificación de la matriz de rigidez estructural y la matriz de masa. Reproducen los datos experimentales exactos y, al no requerir iteraciones, se convierten en algoritmos eficientes computacionalmente (Marwala et al., 2017). Son métodos que se siguen usando actualmente y que, mediante la ayuda de sensores, pueden superar algunas de sus desventajas como la necesidad de grandes cantidades de datos para reproducir con exactitud el comportamiento real de la estructura. Sin embargo, su uso se ha visto reducido a favor de los métodos indirectos, en los que nos centraremos en este trabajo. 2.2.2 Métodos indirectos Estos métodos iterativos han sido desarrollados, evitando las limitaciones del método anterior, para actualizar modelos de elementos finitos complejos. Estos métodos usan ecuaciones no lineales para problemas de optimización no convexa propios de la actualización de sistemas complejos. Son métodos en los que una serie de parámetros son ajustados iterativamente para minimizar una función objetivo donde, en la mayoría de los casos de actualización de modelos, contienen datos modales y/o de funciones respuesta (Marwala et al., 2017). Los métodos directos no son computacionalmente costosos, pero reproducen propiedades modales que no son físicamente reales. Mientras que las técnicas iterativas conciben cambios en los parámetros físicos para actualizar los modelos y, por lo tanto, generar modelos que son físicamente realistas. Estos métodos iterativos son los más usados para el proceso de calibración en las estructuras de ingeniería civil y pueden ser de dos tipos: deterministas y estocásticos. A pesar de que el problema sea fundamentalmente el mismo, aminorar la diferencia entre el comportamiento estructural predicho por el modelo numérico y el comportamiento real, tendrá una naturaleza u otra en función del método que se use. A continuación, se introducirán estas dos familias de métodos puesto que, al ser materia principal del objetivo de este documento, le dedicaré una explicación más desarrollada en los siguientes puntos de este documento. Seguidamente se mencionarán otros métodos dada su relevancia en la actualidad, pero no se entrará en detalle en ellos. 2.2.2.1 Métodos deterministas. Métodos de máxima verosimilitud (MLE) Convirtiendo el problema de actualización en uno de optimización, la actualización se lleva a cabo a través del uso del método de máxima verosimilitud (MLE), en el que se define una función objetivo considerando los valores residuales de distintos tipos de magnitudes. Sin embargo, las más comunes son las propiedades dinámicas de la estructura: frecuencias naturales y modos de vibración. Estos datos son los mejores indicadores del comportamiento real de la estructura pues un cambio en la estructura implica modificaciones en la rigidez estructural, modificando las propiedades dinámicas de ella (Ereiz et al., 2022). 2.2.2.2 Métodos estocásticos Entendemos como estocástico el sistema cuyo comportamiento intrínseco no es determinista, sino que se puede analizar en términos de probabilidad. Por lo tanto, en contraposición al problema de optimización anterior, usar este método para actualizar el modelo de elementos finitos significa realizar un problema estadístico considerando los errores de medida y/o las incertidumbres del modelo. 2.2.3 Otros métodos 2.2.3.1 Métodos de Inteligencia Artificial Las técnicas basadas en la inteligencia artificial son herramientas computacionales que son inspiradas por la forma en la que los sistemas biológicos y naturales funcionan. Desarrollaré aquí dos, de los tantos métodos que componen este grupo, por el papel que juegan en la aplicación práctica de este escrito. 1. Algoritmos genéticos Simula la evolución natural, donde la supervivencia del más fuerte se aplica a una población de individuos. La teoría de la selección natural de Darwin enuncia que los rasgos heredables de un organismo, capaces de aumentar su capacidad de supervivencia y de reproducción, se acaban volviendo más comunes en una población a lo largo del tiempo. Así, el método se basa en esta teoría de la evolución, donde los individuos pasarán cierta función de aptitud y los mejores serán elegidos como una solución eficaz del problema a resolver. Posteriormente se cruzarán unos datos con otros llevando a cabo una recombinación de éstos y lo que se conoce como “mutación” (para evitar bloqueos en la cadena iterativa), volviendo a evaluar la aptitud una vez realizada. El proceso se repetirá una serie de veces hasta encontrar una solución lo suficientemente satisfactoria (o por criterio de parada) logrando optimizar la función objetivo. 2. Optimización mediante enjambre de partículas El método de enjambre de partículas es un algoritmo de optimización que está inspirado en el comportamiento de animales en grupo introducido por Russel Eberhart y James Kennedy en 1995 que destaca por el trasfondo de la cooperación entre las partículas para la ‘’comunicación’’ entre ellas y el aprendizaje de los datos recabados por cada una. Las partículas interactúan y aprenden entre ellas. Además de las variables que se consideran, como son la posición y la velocidad de la partícula, cada una tiene almacenada la información de su mejor experiencia (local) siendo, la mejor experiencia de entre todas las locales de las partículas, la que se denota como mejor experiencia global. De esta manera, para un nuevo movimiento, la partícula usa la experiencia previa de ella misma y la mejor experiencia global. Así, todas las partículas operan para encontrar la mejor solución para el problema de optimización. 2.2.3.2 Métodos basados en la sensibilidad Son métodos que trabajan bajo la premisa de que la respuesta real medida es una perturbación de la respuesta obtenida mediante el modelo de elementos finitos y, por lo tanto, que el primero pueda ser aproximadamente igual a los resultados predichos por el modelo numérico (Marwala et al., 2017). Suele usar propiedades modales como su base para la actualización y calcula, mediante series de Taylor, la diferencia entre los valores predichos. Se ejecutará iteración tras iteración hasta que éstos convergen cumpliendo una tolerancia adecuada. Dada su naturaleza iterativa, podríamos incluir estos métodos en la clasificación de métodos indirectos. 34 3 MÁXIMA VEROSIMILITUD Y MÉTODO BAYESIANO A incertidumbre en ingeniería estructural es un asunto importante para el desarrollo de modelos necesarios para el diseño y control de estructuras complejas. Como hemos comentado en el anterior capítulo, en el campo del modelado de estructuras existen diversas fuentes de incertidumbre que ocasionarán un límite en cuanto a exactitud de datos de salida respecta. Estas incertidumbres han de ser cuantificadas y, para ello, podemos clasificar en dos grupos los tipos de problemas de cuantificación de incertidumbre: 1. Problemas directos Aquellos donde, para un valor de entrada, el modelo y el valor de salida, dado que el valor de entrada es incierto, se estudia cómo se desarrolla esa incertidumbre hasta llegar al valor de salida Consisten en desarrollar un modelo estocástico que prediga, en términos de probabilidad, los posibles valores de los datos de salida y sus estadísticas dado unos parámetros, inciertos, de entrada. Se pueden resolver mediante procesos analíticos o, generalmente, mediante simulaciones Monte Carlo. Las estadísticas de los valores de salida se obtienen mediante, primeramente, generando n muestras de los parámetros a partir de su función de densidad de probabilidad conocida, para luego, para cada valor, evaluar el modelo y obtener su respuesta. Repitiendo este proceso para todas las muestras definidas, acabaremos obteniendo una distribución de los valores de salida. Por lo tanto, este problema acaba siendo de una implementación sencilla una vez definimos la distribución de esos parámetros a actualizar. Sin embargo, es ahí donde se encuentra el verdadero reto, el encontrar la distribución más apropiada para predecir las medidas experimentales (Lye et al., 2021). Serán, este tipo de problemas, el objetivo final de análisis de este trabajo. 2. Problemas inversos Consiste en identificar un modelo adecuado, con sus parámetros variables, para que se relacionen adecuadamente los datos de entrada y salida siendo ambos inciertos. Para cada grupo hay una vertiente principal que los estudia, la bayesiana para el problema directo y la frecuentista para el inverso. El estimador de máxima verosimilitud es uno de los estimadores más usados que consigue encontrar los valores de los parámetros que maximizan la probabilidad de observar los valores experimentales, mientras que el acercamiento Bayesiano trabaja el problema como una inferencia Bayesiana, que es un tipo de inferencia estadística (parte de la estadística que se basa en la inducción para determinar propiedades de una población) en la que las evidencias (observaciones) se emplean para actualizar la probabilidad de que una hipótesis pueda ser cierta. L 2.1 Método de máxima verosimilitud Es el método más popular para aproximar funciones de distribución y, como comentamos en el anterior apartado, es, esencialmente, un problema de optimización, donde el objetivo es reducir la diferencia entre el valor medido y el numérico. De esta manera, la función de probabilidad se define como el error entre ambos tipos de datos, mientras que la solución aproximada se obtiene maximizando su expresión (Lye et al., 2021). 𝜃∗=argmax (log (𝑃(𝐷|𝜃)) Ecuación 1. Resolución del método de máxima verosimilitud Donde D es el vector que contiene las medidas experimentales y 𝜃 los parámetros a actualizar. 2.3. Inferencia bayesiana Las estadísticas Bayesianas se usan para obtener la probabilidad de los parámetros desconocidos de un modelo de elementos finitos dados unos datos experimentales, es decir, lo que llamamos función de densidad de probabilidad posterior. El acercamiento, al problema, bayesiano relaciona esto con la función de probabilidad (que mide la diferencia entre la predicción mediante el modelo numérico y los datos medidos), la probabilidad prior y la evidencia. Incidimos en el método bayesiano pues su principio básico es que esos parámetros inciertos, así como los errores de modelado y medida, se modelan como variables de diseño. La clave de usar esta modalidad reside en su habilidad para combinar información a priori de los datos de interés con los datos de observación, de manera que podamos formar una caracterización estocástica de los datos de estudio. Así, cuando accedemos a una serie de observaciones independientes e idénticamente distribuidas, la información a priori es actualizada mediante el uso del teorema de Bayes, obteniendo la distribución posterior del parámetro de interés. 36 Ilustración 10. Diagrama de flujo de la actualización de un modelo de elementos finitos bajo el método bayesiano Es un modelo estadístico que se basa en la teoría de la probabilidad, donde la incertidumbre se modela usando una función de densidad de probabilidad (PDF). La distribución de esta se puede definir mediante la de Bernoulli, binomial, de Poisson o exponencial, siendo, las más usadas en ingeniería civil y actualización de MEF, las distribuciones uniforme y normal. A cada serie de parámetros se le asigna una función de densidad de probabilidad que incorpora información sobre el valor de los parámetros estructurales. La incertidumbre cuantificada de la PDF prior es actualizada usando el teorema de Bayes, que usa la evidencia (observaciones) para actualizar la probabilidad de que una hipótesis sea verdadera, es decir, de que el valor de un parámetro se de; obteniendo, de esta manera, la distribución posterior del parámetro de interés. Este teorema se expresa matemáticamente, con unas variables de simbología más clarificativa, como: Ecuación 2. Expresión del teorema de Bayes Donde la E representa las evidencias (mediciones/observaciones), que serán las propiedades dinámicas de la estructura en el problema de aplicación de este documento y la H representa las hipótesis (parámetros inciertos). Esta expresión implica que la probabilidad de que la hipótesis sea correcta dada la evidencia (función posterior) depende de la relación entre la probabilidad de los casos donde la hipótesis sea correcta y la probabilidad donde se de la evidencia. No obstante, en los casos en los que los datos obtenidos experimentalmente sean constantes, la probabilidad de la evidencia no depende de los parámetros, por lo que ésta es fija. Siendo de interés la relación entre los parámetros y las observaciones, la probabilidad de la evidencia puede ser ignorada, simplificando la ecuación anterior a: 𝑃 (𝐻|𝐸)∝ 𝑃(𝐻)·𝑃(𝐸|𝐻) Ecuación 3. Simplificación del teorema de Bayes para su aplicación práctica Que, en términos de la simbología que venía usando, quedaría como: 𝑃 (𝜃|𝐷)∝ 𝑃(𝜃)·𝑃(𝐷|𝜃) Ecuación 4. Teorema de Bayes simplificado Donde el primer término del producto, la probabilidad de que se de la hipótesis, se trata de la función de densidad de probabilidad de los parámetros a actualizar, prior, y el segundo término, la probabilidad de que se de la evidencia si se da la hipótesis como correcta, la función de probabilidad de la evidencia dada la existencia de los parámetros inciertos. La distribución posterior representa la distribución actualizada de los parámetros del modelo tras la obtención de mediciones. 3.1.1 Distribución prior Refleja el conocimiento a priori o la hipótesis inicial sobre los parámetros del modelo antes de considerar ninguna medida experimental, pudiendo provenir de experiencia previa, testeo en escala de laboratorio, etc. Teóricamente, se podría usar cualquier tipo de distribución dependiendo de la cantidad de información disponible, sin embargo, en la práctica, las distribuciones uniformes y normales son las más comunes (Lye et al., 2021). Caracterizando las ventajas del uso de estas dos últimas distribuciones: - Uso de una distribución uniforme Puede ser lo más recomendado en caso de que el conocimiento que tengamos de el/los parámetro/s a actualizar sean los valores de los límites superiores e inferiores (con mayor o menor exactitud). La ventaja, además, de tomar esta distribución como prior, la posterior será, simplemente, proporcional a la función de probabilidad. - Uso de una distribución normal Su uso se dedica a los casos donde el conocimiento consiste en los valores de media y error relativo de los parámetros a actualizar. 38 Es por ello que, en los casos prácticos donde usaremos la inferencia bayesiana, reconociendo la característica aleatoria de las variables a actualizar, asumiremos que se distribuyen mediante una tipología normal. 3.1.2 Función de probabilidad Representa el grado de parecido entre los valores experimentales (evidencias) y los resultados del modelo numérico. Para un enfoque determinista, la función de probabilidad evalúa la diferencia entre los valores medidos y la salida del modelo numérico para cada serie de posibles valores que tienen los parámetros a actualizar. Así, la función de probabilidad es una función de θ (parámetros) y no de los datos experimentales (D). Si suponemos que las mediciones Dk son independientes e idénticamente distribuidas, la función de probabilidad se expresa de la siguiente manera (Lye et al., 2021): 𝑃(𝐷|𝜃)= ∏𝑃(𝐷𝑘|θ) 𝑛 𝑘=1 Ecuación 5. Expresión de la función de probabilidad para Dk independientes e idénticamente distribuidas Entonces, si asumimos que la diferencia entre las lecturas experimentales y los resultados del modelo (M(𝜃)) siguen una distribución normal de media cero y una varianza fija, la función de probabilidad que elegiremos será una de distribución también normal, obtendríamos la expresión: 𝑃(𝐷|𝜃)= ∏(1 𝜎𝑘· √2·𝜋) 𝑛 𝑘=1 · 𝑒−∑(𝐷𝑘−𝑀(𝜃))2 2· 𝜎𝑘2 𝑛 𝑘=1 Ecuación 6. Expresión de la función de probabilidad para una distribución normal El valor de σk puede estimarse mediante un ajuste de mínimos cuadrados del término diferencia entre los valores experimentales y los numéricos o, como de costumbre, basar su valor a partir de experiencia previa o buen juicio ingenieril. El acercamiento más común suele ser el de definir su valor como la desviación estándar del error comentado (Lye et al., 2021). Por último, es de vital importancia recalcar el, más apropiado, uso de la función de probabilidad en su formato logarítmico, por su crucial rol en la eficiencia computacional al transformar la expresión gobernada por productos en una predominante de operaciones de suma, más rápidas computacionalmente que las primeras operaciones: log (𝑃(𝐷|𝜃))= ∑ log (𝑃(𝐷𝑘|θ)) 𝑛 𝑘=1 Ecuación 7. Función de probabilidad logarítmica 3.1.3 Distribución posterior Representa la distribución actualizada de los parámetros del modelo basándose en la nueva información obtenida de las observaciones. Usualmente, esta distribución no encaja especialmente con una función de distribución bien, paramétricamente, conocida, por lo que resultaría muy costoso computacionalmente trabajar con ella. De esta manera, la manera de proceder es mediante el uso de una técnica numérica para estimar distribuciones denominada como el método de Monte Carlo que, generalmente, se usa para evaluar las estadísticas del parámetro a actualizar más que la distribución posterior completa. 2.4. Comparativa Para cerrar el capítulo, cabe recalcar que la aproximación obtenida a partir de la resolución del método de máxima verosimilitud, a pesar de sus ventajas como método analítico, es una solución local, mientras que, los resultados que obtenemos aproximando la función de distribución posterior representan una solución global. Sin embargo, cuando la función de probabilidad está bien definida, los resultados obtenidos a partir del método frecuentista son muy parecidos a los obtenidos por el bayesiano (resolviendo la función posterior) (Marwala et al., 2017). Una conclusión que podremos apreciar en el apartado reservado a la aplicación práctica de todo este marco teórico. No obstante, como hecho diferenciador principal, el método de máxima verosimilitud se ajusta expresamente a los datos, por lo que no considera una vista probabilística del problema de actualización, algo que si ofrece el método bayesiano. Esto es una limitación a considerar pues la medición experimental viene acompañada de errores por parte de los sensores; el método determinista, al considerar valores de esas características, lo asume como verdadero y desarrolla el ajuste considerándolo. Sin embargo, contrariamente, el método bayesiano considera probabilidades, por lo que, según la dispersión de los resultados obtenidos, evaluables en función de la distribución hallada, se podrá valorar cuán buena es la estimación lograda. De ahí surgen técnicas de muestreo, como los métodos Monte Carlo, usadas para resolver distribuciones posterior complejas y su ámbito de aplicación reside en aquellas situaciones donde no se puede estimar una solución mediante un método determinista y conviene, precisamente por esa complejidad, el uso de herramientas que generen poblaciones aleatorias en masa que evalúen la solución muchas veces, superando en idoneidad el estudio analítico y exacto de una distribución ardua. Esto significa que no llevarán a una única solución, sino que existirán varios resultados posibles. Estos métodos serán los protagonistas de este trabajo y, por lo tanto, ahondaremos más en ellos en los próximos capítulos. 40 4 MÉTODOS MONTE CARLO Como hemos visto en el capítulo anterior, las estadísticas bayesianas se usan para obtener la probabilidad de los parámetros desconocidos de un modelo de elementos finitos dados unos datos experimentales, es decir, lo que llamamos función de densidad de probabilidad posterior. El acercamiento bayesiano la estima a partir de la información que tenemos de los parámetros desconocidos que gobiernan el problema (prior), la distancia entre el modelo y los datos reales (función de probabilidad), cuyo uso implica considerar la diferencia entre los valores experimentales (realidad) y los resultados numéricos (digital) desde un punto de vista probabilístico, algo que cobra sentido cuando reconocemos que ambas variables son desconocidas, y los datos medidos (evidencia). Sin embargo, conseguir esto es complejo analíticamente, por lo que se recurre a métodos que llevan a cabo simulaciones Monte Carlo (Marwala et al., 2017). 4.1 Simulaciones Monte Carlo Son técnicas cuya aplicación reside en el estudio de la respuesta de modelos a entradas generadas aleatoriamente. Su funcionamiento es simple; consiste en la generación aleatoria de una cantidad N de muestras y ejecutando, para cada una de ellas, una simulación, se analiza las salidas que da como respuesta. Traduciendo esto a nuestro caso de estudio, el proceso sería el mismo, pero extrayendo las muestras a partir de una función de densidad de probabilidad conocida y obteniendo, una vez evaluadas todas, la distribución de ellas. El uso de estas técnicas tiene varias ventajas como: - La convergencia asociada con la simulación Monte Carlo es independiente del número de variables aleatorias - Se pueden paralelizar los cálculos numéricos - Son aplicables a cualquier problema, incluidos los no lineales La problemática reside en que la obtención de muestras independientes de la posterior es difícil al usar únicamente este tipo de simulaciones pues esta distribución es conocida sólo a efecto implícitos (requiere la evaluación de la distribución prior y la de la función de probabilidad), por lo que no es posible generar muestras directamente de la posterior (Lye et al., 2021). Ante esto surge unas técnicas de muestreo para aproximar la inferencia en el contexto de la actualización bayesiana de modelos, entre las cuales destacamos unas de las más conocidas, que serán fruto de estudio en este documento. Nombrándolos mediante la terminología por la que se reconocen universalmente: - Markov Chain Monte Carlo (MCMC) - Transitional Markov Chain Monte Carlo (TMCMC) - Sequential Monte Carlo (SMC) 4.2 Método Monte Carlo basado en cadenas Markov (Markov Chain Monte Carlo - MCMC) El método Monte Carlo basado en cadenas Markov (MCMC) es una técnica Monte Carlo de paso aleatorio que se implementa mediante la producción de una cadena Markov para estimar la función de distribución de probabilidad posterior del problema de actualización de un modelo de elementos finitos. Por lo tanto, veremos cómo es un proceso de paso aleatorio que usa dos técnicas principalmente: una simulación Monte Carlo y el proceso Markov. Una cadena Markov es una secuencia de variables tal que la distribución de 𝜃𝑡 depende únicamente de 𝜃𝑡−1, es decir, es un proceso estocástico donde la probabilidad de un estado depende del estado inmediatamente anterior. Este proceso es el que caracterizará la evolución del sistema en el método MCMC, que estará formado por variables aleatorias que han de ser actualizadas. Cadena en la cual se busca obtener que su distribución sea estacionaria, eliminando los primeros valores simulados (lo que denotaremos como el burn-in, del método, en su aplicación) al no encontrarse aún en ese estado estabilizado. 4.2.1 Algoritmo Metropolis-Hastings El nuevo estado de la cadena Markov será aceptado o rechazado a partir de la aplicación del algoritmo que emplea como método de muestreo, el algoritmo de paso aleatorio Metropolis-Hastings. La idea general de su uso consiste en proponer una función de distribución de probabilidad (que será la que dicte el criterio de aceptación o rechazo del nuevo valor propuesto como el siguiente estado en la cadena Markov) y usarla para generar valores propuestos. Podemos concluir, entonces, que el algoritmo M-H consiste, fundamentalmente, en dos pasos: - Obtención de una muestra del paso de la densidad propuesta - Aceptación o rechazo del paso del estado propuesto Para poder visualizar más claramente el flujo de este algoritmo, cuya distribución estacionaria fuese 𝑃(𝜃|𝐷), podría esquematizarse de la siguiente manera: 1º - Abrimos el proceso con un valor inicial 𝜃0 2º - En la iteración t, obtenemos el valor candidato aleatorio 𝜃∗ de la cadena Markov de una función de densidad de probabilidad propuesta 𝑞(𝜃∗|𝜃𝑡−1) 3º - Actualizamos el modelo de elementos finitos para obtener los nuevos resultados y evaluamos la distribución posterior para la muestra propuesta. Seguidamente, calculamos la aceptación del valor generado mediante: 𝛼(𝜃∗,𝜃𝑡−1)=min{ 1 , 𝑃(𝜃∗|𝐷)·𝑞(𝜃𝑡−1|𝜃∗) 𝑃(𝜃𝑡−1|𝐷)·𝑞(𝜃∗|𝜃𝑡−1) } Ecuación 8. Criterio de aceptación o rechazo del algoritmo M-H Donde 𝑃(𝜃∗|𝐷) es el valor de la función posterior evaluado para la muestra 𝜃∗, análogamente para 𝑃(𝜃𝑡−1|𝐷) pero evaluado para la muestra 𝜃𝑡−1. Mientras que 𝑞(𝜃𝑡−1|𝜃∗) representa la probabilidad de muestrear 𝜃𝑡−1 siendo la muestra actual 𝜃∗ y viceversa para el caso de 𝑞(𝜃∗|𝜃𝑡−1). Considerando que la probabilidad de la evidencia no depende de los parámetros, podremos obviarla, lo que dejaría a la expresión de la siguiente manera: 𝛼(𝜃∗,𝜃𝑡−1)=min{ 1 , 𝑃(𝜃∗|𝐷) 𝑃(𝜃𝑡−1|𝐷) } Ecuación 9. Criterio de aceptación o rechazo del algoritmo M-H simplificado 48 Ilustración 14. Tablero compuesto real De esta manera, elegimos homogeneizar a hormigón por ser el material mayoritario, obteniendo un tablero equivalente de hormigón cuyas propiedades físicas, que influirán en las matrices de masa y rigidez de la estructura, contemplan la verdadera existencia de dos materiales. La forma en la que se lleva a cabo este proceso es manteniendo constantes, tras este proceso, el módulo de deformación y la inercia del tablero (en ambas direcciones; factor que consideraremos en el siguiente subapartado). Primero, como hemos comentado, homogeneizamos a hormigón, consiguiéndolo mediante la modificación de una dimensión física del elemento, esto es, el ancho del acero. Esto consigue alterar el área de las caras de la pieza de estudio, influyendo en su inercia y, consecuentemente, en las matrices, comentadas anteriormente, que gobiernan la ecuación de equilibrio dinámico. Optaremos por la alteración del espesor del tablero por facilidad de resolución al tratarse de la geometría de una chapa grecada y, para lograr esto, se calcula un coeficiente de homogeneización que representa la relación entre el módulo de elasticidad del acero y el del hormigón de la siguiente manera: 𝑛= 𝐸𝑎 𝐸𝑐,28= 210000 𝑀𝑃𝑎 32836,568 𝑀𝑃𝑎=6,395 Ecuación 15. Coeficiente de homogeneización Donde n es el coeficiente de homogeneización, 𝐸𝑎 es el módulo de deformación del acero estructural y 𝐸𝑐,28 es el módulo de deformación longitudinal inicial del hormigón a los 28 días de edad. Así, si pasamos todo el acero de la sección a hormigón, sabiendo que las nuevas dimensiones se calcularían siguiendo la siguiente expresión: 𝐴𝑛𝑐ℎ𝑜𝐻,𝑒𝑞= 𝐴𝑛𝑐ℎ𝑜𝐴𝑐𝑒𝑟𝑜·𝑛 Ecuación 16. Obtención del ancho equivalente en la homogeneización Obtendremos la siguiente geometría: Ilustración 15. Geometría de un nervio homogeneizada a hormigón Una vez tengo la sección completa de hormigón, puedo calcular su inercia. La idea tras el análisis de un único nervio es la de obtener la inercia por metro lineal, hallando que la inercia total, respecto al eje y, por metro lineal es: 𝐼𝑦𝑦,𝑡𝑜𝑡= 7,3337· 10−5 𝑚4 𝑚𝑙 5.2.2 Ortotropía La idea de homogeneizar nace de la necesidad de contemplar un único elemento a la hora de realizar cálculos estructurales. Sin embargo, para mantener la veracidad del modelo respecto a la realidad, las propiedades físicas del módulo de deformación y de inercia modeladas han de ser las adecuadas. Por lo que, daré por satisfecha esa condición, si se cumple la siguiente ecuación, donde participan los parámetros influyentes en la matriz de rigidez: 𝐸· 𝐼 |𝑟𝑒𝑎𝑙= 𝐸· 𝐼 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜 Ecuación 17. Igualdad necesaria para la veracidad del modelado Donde E e I son los módulos de deformación e inercia respectivamente, tanto en el espacio real como el modelado. Sin embargo, el tablero posee una característica peculiar que se trata de la, ya mencionada, ortotropía, al tener dispuesto un sistema de nervios longitudinales. De esta manera, aproximaremos el siguiente comportamiento mecánico: - Para la dirección longitudinal (x): trabajará la losa de hormigón y la chapa de acero colaborante - Para la dirección transversal (y): trabajará, únicamente, la losa de hormigón (la chapa grecada no ofrece inercia en esa dirección; su función se reducirá a la de transmisión de esfuerzos) 50 Así, diseccionamos la expresión anterior en la siguiente: 𝐸𝑥·𝐼𝑦𝑦 |𝑟𝑒𝑎𝑙= 𝐸𝑥·𝐼𝑦𝑦 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜 𝐸𝑦·𝐼𝑥𝑥 |𝑟𝑒𝑎𝑙= 𝐸𝑦·𝐼𝑥𝑥 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜 Ecuación 18. Sistema de ecuaciones definitivo a cumplir para la veracidad del modelado Si modelamos, al igual de como hemos homogeneizado a hormigón, 𝐸𝑥,𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜= 𝐸𝑐,28= 32836,568 𝑀𝑃𝑎. Observamos que la primera ecuación anterior nos deja con: 𝐸𝑥 |𝑟𝑒𝑎𝑙= 𝐸𝑥 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜 Ecuación 19. Primera condición definida para el modelado Igualdad que, irremediablemente, nos impone: 𝐼𝑦𝑦 |𝑟𝑒𝑎𝑙= 𝐼𝑦𝑦 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜 Ecuación 20. Condición consecuente impuesta por la definición anterior Por simplicidad, a todos los niveles, vamos a modelar el tablero como un tablero prismático de canto constante, cuya inercia respecto al eje y será, como sabemos, 𝐼𝑦𝑦 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜= 1 12·𝑏·ℎ3. De esta manera, si mantenemos el ancho real del tablero, la única incógnita será el espesor de este. Recordando que ya teníamos calculado la inercia respecto al eje y, obtenemos: 𝐼𝑦𝑦 |𝑟𝑒𝑎𝑙= 𝐼𝑦𝑦 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜= 7,3337· 10−5 𝑚4 𝑚𝑙→ℎ𝑒𝑞=9,583 𝑐𝑚 Ecuación 21. Obtención del espesor equivalente del tablero modelado Sin embargo, como hemos comentado anteriormente, la rigidez longitudinal y transversal del tablero ha de modelarse acorde a la real. Se verificará, a continuación, la dirección referente al eje y. Sabemos que, modelando el tablero con un canto constante, su inercia, por metro lineal, referente al eje x sería la misma que 𝐼𝑦𝑦. Sin embargo, como mencionamos anteriormente, en realidad, en esta dirección sólo trabajará la losa de hormigón, por lo tanto, el canto efectivo del tablero será su espesor, 4 cm, obteniendo: 𝐼𝑥𝑥 |𝑟𝑒𝑎𝑙= 1 12·𝑏·ℎ3=5,333· 10−6 𝑚4 𝑚𝑙 Esto generará un desajuste en cuanto a la verosimilitud de los pesos modelados, el cual solucionaremos en el siguiente subapartado. Puesto que se tiene que cumplir la condición anteriormente mencionada, podemos obtener: 𝐸𝑦 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜= 𝐸𝑦·𝐼𝑥𝑥 |𝑟𝑒𝑎𝑙 𝐼𝑥𝑥 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜= 32836,568 𝑀𝑃𝑎 · 5,333· 10−6 𝑚4 𝑚𝑙 7,3337· 10−5 𝑚4 𝑚𝑙 ~ 2388 𝑀𝑃 Ecuación 22. Obtención de la última propiedad mecánica del material para su modelado 5.3 Masa La consideración de distintos espesores genera un desajuste en la masa modelada de la estructura respecto a la real, por ello, calcularemos la real: 𝑉ℎ𝑜𝑟𝑚𝑖𝑔ó𝑛/𝑛𝑒𝑟𝑣𝑖𝑜|𝑟𝑒𝑎𝑙= 𝐴𝐻· 1 𝑚=0,04 𝑚· 0,16 𝑚 + 0,058 𝑚 · 0,065 𝑚=0,01017 𝑚3 𝑚𝑙 1 𝑛𝑒𝑟𝑣𝑖𝑜 𝑚ℎ𝑜𝑟𝑚𝑖𝑔ó𝑛/𝑛𝑒𝑟𝑣𝑖𝑜|𝑟𝑒𝑎𝑙= 𝜌𝐻· 𝑉𝐻=2400 𝑘𝑔 𝑚3·0,01017 𝑚2=24,408 𝑘𝑔 𝑚𝑙 1 𝑛𝑒𝑟𝑣𝑖𝑜 𝑉𝑎𝑐𝑒𝑟𝑜/𝑛𝑒𝑟𝑣𝑖𝑜|𝑟𝑒𝑎𝑙= 𝐴𝑎𝑐𝑒𝑟𝑜· 1 𝑚=0,001 𝑚 ·0,095 𝑚 + 0,001 𝑚 · 0,065 𝑚=0,00016 𝑚3 𝑚𝑙 1 𝑛𝑒𝑟𝑣𝑖𝑜 𝑚𝑎𝑐𝑒𝑟𝑜/𝑛𝑒𝑟𝑣𝑖𝑜|𝑟𝑒𝑎𝑙= 𝜌𝑎𝑐𝑒𝑟𝑜· 𝑉𝑎𝑐𝑒𝑟𝑜=8050 𝑘𝑔 𝑚3· 0,00016 𝑚2=1,288 𝑘𝑔 𝑚𝑙 1 𝑛𝑒𝑟𝑣𝑖𝑜 𝑚𝑡𝑜𝑡𝑎𝑙|𝑟𝑒𝑎𝑙=(𝑚ℎ𝑜𝑟𝑚𝑖𝑔ó𝑛 𝑛𝑒𝑟𝑣𝑖𝑜 |𝑟𝑒𝑎𝑙+𝑚𝑎𝑐𝑒𝑟𝑜 𝑛𝑒𝑟𝑣𝑖𝑜|𝑟𝑒𝑎𝑙)·19 𝑛𝑒𝑟𝑣𝑖𝑜𝑠=488,224 𝑘𝑔 𝑚𝑙 Ecuación 23. Cálculo de volúmenes y masas reales de cada material que compone el tablero compuesto Mientras que el modelado consiste en: 𝑚𝑡𝑎𝑏𝑙𝑒𝑟𝑜|𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜= 𝜌𝐻· 𝑉𝐻=2400 𝑘𝑔 𝑚3· 3 𝑚 · ℎ𝑒𝑞· 1𝑚 =689,976 𝑘𝑔 𝑚𝑙 Ecuación 24. Cálculo de la masa del tablero modelado Como se puede apreciar, ambos valores no coinciden. Esto es fácilmente corregible modificando la densidad con la que modelaremos el tablero. 𝑚𝑡𝑜𝑡𝑎𝑙|𝑟𝑒𝑎𝑙=488,224 𝑘𝑔 𝑚𝑙= 𝑚𝑡𝑎𝑏𝑙𝑒𝑟𝑜|𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜= 𝜌𝑡𝑎𝑏· 3 𝑚 · ℎ𝑒𝑞· 1𝑚→ 𝜌𝑡𝑎𝑏=1698,23 𝑘𝑔 𝑚3 Ecuación 25. Obtención de la densidad equivalente modelada del tablero para el ajuste de masas Es este el proceso que ha permitido corregir la masa del tablero y que permitirá el uso de esta nueva densidad para modelar el material del tablero en lugar de la tradicional del hormigón. 5.4 Otras propiedades Este punto se centrará en el análisis del resto de propiedades que están bajo consideración en un cálculo estructural de la índole que respecta a este documento. 52 5.4.1 Coeficiente de Poisson Definimos la siguiente propiedad: 𝜈𝑥𝑦=𝜈ℎ𝑜𝑟𝑚𝑖𝑔ó𝑛=0,2 y, sabiendo que se tiene que cumplir la siguiente igualdad, obtenemos: 𝜈𝑦𝑥 𝐸𝑦 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜= 𝜈𝑥𝑦 𝐸𝑥 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜→ 𝜈𝑦𝑥=𝜈𝑥𝑦· 𝐸𝑦 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜 𝐸𝑥 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜 =0,0145447 Ecuación 26. Condiciones de igualdad a cumplir referentes a los coeficientes de Poisson y módulos de deformación modelados y obtención de los primeros 5.4.2 Módulo de torsión Sabiendo que el módulo de torsión del material principal es 𝐺𝑥𝑦, lo podremos calcular como sigue: 𝐺𝑥𝑦=𝐸𝑥 |𝑚𝑜𝑑𝑒𝑙𝑎𝑑𝑜 2·(1 + 𝜈𝑥𝑦)=13681,9 𝑀𝑃𝑎 Ecuación 27. Obtención del módulo de torsión del material modelado 5.5 Modelado Una vez hemos conseguido definir todas las propiedades físicas y mecánicas que dominan el problema, así como su geometría, ya tenemos lo necesario para modelar la estructura al completo. El modelo de elementos finitos de la pasarela a analizar será el siguiente, el cual podremos descomponer en los siguientes elementos: - Unas vigas cuadrangulares longitudinales de 400 x 200 x 10 mm Ilustración 16. Modelado - Vigas cuadrangulares longitudinales - Unas vigas cuadrangulares transversales y verticales que conforman los módulos de la estructura de 300 x 200 x 12 mm Ilustración 17. Modelado - Vigas cuadrangulares transversales y verticales - Unas vigas transversales en dirección diagonal con sección en L de 150 x 150 x 10 mm las inferiores y de 120 x 120 x 15 mm las superiores: Ilustración 18. Modelado - Vigas transversales en dirección diagonal con sección en L – Visión global Ilustración 19. Modelado - Vigas transversales en dirección diagonal con sección en L – Visión local - Y, finalmente, un tablero cuyas propiedades hemos definido en el apartado anterior: Ilustración 20. Modelado – Tablero – Visión global Ilustración 21. Modelado – Tablero – Visión local 54 Obteniendo el conjunto estructural al completo: Ilustración 22. Modelo de elementos finitos completo Habiendo hablado del material, sus propiedades, definido la geometría de la estructura y los elementos que la componen con sus respectivas secciones, sólo falta por comentar las condiciones de contorno a las que está sometida. Supondremos dispuestos unos aparatos de apoyo en las esquinas del tablero, es decir, en los extremos de este, uno para cada rama longitudinal de vigas, en los que impediremos de manera rígida el movimiento vertical. Impidiendo, en los nodos correspondientes a los apoyos, la libertad de movimiento en la dirección vertical (eje z), nos planteamos modelar las condiciones de contorno, en dichos puntos, en las direcciones x e y mediante unos apoyos flexibles. Sin embargo, debido a la poca influencia que ha de tener la rigidez al desplazamiento transversal en el estudio de los modos de vibración vertical, he decidido definir unas condiciones rígidas, en los nodos comentados, para la dirección transversal y también. Ilustración 23. Modelado – Condiciones de contorno en direcciones transversal y vertical Finalmente, sólo queda definir las condiciones de contorno en la dirección longitudinal que, para una correcta consideración del comportamiento estructural, debo disponerlo como un apoyo flexible, es decir, construir un sistema en el que un muelle domine las condiciones de sustentación. De esta manera, unimos los nodos correspondientes a los aparatos de apoyo con un muelle en dirección longitudinal, que será el responsable de modelar la rigidez longitudinal de la estructura, que, a su vez, irá unido a un empotramiento necesario por cuestiones de computación. Quedando como: Ilustración 24. Modelado - Condiciones de contorno en dirección longitudinal - Elemento muelle unido a los extremos del tablero y a un empotramiento 56 5.6 Análisis modal Con anterioridad, he comentado que las propiedades más usadas a la hora de actualizar un modelo numérico son las dinámicas. De esta manera, una estructura, en términos dinámicos, se caracteriza por sus frecuencias naturales y por la forma de los modos de vibración, es decir, las formas básicas en las que vibra una estructura, que cuantificaremos en función de la amplitud del movimiento que se de en ciertos puntos de la estructura, y los periodos asociados (Javier et al., 2020). El estudio de estas dos propiedades y su determinación construyen lo denominado como análisis modal, que dependerá de las propiedades y condiciones de contorno de la estructura. Así, someteremos a nuestra pasarela a este tipo de análisis. Habiendo escogido como objetivo de estudio los cuatro primeros modos de vibración verticales y sus frecuencias correspondientes: Ilustración 25. Análisis modal - Primer modo de vibración vertical – Correspondiente a una frecuencia de 3,78904 Hz Ilustración 26. Análisis modal - Segundo modo de vibración vertical – Correspondiente a una frecuencia de 12,6789 Hz Ilustración 27. Análisis modal - Tercer modo de vibración vertical – Correspondiente a una frecuencia de 21,9883 Hz Ilustración 28. Análisis modal - Cuarto modo de vibración vertical – Correspondiente a una frecuencia de 53,3202 Hz Esta amplitud del movimiento que producen los modos de vibración de la estructura se cuantifica en función del nodo que se esté considerando. Para ello, elegiremos una serie de nodos característicos que serán los encargados de representar estos modos adecuadamente puesto que, para computar esta propiedad dinámica, necesito valores numéricos concretos y resulta innecesario considerar todos los nodos de los que se compone el modelo. De esta manera, escogemos 30 nodos correspondientes al tablero de la estructura, que siguen la línea que forman las vigas longitudinales de éste. Serán 15 nodos en cada rama longitudinal, correspondientes a aquellos que separan los módulos de la estructura: Ilustración 29. Análisis modal - Elección de nodos representativos para la determinación de la forma de los modos de vibración Así, los modos de vibración podrán ser comparados a través de: - El desplazamiento vertical, mediante un análisis nodo a nodo, de los 30 nodos escogidos pertenecientes al tablero que experimentarán las deformaciones asociadas a los modos de vibración verticales y que, en su totalidad, representarán con suficiente exactitud, el respectivo modo y las diferencias respecto al resto de ellos. - El MAC (Modal Assurance Criterion), que es un indicador estadístico que sirve para medir el grado de similitud que tienen dos modos de vibración distintos en términos de sus movimientos nodales. Es una herramienta que implicará un análisis modo a modo (aunque, para su determinación, requiere del estudio de la amplitud de cada nodo). 𝑀𝐴𝐶𝑖(𝜃)= (𝜙𝑛𝑢𝑚,𝑖 𝑇(𝜃)· 𝜙𝑒𝑥𝑝,𝑖)2 (𝜙𝑛𝑢𝑚,𝑖 𝑇(𝜃)· 𝜙𝑛𝑢𝑚,𝑖(𝜃))·(𝜙𝑒𝑥𝑝,𝑖 𝑇· 𝜙𝑒𝑥𝑝,𝑖) Ecuación 28. Expresión del indicador estadístico MACG Donde 𝜙𝑛𝑢𝑚,𝑖(𝜃) es el modo de vibración numérico i que depende de los parámetros 𝜃 y 𝜙𝑒𝑥𝑝,𝑖 representa el modo de vibración experimental i cuyos valores no dependen de ninguna variable pues es un dato conocido y fijo. El MAC es un indicador usado típicamente para la validación de modelos experimentales, así como el estudio de correlación entre modelos analíticos, entre otros. De igual manera, para ambos casos, las amplitudes de los nodos serán normalizadas de manera relativa al máximo desplazamiento de cada modo, de manera que el nodo que experimente la mayor amplitud del modo pasará a tener un desplazamiento del valor unidad y, el resto de los nodos, un valor proporcional a este. 64 5.8 Método determinista Ya es bien sabido que la resolución consistirá en un ejercicio de optimización, donde minimizaremos la diferencia entre el modelo numérico y la realidad. La falta de datos experimentales no tiene por qué ser crucial en un caso puramente académico pues, grosso modo, conocemos los valores por los que se tienen que mover las variables. - Para la definición del modelo numérico se asumirán unos valores típicos de los parámetros escogidos para calibrar - Para simular una medición real de la estructura, primero, se proponen unos valores que podrían representar la realidad y que distan lo suficiente de los tradicionales para, así, poder desarrollar el ejercicio cómodamente y obtener, finalmente, unas propiedades dinámicas. Estas propiedades representarán los datos que podría haber obtenido in-situ y serán los datos objetivo a los que el modelo numérico tratará de aproximarse en la actualización de éste. Por simplicidad de resolución, en lugar de obtener el parámetro actualizado, desarrollo el código para que el programa me actualice el número por el cual he de multiplicar el valor numérico para hallar el experimental. Es decir: 𝑉𝑎𝑙𝑜𝑟 𝑛𝑢𝑚é𝑟𝑖𝑐𝑜 |𝑐𝑜𝑛𝑜𝑐𝑖𝑑𝑜 · 𝑐𝑜𝑒𝑓𝑖𝑐𝑖𝑒𝑛𝑡𝑒 = 𝑉𝑎𝑙𝑜𝑟 𝑟𝑒𝑎𝑙 |𝑑𝑒𝑠𝑐𝑜𝑛𝑜𝑐𝑖𝑑𝑜 Ecuación 30. Igualdad que seguirá el resultado obtenido por el método determinista Sin embargo, considerando que conozco tanto los valores numéricos de partida como los valores reales que me darán las propiedades dinámicas a las que aspiro (desconocidos en un problema real y objetivo del procedimiento), puedo saber, de primera mano, el coeficiente exacto que me debería de dar para cada parámetro. Estableciendo, así, con los siguientes valores: Ilustración 36. Valores numéricos y reales de los parámetros a actualizar El método requerirá, como datos de partida para realizar la actualización, unos rangos de variación entre los que consideramos que se moverán los parámetros. De esta manera, al conocer ambas familias de datos por ser un ejemplo académico, conseguiré mejorar la eficiencia computacional otorgándole un rango de búsqueda más acertado desde el primer momento. Los rangos de variación suelen concebirse entre 0 y 1. Sin embargo, en esta aplicación, una variación de 0,5 implica reducir un módulo de deformación a la mitad, un valor totalmente fuera del orden de magnitud de este campo. De esta manera, observando los valores reales que hemos asumido, que me darán las propiedades dinámicas que tomaríamos como medias experimentales, vemos que rondan un 5-10% de variación respecto a los valores típicos de los parámetros escogidos. Así, definiremos los rangos de variación considerando una variación de la mitad a la comentada, generando un rango de 0,5 en su totalidad, es decir, una variación de 0,25 para el límite superior e inferior a partir del coeficiente exacto que he calculado. El flujo de trabajo consiste en la ejecución de un código, en el programa de cálculo numérico Matlab, mediante el cual realizamos un proceso iterativo donde minimizamos una función, con unas opciones de optimización previamente definidas, mediante el uso de un comando que usa un algoritmo genético. El código seguirá su ejecución hasta que haya testeado, para toda la población, el número de iteraciones indicadas. Particularizando, el algoritmo genético se ejecutará a partir de la función ga, predefinida en el software, que consiste en un algoritmo genético de optimización que trata de encontrar el mínimo de una función objetivo dada gracias al uso de unas poblaciones y generaciones de individuos. Estos, llamados, individuos, son cualquier punto a los que se le puede aplicar la función objetivo siendo, el valor de ésta para un individuo, su puntuación. El código desarrollará su ejecución hasta que haya testeado toda la población de individuos el número de iteraciones indicadas, puesto que, en cada iteración, el algoritmo realiza una serie de cálculos en la población actual para producir una nueva. Este algoritmo genético permite configurar su funcionamiento mediante la definición de algunos parámetros. Un ejemplo de ello resulta ser las opciones de reproducción, que especifica cómo el algoritmo genera descendientes para la siguiente generación, aunque nos centraremos, especialmente, en la variación del tamaño de la población y el máximo número de generaciones como consideración de precisión de la solución y, a su vez, criterio de parada. De esta manera, el algoritmo genético lanzará tantos puntos (individuos) como tamaño de la población indique, tantas veces como número de generaciones defina donde, cada punto, implicará una entrada de datos en el sistema. Esto significa que, para cada iteración, el algoritmo genera un vector de cuatro valores (al tener cuatro parámetros) que serán los cuatro coeficientes, comprendidos entre los rangos superior e inferior, que multiplican los valores numéricos de los parámetros a actualizar, con la intención de aproximar los valores reales. El software, para cada generación, calculará el valor medio de la función objetivo y su mejor valor o, lo que es lo mismo, su valor mínimo. Una vez entendemos el proceso de manera cualitativa, podemos especificar para nuestro caso de aplicación. En un análisis más profundo del proceso, se lleva a cabo una llamada a Ansys, con previa sobreescritura del vector que contiene los coeficientes, el cual realiza un análisis modal del modelo de elementos finitos de la estructura, obteniendo las frecuencias naturales y modos de vibración de esta. Una vez obtenidos, se compara con las propiedades dinámicas de las mediciones experimentales (introducidas, previamente, a mano). Este método de comparación resulta ser la definición de una función de residuos, es decir, la diferencia entre los valores numéricos recién calculados y los experimentales a los que deseamos llegar (como ya comentamos; la función objetivo resulta ser la diferencia entre ambos modelos, la cual intentaremos minimizar). Puesto que estamos considerando tanto frecuencias naturales como los modos de vibración, considerando estos últimos mediante los desplazamientos dados en los nudos representativos del tablero normalizados a la unidad, la función residual consistirá en una media aritmética de los sumatorios de los residuos de ambas variables, ponderando con mayor trascendencia las frecuencias naturales por ser más representativas que los modos. La expresión relativa a las frecuencias de la estructura quedaría como: 𝑓𝑛𝑖𝑟=0,9 · 𝑓𝑛𝑖𝑛𝑢𝑚− 𝑓𝑛𝑖𝑒𝑥𝑝 𝑓𝑛𝑖𝑒𝑥𝑝 Ecuación 31. Residuo de las frecuencias í-esimas numéricas y experimentales ponderadas al 90% Mientras que la diferencia de amplitud normalizada de desplazamiento del nodo i para el modo j ponderada al 10% quedaría tal que: 𝑚𝑜𝑑𝑖,𝑗 𝑟=0,1 · 𝑚𝑜𝑑𝑖,𝑗 𝑛𝑢𝑚− 𝑚𝑜𝑑𝑖,𝑗 𝑒𝑥𝑝 𝑚𝑜𝑑𝑖,𝑗 𝑒𝑥𝑝 Ecuación 32. Residuo de los modos de vibración i-ésimos numéricos y experimentales ponderados al 10% De esta manera, la función objetivo a minimizar quedaría finalmente como: 66 𝑓= 12 · √∑𝑓𝑛𝑖𝑟+ ∑ 𝑚𝑜𝑑𝑖𝑟 𝑛 𝑖=1 𝑛 𝑖=1 Ecuación 33. Función objetivo que minimizar Así el programa buscará converger la solución, iteración tras iteración, hasta conseguir minimizar lo máximo posible la diferencia entre el modelo numérico y el real, es decir, la función objetivo, siguiendo los anteriores criterios de parada de máximas generaciones definidas. En el caso que nos ocupa, para una población de 10 individuos y 40 generaciones, se obtuvieron los siguientes resultados: Ilustración 37. Evolución de las soluciones mediante el algoritmo genético en el método determinista Donde podemos apreciar que el sistema había encontrado una solución aparentemente satisfactoria para el transcurso de 25 generaciones. Nótese cómo el mejor valor de la función objetivo se mantiene constante prácticamente desde el comienzo y la media de ésta, sacada de los 10 individuos, va convergiendo hasta el valor mínimo. La solución optimizada, que minimiza la diferencia entre el modelo numérico y el real, para esta configuración del algoritmo genético, aporta los siguientes valores de los coeficientes: Ilustración 38. Coeficientes resultantes – Método determinista Comparando con el coeficiente exacto que podíamos esperar: Ilustración 39. Comparativa coeficientes teóricos y numéricos – Método determinista Observamos que los valores son considerablemente próximos a los teóricos, dando una idea de la suficiente precisión que hemos obtenido con el método. Comparando las frecuencias para estos parámetros aproximados respecto a las reales: Ilustración 40. Comparativa frecuencias numéricas y experimentales – Método determinista El archivo de texto que contiene las frecuencias naturales de la estructura correspondiente a los modos de vibración verticales, obtenido a partir del modelo numérico de elementos finitos con unos valores de los parámetros estimados pero resultado de la optimización que hemos comentado, da unos valores prácticamente equivalentes a los reales (dando por bueno valores por debajo del 5%), existiendo una diferencia más que despreciable y concluyendo con que los resultados son clavados a los esperados. Por otra parte, considerando, también, los modos de vibración, procedemos a estudiar el MAC. Organizando las coordenadas de los nodos seleccionados correspondientes a la posición deformada de los modos de vibración de interés ponderados a la unidad, podremos evaluar dicho indicador: Ilustración 41. Resultado de los indicadores MAC – Método determinista Apreciando que se consigue obtener un valor muy cercano a la unidad para los modos de vibración verticales escogidos, es decir, que estamos ante modos de vibración numéricos y experimentales casi idénticos. Finalmente, considerando esto y la positiva conclusión hallada en el cálculo de las frecuencias naturales, determinamos que la aproximación realizada con el método determinista es completamente satisfactoria para las características ya mencionadas. En un problema real, el objetivo sería saber los valores de los parámetros escogidos para actualizar. De esta manera, para el modelo que simula la pasarela entre Olivares y Villanueva del Ariscal, los resultados, plenamente satisfactorios, de la actualización mediante el método determinista, con los que obtenemos los parámetros actualizados de los módulos de deformación del acero y del tablero, en la zona central y lateral, en ambas 68 direcciones, son: Ilustración 42. Resultados – Variables actualizadas – Método determinista 5.9 Método bayesiano Como se comentó en el marco teórico de este documento, aunque todos los parámetros que definen este problema son variables aleatorias, se asumirá que siguen una distribución normal considerando que los conocimientos que tenemos sobre ellos son su valor medio aproximado junto con un margen de tolerancia de éste. Por otro lado, como se verá a continuación en los dos ejemplos de aplicación a desarrollar, la generación de los números aleatorios por parte de las simulaciones Monte Carlo se harán a partir de una distribución uniforme con la que se consigue otorgar la misma probabilidad de ocurrencia a cada valor propuesto. Una vez extraídos los valores de la funciones uniformes, se les aplicarán los algoritmos correspondientes a estas técnicas numéricas con previa definición de las funciones participantes en el teorema de Bayes. Obteniendo, según este teorema y la distribución uniforme prior, una distribución resultante de la posterior de tipo normal, al comparar con las series de datos experimentales normalmente distribuidas. Por otra parte, en el caso de los valores experimentales, si asumo que las variables aleatorias a considerar siguen cierta distribución, la respuesta estructural en la que base las mediciones también lo harán. Esto significa que al definir la distribución normal como aquella por la que se rigen los parámetros a calibrar, los que reflejan la respuesta con la que estudiamos el modelo, es decir, los modos de vibración y sus frecuencias naturales asociadas, también lo harán. Así, a modo de recapitulación, se enumera: - El conjunto estructural parte de unos parámetros de diseño que, aún siendo variables aleatorias, asumimos que se distribuyen de manera normal como método de aproximación para estudiar el problema. Declaración que implica que las propiedades dinámicas que gobiernan la respuesta de la estructura se comporten de la misma manera. - Las mediciones experimentales, de esta forma, consistirán en una serie de datos relativas a estas dos propiedades que seguirán una distribución normal. - Por parte del análisis numérico, las simulaciones Monte Carlo se encargarán de llevar a cabo la actualización del modelo de elementos finitos siguiendo el teorema de Bayes en el que: o Se partirá de una generación aleatoria de valores que se distribuirán uniformemente para otorgar esa equidad de probabilidad de ocurrencia. Esto será la función prior. o Con esos valores generados se llevará a cabo una simulación numérica con el que se obtendrán unos resultados que, posteriormente, se compararán con los experimentales. Esto será la función de probabilidad. o Se obtendrá la función posterior mediante la relación entre la función prior y la de probabilidad que, debido a sus naturalezas comentadas, otorgarán un resultado que se distribuye normalmente. Una vez clarificado el entorno estadístico que domina los problemas de aplicación, se presenta, a continuación, la puesta en práctica de éstos conceptos probabilísticos. 5.9.1 Edificio a cortante Mientras la estructura que aspiramos estudiar se trata de una pasarela ciclista, la implementación de los métodos estadísticos era un asunto complejo en términos de coste computacional considerando la magnitud del modelo numérico, así como mi carácter novel en la materia. De esta manera, se decidió llevar a cabo la realización de una primera aplicación práctica más sencilla numérica y conceptualmente hablando que se trató de un modelo, clásico, de un edificio a cortante. Las ventajas que tiene este ejemplo son, además de ser definido por un modelo numérico mucho más reducido y simplificado, la capacidad de obtener una solución analítica en términos de las propiedades dinámicas de la estructura, lo que permite la obtención de los resultados numéricos mediante su programación en el software Matlab. Así conseguí aliviar el flujo de trabajo al no necesitar la interacción Ansys- Matlab, algo crucial en términos de ahorrar tiempo de computación. Se define entonces el siguiente problema que representará un acercamiento, de resolución simplificada respecto al ejemplo de la pasarela, a los métodos bayesianos al ser más complejos que los métodos frecuentistas. El modelo de edificio a cortante viene definido por la siguiente configuración: Ilustración 43. Modelo de edificio a cortante bajo estudio Se trata de un modelo simplificado de la estructura de un edificio de cuatro plantas que modelaremos en un espacio bidimensional por simplificar el análisis y la gestión de datos. Siguiendo esta rama de enfocarnos en evitar la complejidad del problema, dotaremos al edificio de pilares de sección rectangular maciza puesto que el fin del ejercicio es puramente de testeo, por lo que unas propiedades geométricas que favorezcan la fácil resolución serán prioritarias. Por último, estimaremos que el peso de cada forjado es de 10000 kg y, conociendo el comportamiento mecánico del conjunto estructural para la configuración dada, podremos calcular la respuesta 70 dinámica del edificio de una manera analítica. Así, identificamos que el parámetro a calibrar en este ejercicio sería el ancho de los pilares rectangulares de cada planta, a los que llamaremos d1, d2, d3 y d4, refiriéndonos al ancho del pilar de la planta 1 a la 4 respectivamente. Con el fin de mejorar la comprensión de la resolución, describiré el proceso inverso de éste, empezando por el final: 1. Sabemos que las matrices de masa y rigidez de un edificio de cuatro plantas se ven de la siguiente forma: Ilustración 44. Matriz de masa - Edificio a cortante Ilustración 45. Matriz de rigidez - Edificio a cortante Donde m es la masa de los forjados y k es la rigidez de cada planta (marcando, el subíndice, la planta de la que se trata) 2. Las componentes de la rigidez se obtienen según, como pude comprobar según los prontuarios de matrices de rigidez, el valor 12·𝐸·𝐼 𝐿3. Especificando para cada planta, al tener cuatro pilares, obtendremos: 𝐾𝑖 =4 · 12·𝐸·𝐼𝑖 𝐿3 Ecuación 34. Rigidez de la planta i-ésima 3. Sabiendo que el momento de inercia de área de la sección del pilar, un cuadrado, es: 𝐼𝑖= 1 12· 𝑑4 Ecuación 35. Momento de inercia de un pilar de sección cuadrada 4. Así, suponiendo un módulo de elasticidad del acero, E = 210e9 Pa, sólo falta por saber el ancho del pilar de cada planta para poder obtener las matrices de rigidez y de masa. De esta manera, será este último parámetro, del que depende la resolución del problema, la variable con la que calibraremos el modelo (4 parámetros, 1 por cada planta). Por consiguiente, para obtener unos datos experimentales ficticios, supondré, al igual que en el análisis de sensibilidad, que los parámetros se distribuyen mediante una función normal y generaré una población de datos que sigan un muestreo, también, de hipercubo latino. Sabiendo que todo el problema depende de los cuatro parámetros correspondientes a los anchos de los pilares de cada planta, obtendré tantas matrices de masa y rigidez como simulaciones haga. Una vez obtenga la matriz de masa (M) y rigidez (K) para cada modelo distinto, podré calcular de manera analítica sus: - Frecuencias naturales: 𝑤= √𝑲 𝑴 [𝑟𝑎𝑑 𝑠] Ecuación 36. Expresión de las frecuencias naturales de un sistema 𝑓= 𝑤 2 · 𝜋 [𝐻𝑧] Ecuación 37. Frecuencias naturales de un sistema en hercios - Modos de vibración: (−𝑴𝑤𝑖2+𝑲)· 𝒙𝟎𝒊=𝟎 Ecuación 38. Ecuación para determinar los autovectores del sistema - Modos de vibración asociados a la frecuencia wi Así, para los valores experimentales que mediríamos in-situ pero que, al tratarse de un trabajo académico, podemos simular su obtención a partir de una simulación, que repetiremos 100 veces, con valores de entrada aleatorios que se muevan en el entorno de los exactos, obtendremos las frecuencias naturales y modos de vibración del edificio. De esta manera, obtendremos 100 muestras del ancho de cada pilar a partir de un muestreo, ya conocido, de hipercubo latino asumiendo que las cuatro variables se distribuyen mediante una función normal. Llevando a cabo el cálculo de las propiedades dinámicas mediante la formulación anteriormente definida, optaré por sacar la media y desviación típica de todos los valores hallados para las cuatro frecuencias y modos de vibración verticales, generando las medidas experimentales a partir de la media y desviación típica de cada propiedad dinámica y multiplicando por un valor aleatorio obtenido a partir de una distribución normal estándar: 𝑉𝑎𝑙𝑜𝑟 𝑒𝑥𝑝𝑒𝑟𝑖𝑚𝑒𝑛𝑡𝑎𝑙𝑖= 𝑚𝑒𝑑𝑖𝑎𝑖+ 𝑑𝑒𝑠𝑣𝑖𝑎𝑐𝑖ó𝑛 𝑡í𝑝𝑖𝑐𝑎𝑖· 𝑟𝑎𝑛𝑑𝑜𝑚 𝑛𝑢𝑚𝑏𝑒𝑟∈[0,1] Ecuación 39. Ecuación para la obtención de las medidas experimentales para la propiedad dinámica i – Edificio a cortante Así, obtenemos los valores experimentales correspondientes a los desplazamientos normalizados de los primeros cuatro modos de vibración y sus frecuencias asociadas. 72 Ilustración 46. Edificio a cortante - Valores experimentales de las cuatro primeras frecuencias Ilustración 47. Edificio a cortante - Histogramas de las cuatro primeras frecuencias Por consiguiente, quedará evaluar los modelos numéricos hasta obtener una convergencia de similitud con los datos experimentales obtenidos. Para ello, resolveremos el problema de solución analítica en el que obtendremos las frecuencias naturales de la estructura a partir de las ecuaciones Ecuación 36. Expresión de las frecuencias naturales de un sistema y Ecuación 37. Frecuencias naturales de un sistema en hercios. Mientras que los modos de vibración asociados los obtendremos hallando los autovectores del sistema de la Ecuación 38. Ecuación para determinar los autovectores del sistema - Modos de vibración asociados a la frecuencia wi. El desarrollo optimizador resolverá el problema analítico tantas veces como se indique al método. 5.9.1.1 Desarrollo del programa en Matlab En este apartado se redactarán los aspectos comunes a los tres métodos estadísticos comentados con el objetivo de clarificar el procedimiento, así como las diferencias entre ellos. Los programas tendrán el siguiente proceso de desarrollo: 1. Definir los valores inferiores y superiores de la función uniforme que conformará la función prior. Esto se llevará a cabo manualmente en base a la siguiente expresión: 𝑑𝑖= 𝑑𝑖,𝑚𝑒𝑑𝑖𝑜± 0,1· 𝑑𝑖,𝑚𝑒𝑑𝑖𝑜 5 Ecuación 40. Obtención del valor superior o inferior de los parámetros a actualizar Así, podemos definir una función de densidad de probabilidad uniforme y continua para cada ancho de los pilares de cada planta. Quedando, la PDF de la prior conjunta, como el producto de las cuatro anteriores. Generaremos un número aleatorio de la PDF prior a partir de esos límites superiores e inferiores de los valores, que será, ese primer valor, la muestra inicial aleatoria que dará comienzo a la cadena Markov. 2. Se definirá la función de probabilidad en su versión logarítmica que, al igual que el paso anterior y las definiciones posteriores, se harán de manera anónima para indicar que los valores de entrada se darán una vez se lleve a cabo el algoritmo correspondiente. Esta ejecutará las siguientes fases del cálculo: o Resolución numérica del problema de aplicación o Cargado de los resultados experimentales previamente obtenidos o Realización de un ajuste de mínimos cuadrados entre ambas familias de resultados en términos de modos de vibración y frecuencias asociadas a ellos, ponderando en 0.1 los primeros y en 0.9 en los segundos por la calidad de la información obtenida de cada uno. 3. De esta manera, la función posterior se define como la suma entre el logaritmo de la función PDF prior y la función logarítmica de máxima verosimilitud gracias a esta herramienta matemática. 4. Entramos en el apartado de configuración del muestreador Metropolis-Hastings y definimos: o La matriz de covarianza como la matriz identidad multiplicada por 𝜎𝑃, parámetro que veremos posteriormente. o La función PDF propuesta, que será una función de densidad de probabilidad normal multivariable cuya media será la muestra actual y la matriz de identidad de covarianza comentada anteriormente será la desviación, función que se evaluará para el valor de la muestra candidata. o Generación de unos números aleatorios a partir de una distribución normal cuya media será la muestra actual y desviación la matriz de covarianza. 5. Finalmente, el paso restante consiste en atribuir los valores necesarios a los algoritmos en materia de número de muestras y burn-in . Por lo que ya podemos analizar, individualmente, cada uno de ellos. 5.9.1.2 Método MCMC Además del procedimiento recientemente visto, la característica destacable de este método es el parámetro 𝜎𝑃, que consiste en la desviación estándar de la distribución normal propuesta, denominado como parámetro de tuning del algoritmo MH, sobre el cual profundizamos un poco sobre su significado en el apartado 4.2.2. La elección de este parámetro, como ya hemos comentado, es compleja y crucial. Tanto es así, que se ha llevado a cabo un programa que optimiza el valor de éste, minimizando la diferencia de la aceptación obtenida con el valor ideal que presenta la literatura científica de 0,234. Dicha optimización se ha llevado a cabo mediante un método de enjambre de partículas, cuyo método de funcionamiento se detalló en el apartado 2.2.3.1. Sin 80 Obteniendo los siguientes coeficientes de variación e histogramas más completos de entre los tres métodos, apreciando la mejor distribución normal hasta ahora: Ilustración 62. Edificio a cortante – SMC - COV Ilustración 63. Edificio a cortante – SMC - Histogramas Hallando los siguientes valores medios de los anchos de los pilares para cada planta: Ilustración 64. Edificio a cortante – SMC – Medias de las iteraciones de los anchos de los pilares para cada planta Para un tiempo total de computación bajo de 7 segundos, parecido al MCMC y bastante reducido respecto al TMCMC: Ilustración 65. Edificio a cortante – SMC – Tiempo invertido 5.9.1.5 Comparación de métodos De acuerdo con los apartados anteriores, en este se llevará a cabo la elección del método que será de aplicación en el ejemplo real de la pasarela ciclista y ámbito de aplicación de este trabajo. Para ello, se ha esquematizado, en una tabla, los resultados de cada método: Ilustración 66. Edificio a cortante - Comparativa métodos Como podemos apreciar, el método TMCMC es el que logra una mejor solución a efectos globales, llevándose la mejor o segunda mejor puntuación en cada uno de los aspectos. Sin embargo y como destaca por encima de todo, este mismo método tarda unas 4-5 veces lo que tardan los demás debido al mayor número de evaluaciones que realiza el algoritmo debido a la necesidad de generar muestras no sólo de la función posterior sino de las distribuciones transicionales. Por otra parte, el SMC carece de consideración del burn-in debido al proceso inicial del algoritmo que obtiene una serie inicial de muestras directamente de la distribución prior y, por lo tanto, asegurando que las muestras obtenidas al final se encuentran en la posterior (Lye et al., 2021). Cabe a destacar que los tres métodos logran una estimación de las variables a actualizar muy buena como hemos podido apreciar mediante los resultados estadísticos y medidas de aceptación. Sin embargo, de elegir uno para aplicarlo al modelo numérico complejo, se ha escogido el método TMCMC porque, a pesar de su considerable mayor coste computacional, la exactitud de sus resultados es innegable y presenta una gran ventaja respecto a la competencia y es la no necesidad de definición del parámetro de tuning. Parámetro que, una vez determinado, haría mucho más eficiente el proceso de calibración en términos de tiempo invertido por el algoritmo, pero que su búsqueda y la incertidumbre que gira en torno a él, reducen mucho su idoneidad. 5.9.2 Método bayesiano – Pasarela ciclista Se ha abarcado toda la toma de decisiones necesaria para abordar el problema protagonista de este documento, tanto la elección del método a ejecutar como las variables y el modelo de elementos finitos a actualizar. Habiendo hecho un repaso sobre las tres técnicas de muestreo que eran ámbito de estudio de este documento y habiendo determinado el óptimo para la calibración que se acontece, tenemos el trasfondo teórico y práctico necesario para abarcar el objetivo final de este trabajo: la actualización de la pasarela de Olivares que, como ha sido previamente detallada, en este apartado analizaremos únicamente la aplicación del método. 5.9.2.1 Medidas experimentales Como determinamos en el apartado del coeficiente de Spearman, los parámetros a actualizar serán los módulos 82 de deformación del: - Acero que compone los perfiles: 𝐸𝑎𝑐𝑒𝑟𝑜 - Tablero en dirección longitudinal de la zona lateral: 𝐸𝑡𝑎𝑏𝑙𝑒𝑟𝑜,𝑙𝑎𝑡𝑒𝑟𝑎𝑙,𝑙𝑜𝑛𝑔 - Tablero en dirección transversal de la zona central: 𝐸𝑡𝑎𝑏𝑙𝑒𝑟𝑜,𝑐𝑒𝑛𝑡𝑟𝑎𝑙,𝑡𝑟𝑎𝑛𝑠𝑣 - Tablero en dirección transversal de la zona lateral: 𝐸𝑡𝑎𝑏𝑙𝑒𝑟𝑜,𝑙𝑎𝑡𝑒𝑟𝑎𝑙,𝑡𝑟𝑎𝑛𝑠𝑣 A su vez, como se realizó en ese mismo apartado, se comenzará con la generación artificial de unas mediciones experimentales a partir de las medias asumidas de los valores de estos parámetros, con las que, atribuyéndoles una variación, asignaremos los valores de entrada al modelo de elementos finitos. De la misma manera que se llevó a cabo, se ejecutan tantas combinaciones de inputs como valores experimentales queramos, llamando al programa Ansys para que realice el análisis modal con esos parámetros, obteniendo las frecuencias experimentales y modos de vibración. Estas propiedades dinámicas serán comparadas con unos resultados previamente obtenidos manualmente que corresponden con los cuatro primeros modos de vibración verticales de la estructura. De esta manera, podemos, mediante el MAC, obtener, para cada nueva configuración del modelo numérico, los modos que correspondan con los verticales. Proceso que implicará, como comentamos en el 5.7.1.2, comparar los 150 modos de vibración con los modos previamente introducidos manualmente, para cada una de las 100 configuraciones del modelo computacional que se han creido adecuadas como serie de datos experimentales. Así, obtenemos tales mediciones que nos presentan los siguientes histogramas: Ilustración 67. Pasarela - Histogramas - Mediciones experimentales de las propiedades dinámicas 5.9.2.2 Resultados numéricos Se obtendrán de la misma manera que se obtuvieron en el ejemplo del edificio a cortante. Se definirán unos valores límite superiores e inferiores de una función de densidad de probabilidad uniforme y continua para cada parámetro a actualizar a partir de la media de cada uno y una pequeña variación. Posteriormente se definirá la función prior como producto de estas cuatro funciones de densidad así como la función logarítmica posterior mediante la suma, en sus versiones logarítmicas, de esta función prior y la función de máxima verosimilitud. La función de máxima verosimilitud variará, en este caso de aplicación, al tratarse de un problema cuya solución depende del uso de un modelo numérico. Para cada iteración, el programa llamará al software Ansys, que realizará el análisis modal de la pasarela ciclista con la configuración de valores de los parámetros que apliquen en esa iteración. El requisito de esta interacción de programas será el culpable del drástico aumento del coste computacional, pero será la forma de obtener unas frecuencias y modos de vibración numéricos a partir del cual, de los últimos, obtendremos los indicadores MAC asociados a los modos verticales con previa comparación con unos modos verticales introducidos manualmente. 5.9.2.3 TMCMC aplicado a la pasarela De esta manera, cargando los datos experimentales y estos últimos resultados numéricos, realizamos un ajuste de mínimos cuadrados entre ambas series de datos como realizamos en el ejemplo del edificio a cortante. Una vez definida la función logarítmica de la posterior, entramos en el apartado de configuración de la técnica de muestreo TMCMC. Apartado que, como hemos visto en el ejemplo simplificado anterior, basta con definir una serie de iteraciones y un burn-in que, como proceso iterativo, estableceremos cero y, una vez visualizando el proceso del algoritmo, hallaremos la solución final. Configurando el sampleador para 1000 iteraciones y 0 burn-ins realizará el flujo de operaciones ya definidas, obteniendo: Ilustración 68. Pasarela - 1000 iteraciones; 0 burn-ins Pudiendo apreciar un adecuado descarte de los primeros 100 resultados, para 1000 iteraciones y 100 burn-ins hallamos: 84 Ilustración 69. Pasarela – 1000 iteraciones; 100 burn-ins Con dicha traza obtenemos los siguientes coeficientes de variación: Ilustración 70. Pasarela – COV Que son algo superiores a los del ejemplo simplificado del edificio a cortante. Sin embargo, si analizamos las medias obtenidas: Ilustración 71. Pasarela – Medias de los parámetros estimadas Podemos apreciar que el programa es más que capaz de replicar la media de la que partían los datos experimentales. No obstante, esto no dejan de ser valores deterministas, algo que no concierne a estos métodos de naturaleza estadística. Observando los histogramas: Ilustración 72. Pasarela – Histogramas de los valores de los parámetros actualizados Observamos que, si bien es cierto, la asunción de la distribución de los parámetros mediante una normal requeriría una mayor semejanza de los histogramas con la famosa campana de este tipo de funciones, la desviación estándar/COV que caracteriza a estas variables y, especialmente, el relativamente reducido número de muestras aportadas al algoritmo provoca una apariencia menos ideal de lo esperado. 86 6 CONCLUSIONES N este trabajo se ha llevado a cabo un problema de actualización de un modelo de elementos finitos que simula el comportamiento de una pasarela ciclista. Documento en el que, tras un recorrido por el estado del arte de los distintos métodos existentes, sus trasfondos teóricos y distintas aplicaciones, realiza una calibración de dicho modelo mediante un método determinista y bayesiano, con el fin de estudiar e indagar en estos últimos al hacer un análisis previo de varias técnicas de muestreo en un ejemplo simplificado de resolución analítica. La principal problemática que se presenta al plantearse el desarrollo del modelo reside en las características que presenta el tablero de la pasarela, formada por un material compuesto y ortótropo. Donde la primera cuestión se ha tratado mediante un problema de homogeneización, que conlleva la modelización de un único material con propiedades geométricas mecánicamente equivalentes a la composición real, mientras que la segunda ha sido solucionada a través del mismo principio físico que consiste en el mantenimiento de las propiedades mecánicas del elemento. Se lleva a cabo, a nivel de modelado, una variación de la geometría modelada respecto a la real para solucionar el problema de varios materiales que componen el mismo elemento para equivaler las propiedades reales y modeladas. No obstante, la ortotropía implica que estas propiedades varían en función de la dirección de estudio, por lo que fue necesario un ajuste de ciertos parámetros para mantener esas propiedades equivalentes en todas las direcciones del espacio. Una vez analizado estas características, se estuvo en posición de llevar a cabo el modelo de elementos finitos, cuyo análisis modal fue realizado con éxito y comprobado con un modelo analítico simplificado respecto al que experimentaba ciertas variaciones pero que se asumieron como despreciables dadas las diferencias claras de geometría, tipología estructural y cuestiones propias de un modelo extenso. Con el modelo de elementos finitos correctamente definido, la siguiente cuestión pasó a ser la elección de los parámetros que podían ser susceptibles de actualizar entre todos aquellos que gobernaban el comportamiento de la estructura. En un primer acercamiento, y para poner en práctica los programas y técnicas que posteriormente se explicarían, se trató de estimar estas variables basando, esta decisión, en el juicio ingenieril. Sin embargo, tras una serie de pruebas y no siendo, los resultados, plenamente satisfactorios, inexorablemente se estipuló como necesaria la correcta justificación de esta elección con el objeto de lograr una respuesta más satisfactoria. Para realizar esto se llevó a cabo un análisis de sensibilidad que, tras varias interpretaciones, se definió como definitivo, por la precisión del análisis que conllevaba, lo obtenido mediante el estudio del coeficiente de Spearman, materia en la que la aplicación también resultaba novedosa para mi. Habiendo estudiado las propiedades mecánicas de la estructura y adaptadas para mantener constante la respuesta de ésta, así como habiéndola modelado y llevado a cabo un análisis de sensibilidad para obtener los parámetros ideales para la calibración del modelo de elementos finitos, sólo queda la aplicación de los métodos de actualización que son el propio objeto de estudio de este documento. Los métodos de actualización a emplear han sido de dos naturalezas: deterministas y estocásticas. En la primera, se llevó a cabo una optimización del problema cuyo objetivo residía en la minimización de la diferencia de los valores de las frecuencias y modos de vibración de la estructura numéricos y aquellos tomados como experimentales gracias al uso de un algoritmo genético que generaba poblaciones de individuos encargados de realizar la evaluación de la función y convergiendo, cada iteración, en soluciones más precisas. Mientras que el E segundo caso resultó más complejo en lo que respecta tanto a la teoría como a la práctica, por lo que resultó conveniente recurrir a ejemplos más abarcables para ese primer acercamiento. De esta manera, para el método bayesiano se recurrió, como introducción, a un ejemplo más sencillo en términos de computación: el modelo de edificio a cortante. Gracias a su solución analítica, su resolución podía llevarse a cabo íntegramente en el programa de cálculo numérico Matlab, no precisando la interacción entre dos programas, la cual ha sido la principal fuente de aumento del coste computacional. Esto permitió un análisis de los tres métodos bayesianos estudiados en el documento, lo que serviría como preámbulo al problema de aplicación inicial del trabajo: la pasarela ciclista, con el que se obtuvieron las siguientes conclusiones: - El método MCMC resulta ser el más rápido en ejecutarse, seguido por el SMC y dejando, en última posición, el TMCMC, que tarda considerablemente más que éstos primeros. Esto último se debe al mayor número de evaluaciones que realiza debido a la generación de una mayor cantidad de muestras basadas en su algoritmo. - A costa de esa primera conclusión, se pudo observar que los mejores valores del coeficiente de variación se hallaron con el método TMCMC. No alejándose mucho de ellos el resto de métodos. - De esta misma manera, las medias obtenidas de los valores de los parámetros a actualizar no difieren de un método a otro, otorgando buenos resultados los 3 métodos analizados. - Uno de los métodos de comprobación de la correcta implantación de estas técnicas estadísticas de muestreo consiste en la asimilación de los resultados obtenidos a una distribución normal como la que se ha asumido que se distribuyen los parámetros a actualizar. En esta materia, el MCMC otorga unos resultados más pobres en favor a los restantes. Sin embargo, aunque se podría esperar un histograma, por parte del TMCMC, mucho más acorde a lo esperado, por esa numerosa cantidad de evaluaciones que hace de más en comparación al resto de técnicas, acabó graficando uno a la par del método SMC, el cual era conseguido con un tiempo de ejecución 4 veces inferior. Si bien, en un primer acercamiento, pudiera parecer que el TMCMC y el SMC generaron unos mejores resultados y que, en base al notablemente menor tiempo de ejecución del segundo, fuera éste el elegido como el ganador, existió un gran factor condicionante que limitaba la idoneidad tanto de ésta técnica como de la principal, el MCMC, éste es el parámetro de tuning. A lo largo del trabajo se menciona la dificultad de su definición por lo que, en base a la exactitud de los resultados obtenidos con el método TMCMC y su capacidad de automatización, se decantó por su elección para el problema real. Pudiendo concluir con que la idoneidad de los métodos MCMC y SMC reside en el conocimiento de este parámetro, diferenciando entre ambos que el primero es algo más rápido de computar, aunque algo menos preciso mientras que, en casos en los que no se disponga de este valor, la técnica TMCMC gana mucha ventaja. Así, el sampleador transicional elegido, junto al algoritmo de aceptación o rechazo Metropolis-Hastings, llevaron a cabo un total de, para 1000 muestras, más de 17000 análisis modales de la pasarela, obteniendo unos histogramas que, a pesar de poder presentar mejor geometría, aciertan con suficiente proximidad las variables a estimar. Concluyendo una exitosa aplicación de esta técnica de muestreo sobre el modelo de elementos finitos al obtener un, ahora sí, gemelo digital de la pasarela ciclista construida entre ambos municipios que, de disponer una serie de aparatos de medida en ella y realizar una medición continuada en el tiempo, podría analizarse, a “tiempo real” el estado de la estructura, logrando su monitorización constante para trabajos de mantenimiento, conservación, estudios de viabilidad,... entre un sinfín de posibilidades al poseer un control absoluto sobre la respuesta y estado de la infraestructura. 88 7 FUTURAS LÍNEAS DE INVESTIGACIÓN El desarrollo de este problema se ha realizado de la manera más veraz, dentro de ciertas limitaciones, que se ha podido plantear. Un ejemplo de esto ha sido el modelo de elementos finitos en los que se ha querido definir con precisión sus características, algo que ha llevado a un modelo muy extenso y que ha sido el motivo del inconveniente protagonista de este trabajo, el coste computacional. El hecho de tener un modelo de elementos finitos de este tamaño implica una mayor tardanza en el cálculo del análisis modal, pero, además de eso, también, en este caso, ha implicado una gran presencia de modos de vibración locales que han colocado a los globales en números tan altos como el 150. Una computación de un modelo extenso para un número tan alto de modos de vibración ha derivado en que el último ejercicio y objetivo de este documento tenga un coste de computación del orden de numerosas horas pues ha llegado a implicar un intervalo de tiempo de 75 segundos cada uno de los análisis modales que han sido evaluados. De esta manera, como futuras líneas de investigación, se podría estudiar la paralelización de los algoritmos pues el cálculo numérico de Matlab lo permite y el programa de elementos finitos Ansys lo concibe. Sin embargo, requiere un estudio previo, especialmente en este segundo al presentar menos información y es posible que una licencia educativa no presente ese servicio. Una modelización más fidedigna de la estructura, en lo que respecta a la geometría, fue perseguida con el fin de desarrollar un acercamiento más visual e intuitivo a todos los niveles. Por lo que, por otro lado y como consideraciones a esos posibles trabajos, algunas simplificaciones que podrían ser materia de estudio para ese ahorro computacional, en trabajos con otros alcances, podrían ser la de análisis de un modelo simplificado equivalente mecánicamente. Sin olvidar que la existencia de tantos modos de vibración locales fue debido a esa extensión del modelo que derivó en un análisis modal más costoso. De esta manera, la consideración de las tres primeras frecuencias asociadas a los tres primeros modos de vibración verticales puede que cumplan con los requisitos de ese posible futuro documento en cuanto a actualización del modelo respecta, siendo consciente de que, en el caso de este trabajo, reduce el tiempo de computación de cada análisis modal de 75 segundos a sólo 5. Sin embargo, esto es, de nuevo, un factor a considerar en modelos extensos, en caso de haber optado por la vía de simplificación del modelo, seguramente no resulte necesaria la no consideración de un modo de vibración por completo que no deja de ser la fuente principal de información, de entrada, de estas técnicas estadísticas de muestreo. Por último, queda por mencionar la inolvidable importancia que tienen los resultados hallados gracias a los métodos probabilísticos de actualización de modelos, pues son éstos quienes sientan las bases de lo conocido como análisis de fiabilidad: una técnica que permite, una vez considerada la probabilidad de fallo exigida a la estructura, no sobredimensionar la misma al no requerir la utilización de coeficientes de seguridad en su diseño. Contemplando, como buena continuación a este trabajo, el estudio de fiabilidad de una estructura partiendo de unos resultados como los hallados al final de éste. 8 REFERENCIAS Ansys® Academic Research Mechanical, Release 18.1 Bado, M. F., Tonelli, D., Poli, F., Zonta, D., & Casas, J. R. (2022). Digital Twin for Civil Engineering Systems: An Exploratory Review for Distributed Sensing Updating. In Sensors (Vol. 22, Issue 9). MDPI. https://doi.org/10.3390/s22093168 Barišić, L. (2020). A Digital Twin of an Asphalt Road. www.byg.dtu.dk Brocca, L., Barbetta, S., Camici, S., Ciabatta, L., Dari, J., Filippucci, P., Massari, C., Modanesi, S., Tarpanelli, A., Bonaccorsi, B., Mosaffa, H., Wagner, W., Vreugdenhil, M., Quast, R., Alfieri, L., Gabellani, S., Avanzi, F., Rains, D., Miralles, D. G., … Fernandez, D. (2024). A Digital Twin of the terrestrial water cycle: a glimpse into the future through high-resolution Earth observations. Frontiers in Science, 1. https://doi.org/10.3389/fsci.2023.1190191 Ching, J., & Chen, Y.-C. (2007). Transitional Markov Chain Monte Carlo Method for Bayesian Model Updating, Model Class Selection, and Model Averaging. Journal of Engineering Mechanics, 133(7), 816–832. https://doi.org/10.1061/(ASCE)0733-9399(2007)133:7(816) Ereiz, S., Duvnjak, I., & Fernando Jiménez-Alonso, J. (2022). Review of finite element model updating methods for structural applications. Structures, 41, 684–723. https://doi.org/10.1016/J.ISTRUC.2022.05.041 Friederich, J., Francis, D. P., Lazarova-Molnar, S., & Mohamed, N. (2022). A framework for data-driven digital twins of smart manufacturing systems. Computers in Industry, 136, 103586. https://doi.org/10.1016/j.compind.2021.103586 Gómez-Coronel, L., Santos-Ruiz, I., Torres, L., López-Estrada, F. R., Gómez-Peñate, S., & Escobar- Gómez, E. (2023). Digital Twin of a Hydraulic System with Leak Diagnosis Applications. Processes, 11(10). https://doi.org/10.3390/pr11103009 Grieves, M. (2016). Origins of the Digital Twin Concept. https://doi.org/10.13140 Javier, F., Somet, P., José, J., Fuentes, M., Romero Ordóñez, A., Barrera, P. G., & Solís Muñiz, M. (2020). Apuntes para el Análisis de Estructuras Universidad de Sevilla. Liu, C., Zhang, P., & Xu, X. (2023). Literature review of digital twin technologies for civil infrastructure. Journal of Infrastructure Intelligence and Resilience, 2(3), 100050. https://doi.org/10.1016/j.iintel.2023.100050 Lye, A., Cicirello, A., & Patelli, E. (2021). Sampling methods for solving Bayesian model updating problems: A tutorial. Mechanical Systems and Signal Processing, 159. https://doi.org/10.1016/j.ymssp.2021.107760 Marwala, T., Boulkaibet, I., & Adhikari, S. (2017). PROBABILISTIC FINITE ELEMENT MODEL UPDATING USING BAYESIAN STATISTICS. MATLAB. (2010). version 9.14.0.2286388 (R2023a) Update 3. Natick, Massachusetts: The MathWorks In