Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Autor: Marta Montero Fernández-Cavada Tutor: Luis Rodríguez de Tembleque Solano Dpto. de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 iii Trabajo de Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Autor: Marta Montero Fernández-Cavada Tutor: Luis Rodríguez de Tembleque Solano Profesor Titular de Universidad Dpto. de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 v Trabajo de Fin de Grado: Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Autor: Marta Montero Fernández-Cavada Tutor: Luis Rodríguez de Tembleque Solano El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2023 El Secretario del Tribunal vii Agradecimientos A mis padres, porque ellos siempre creyeron en mí. A mi tutor, Luis Rodríguez de Tembleque, por su dedicación y su paciencia con este trabajo. Marta Montero Fernández-Cavada 2023 ix Índice Agradecimientos vii Índice ix Índice de Tablas xi Índice de Figuras xiii Notación xv 1 Introducción 3 1.1 Introducción 3 1.2 Motivación y Objetivos 7 1.3 Organización del documento 7 2 Modelización teórica de las estructuras tipo honeycomb 13 2.1 Leyes de comportamiento 13 2.1.1 Materiales monoclínicos 16 2.1.2 Materiales ortótropos 16 2.1.3 Materiales isótropos 17 2.2 Modelo de estructura honeycomb 18 2.3 Formulación teórica de las propiedades elásticas efectivas 20 2.3.1 Formulación teórica Gibson 1982 22 2.3.2 Formulación teórica Gibson y Malek 2015 23 2.3.3 Formulación teórica Sorohan 2015 25 2.3.4 Formulación teórica para la determinación de las propiedades elásticas efectivas 26 3 Análisis por volumen representativo basado en el MEF 29 3.1 Método del elemento de volumen representativo (RVE) 29 3.1.1 Concepto de volumen representativo 29 3.1.2 Análisis multiescala 31 3.1.3 Condiciones de contorno y obtención de las propiedades efectivas 33 3.1.4 Implementación numérica 34 3.2 Implementación en ANSYS 35 3.2.1 Modelización en ANSYS 36 3.2.2 Elección del mallado 37 3.2.3 Obtención de las propiedades efectivas 39 3.2.4 Obtención de las constantes ingenieriles 41 4 Resultados 43 4.1 Comparativa entre los modelos teóricos y el MEF (RVE) 43 σij σi 𝜎𝑖  Tensor de tensiones Pseudovector de tensiones Pseudovector de tensión media𝐶𝑖𝑗     𝜎𝑖𝑝 𝜎𝑖𝑒 Pseudovector de tensión periódico 𝐶𝑖𝑗     Pseudovector de tensión de cada elemento 1 Índice de Figuras 2 3 3 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 1 INTRODUCCIÓN 1.1 Introducción lo largo de la historia, la Ingeniería, al igual que otras ramas del conocimiento como la Medicina, ha aportado numerosas soluciones inspiradas en la propia naturaleza. En general, el ser humano ha buscado estudiar y reproducir los comportamientos de la naturaleza para mejorar su calidad de vida. Es curioso advertir que, todos los sistemas de alta ingeniería que hoy conocemos, la naturaleza los hizo primero. Desde que el hombre es hombre, este empezó a mostrar una gran curiosidad por las enrevesadas estructuras que constituían los panales de las abejas, observando las ventajas naturales de estas estructuras hexagonales, así como su eficiencia para el almacenaje y su sorprendente resistencia. No es coincidencia la tendencia de importantes civilizaciones, como es la romana o la griega, de usar esta geometría es su arquitectura para ahorrar espacio y material. Los paneles hexagonales periódicos, o más comúnmente conocidos como estructuras honeycombs -de panal de abeja-, consisten en una disposición uniforme de celdas hexagonales conectadas entre sí, formando una estructura que reproduce a la de los panales de las abejas como se ve en la Figura 1-1. Esta estructura de red tridimensional es particularmente eficaz a la hora de distribuir de forma uniforme la carga haciendo que esta sea fácil de absorber sin sufrir deformaciones significativas y evitando la ruptura de la estructura [1]. En lo que sigue se hará referencia a estas estructuras como estructuras honeycombs. Figura 1-1. Red de paneles hexagonales conectados entre sí. Su configuración geométrica supone que estas celdas tengan una muy buena relación entre su resistencia y su peso, así como unas propiedades mecánicas muy interesantes, tales como alta resistencia a la compresión y a la flexión, principalmente en la dirección perpendicular al plano de la estructura, además de una rigidez muy favorable también en relación con la ligereza [2], [3]. En la A Introducción 4 dirección paralela al plano de la estructura la carga se distribuye mediante las paredes de las celdas lo que supone una gran resistencia lateral. Las estructuras honeycombs han sido objeto de numerosos estudios con el fin de profundizar en estas propiedades [2-5]. Gracias a los mismos, se llega a la conclusión de que las estructuras hexagonales aportan ligereza y mejoran las propiedades mecánicas de los materiales en los que se implementan [6]. Es por esto por lo que la configuración más empleada de las estructuras honeycomb es en paneles sándwich como parte del core. Estos paneles compuestos se componen de dos placas frontales separadas por un material central para aumentar la rigidez a la flexión y proporcionar una estructura ligera y rígida (ver Figura 1-2). Figura 1-2. Esquema de un panel sándwich con core o núcleo de panel de abeja. Figura extraída de [7]. El núcleo de panel de abeja constituye el core ideal para multitud de aplicaciones [8], [9], [7]. No es hasta el siglo pasado que se empiezan a utilizar estos paneles en el ámbito de la ingeniería aeroespacial, de la mano de una creciente necesidad de materiales con unas propiedades capaces de satisfacer requisitos cada vez más exigentes [10-13]. Y es desde este momento que se implementan estas estructuras hexagonales en materiales destinados a aeronaves, con el objetivo de aumentar su resistencia y conseguir estructuras más ligeras sin perder rigidez. El material más empleado en diferentes elementos de aeronaves y aviones se trata del Nomex [14-17]. Ejemplos de esto son el avión comercial Boeing 747 (primer vuelo el 9 de febrero de 1969) y el ATR 72 (primer vuelo el 27 de octubre de 1988), que están, en mayor parte, formados por paneles sándwich con núcleo de estructura honeycomb Nomex [10] como puede observarse en la Figura 1-3.a que muestra las zonas concretas de este avión comercial donde se emplea este material. El creciente interés en estos paneles hexagonales supone que estos empiecen a aparecer cada vez más en la ingeniería moderna abarcando ámbitos tan diversos como son el campo de la medicina [12] o el del embalaje y empacado [19], [20] (ver Figura 1-3.b). Las estructuras honeycombs son estructuras porosas ampliamente usadas como injertos óseos por conseguir la ligereza y la similitud con las propiedades propias de un verdadero hueso humano [18], [21], [22] como se observa en la Figura 1- 3.c donde se esquematiza como sería la utilización de estas estructuras en implantes óseos en animales. Además, son popularmente empleadas para utilización militar debido a su capacidad de absorción tanto acústica como de energía [23], su eficiencia mecánica y su resistencia para soportar impactos de 5 5 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos balas [24-26]. (a) (b) (c) Figura 1-3. Ilustraciones de las aplicaciones de las estructuras honeycombs (a) Materiales compuestos presentes en el modelo de avión comercial ATR 71. Figura extraída de [10]. (b) Estructuras honeycomb de espuma para empacaje. Figura extraída de [27]. (c) Ejemplos de implantes óseos mediante estructuras honeycomb. Figura extraída de [22] Introducción 6 Las estructuras honeycombs también han ganado cierta importancia en el ámbito de la ingeniería de la energía. Pueden ser encontrados en paneles solares y recolectores de energía [28-31]. Numerosos estudios son realizados con paneles sándwich con núcleo de panel de abeja para reducir el consumo de energía en edificios como se explica en [32]. Cabe destacar que, debido a su buena relación entre resistencia y peso permite utilizar menos material consiguiendo estructuras que presentan buenas propiedades mecánicas. Esto supone una mejora de la eficiencia energética de los materiales en los que se implementa [6] siendo la problemática de la sostenibilidad una de las cuestiones que condicionan todos los proyectos en la ingeniería moderna. Finalmente, cabe destacar la presencia en los últimos años de estas estructuras en el campo de la impresión 3D. La impresión 3D se refiere a la impresión de objetos físicos mediante capas de material desde un ordenar [33]. Esta tecnología ha realizado grandes conquistas recientemente en cuanto a avances se refiere. Permite diseñar e imprimir geometrías más complejas que las conseguidas en los procesos tradicionales. Los rellenos de estas impresiones tienen diferentes geometrías. Una muy utilizada para el relleno es la geometría honeycomb [34], [35]. Este tipo de relleno y un ejemplo de impresora 3D que lleva a cabo estos procedimientos se puede ver en la Figura 1-4. Este es el tipo de patrón para el relleno más eficiente al usar una menor cantidad de material y proporcionar mayor resistencia que otros patrones como pueden ser el rectangular o el triangular. Se han utilizado los paneles hexagonales para técnicas de impresión 3D como es el FDM [36] y con diferentes materiales como los polímeros o las fibras de carbono [36], [37]. Las estructuras honeycombs se estudian en la actualidad para aplicar sus propiedades de la manera más eficiente. En [38] estudian la buena respuesta de estas estructuras a largas deformaciones de compresión y en [39] se estudia la facilidad de recuperación de la forma inicial de los paneles hexagonales. Figura 1-4. (a) Impresora 3D (b) Panel con relleno de estructura honeycomb. Figura extraída de [36]. Por todo lo anterior, es fácil comprender el interés que tienen estas estructuras en ámbitos tan diversos como son la industria aeroespacial, la militar o la medicina. 7 7 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 1.2 Motivación y Objetivos La gran diversidad de aplicaciones industriales (de carácter mecánico), viene motivando el desarrollo de numerosos trabajos que tienen por objeto modelizar el comportamiento mecánico de los sistemas estructurales honeycomb. Desde los pioneros trabajos de Gibson [40], [41], en los últimos años se vienen desarrollando nuevos modelos teóricos [42] y nuevos modelos numéricos basados en el Modelo de Elementos Finitos [7], que tratan de acercarse más al comportamiento real de dichas estructuras. En este contexto se presenta este trabajo, donde se presenta un marco computacional basado en las técnicas de Elementos de Volumen Representativo (RVE) con el MEF. La idea es desarrollar dicho esquema de modelización numérica para posteriormente comparar sus estimaciones con los principales modelos teóricos que encontramos en la literatura. Además, el modelo RVE permitirá estudiar la influencia que tienen diversos parámetros geométricos sobre la respuesta mecánica estas estructuras honeycomb. De esta manera, los objetivos concretos que persigue este proyecto son: 1. Se busca implementar un modelo de cálculo para obtener las constantes ingenieriles del material de las estructuras de panal de abeja o estructuras honeycombs así como los coeficientes efectivos mediante la aplicación de diferentes condiciones de contorno. Esto se llevará a cabo mediante la modelización de dicha estructura haciendo uso del software de elementos finitos ANSYS junto con el desarrollo de ficheros de cálculo en Matlab. 2. Es competencia también de este proyecto validar las expresiones analíticas utilizadas para dicho desarrollo analítico, así como lo es validar el modelo RVE implementado en ANSYS. 3. Una vez desarrollado el modelo de la estructura honeycomb y obtenidas sus propiedades efectivas se llevarán a cabo una serie de estudios sobre el modelo. En primer lugar, se estudiará el impacto en los resultados obtenidos de variar ciertos parámetros geométricos del modelo. Así mismo, se estudiará influencia del ángulo interno de las células hexagonales periódicas en las propiedades efectivas de la misma además de su efecto en las tensiones y desplazamientos una vez aplicadas diferentes condiciones de contorno al modelo de la estructura. 1.3 Organización del documento Este documento va a estar dividido en cinco capítulos además de los anexos donde se desarrollan los contenidos del presente proyecto. A continuación, se describen brevemente el contenido de cada uno de ellos. • Capítulo 2 – Modelización teórica de estructuras tipo honeycomb En este capítulo se describirán las leyes de comportamiento de la Teoría de la Elasticidad para clasificar los materiales compuestos haciendo hincapié en el caso ortótropo que será Introducción 8 objeto de este trabajo. Además, se presentará la geometría que se ha utilizado para modelar y estudiar las estructuras hexagonales. Una vez presentado el modelo de la estructura se exponen las propiedades efectivas que definen el problema y se presentan las tres formulaciones que se utilizarán para calcularlas. Finalmente, se definirán las expresiones que relacionan los coeficientes efectivos con las constantes ingenieriles. • Capítulo 3 – Análisis por volumen representativo basado en el MEF En este capítulo se desarrollarán las ideas principales del método del elemento de volumen representativo, en el que está basado el modelo realizado en este trabajo para la estructura honeycomb. Además, se describe el modelo realizado con el software de elementos finitos ANSYS y los pasos que se han seguido para modelizarlo. Se desarrollará la formulación y las condiciones de contorno necesarias para obtener todos los coeficientes efectivos. Finalmente, se desarrollarán las expresiones para relacionar las constantes ingenieriles a partir de los coeficientes efectivos. • Capítulo 4 – Resultados En este capítulo se realizará la validación de las expresiones analíticas que se utilizarán para calcular las propiedades efectivas de las estructuras honeycombs mencionadas en capítulos anteriores. Una vez validadas las expresiones se realizará un estudio de la influencia de la variación del ángulo interno de una célula unitaria de una estructura periódica hexagonal. Se mostrarán también las figuras más representativas de la variación de los desplazamientos y las tensiones en función del ángulo interno para cada condición de contorno. • Capítulo 5 – Conclusiones En este capítulo se realizará un resumen del trabajo realizado junto con una interpretación de los resultados obtenidos y las conclusiones más importantes a las que se han llegado. Igualmente, para acabar se propondrán las aplicaciones del propio proyecto para la realización de estudios futuros. • Anexos En los anexos se recogerán las rutinas de MATLAB utilizadas para el cálculo y modelización del problema. 11 Modelización teórica de las estructuras tipo honeycomb 18 18 De igual manera, el tensor de flexibilidad que se obtiene tiene la siguiente forma: [ 𝜀1 𝜀2 𝜀3 𝛾23 𝛾13 𝛾12 ] = [ 𝑆11 𝑆12 𝑆12 𝑆12 𝑆11 𝑆12 𝑆12 𝑆12 𝑆11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2(𝑆11−𝑆12) 0 0 0 2(𝑆11−𝑆12) 0 0 0 2(𝑆11−𝑆12) ] [ 𝜎1 𝜎2 𝜎3 𝜎23 𝜎13 𝜎12 ] (2.18) Por tanto, un material isótropo queda definido mediante 2 constantes elásticas independientes. 2.2 Modelo de estructura honeycomb Como se ha descrito previamente, el objetivo de este estudio es modelizar una célula unitaria hexagonal como se muestra en la Figura 2-1. En este apartado se van a definir los parámetros que parametrizan el modelo que se ha utilizado en este proyecto para realizar los diferentes estudios que se exponen en la presente memoria. Figura 2-1. Celda unitaria. Figura extraída de [42] El modelo de la celda unitaria quedaría definido con cuatro dimensiones geométricas de diseño como son θ, t, h y l. Se define el ángulo θ como el ángulo interno de la célula hexagonal que se trata del ángulo que forma la pared vertical con la horizontal. Cabe destacar que, debido a que la dimensión del espesor de la pared del núcleo es mucho más pequeña que las dimensiones h y l se ha planteado un modelo bidimensional como se muestra en la Figura 2-2, donde a1 en la semilongitud en la dirección principal x1, y a2 será la semilongitud en la dirección principal x2. La profundidad de la celda quedará definida por la semilongitud a3, a la cual, en lo que sigue, se hará referencia como distancia b (ver Figura 2-3. b) que es equivalente al doble de la semilongitud a3. 19 19 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Figura 2-2. Modelo de célula unitaria hexagonal con los parámetros geométricos que la definen. En la Figura 2-3 se muestra el modelo del núcleo de la celda unitaria para un modelo tridimensional con una profundidad b. Se va a estudiar, como se muestra en la Figura 2-2, un núcleo con un espesor t1 = t2 = t. (a) (b) Figura 2-3. (a) Parámetros que definen el modelo de célula hexagonal en un modelo tridimensional. Figura extraída de [45]. (b) Dimensión del ancho de la celda unitaria. Figura extraída de [7] Modelización teórica de las estructuras tipo honeycomb 20 20 2.3 Formulación teórica de las propiedades elásticas efectivas El conocimiento de las constantes ingenieriles de un material es indispensable para el diseño y cálculo de estructuras en el campo de la ingeniería. Estas proporcionan información sobre sobre las propiedades físicas y mecánicas de los materiales con los que se trabaja. Previamente a desarrollar la formulación teórica utilizada para calcular las constantes ingenieriles de las estructuras hexagonales periódicas se van a definir las constantes que serán objetivo de este trabajo. • Módulos de Young (E1, E2 y E3): se trata de las rigideces de un material en las direcciones principales 1, 2 y 3 respectivamente. Se define como el esfuerzo por unidad de deformación unitaria de un material cuando este está sometido a tracción en una de estas direcciones. • Módulo de cizalladura (Gij): también llamado módulo de elasticidad transversal es una constante elástica que define la respuesta de un material ante la acción de un esfuerzo cortante que lo deforma, en este caso, en el plano ij. • Coeficiente de Poisson (νij): define la relación entre las deformaciones longitudinales que sufre un material en dirección j y las deformaciones longitudinales en dirección i dada una tensión en la dirección i. • Densidad (ρ): se trata de una propiedad intrínseca del material que se define como la masa por unidad de volumen. • Módulo de rigidez volumétrica (K): en los materiales de tipo isótropo representa la presión hidrostática que es necesaria para obtener una deformada volumétrica unitaria como se expresa a continuación: 𝑝= −𝐾∆𝑉 𝑉 (2.19) Las constantes ingenieriles descritas anteriormente, se relacionan con la ley de comportamiento de forma directa mediante el tensor de flexibilidad. A pesar de ello, es posible también obtener dichas constantes a partir de las componentes del tensor de rigidez teniendo en cuenta su relación inversa con el tensor de flexibilidad (2.8). En la Tabla 2-1 se muestra un resumen de las expresiones para calcular el tensor de rigidez a partir de las constantes ingenieriles en función del tipo de material. 21 21 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Tabla 2-1. Resumen de las constantes ingenieriles en función del tipo de material Tipo de Material Ortótropo Isótropo Constante ingenieriles independientes E1, E2, E3, 𝜈23, 𝜈13, 𝜈12, G23, G13 y G12 K y G Constantes ingenieriles dependientes 𝜈32= 𝜈23𝐸3 𝐸2 𝜈31= 𝜈13𝐸3 𝐸1 𝜈21= 𝜈12𝐸2 𝐸1 𝐸= 9𝐺𝐾 3𝐾+𝐺 𝜈= 𝐸−2𝐺 2𝐺 𝐸1= 𝐸2= 𝐸3=𝐸 𝜈23= 𝜈13= 𝜈21= 𝜈 𝐺23= 𝐺13= 𝐺12=𝐺 Tensor de Flexibilidad [ 1 𝐸1−𝜈21 𝐸2−𝜈31 𝐸3 −𝜈12 𝐸11 𝐸2−𝜈32 𝐸3 −𝜈13 𝐸1−𝜈23 𝐸21 𝐸3 0 00 0 01 𝐸10 0 0 00 0 01 𝐸1 0 0 01 𝐸10 0 01 𝐸10 0 1 𝐺23 0 0 01 𝐺13 0 0 0 1 𝐺12 ] [ 1𝐸−𝜈 𝐸−𝜈 𝐸 −𝜈 𝐸1𝐸−𝜈 𝐸 −𝜈 𝐸−𝜈 𝐸1𝐸 0 00 0 01 𝐸10 0 0 00 0 01 𝐸1 0 0 01 𝐸10 0 01 𝐸10 0 1𝐺 0 0 01𝐺0 0 0 1𝐺 ] Obtención a partir del tensor de rigidez 𝐶11= 1− 𝜈23𝜈32 𝐸2 𝐸3∆ 𝐶22= 1− 𝜈13𝜈31 𝐸1 𝐸3∆ 𝐶33= 1− 𝜈12𝜈21 𝐸1 𝐸3∆ 𝐶23= 𝜈23− 𝜈21𝜈13 𝐸1 𝐸2∆ 𝐶13= 𝜈13− 𝜈12𝜈23 𝐸1 𝐸2∆ 𝐶12= 𝜈12− 𝜈32𝜈13 𝐸1 𝐸3∆ 𝐶44=𝐺23 𝐶55=𝐺13 𝐶66=𝐺12 𝐸= 𝐶11− 2𝐶12 2 𝐶11+ 𝐶12 𝐺= 𝐶44 Modelización teórica de las estructuras tipo honeycomb 22 22 Donde ∆ = (1 – 𝜈12 𝜈21 – 𝜈23 𝜈32 – 𝜈31 𝜈13 – 2 𝜈21 𝜈32 𝜈13) / (E1 E2 E3). Para obtener las constantes ingenieriles de las estructuras honeycombs se utilizarán tres formulaciones diferentes con la finalidad de validar y comparar los resultados obtenidos. Las distintas expresiones teóricas han sido extraídas de [40], [42] y [45]. 2.3.1 Formulación teórica Gibson 1982 A continuación, se desarrollará la formulación teórica expuesta en el artículo de Gibson [40] para obtener las constantes ingenieriles de una célula unitaria hexagonal. La densidad del material se obtiene con la siguiente expresión: 𝜌𝑐 𝜌 = 𝑡𝑙 ⁄(ℎ𝑙 ⁄ + 2) 2(ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑐𝑜𝑠𝜃 (2.20) Por otro lado, los módulos de Young en las diferentes direcciones según la formulación del artículo de Gibson se obtienen mediante las siguientes expresiones: 𝐸1𝑐 = 𝐸 (𝑡𝑙 ⁄)3𝑐𝑜𝑠𝜃 (ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑠𝑖𝑛2𝜃 (2.21) 𝐸2𝑐 = 𝐸 (𝑡𝑙 ⁄)3 (ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃) 𝑐𝑜𝑠3𝜃 (2.22) El módulo de Young en la dirección X3 se puede obtener asumiendo que para la carga normal en dicha dirección la deformación es axial. Por lo tanto, el módulo de Young longitudinal de la estructura honeycomb se obtiene como: 𝐸3𝑐 = 𝐸 𝜌𝑐 𝜌 (2.23) 23 23 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Los coeficientes de Poisson se obtienen mediante las siguientes expresiones: 𝜈12 𝑐 = 𝑐𝑜𝑠2𝜃 (ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑠𝑖𝑛 𝜃 ; 𝜈21 𝑐 = (ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑠𝑖𝑛 𝜃 𝑐𝑜𝑠2𝜃 (2.24) Una vez obtenidos los módulos de Young en las diferentes direcciones principales (E1, E2 y E3), los coeficientes de Poisson restantes pueden calcularse como: 𝜈13 𝑐 = 𝜈 𝐸1 𝐸3; 𝜈23 𝑐 =𝜈 𝐸2 𝐸3 (2.25) 2.3.2 Formulación teórica Gibson y Malek 2015 En el siguiente apartado, se desarrollará la formulación teórica extraída del artículo de Gibson y Malek [42] para obtener las constantes ingenieriles de una célula unitaria hexagonal. En esta formulación se introducen dos nuevos parámetros de corrección de los parámetros anteriormente definidos h y l que se obtienen según las siguientes expresiones: 𝑙𝑏 = 𝑙 − 𝑡 (2 𝑐𝑜𝑠𝜃) ; ℎ𝑏 = ℎ − 𝑡 ( 1 − 𝑠𝑖𝑛𝜃) 𝑐𝑜𝑠𝜃 (2.26) Donde t es el espesor de la pared de la célula y θ es el ángulo interno entre la horizontal y la pared inclinada como se mostró en el apartado 2.2 de esta memoria. La densidad se obtiene como: 𝜌𝑐 𝜌 =1− 𝑙𝑏 ( ℎ𝑏 + 𝑙𝑏 𝑠𝑖𝑛𝜃 ) 𝑙 ( ℎ + 𝑙 𝑠𝑖𝑛𝜃) (2.27) Los módulos de Young del núcleo en las tres direcciones principales se obtienen mediante las siguientes expresiones: 𝐸1𝑐 = 𝐸 (𝑡𝑙𝑏 ⁄)3𝑐𝑜𝑠𝜃 (ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑠𝑖𝑛2𝜃 [ 1 1 + (2.4 + 1.5𝜈 + 𝑐𝑜𝑡2𝜃) (𝑡𝑙𝑏 ⁄)2 ] (2.28) Modelización teórica de las estructuras tipo honeycomb 24 24 𝐸2𝑐 =𝐸(𝑡𝑙𝑏 ⁄)3 ℎ𝑙 ⁄+𝑠𝑖𝑛𝜃 𝑐𝑜𝑠3𝜃 [ 1 1+ (𝑡𝑙𝑏 ⁄)2(2.4+1.5𝜈+𝑡𝑎𝑛2𝜃+2 ℎ𝑏 𝑙𝑏 𝑐𝑜𝑠2𝜃 ) ] (2.29) Como se mencionaba en el desarrollo de la anterior formulación, el módulo de Young en la dirección que sale del plano de la célula o la dirección X3 se estima asumiendo que la deformación provocada es axial, como se expresa a continuación: 𝐸3𝑐 = 𝐸 𝜌𝑐 𝜌 (2.30) Por último, se obtienen los coeficientes de Poisson con las siguientes expresiones: 𝜈12 𝑐 = 𝑐𝑜𝑠2𝜃 (ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑠𝑖𝑛 𝜃 [ 1+(𝑡𝑙𝑏 ⁄)2(1.4 + 1.5𝜈) 1+(𝑡𝑙𝑏 ⁄)2(2.4+1.5𝜈+𝑐𝑜𝑡2𝜃)] (2.31) 𝜈21 𝑐 = (ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃) 𝑠𝑖𝑛 𝜃 𝑐𝑜𝑠2𝜃 [ 1+ (𝑡𝑙𝑏 ⁄)2(1.4 + 1.5𝜈) 1+(𝑡𝑙𝑏 ⁄)2(2.4+1.5𝜈+𝑡𝑎𝑛2𝜃 + 2 ℎ𝑏 𝑙𝑏 𝑐𝑜𝑠2𝜃) ] (2.32) De forma similar a como se expresa en [40], una vez definidos los módulos de Young en las tres direcciones principales (E1, E2 y E3), se pueden obtener el resto de los coeficientes de Poisson con las siguientes relaciones: 𝜈31 𝑐 = 𝜈32 𝑐 = 𝜈 (2.33) 𝜈13 𝑐 = 𝜈31 𝑐 𝐸1 𝐸3; 𝜈23 𝑐 = 𝜈32 𝑐 𝐸2 𝐸3 (2.34) 25 25 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 2.3.3 Formulación teórica Sorohan 2015 Finalmente, en este apartado se desarrollará la última formulación teórica utilizada en este proyecto para el cálculo de las constantes ingenieriles del núcleo de la célula unitaria hexagonal que corresponde con la formulación de Sorohan expuesta en el artículo [45]. La densidad del material que forma el núcleo de la célula hexagonal se calculará con la siguiente expresión: 𝜌𝑐 𝜌 = 𝑡𝑙 ⁄(ℎ𝑙 ⁄ + 2) 2(ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑐𝑜𝑠𝜃 (2.35) Por otro lado, las expresiones para obtener los módulos de Young en las diferentes direcciones principales son las siguientes: 𝐸1𝑐 = 𝑘1 𝐸 (𝑡𝑙 ⁄)3𝑐𝑜𝑠𝜃 (ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑠𝑖𝑛2𝜃 (2.36) 𝐸2𝑐 = 𝑘2𝐸 (𝑡𝑙 ⁄)3 ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃 𝑐𝑜𝑠3𝜃 (2.37) 𝐸3𝑐 = 𝐸 𝜌𝑐 𝜌 (2.38) Donde: 𝑘1 = 1 1+ (𝑡𝑙 ⁄)2 (2.4 + 1.5𝜈 + 𝑐𝑜𝑡2𝜃) (2.39) 𝑘2 = 1 1+ (𝑡𝑙 ⁄)2 (2.4 + 1.5𝜈 + 𝑡𝑎𝑛2𝜃 + 2 ℎ 𝑙 𝑐𝑜𝑠2𝜃 ) (2.40) Modelización teórica de las estructuras tipo honeycomb 26 26 Los coeficientes de Poisson del núcleo de la célula hexagonal pueden ser calculados como: 𝜈12 𝑐 = 𝑐12 𝑐𝑜𝑠2𝜃 (ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑠𝑖𝑛 𝜃; 𝜈21 𝑐 = 𝑐21(ℎ𝑙 ⁄ + 𝑠𝑖𝑛𝜃)𝑠𝑖𝑛 𝜃 𝑐𝑜𝑠2𝜃 (2.41) Donde: 𝑐12= 1+(𝑡𝑙 ⁄)2(1.4 + 1.5𝜈) 1+(𝑡𝑙 ⁄)2(2.4+1.5𝜈+𝑐𝑜𝑡2𝜃) (2.42) 𝑐21= 1+ (𝑡𝑙 ⁄)2(1.4 + 1.5𝜈) 1+(𝑡𝑙 ⁄)2(2.4+1.5𝜈+ 𝑡𝑎𝑛2𝜃 + 2 ℎ 𝑙 𝑐𝑜𝑠2𝜃) (2.43) Igualmente, el resto de los coeficientes de Poisson se calculan a partir de los módulos de Young previamente calculados: 𝜈13 𝑐 = 𝜈 𝐸1 𝐸3; 𝜈23 𝑐 =𝜈 𝐸2 𝐸3 (2.44) 2.3.4 Formulación teórica para la determinación de las propiedades elásticas efectivas A partir de las constantes ingenieriles obtenidas haciendo uso de las ecuaciones (2.20 - 2.25), (2.26 - 2.34) o (2.35 - 2.44) se puede obtener la matriz de rigidez para materiales ortótropos. En primer lugar, se obtiene la matriz de flexibilidad. Como se ha mencionado en el apartado 2 de esta memoria la matriz de flexibilidad y la matriz de rigidez guardan una relación inversa. De esta forma, se obtienen los coeficientes efectivos con las expresiones expuestas a continuación. 27 27 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 𝑆= [ 1 𝐸1𝑐−𝜈21 𝑐 𝐸2𝑐−𝜈31 𝑐 𝐸3𝑐 −𝜈12 𝑐 𝐸1𝑐1 𝐸2𝑐−𝜈32 𝑐 𝐸3𝑐 −𝜈13 𝑐 𝐸1𝑐−𝜈23 𝑐 𝐸2𝑐1 𝐸3𝑐 0 00 0 01 𝐸10 0 0 00 0 01 𝐸1 0 0 01 𝐸10 0 01 𝐸10 0 1 𝐺23 𝑐 0 0 01 𝐺13 𝑐0 0 0 1 𝐺12 𝑐 ] (2.45) Donde: 𝑪 = 𝑺−𝟏 (2.46) Con estas expresiones, más adelante en esta memoria, se relacionarán las constantes ingenieriles con los coeficientes efectivos, y viceversa. Con el objetivo de aplicar mayor claridad al documento, a continuación, se añaden las expresiones explícitas que calculan los coeficientes efectivos. (2.47) (2.48) (2.49) (2.50) (2.51) Análisis por volumen representativo basado en el MEF 34 34 (3.7) Atendiendo a la expresión (3.6) junto con estas condiciones de contorno (3.7) se obtienen las siguientes condiciones: (3.8) Con estas ecuaciones restrictivas en función de las deformaciones medias impuestas el problema queda completamente definido al fijar un punto cualquiera del volumen, impidiendo su desplazamiento como sólido rígido. 3.1.4 Implementación numérica Para la resolución de las ecuaciones de la elasticidad a las que se les ha impuesto las condiciones de contorno (3.8) se hará uso del Método de los Elementos Finitos (Finite Element Method) o MEF. Este método realiza una discretización del volumen, que pasa a estar constituido por elementos continuos acoplados mediante nodos. El desarrollo de la metodología de cálculo del MEF se encuentra descrita en [43]. Para obtener las propiedades efectivas haciendo uso de estas condiciones de contorno es necesario tomar valores medios de las tensiones y las deformaciones encontrando una relación como (2.1) y (2.9) solo que empleando términos medios: 𝜎𝑖 = 𝐶𝑖𝑗     𝜀𝑗  (3.9) 𝜀𝑖  = 𝑆𝑖𝑗     𝜎𝑗  (3.10) 35 35 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Se busca, por tanto, obtener el tensor de rigidez efectivo que aparece en la expresión (3.9) que de forma matricial se escribe como: [ 𝜎1    𝜎2    𝜎3    𝜎4    𝜎5    𝜎6    ] = [ 𝐶11     𝐶12     𝐶13     0𝐶22     𝐶23     0 0 𝐶33     𝐶14     𝐶15     𝐶16     𝐶24     𝐶25     𝐶26     𝐶34     𝐶35     𝐶36     0 0 0 0𝑆𝑖𝑚 0 0 0 0 𝐶44     𝐶45     𝐶46     0𝐶55     𝐶56     0 0 𝐶66     ] [ 𝜀1  𝜀2  𝜀3  𝜀4  𝜀5  𝜀6  ] (3.11) Para obtener este tensor es necesario resolver 6 problemas independientes, uno por cada deformación media que se va a imponer (𝜀𝑖 ). Para ello se cuenta con las condiciones de contorno desarrolladas previamente y se obtendrán las tensiones medias (𝜎𝑖 ) del volumen representativo. Es fácil observar que para poder obtener las componentes efectivas del tensor de rigidez (𝐶𝑖𝑗     ) bastará con imponer deformaciones medias unitarias. El cálculo de las tensiones medias se trata del sumatorio de los valores de cada uno de los elementos del mallado, por lo que la expresión (3.1) se puede expresar como: 𝜎𝑖 = 1 𝑉𝑅𝑉𝐸∑𝑉𝑒𝜎𝑖𝑒 𝑁 𝑖=1 (3.12) Se utiliza un razonamiento análogo para las deformaciones. Se expresa la relación (3.2) como el sumatorio: 𝜀𝑖  = 1 𝑉𝑅𝑉𝐸∑𝑉𝑒𝜀𝑖𝑒 𝑁 𝑖=1 (3.13) 3.2 Implementación en ANSYS En este apartado se va a presentar el modelo de célula honeycomb implementado en el software ANSYS, haciendo inciso en el mallado elegido para obtener los resultados. Así mismo, se presentan las condiciones de contorno en deformaciones que es necesario aplicar al modelo para obtener las propiedades efectivas. Análisis por volumen representativo basado en el MEF 36 36 3.2.1 Modelización en ANSYS El cálculo de las constantes ingenieriles (2.20) – (2.44) y de los coeficientes efectivos (2.45) – (2.46) de la estructura, así como el consiguiente estudio de la influencia del ángulo interno de la estructura se va a abordar mediante el Método de Elementos Finitos (MEF) en el software de elementos finitos ANSYS. A continuación, se define la metodología que se ha llevado a cabo para realizar la modelización del RVE usando el software ANSYS. En primer lugar, para trabajar con el RVE en ANSYS se realizará un modelo 3D del volumen representativo en el software de elementos finitos representando la geometría que se ha definido en el apartado 2.2. Esta geometría dependerá de los parámetros definidos en el apartado 2.2 que nos permitirán estudiar el RVE para diferentes geometrías. Se procederá indicando el tipo de elemento con el que se va a modelar el RVE que será un elemento SOLID 186. Además, se definirán las propiedades del material E, ν y ρ. Se estudiará un RVE en el que el núcleo de la matriz estará constituido por un material con unos valores de E, ν y ρ definidos mientras que los huecos del RVE estarán vacíos. Esto se modelará dándole a los huecos unos valores de E, ν y ρ infinitamente menores que los valores dados al núcleo. Una vez definidas las propiedades del material, se procederá a mallar el modelo del RVE. Este mallado se realizará realizando divisiones en cada línea que conforma la celda hexagonal. Se busca un mallado lo suficientemente fino para que los resultados sean precisos, pero suficientemente grande para que no requiera unos costes computacionales demasiado elevados. Tras realizar el mallado del modelo RVE, se aplican las condiciones de contorno CC1, CC2 y CC3 que además imponen cada una de ellas una condición de deformación para obtener los coeficientes efectivos (2.45)- (2.46). Estas condiciones de contorno están definidas en el apartado 3.2.3.1. Una vez aplicadas se pedirá a ANSYS que abra un documento de texto donde almacenará estos coeficientes. Una vez realizadas todas las validaciones necesarias para verificar la eficiencia de este modelo se realizará el estudio de la influencia del ángulo interno de la matriz hexagonal en el RVE donde se extraerán los valores de los desplazamientos y las tensiones más representativos de cada condición de contorno. Por tanto, la metodología de cálculo en ANSYS de los coeficientes efectivos se describirá en el apartado 3.2.3. El proceso descrito en este apartado se realizará idéntico para dos modelos distintos de la estructura, uno para los valores del ángulo interno θ menores que 0º y otro para los valores positivos de θ. Esto se hace así para garantizar la exactitud de los resultados para altos rangos de valores del ángulo θ, ya que conforme en ángulo cambia de signo se pierde precisión de la geometría del modelo. Los dos modelos realizados son los que se observan en la Figura 3-6. 37 37 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos (a) (b) Figura 3-6. (a) Célula hexagonal para valores de θ > 0º (b) Célula hexagonal para valores de θ < 0º. 3.2.2 Elección del mallado Para la elección del espesor del mallado se lleva a cabo un proceso iterativo en el que se busca un mallado suficientemente fino para obtener resultados precisos sin que resulte un coste computacional demasiado elevado. Para mejorar la malla predefinida de ANSYS, se mallan manualmente, una a una, las líneas de la estructura. Iremos iterando con el número de divisiones de las líneas más grandes, ya que, estas tendrán un mayor número de divisiones y, por tanto, un mayor peso en la precisión que va a tener la malla. Este estudio de sensibilidad se llevará a cabo con mallas de elementos triangulares ya que este permite una mejor transición entre elementos de distintos tamaños. El mallado finalmente elegido que proporciona resultados con un error suficientemente pequeño que permite abordar el problema del estudio de la estructura honeycomb sin tener un tiempo computacional demasiado elevado es el que se muestra en la Figura 3-7 y la Figura 3-8, una por cada modelo de ANSYS que se ha implementado. Análisis por volumen representativo basado en el MEF 38 38 Figura 3-7. Mallado elegido para modelar la célula hexagonal para θ = 30º. Figura 3-8. Mallado elegido para modelar la célula hexagonal para θ = -30º 39 39 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 3.2.3 Obtención de las propiedades efectivas En este trabajo se ha considerado que el material del que está formado el núcleo de la estructura honeycomb es ortótropa y que los huecos están vacíos. Por ser ortótropo la ley de comportamiento queda definida mediante 9 constantes. Las ecuaciones constitutivas se expresan de forma matricial como se muestra a continuación: [ 𝜎1 𝜎2 𝜎3 𝜎23 𝜎13 𝜎12 ] = [ 𝐶11 𝐶12 𝐶13 𝐶21 𝐶22 𝐶23 𝐶31 𝐶32 𝐶33 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 𝐶44 0 0 0 𝐶55 0 0 0 𝐶66 ] [ 𝜀1 𝜀2 𝜀3 𝛾23 𝛾13 𝛾12 ] (3.14) Los coeficientes efectivos se obtendrán, como se ha mencionado anteriormente, imponiendo valores constantes de las deformaciones. Para ello, se han realizado 3 grupos de condiciones de contorno que permiten obtener las tres primeras columnas de la matriz de rigidez a partir de las ecuaciones (3.14) imponiendo un valor unitario de las deformaciones. La metodología que se ha seguido se trata de aplicar las condiciones de contorno al modelo RVE de elementos finitos, se obtienen los valores medios de las tensiones y a partir de ellos se obtienen los valores de los coeficientes efectivos. Seguidamente, se definen cada una de estas condiciones de contorno y los coeficientes efectivos que se obtienen con cada una de ellas. 3.2.3.1 Condiciones de contorno • CC1 Para obtener la primera columna de la matriz de rigidez se impone unitaria la deformación en la dirección principal x1 (ε1 = 1). Además, se impone que el resto de las componentes de la deformación sean nulas (εij = γij = 0). Mediante esta condición de contorno se calculan los coeficientes efectivos 𝐶11 𝑒𝑓𝑓, 𝐶21 𝑒𝑓𝑓y 𝐶31 𝑒𝑓𝑓. 𝐶11 𝑒𝑓𝑓=𝜎11     𝜀11     ⁄ (3.15) 𝐶21 𝑒𝑓𝑓=𝜎22     𝜀11     ⁄ (3.16) Análisis por volumen representativo basado en el MEF 40 40 𝐶31 𝑒𝑓𝑓=𝜎33     𝜀11     ⁄ (3.17) • CC2 Para obtener la segunda columna de la matriz de rigidez se impone unitaria la deformación en la dirección principal x2 (ε2 = 1). Además, se impone que el resto de las componentes de la deformación sean nulas (εij = γij = 0). Mediante esta condición de contorno se calculan los coeficientes efectivos 𝐶12 𝑒𝑓𝑓, 𝐶22 𝑒𝑓𝑓y 𝐶32 𝑒𝑓𝑓. 𝐶12 𝑒𝑓𝑓=𝜎11     𝜀22     ⁄ (3.18) 𝐶22 𝑒𝑓𝑓=𝜎22     𝜀22     ⁄ (3.19) 𝐶32 𝑒𝑓𝑓=𝜎33     𝜀22     ⁄ (3.20) • CC3 Para obtener la tercera columna de la matriz de rigidez se impone unitaria la deformación en la dirección principal x3 (ε3 = 1). Además, se impone que el resto de las componentes de la deformación sean nulas (εij = γij = 0). Mediante esta condición de contorno se calculan los coeficientes efectivos 𝐶13 𝑒𝑓𝑓, 𝐶23 𝑒𝑓𝑓y 𝐶33 𝑒𝑓𝑓. 𝐶13 𝑒𝑓𝑓=𝜎11     𝜀33     ⁄ (3.21) 𝐶23 𝑒𝑓𝑓=𝜎22     𝜀33     ⁄ (3.22) 𝐶33 𝑒𝑓𝑓=𝜎33     𝜀33     ⁄ (3.23) 41 41 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 3.2.4 Obtención de las constantes ingenieriles Anteriormente se ha mencionado como se calculan los coeficientes efectivos a partir de las constante ingenieriles haciendo uso de la formulación teórica. En este apartado se explica cómo se van a calcular las constantes ingenieriles realizando el proceso inverso. Una vez calculadas en ANSYS los coeficientes efectivos imponiendo las condiciones de contorno mencionadas en el apartado anterior (Apartado 3.2.3.1), se pueden obtener las constantes ingenieriles mediante las relaciones (2.45) y (2.46). Con el objetivo de aclarar el proceso de cálculo que se ha seguido se añaden seguidamente las expresiones explícitas para calcular las constantes ingenieriles a partir del tensor de flexibilidad, que se ha obtenido aprovechando la relación inversa que guarda con el tensor de rigidez obtenido en ANSYS. 𝐸1    = 1 𝑆11 𝑒𝑓𝑓 (3.24) 𝐸2    = 1 𝑆22 𝑒𝑓𝑓 (3.25) 𝐸3    = 1 𝑆33 𝑒𝑓𝑓 (3.26) 𝜈12     = –𝑆12 𝑒𝑓𝑓𝐸1    𝜈21     = –𝑆21 𝑒𝑓𝑓𝐸2    (3.27) 𝜈13     = –𝑆13 𝑒𝑓𝑓𝐸1    𝜈31     = –𝑆31 𝑒𝑓𝑓𝐸3    (3.28) 𝜈23     = –𝑆23 𝑒𝑓𝑓𝐸2    𝜈32     = –𝑆32 𝑒𝑓𝑓𝐸3    (3.29) 𝐺23      = 1 𝑆44 𝑒𝑓𝑓 (3.30) 𝐺13      = 1 𝑆55 𝑒𝑓𝑓 (3.31) 𝐺12      = 1 𝑆66 𝑒𝑓𝑓 (3.32) Análisis por volumen representativo basado en el MEF 42 42 Para obtener la densidad del núcleo se utiliza la expresión que la relaciona el módulo de Young del núcleo en la dirección principal X3. 𝜌= 𝜌𝐸3 𝐸 (3.33) Por tanto, se ilustra en la Figura 3-9 el proceso ha seguir partiendo del RVE en ANSYS para finalmente obtener las constantes ingenieriles del panel hexagonal. Figura 3-9. Diagrama de flujo de la metodología empleada para obtener las constantes ingenieriles del modelo RVE en ANSYS. 43 43 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 4 RESULTADOS n este capítulo se desarrollará la validación de la formulación teórica presentada en los capítulos 2 y 3 para el cálculo de las constantes ingenieriles y los coeficientes efectivos de las estructuras honeycombs. Asimismo, se validará el modelo de elementos finitos RVE en ANSYS comparando los resultados obtenidos con los resultados de las expresiones analíticas. Para ello se implementarán los archivos de Matlab y ANSYS usando la configuración de la celda utilizada en [45]. Además, se realizará un ejemplo para un material cualquiera, en este caso un aluminio. Una vez validadas estas expresiones y el modelo RVE se realizará un estudio de la influencia de la relación t/h para la misma configuración de la celda unitaria hexagonal. Por último, se realizará un estudio de la influencia del ángulo interno en el modelo RVE según el caso 6 resuelto en [45]. 4.1 Comparativa entre los modelos teóricos y el MEF (RVE) En este apartado se va a realizar la validación de la implementación numérica de los modelos teóricos (2.20) – (2.44) y el modelo MEF. Para ello se realizará la comparación de los valores obtenidos con estas expresiones y los valores obtenidos en el artículo de Sorohan [45] para las constantes ingenieriles. Una vez validadas estas expresiones se utilizarán las relaciones (2.45) y (2.46) para obtener los coeficientes efectivos de las tres primeras columnas de la matriz de rigidez. Por tanto, se tomarán los mismos valores para la geometría de la célula hexagonal que los utilizados en el artículo mencionado anteriormente. Esta configuración está definida por los siguientes valores: t = 1 mm l = h = 10 mm b = 15 mm θ = 30º Además, para las propiedades del material se tomarán los siguientes valores también coincidentes con los valores del artículo: E = 1000 MPa ν = 0,3 ρ = 1000 kg/m3 En la siguiente tabla se encuentran recogidos los valores obtenidos mediante las expresiones teóricas para los datos de este ejemplo. La primera columna corresponde a los resultados obtenidos haciendo E Resultados 50 50 Tabla 4-6. Resultados de los coeficientes efectivos para las diferentes formulaciones para el caso 2. Gibson 1982 Malek & Gibson 2015 RVE-MEF C11 7,9396 9,6045 8,1047 C12 7,9216 9,5898 7,9697 C13 4,7584 5,7583 4,8223 C21 7,9216 9,5898 7,9697 C22 7,9396 9,6127 7,8756 C23 4,7584 5,7607 4,7536 C31 4,7584 5,7583 4,8223 C32 4,7584 5,7607 4,7536 C33 31,7225 40,8658 31,3300 4.3 Estudio de la influencia del ángulo interno θ Para este estudio se va a realizar partiendo del caso 6.1 resuelto en [45] de donde se tomarán los datos de partida y las expresiones teóricas para resolver el problema. Se va a considerar una estructura del tipo honeycomb que va a presentar los siguientes valores para su geometría: t = 1 mm l = 10 mm h = 20 mm b = 10 mm El ángulo interno de la célula hexagonal θ irá tomando valores entre -45º y 45º. En este caso, para poder modelar la estructura en ANSYS es necesario que las longitudes l y h no sean iguales, sino que h. Esto se debe a que en caso de que fueran iguales existirían ciertos ángulos en los que la geometría en ANSYS se sobrepondría y no se podría calcular. En este apartado, el material que se va a emplear es el mismo que en los ejercicios realizados anteriormente y, por tanto, presenta las siguientes propiedades: E = 1000 MPa ν = 0.3 ρ = 1000 kg/m3 La metodología que se va a seguir consiste en obtener las constantes ingenieriles para la estructura definida por los parámetros de geometría y propiedades del material, previamente descritos en este apartado, para ver como varían estas constantes. A su vez, se van a aplicar las condiciones de contorno definidas en el apartado 3.2.3.1 de esta memoria, al modelo RVE de elementos finitos para valores del 51 51 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos ángulo θ también variando entre -45º y 45º. De esta forma se obtendrán los coeficientes efectivos de la matriz de rigidez. A partir de las mismas, se calcularán las constantes ingenieriles del material haciendo uso de las relaciones (2.45) – (2.46). Se podrá entonces comparar la variación de las propiedades efectivas con el ángulo por estos dos métodos. Para ello se implementa una rutina de MATLAB para graficar estas constantes ingenieriles al variar con el ángulo interno θ usando como formulación teórica las expresiones (2.35) – (2.44) frente a los resultados obtenidos mediante el MEF en ANSYS para los tres grupos de condiciones de contorno. Estas gráficas se encuentran en [45] de forma que se pueden validar los resultados. Una vez graficadas estas propiedades y validada la exactitud del modelo RVE en ANSYS se estudiarán las distribuciones de desplazamientos y tensiones más significativas para cada condición de contorno impuesta. 4.3.1 Variación de las constantes ingenieriles con el ángulo θ Se puede observar que la densidad del núcleo tiene un mínimo cuando el ángulo interno θ es aproximadamente 20º. Las estimaciones analíticas de la densidad son mayores que los resultados obtenidos para el modelo de elementos finitos del RVE (ver Figura 4-1). Se distingue una buena precisión en los resultados del modelo de ANSYS con respecto a los resultados obtenidos de la formulación teórica de Sorohan [45]. Figura 4-4. Variación de la densidad frente al ángulo θ. Comparativa de un modelo teórico y un modelo resuelto mediante el MEF en ANSYS. Resultados 52 52 En cuanto a los módulos de Young, en la Figura 4-5 se puede observar una buena coincidencia entre los valores analíticos y los que han sido obtenidos por el MEF. Se puede distinguir que los valores teóricos de E2 están por debajo de los calculados en ANSYS cuando el ángulo θ es menor que 0, mientras que para valores positivos del mismo estos se encuentran por encima. Además, el error disminuye cuando el ángulo interno θ aumenta. El valor de E1 tiene una discontinuidad cuando θ = 0º. Esto se debe a la presencia de un coseno en el numerador de la expresión utilizada para calcular esta constante (ver expresión 2.36). ANSYS por el contrario si calcula para este valor del ángulo, pero se observa que ambos tienen la misma tendencia. Figura 4-5. Variación del Módulo de Young frente al ángulo θ. Comparativa de un modelo teórico y un modelo resuelto mediante el MEF en ANSYS. Se observa en la Figura 4-6, que muestra los coeficientes de Poisson, una precisión suficiente que asegura una buena solución del modelo de elementos finitos por el método del RVE en ANSYS. Vuelve a ocurrir una discontinuidad, en este caso en la gráfica de 𝜈12 𝑐, debido igualmente a la presencia de un coseno en el numerador de la expresión que calcula este coeficiente (ver expresión 2.41). 53 53 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Figura 4-6. Variación del Coeficiente de Poisson frente al ángulo θ. Comparativa de un modelo teórico y un modelo resuelto mediante el MEF en ANSYS. Resultados 54 54 4.3.2 Condición de contorno 1 θ = -30º θ = -30º θ = 0º θ = 0º θ = 30º θ = 30º (a) (b) Figura 4-7. (a) Distribución del desplazamiento 𝑢𝑥 para los ángulos θ = {-30º, 0º, 30} (b) Distribución de tensiones 𝜎1 para los ángulos θ = {-30º, 0º, 30º} 55 55 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 4.3.3 Condición de contorno 2 θ = -30º θ = -30º θ = 0º θ = 0º θ = 30º θ = 30º Figura 4-8. (a) Distribución del desplazamiento 𝑢𝑦 para los ángulos θ = {-30º, 0º, 30} (b) Distribución de tensiones 𝜎2 para los ángulos θ = {-30º, 0º, 30º} Resultados 56 56 4.3.4 Condición de contorno 3 θ = -30º θ = -30º θ = 0º θ = 0º θ = 30º θ = 30º Figura 4-9. (a) Distribución del desplazamiento 𝑢𝑧 para los ángulos θ = {-30º, 0º, 30} (b) Distribución de tensiones 𝜎3 para los ángulos θ = {-30º, 0º, 30º} 57 57 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 5 RESUMEN Y CONCLUSIONES inalmente, en este capítulo se realizará un resumen del trabajo que se ha llevado a cabo y se enumerarán las principales conclusiones a las que se ha llegado. Además, se plantearán posibles trabajos futuros que pueden partir del estudio realizado en este proyecto o la metodología empleada. 5.1. Resumen y conclusiones En la presente memoria se ha expuesto un modelo de estructura honeycomb en función de ciertos parámetros variables que definen la geometría de la celda. Este modelo se ha implementado en ANSYS y se ha comparado con los resultados obtenidos de distintas rutinas de MATLAB. Los archivos desarrollados en MATLAB obtienen las constantes ingenieriles y los coeficientes efectivos de una estructura honeycomb para un material y una geometría determinados por diferentes métodos. A este modelo se le han aplicado una serie de condiciones de contorno, en las que se imponen deformaciones unitarias, con el fin de obtener la matriz de rigidez y con ella las constantes ingenieriles. Se han realizado todos los estudios para una celda con un núcleo constituido por un cierto material y unos huecos vacíos. En primer lugar, se ha realizado la validación de las ecuaciones implementadas en ANSYS comparándolas con el estudio realizado en [45]. Se ha comprobado la validez de los métodos de cálculo implementados en las rutinas de MATLAB para obtener las constantes ingenieriles y los coeficientes efectivos. Esta validación se ha realizado comparando los resultados analíticos con el estudio de una estructura honeycomb en [45]. Una vez validados los métodos de obtención de resultados, se ha disminuido el valor del espesor de la pared del núcleo para observar la influencia de la relación de apariencia t/h manteniendo el resto de los parámetros geométricos idénticos y utilizando el mismo material. Por último, se ha estudiado la influencia del ángulo interno de la pared del núcleo en los valores de las constantes ingenieriles. A modo ilustrativo, se han añadido las representaciones de las tensiones y los desplazamientos más significativos en cada condición de contorno (ver Figura 4-4, 4-5 y 4-6). De los estudios realizados y a la vista de los resultados obtenidos se puede concluir: • Se ha comprobado la exactitud y validez del modelo de celda hexagonal implementado en ANSYS, comparando los resultados obtenidos imponiendo condiciones de contorno en deformaciones con los resultados del estudio realizado en [45]. • Cabe, por tanto, destacar el grado de concordancia tan alto resultante entre los resultados analíticos y los resultados del RVE mediante el método de elementos finitos en ANSYS. • Las condiciones de contorno (3.15) – (3.23) permiten obtener con gran exactitud las F 58 58 componentes de la matriz de rigidez y posteriormente las constantes ingenieriles. • Los archivos de ANSYS y MATLAB pueden ser fácilmente adaptados para calcular otros materiales y otras geometrías de celdas de estructuras honeycombs. • El espesor del mallado es suficiente para obtener resultados con un alto grado de acuerdo con los valores teóricos y no supone un coste computacional elevado. • Del estudio de la influencia en los resultados de la variación de la relación de aspecto t/h (caso 1 y caso 2) se obtienen unos valores de los resultados muy inferiores cuando el espesor tiene una influencia menor que el espesor utilizado en el caso anterior. Esto se explica observando las ecuaciones (2.20) – (2.44) donde el espesor se encuentra multiplicando en la mayoría de las ecuaciones. • Se concluye que los modelos RVE funcionan bien para relaciones t/h grandes y pequeñas. Se ha trabajado con los valores de la relación [0.1 0.05 0.025] y se observa que los acuerdos con los modelos teóricos son aceptables. • El RVE trabaja correctamente bajo condiciones auxéticas, ya que observamos un coeficiente de Poisson negativo para θ < 0º en la Figura 4-6. • De las Figuras 4-7, 4-8 y 4-9 también puede observarse que las tensiones se concentran en el núcleo de la celda. Esto es debido a que los huecos de la célula hexagonal se encuentran vacíos de material. 5.2. Aplicación para trabajos futuros Para trabajos futuros, sería interesante estudiar la influencia del ángulo interior cuando el material existente en los huecos de la matriz tiene valores de E, ν y ρ cada vez más cercanos a los de la matriz. En este trabajo se ha supuesto que estos huecos están vacíos implementándolo en ANSYS como un material con propiedades E, ν y ρ infinitamente más pequeñas que las de la matriz y, por tanto, despreciables. Otro estudio que podría estar basado en la metodología seguida en este proyecto es el estudio del efecto que tendría la variación de la profundidad b en el valor de las constantes ingenieriles. Podría validarse este estudio con uno similar realizado en el apartado 6.2 del artículo [45]. Para validar los resultados también se puede recurrir al estudio de la influencia del espesor del núcleo realizado en [3] adaptando la geometría de la celda a la utilizada en dicho artículo. 59 59 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos REFERENCIAS [1] L. Zhang, B. Liu, Y. Gu y X. H. Xu, ‹‹Modelling and characterization of mechanical properties of optimized honeycomb structure››, International Journal of Mechanics and Materials in Design, vol. 16, pp. 155-166, 2019. [2] Z. Wang, J. Liu, Z. Lu y D. Hui, ‹‹Mechanical behavior of composited structure filled with tandem honeycombs››, Composites Part B: Engineering, vol. 114, pp. 128-138, 2017. [3] W. Becker, ‹‹Closed-form analysis of the thickness effect of regular honeycomb core material››, Composite Structures, vol. 48, nº 1-3, pp. 67-70, 2000. [4] M. He y W. Hu, ‹‹A study of composite honeycomb sándwich panel structure››, Materials & Design, vol. 29, nº 3, pp. 709-713, 2008. [5] Z. Wang, J. Liu y D. Hui, ‹‹Mechanical behaviors of inclined cell honeycomb structure subjected to compression. ››, vol. 110, pp. 307-314,2017. [6] C. Qi, F. Jiang y S. Yang, ‹‹Advanced honeycomb designs for improving mechanical properties: A review››, Composites Part B: Engineering, vol. 227, 109393, 2021. [7] E. Sather y T. Krishnamurthy, ‹‹An Analytical Method to Calculate Effective Elastic Properties of General Multifunctional Honeycomb Cores in Sandwich Composites››, NASA/TM–2019–220275, 2019. [8] N.D. Duc y P.M. Cong, ‹‹Nonlinear dynamic response and vibration of sandwich composite plates with negative Poisson’s ratio in auxetic honeycombs››, Journal of Sandwich Structures and Materials, vol. 20, nº 6, pp 692-717 ,2018. [9] R. Moradi-Dastjerdi y K. Behdinan, ‹‹Electromechanical Natural Frequency Analysis of an Eco- Friendly Active Sandwich Plate››, Actuators, vol. 11, nº 261, 2022. [10] B. Castanie, C. Bouvet y M. Ginot, ‹‹Review of composite sandwich structure in aeronautic applications››, Composites Part C: Open Access, vol. 1, 2020 [11] S.E. Sadiq, M. J. Jweeg y S. H. Bakhy, ‹‹Strength Analysis of an Aircraft Sandwich Structure with a Honeycomb Core: Theoretical and Experimental Approaches››, Engineering and Technology Journal, vol. 39, nº 1, pp. 153-166, 2021. Anexos 66 66 % le damos valores unitarios mietras tanto: C44 = 1; C55 = 1; C66 = 1; %------------------------------------------------------- % Matiz C %------------------------------------------------------- C=[C11 C12 C13 0 0 0; C12 C22 C23 0 0 0; C13 C23 C33 0 0 0; 0 0 0 C44 0 0; 0 0 0 0 C55 0; 0 0 0 0 0 C66]; %------------------------------------------------------- % Cálculo de la constantes ingenieriles %------------------------------------------------------- S = inv(C(1:3,1:3)); S = inv(C) Ec11 = 1/S(1,1); Ec22 = 1/S(2,2); Ec33 = 1/S(3,3); NUc12=-S(1,2)*Ec11; NUc13=-S(1,3)*Ec11; NUc23=-S(2,3)*Ec22; NUc21=-S(2,1)*Ec22; NUc31=-S(3,1)*Ec33; NUc32=-S(3,2)*Ec33; Gc23 = 1/S(4,4); Gc13 = 1/S(5,5); 67 67 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Gc12 = 1/S(6,6); rho_c=Ec33*rho/E; [rho_c Ec11 Ec22 Ec33 Gc12 Gc23 Gc13 NUc12 NUc21 NUc23 NUc32 NUc13 NUc31] • Ecuaciones_gibson1982_1 Archivo que calcula la matriz C a partir de la formulación [40] para el apartado 4.1. % ------------------------------------------ % Material Properties %----------------------------------------- E = 1000; % [MPa] nu = 0.3; rho = 1000; % [kg/m3] G = E/(2*(1+nu)); % [Pa] %------------------------------------ % Geometría de un panel %------------------------------------ h = 10; l = 10; t=1; % espesor Anexos 68 68 theta= 30*pi/180; % ANGULO b=15; a1 = h + l*sin(theta); a2 = t/2 + l*cos(theta); a3 = b/2; alpha = pi/2 - theta; % Constantes Ingenieriles Ec11=[]; Ec22=[]; Ec33=[]; Gc12=[]; Gc13=[]; Gc23=[]; NUc12=[]; NUc13=[]; NUc23=[]; rho_c=[]; for i=1:max(size(theta)) alpha_x = h/l; %ratio of cell wall length beta_x = t/l; %slenderness ratio of the unit of cell wall % a1=l*( beta_x/2 + cos(theta(i)) ); % valores de aproximadamente=1 % a2=l*( alpha_x + sin(theta(i)) ); % classical Honeycomb rho_c_aux = rho*beta_x*( alpha_x + 2 )/( 2*cos(theta(i))*( alpha_x + sin(theta(i)) ) ); k1 = 1/(1 + (beta_x^2)*(2.4 + 1.5*nu + (1/tan(theta(i)))^2) ); k2 = 1/(1 + (beta_x^2)*(2.4 + 1.5*nu + tan(theta(i))^2 + 2*alpha_x*(1/cos(theta(i)))^2) ); 69 69 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos c12= k1*(1 + (beta_x^2)*(1.4 + 1.5*nu)); c21= k2*(1 + (beta_x^2)*(1.4 + 1.5*nu)); k12=(1+2*h/l)/( 1+2*h/l+(t/l)^2*(l/h*(2.4+1.5*nu)*(2+alpha_x+sin(theta(i)))+ ((alpha_x+sin(theta(i)))/alpha_x^2 ))*((alpha_x+sin(theta(i)))*tan(theta(i))^2+sin(theta(i))) ); Ec11_aux = k1*E*( (beta_x)^3 )*cos(theta(i))/( (alpha_x + sin(theta(i)))*sin(theta(i))^2 ); Ec22_aux = k2*E*( (beta_x)^3 )*( alpha_x + sin(theta(i)) )/( cos(theta(i))^3 ); Ec33_aux = E*rho_c_aux/rho; rho_c=[rho_c rho_c_aux]; Ec11= [Ec11 Ec11_aux]; Ec22= [Ec22 Ec22_aux]; Ec33= [Ec33 Ec33_aux]; Gc12_aux = k12*E*(beta_x^3)*( alpha_x + sin(theta(i)) )/( (alpha_x^2)*( 1 + 2*alpha_x )*cos(theta(i)) ); Gc13_aux = G*beta_x*cos(theta(i))/( alpha_x + sin(theta(i)) ); Gc23_aux = G*beta_x*( 1 + 2*sin(theta(i))^2 )/( 2*cos(theta(i))*( alpha_x + sin(theta(i)) ) ); Gc12= [Gc12 Gc12_aux]; Gc13= [Gc13 Gc13_aux]; Gc23= [Gc23 Gc23_aux]; NUc12_aux = c12*( cos(theta(i))^2 )/( (alpha_x+sin(theta(i)))*sin(theta(i)) ); NUc13_aux = nu*Ec11_aux/Ec33_aux; NUc23_aux = nu*Ec22_aux/Ec33_aux; NUc12= [NUc12 NUc12_aux]; NUc13= [NUc13 NUc13_aux]; NUc23= [NUc23 NUc23_aux]; end NUc32=NUc23.*Ec33./Ec22; Anexos 70 70 NUc31=NUc13.*Ec33./Ec11; NUc21=c21*( (alpha_x+sin(theta(i)))*sin(theta(i)) )/( cos(theta(i))^2 ); %NUc21=NUc12.*Ec22./Ec11; [rho_c Ec11 Ec22 Ec33 Gc12 Gc23 Gc13 NUc12 NUc21 NUc23 NUc32 NUc13 NUc31] %------------------------------------------------------- % Cálculo de la matiz C %------------------------------------------------------- % Matriz flexibilidad S11=1/Ec11; S12=-NUc12/Ec11; S13=-NUc13/Ec11; S22=1/Ec22; S23=-NUc23/Ec22; S33=1/Ec33; S44=1/Gc23; S55=1/Gc13; S66=1/Gc12; S=[S11 S12 S13 0 0 0; 71 71 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos S12 S22 S23 0 0 0; S13 S23 S33 0 0 0; 0 0 0 S44 0 0; 0 0 0 0 S55 0; 0 0 0 0 0 S66] % Matriz C C=inv(S) • Ecuaciones_gibson1982_2 Archivo que calcula la matriz C a partir de la formulación [40] para el Apartado 4.2. % ------------------------------------------ % Material Properties %----------------------------------------- E = 1000; % [MPa] nu = 0.3; rho = 1000; % [kg/m3] G = E/(2*(1+nu)); % [Pa] %------------------------------------ % Geometría de un panel %------------------------------------ % Caso 1 h = 20; l = 20; % Caso 2 % h = 40; % l= 40; t = 1; % espesor theta= 30*pi/180; % ANGULO b=15; Anexos 72 72 a1 = h + l*sin(theta); a2 = t/2 + l*cos(theta); a3 = b/2; alpha = pi/2 - theta; % Constantes Ingenieriles Ec11=[]; Ec22=[]; Ec33=[]; Gc12=[]; Gc13=[]; Gc23=[]; NUc12=[]; NUc13=[]; NUc23=[]; rho_c=[]; for i=1:max(size(theta)) alpha_x = h/l; %ratio of cell wall length beta_x = t/l; %slenderness ratio of the unit of cell wall % a1=l*( beta_x/2 + cos(theta(i)) ); % valores de aproximadamente=1 % a2=l*( alpha_x + sin(theta(i)) ); % classical Honeycomb rho_c_aux = rho*beta_x*( alpha_x + 2 )/( 2*cos(theta(i))*( alpha_x + sin(theta(i)) ) ); k1 = 1/(1 + (beta_x^2)*(2.4 + 1.5*nu + (1/tan(theta(i)))^2) ); k2 = 1/(1 + (beta_x^2)*(2.4 + 1.5*nu + tan(theta(i))^2 + 2*alpha_x*(1/cos(theta(i)))^2) ); c12= k1*(1 + (beta_x^2)*(1.4 + 1.5*nu)); c21= k2*(1 + (beta_x^2)*(1.4 + 1.5*nu)); 73 73 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos k12=(1+2*h/l)/( 1+2*h/l+(t/l)^2*(l/h*(2.4+1.5*nu)*(2+alpha_x+sin(theta(i)))+ ((alpha_x+sin(theta(i)))/alpha_x^2 ))*((alpha_x+sin(theta(i)))*tan(theta(i))^2+sin(theta(i))) ); Ec11_aux = k1*E*( (beta_x)^3 )*cos(theta(i))/( (alpha_x + sin(theta(i)))*sin(theta(i))^2 ); Ec22_aux = k2*E*( (beta_x)^3 )*( alpha_x + sin(theta(i)) )/( cos(theta(i))^3 ); Ec33_aux = E*rho_c_aux/rho; rho_c=[rho_c rho_c_aux]; Ec11= [Ec11 Ec11_aux]; Ec22= [Ec22 Ec22_aux]; Ec33= [Ec33 Ec33_aux]; Gc12_aux = k12*E*(beta_x^3)*( alpha_x + sin(theta(i)) )/( (alpha_x^2)*( 1 + 2*alpha_x )*cos(theta(i)) ); Gc13_aux = G*beta_x*cos(theta(i))/( alpha_x + sin(theta(i)) ); Gc23_aux = G*beta_x*( 1 + 2*sin(theta(i))^2 )/( 2*cos(theta(i))*( alpha_x + sin(theta(i)) ) ); Gc12= [Gc12 Gc12_aux]; Gc13= [Gc13 Gc13_aux]; Gc23= [Gc23 Gc23_aux]; NUc12_aux = c12*( cos(theta(i))^2 )/( (alpha_x+sin(theta(i)))*sin(theta(i)) ); NUc13_aux = nu*Ec11_aux/Ec33_aux; NUc23_aux = nu*Ec22_aux/Ec33_aux; NUc12= [NUc12 NUc12_aux]; NUc13= [NUc13 NUc13_aux]; NUc23= [NUc23 NUc23_aux]; end NUc32=NUc23.*Ec33./Ec22; NUc31=NUc13.*Ec33./Ec11; NUc21=c21*( (alpha_x+sin(theta(i)))*sin(theta(i)) )/( cos(theta(i))^2 ); %NUc21=NUc12.*Ec22./Ec11; Anexos 74 74 [rho_c Ec11 Ec22 Ec33 Gc12 Gc23 Gc13 NUc12 NUc21 NUc23 NUc32 NUc13 NUc31] %------------------------------------------------------- % Cálculo de la matiz C %------------------------------------------------------- % Matriz flexibilidad S11=1/Ec11; S12=-NUc12/Ec11; S13=-NUc13/Ec11; S22=1/Ec22; S23=-NUc23/Ec22; S33=1/Ec33; S44=1/Gc23; S55=1/Gc13; S66=1/Gc12; S=[S11 S12 S13 0 0 0; S12 S22 S23 0 0 0; S13 S23 S33 0 0 0; 0 0 0 S44 0 0; 75 75 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 0 0 0 0 S55 0; 0 0 0 0 0 S66] % Matriz C C=inv(S) • Ecuaciones_gibson2015 Archivo que calcula la matriz C mediante la formulación [42] para el Apartado 4.1. %---------------------------------------------------------------- % Archivo Definitivo % Calculo de la matriz C a partir de las ctes ingenieriles %---------------------------------------------------------------- % Propiedades del material E = 1000; % [MPa] nu = 0.3; rho = 1000; % [kg/m3] G = E/(2*(1+nu)); % [Pa] % Geometría del material h=10; l=10; t = 1; theta= 30*pi/180; % angulo b=15; alpha = pi/2 - theta; % Constantes Ingenieriles Ec11=[]; Ec22=[]; Ec33=[]; Gc12=[]; Anexos 82 82 Ec11 Ec22 Ec33 Gc12 Gc23 Gc13 NUc12 NUc21 NUc23 NUc32 NUc13 NUc31] %------------------------------------------------------ % Cálculo de la matiz C %------------------------------------------------------ % Matriz flexibilidad S11=1/Ec11; S12=-NUc12/Ec11; S13=-NUc13/Ec11; S22=1/Ec22; S23=-NUc23/Ec22; S33=1/Ec33; S44=1/Gc23; S55=1/Gc13; S66=1/Gc12; S=[S11 S12 S13 0 0 0; S12 S22 S23 0 0 0; S13 S23 S33 0 0 0; 0 0 0 S44 0 0; 83 83 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos 0 0 0 0 S55 0; 0 0 0 0 0 S66] % Matriz C C=inv(S) • Mat_properties_MEF_1 Rutina para calcular la matriz de rigidez y las constantes ingenieriles a partir de los resultados obtenidos del modelo de elementos finitos del RVE en ANSYS para el Apartado 4.1. %-------------------------------------- % Material Properties %-------------------------------------- E = 1000; % [MPa] nu = 0.3; rho = 1000; % [kg/m3] G = E/(2*(1+nu)); % [Pa] %-------------------------------------- % Resultados de MEF (ANSYS) %-------------------------------------- Result_CC1 =[ 0.34845E+02 0.31803E+02 0.19994E+02]; Result_CC2 =[ 0.31803E+02 0.31190E+02 0.18898E+02]; Result_CC3 =[ 0.19994E+02 0.18898E+02 0.12083E+03]; C11 = Result_CC1(1,1); C12 = Result_CC1(1,2); C13 = Result_CC1(1,3); C22 = Result_CC2(1,2); C23 = Result_CC2(1,3); C33 = Result_CC3(1,3); % imponemos valores unitarios: C44 = 1; C55 = 1; C66 = 1; %------------------------------------------------------- % Matiz C %------------------------------------------------------- Anexos 84 84 C=[C11 C12 C13 0 0 0; C12 C22 C23 0 0 0; C13 C23 C33 0 0 0; 0 0 0 C44 0 0; 0 0 0 0 C55 0; 0 0 0 0 0 C66]; %------------------------------------------------------- % Cálculo de la constantes ingenieriles %------------------------------------------------------- S = inv(C(1:3,1:3)); S = inv(C) Ec11 = 1/S(1,1); Ec22 = 1/S(2,2); Ec33 = 1/S(3,3); NUc12=-S(1,2)*Ec11; NUc13=-S(1,3)*Ec11; NUc23=-S(2,3)*Ec22; NUc21=-S(2,1)*Ec22; NUc31=-S(3,1)*Ec33; NUc32=-S(3,2)*Ec33; Gc23 = 1/S(4,4); Gc13 = 1/S(5,5); Gc12 = 1/S(6,6); rho_c=Ec33*rho/E; [rho_c Ec11 Ec22 Ec33 Gc12 Gc23 Gc13 NUc12 NUc21 NUc23 NUc32 NUc13 NUc31] 85 85 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos • Mat_properties_MEF_2 Rutina para calcular la matriz de rigidez y las constantes ingenieriles a partir de los resultados obtenidos del modelo de elementos finitos del RVE en ANSYS para el caso Apartado 4.2. %-------------------------------------- % Material Properties %-------------------------------------- E = 1000; % [MPa] nu = 0.3; rho = 1000; % [kg/m3] G = E/(2*(1+nu)); % [Pa] %-------------------------------------- % Resultados de MEF (ANSYS) %-------------------------------------- % Caso 1: h=l=20 Result_CC1 = [ 0.16592E+02 0.15970E+02 0.97685E+01 ]; Result_CC2 = [ 0.15970E+02 0.15670E+02 0.94919E+01 ]; Result_CC3 = [ 0.97685E+01 0.94919E+01 0.61893E+02 ]; % Caso 2: h=l=40 % Result_CC1 = [ 0.81047E+01 0.79697E+01 0.48223E+01 ]; % Result_CC2 = [ 0.79697E+01 0.78756E+01 0.47536E+01 ]; % Result_CC3 = [ 0.48223E+01 0.47536E+01 0.31330E+02 ]; C11 = Result_CC1(1,1); C12 = Result_CC1(1,2); C13 = Result_CC1(1,3); C22 = Result_CC2(1,2); C23 = Result_CC2(1,3); C33 = Result_CC3(1,3); % imponemos valores unitarios: C44 = 1; C55 = 1; C66 = 1; %------------------------------------------------------- Anexos 86 86 % Matiz C %------------------------------------------------------- C=[C11 C12 C13 0 0 0; C12 C22 C23 0 0 0; C13 C23 C33 0 0 0; 0 0 0 C44 0 0; 0 0 0 0 C55 0; 0 0 0 0 0 C66] %------------------------------------------------------- % Cálculo de la constantes ingenieriles %------------------------------------------------------- S = inv(C(1:3,1:3)); S = inv(C); Ec11 = 1/S(1,1); Ec22 = 1/S(2,2); Ec33 = 1/S(3,3); NUc12=-S(1,2)*Ec11; NUc13=-S(1,3)*Ec11; NUc23=-S(2,3)*Ec22; NUc21=-S(2,1)*Ec22; NUc31=-S(3,1)*Ec33; NUc32=-S(3,2)*Ec33; Gc23 = 1/S(4,4); Gc13 = 1/S(5,5); Gc12 = 1/S(6,6); rho_c=Ec33*rho/E; [rho_c Ec11 Ec22 Ec33 NUc12 NUc13 NUc23 NUc21 NUc31 NUc32] 87 87 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos • Graficas_definitivo Archivo para dibujar las gráficas de la influencia del ángulo en E, ρ y ν haciendo uso de las ecuaciones de Sorohan [8] para el estudio de la influencia del ángulo en el caso 3. %---------------------------------------------- % Simulación Numérica % Influencia del ángulo interno %---------------------------------------------- % --------------------------------------------- % Material Properties %---------------------------------------------- E = 1000; % [MPa] nu = 0.3; rho = 1000; % [kg/m3] G = E/(2*(1+nu)); % [Pa] %--------------------------------------------- % Geometría de un panel %--------------------------------------------- h=20; % [mm] l=10; % [mm] t=1; % espesor [mm] b=10; theta=[-45:0.33:45]*pi/180; theta1=[-45:0.33:45]*pi/180; % Constantes Ingenieriles: Formulación Sorohan 2015 Ec11=[]; Ec22=[]; Anexos 88 88 Ec33=[]; Gc12=[]; Gc13=[]; Gc23=[]; NUc12=[]; NUc13=[]; NUc23=[]; rho_c=[]; for i=1:max(size(theta)) alpha_x = h/l; %ratio of cell wall length beta_x = t/l; %slenderness ratio of the unit of cell wall alpha = pi/2 - theta(i); a1=( t/2 + l*cos(theta(i)) ); % valores de aproximadamente=1 a2=( h + l*sin(theta(i)) ); a3= b/2; % classical Honeycomb rho_c_aux = rho*beta_x*( alpha_x + 2 )/( 2*cos(theta(i))*( alpha_x + sin(theta(i)) ) ); k1 = 1/(1 + (beta_x^2)*(2.4 + 1.5*nu + (1/tan(theta(i)))^2) ); k2 = 1/(1 + (beta_x^2)*(2.4 + 1.5*nu + tan(theta(i))^2 + 2*alpha_x*(1/cos(theta(i)))^2) ); c12= k1*(1 + (beta_x^2)*(1.4 + 1.5*nu)); c21= k2*(1 + (beta_x^2)*(1.4 + 1.5*nu)); k12=(1+2*h/l)/( 1+2*h/l+(t/l)^2*(l/h*(2.4+1.5*nu)*(2+alpha_x+sin(theta(i)))+ ((alpha_x+sin(theta(i)))/alpha_x^2 ))*((alpha_x+sin(theta(i)))*tan(theta(i))^2+sin(theta(i))) ); Ec11_aux = k1*E*( (beta_x)^3 )*cos(theta(i))/( (alpha_x + sin(theta(i)))*sin(theta(i))^2 ); Ec22_aux = k2*E*( (beta_x)^3 )*( alpha_x + sin(theta(i)) )/( cos(theta(i))^3 ); Ec33_aux = E*rho_c_aux/rho; rho_c = [rho_c rho_c_aux]; 89 89 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos Ec11 = [Ec11 Ec11_aux]; Ec22 = [Ec22 Ec22_aux]; Ec33 = [Ec33 Ec33_aux]; Gc12_aux = k12*E*(beta_x^3)*( alpha_x + sin(theta(i)) )/( (alpha_x^2)*( 1 + 2*alpha_x )*cos(theta(i)) ); Gc13_aux = G*beta_x*cos(theta(i))/( alpha_x + sin(theta(i)) ); Gc23_aux = G*beta_x*( 1 + 2*sin(theta(i))^2 )/( 2*cos(theta(i))*( alpha_x + sin(theta(i)) ) ); Gc12 = [Gc12 Gc12_aux]; Gc13 = [Gc13 Gc13_aux]; Gc23 = [Gc23 Gc23_aux]; NUc12_aux = c12*( cos(theta(i))^2 )/( (alpha_x+sin(theta(i)))*sin(theta(i)) ); NUc13_aux = nu*Ec11_aux/Ec33_aux; NUc23_aux = nu*Ec22_aux/Ec33_aux; NUc12 = [NUc12 NUc12_aux]; NUc13 = [NUc13 NUc13_aux]; NUc23 = [NUc23 NUc23_aux]; end %------------------------------------------------------------- % Constantes ingenieriles MEF %------------------------------------------------------------- theta_1=[-45,-30,-15,-10,0,10,15,30,45]*pi/180; rho_m = [ 200.8357, 145.5593, 111.1358, 102.6313, 95.2578, 90.0660, 88.8581, 88.5078, 87.3019, 86.2352, 85.6685, 87.3335, 95.9164 ]; Ec11_m = [ 1.5191, 2.9826, 7.9444, 13.5514, 27.9023, 45.5562, 45.8393, 43.7159, 23.6790, 10.2915, 5.4068, 1.4428, 0.51 ]; Ec22_m = [ 3.8860, 2.5060, 1.7802, 1.6002, 1.435, 1.3139, 1.9217, 1.3288, 1.5201, 1.7951, 2.1135, 3.4964, 6.0783 ]; Ec33_m = [ 200.8357, 145.5593, 111.1358, 102.6313, 95.2578, 90.0660, 95.2363, 88.5078, 87.3019, 86.2352, 85.6685, 87.3335, 95.9164 ]; NUc12_m = [ -0.5863, -1.0288, -1.9163, -2.4703, -2.8937, -0.9484, 0.1553, 1.2058,2.6766, 2.0536, 1.4599, 0.6094, 0.2746 ]; NUc23_m = [ 0.0058, 0.0052, 0.0048, 0.0047, 0.0045, 0.0044, 0.0043, 0.0045, 0.0052, 0.0062, 0.0074, Anexos 90 90 0.012, 0.019 ]; NUc13_m = [ 0.0023, 0.0061, 0.0214, 0.0396, 0.0879, 0.1517, 0.1548, 0.1482, 0.0814, 0.0358, 0.0189, 0.005, 0.0016 ]; theta_m=[-45,-30,-15,-10,-5,-1,0,1,5,10,15,30,45]*pi/180; indices=[]; for i=1:max(size(theta_1)) indices = [indices find(theta_m==theta_1(i))]; end figure(1) plot(theta*180/pi,10*Ec11,"r",theta*180/pi,100*Ec22,"b",theta*180/pi,Ec33,"k",... theta_m(indices)*180/pi,10*Ec11_m(indices),"sqr",theta_m(indices)*180/pi,100*Ec22_m(indices)," ob",theta_m(indices)*180/pi,Ec33_m(indices),"dk"); xlabel('\theta [deg.]'); ylabel('Young´s Moduli [MPa]'); axis square legend('10 \cdot E_{1} (Teo.)',... '10^2\cdotE_{2} (Teo.)',... ' E_{3} (Teo.)',... '10 \cdot E_{1} (RVE-MEF)',... '10^2\cdotE_{2} (RVE-MEF)',... ' E_{3} (RVE-MEF)'); figure(2) plot(theta*180/pi,NUc12,"r",theta*180/pi,10*NUc13,"b",theta*180/pi,100*NUc23,"k",... theta_m(indices)*180/pi,NUc12_m(indices),"sqr",theta_m(indices)*180/pi,10*NUc13_m(indices)," ob",theta_m(indices)*180/pi,100*NUc23_m(indices),"dk"); xlabel('\theta [deg.]'); ylabel('Poisson´s Ratios'); legend('\nu_{12} (Teo.)'); axis square legend('10 \cdot \nu_{12} (Teo.)',... '10^2\cdot \nu_{13} (Teo.)',... 91 91 Propiedades Elásticas Efectivas de Estructuras de Panales Hexagonales Periódicos ' \nu_{23} (Teo.)',... '10 \cdot \nu_{12} (RVE-MEF)',... '10^2\cdot \nu_{13} (RVE-MEF)',... ' \nu_{23} (RVE-MEF)'); figure(3) plot(theta*180/pi,rho_c,"m",theta_m(indices)*180/pi,rho_m(indices),"sqm"); xlabel('\theta [deg.]'); ylabel('Density [kg / m^{3}]'); axis square legend('\rho_{c} (Teo.)','\rho_{c} (MEF)');