Modelización micromecánica de materiales compuestos: comparativa entre modelos analíticos y numéricos (MEF)

i Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado Ingeniería Aeroespacial Modelización micromecánica de materiales compuestos: comparativa entre modelos analíticos y numéricos (MEF) Autor: Antonio Castro Sánchez Tutor: Luis Rodríguez de Tembleque Solano Cotutor: Enrique García Macías Dep. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 ii iii Trabajo Fin de Grado Grado Ingeniería Aeroespacial Modelización micromecánica de materiales compuestos: comparativa entre modelos analíticos y numéricos (MEF) Autor: Antonio Castro Sánchez Tutor: Luis Rodríguez de Tembleque Solano Profesor Titular de Universidad Cotutor: Enrique García Macías Doctor Ingeniero de Caminos, Canales y Puertos Dep. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 v Trabajo Fin de Grado: Modelización micromecánica de materiales compuestos: comparativa entre modelos analíticos y numéricos (MEF): comparativa entre modelos analíticos y numéricos (MEF) Autor: Antonio Castro Sánchez Tutor: Luis Rodríguez de Tembleque Solano Cotutor: Enrique García Macías El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2018 El Secretario del Tribunal vii Agradecimientos A mi familia, por apoyarme incondicionalmente en todo momento. A mis amigos, a aquellos que se han hecho por el camino y a aquellos que nunca se han ido. Por todos esos ratos que hemos pasado juntos y los que aún están por llegar. Y a mis profesores, especialmente a aquellos que se esfuerzan día a día en la docencia para conseguir una escuela mejor, y a mi tutor por su gran ayuda en este proyecto. Antonio Castro Sánchez Sevilla, 2018 ix Índice Agradecimientos vii Índice ix Índice de Tablas xi Índice de Figuras xiii Notación xvii 1 Introducción 19 2 Modelos micromecánicos de los materiales compuestos 23 2.1 Leyes de comportamiento 23 2.1.1 Materiales monoclínicos 25 2.1.2 Materiales ortótropos 25 2.1.3 Materiales transversalmente isótropos 26 2.1.4 Materiales isótropos 26 2.2 Constantes ingenieriles 27 2.3 Modelos micromecánicos 29 2.3.1 Ley de mezclas 30 2.3.2 Ley de mezclas modificada: Modelo de Hopkins y Chamis 34 2.3.3 Modelo de Hashin y Shtrikman 35 2.3.4 Ecuaciones de Halpin-Tsai 36 3 RVE para la implementación numérica de un modelo micromecánico 39 3.1 Concepto de volumen representativo 39 3.2 Análisis local-global 40 3.3 Condiciones de contorno periódicas 41 3.4 Implementación numérica 42 4 Ejemplos de aplicaciones 45 4.1 Implementación de un modelo de fibra continua 45 4.1.1 Modelo y mallado 46 4.1.2 Obtención de la matriz de rigidez 50 4.1.3 Obtención de las constantes del material 66 4.2 Implementación de un modelo de inclusiones esféricas 71 4.2.1 Modelo y mallado 72 4.2.2 Obtención de la matriz de rigidez 74 4.2.3 Obtención de las constantes del material 83 5 Resumen y conclusiones 87 5.1 Resumen 87 5.2 Conclusiones 87 5.3 Trabajos futuros 88 Referencias 89 Anexos 91 A. Rutinas de ANSYS 91 FibraContinua.mac 91 xvi refinamientos del mallado y los modelos micromecánicos analíticos especificados. 85 Figura 5-1. Caracterización de la interfase entre fibra y matriz. 88 Figura 5-2. Ejemplos de RVE para compuestos trenzados. 88 xvii Notación 𝜎𝑖𝑗 Tensor de tensiones 𝜖𝑖𝑗 Tensos de deformaciones 𝜎𝑖 Pseudovector de tensiones 𝜖𝑖 Pseudovector de deformaciones 𝛾𝑖𝑗 Deformación tangencial en el plano ij 𝐶𝑖𝑗 Tensor de rigidez 𝑆𝑖𝑗 Tensor de flexibilidad 𝐸𝑖 Módulo de Young en la dirección i 𝜈 Coeficiente de Poisson 𝐺 Módulo de cizalladura 𝐾 Módulo de rigidez volumétrica 𝑘 Módulo de rigidez volumétrica en deformación plana 𝑓 Fibra/Refuerzo 𝑚 Matriz 𝑣𝑓 Fracción volumétrica de fibra/refuerzo 𝑎𝑖 Semilongitud del volumen representativo en dirección i 𝑑𝑓 Diámetro de fibra/Diámetro de la inclusión 𝜎𝑖  Pseudovector de tensión media 𝜖𝑖  Pseudovector de deformación media 𝐶𝑖𝑗     Tensor de rigidez media 𝑆𝑖𝑗     Tensor de flexibilidad media 𝜎𝑖 𝑝 Pseudovector de tensión periódico 𝜖𝑖 𝑝 Pseudovector de deformación periódico 𝑉 𝑒 Volumen de cada elemento 𝜎𝑖 𝑒 Pseudovector de tensión de cada elemento 𝜖𝑖 𝑒 Pseudovector de deformación de cada elemento 𝐂 𝐅𝐄𝐌 Tensor de rigidez calculado mediante elementos finitos 𝐂 𝐓𝐄𝐎 Tensor de rigidez teórico para una cierta ley de comportamiento 19 1 INTRODUCCIÓN a evolución tecnológica ha estado siempre fuertemente ligada al desarrollo de nuevos materiales cuyas propiedades sean capaces de satisfacer unos requisitos cada vez más exigentes. La idea de combinar dos tipos de materiales para conseguir uno con mejores características ha estado presente desde la antigüedad a través de las aleaciones. Estas han mostrado ser lo suficientemente resistentes y maleables para que su uso se haya extendido enormemente en todo tipo de aplicaciones, incluidas la de barcos y vehículos. No obstante, el desarrollo de la industria aeronáutica durante el siglo XX trajo consigo la necesidad de utilizar materiales que, siendo cada vez más livianos, fueran capaces de soportar cargas cada vez mayores. Es por ello que la investigación sobre materiales compuestos ha prosperado, y su uso se ha ido extendiendo cada vez más. La importancia de los materiales compuestos en la industria aeronáutica se ha ido consagrando a lo largo de los años. Los modelos más avanzados como el B787 y el A350 ya presentan más del 50% de su peso asociado a materiales compuestos. En una industria en la que el peso penaliza enormemente el consumo de combustible, y por tanto los gastos para la compañía, el uso de materiales compuestos en los nuevos diseños se ha vuelto indispensable para fabricar aviones cada vez más eficientes. Esto se ha visto favorecido por la mejora en las técnicas de fabricación de materiales compuestos, que nos permiten actualmente fabricar grandes laminados de geometría compleja a menor coste. Figura 1-1. Relación entre el año del primer vuelo y la fracción de peso de materiales compuestos para distintos aviones. L 20 Figura 1-2. Uso de materiales compuestos en el B787. La necesidad de fabricar compuestos de diversas características estructurales ha permitido que se desarrollen técnicas de diseño de compuestos. Controlando la forma en que se dispone la matriz y el refuerzo durante su fabricación, así como la cantidad de material de cada una de las fases, somos capaces de obtener compuestos de distintas propiedades mecánicas. Ejemplo de esto es el apilado de láminas en distintas direcciones, lo que nos permite conseguir un laminado con una dirección favorable de carga. Es por todo esto que la predicción de las propiedades de este tipo de materiales ha pasado a ser un tema de gran interés. Para ello se han propuesto diversas técnicas basadas en distintos modelos micromecánicos. Las primeras estimaciones de las propiedades mecánicas de compuestos fueron realizadas mediante las ecuaciones de Voigt [1] y Reuss [2], más conocidas como leyes de mezclas, y simplificaban el problema suponiendo la misma deformación o tensión en todas las fases. Más adelante se demostró que estas ecuaciones representaban los límites superior e inferior de las propiedades del compuesto, y se establecieron otros límites basados en principios energéticos (Hill [3]) y en principios variacionales (Hashin y Shtrikman [4]). Los métodos basados en establecer límites dieron paso a los modelos auto suficientes. Estos modelos suponían una microestructura concreta para resolver el problema elástico y obtener las propiedades efectivas del material. Ejemplos son el modelo de fibra continua de Hill [5] (aunque resuelto por Hermans [6]), el modelo de inclusiones esféricas de Kerner [7] y el modelo de fibra corta de Cox [8]. Con el objetivo la condensar los resultados de distintos modelos, Halpin y Tsai revisaron los resultados obtenidos hasta la fecha y demostraron que era posible recogerlos en una única formulación que dependía de un cierto parámetro de refuerzo. En [9] se hizo una revisión de estas ecuaciones. Los grandes avances en computación permitieron emplear el Método de los Elementos Finitos (MEF) para resolver numéricamente los modelos. Esto ha conducido a los más recientes modelos numéricos en los que mediante una técnica de homogenización somos capaces de obtener las propiedades del material a partir de una celda o volumen representativo. Un precursor de este método fue Aboudi con su trabajo [10]. Implementaciones más recientes de este método pueden encontrarse en [11], [12] y [13], donde se aplican en modelos de fibra continua y fibra corta tanto alineadas como distribuidas aleatoriamente. En este trabajo se pretende hacer una revisión de los modelos micromecánicos que han sido propuestos a lo largo de los años y comparar sus resultados con los obtenidos a través de modelos basados en el MEF. Esto se hará para el caso de compuestos con fibra circular, continua y alineada y en el caso de inclusiones de tipo esférico. Una vez definidos los modelos que se emplearán para la resolución del problema elástico en el software de elementos finitos ANSYS, podremos extraer y analizar los resultados obtenidos mediante MATLAB y comprobar la validez de los modelos comúnmente aceptados mediante gráficas comparativas. De esta forma, 21 dispondremos de una potente herramienta que nos permitirá estudiar compuestos con cualquier otro tipo de estructura micromecánica interna. A continuación, se detallará la estructura del documento y el contenido de cada capítulo: • Capítulo 2: Micromecánica de los materiales compuestos Aquí se pretende dar a conocer los modelos micromecánicos más importantes, haciendo especial hincapié en los modelos empleados para la comparación con los resultados obtenidos en el Capítulo 4. Comenzaremos con un somero resumen de las leyes de comportamiento de la Teoría de la Elasticidad que nos permiten caracterizar los distintos tipos de compuestos. También se explicarán las constantes ingenieriles más importantes debido al gran significado físico que guardan. Finalmente, se recogen los modelos micromecánicos con sus consideraciones y formulación matemática. • Capítulo 3: RVE para la implementación numérica de un modelo micromecánico El objetivo de este capítulo es profundizar sobre el concepto de volumen representativo (Representative Volumen Element, RVE), usado en los distintos modelos del Capítulo 2. Desarrollar la idea de los análisis local y global es fundamental para entender como somos capaces de obtener las propiedades del compuesto a partir de un volumen que representa una pequeña parte de este. En este apartado también se pretende cerrar el problema elástico dando las condiciones de contorno, que surgen de una técnica de homogenización que se basa en la periodicidad del compuesto. Para acabar, se explicará brevemente en qué consiste el MEF y cuáles son las ecuaciones que nos permiten obtener las propiedades efectivas a partir de la solución del problema elástico. • Capítulo 4: Ejemplos de aplicaciones En este capítulo se recogen los resultados numéricos obtenidos mediante el método descrito en el capítulo anterior y constituye la parte principal de este trabajo. En este apartado se hace hincapié en las propuestas de solución elegidas, refiriéndonos al software y a los códigos empleados. Comenzaremos describiendo los modelos usados, que serán un modelo de fibra continua en disposición hexagonal, y un modelo de inclusiones esféricas en disposición FCC. Para ambos modelos se mostrarán numérica y gráficamente las soluciones del problema elástico. También se muestran las gráficas comparativas y se hacen valoraciones de los resultados obtenidos. • Capítulo 5: Conclusiones Finalmente se hará un resumen del trabajo realizado y se enumerarán las conclusiones más importantes a las que se ha llegado. Igualmente se propondrán trabajos futuros íntimamente relacionados con la metodología aquí utilizada. • Anexos En los anexos se recogen el conjunto de códigos de ANSYS y MATLAB usados para la resolución del problema. 23 2 MODELOS MICROMECÁNICOS DE LOS MATERIALES COMPUESTOS E recomienda consultar la norma UNE 50136:1997: Documentación, presentación de tesis y documentos similares. Por otro lado, existen recursos en línea que ayudan a organizar todo el proceso de elaboración del Trabajo. El personal de la Biblioteca imparte formación y asesoramiento sobre su uso. Para establecer la estructura y orden de los datos y partes del Trabajo Fin de Grado nos basamos en la Normativa de los Trabajos Fin de Grado de la Escuela Técnica Superior de Ingeniería de Sevilla (Aprobada en la sesión de la Junta de Escuela de 12 de julio de 2013, modificada en Junta de Escuela de 05 de febrero de 2014) y en la norma UNE 50136:1997: Documentación, presentación de tesis y documentos similares. 2.1 Leyes de comportamiento Antes de presentar los modelos micromecánicos analíticos se procede a hacer una breve descripción de las leyes de comportamiento más importantes, haciendo especial hincapié en aquellas que describen con mayor precisión la estructura interna de los materiales compuestos. Puesto que nos estamos moviendo dentro de la teoría de la elasticidad lineal 1 podemos establecer una relación biunívoca entre tensiones y deformaciones, relación que se caracteriza de forma general a través del tensor de rigidez: ij ijkl kl C  = (2–1) A partir de ahora aprovecharemos la simetría de los tensores de tensión y deformación y se usará la formulación compacta para definir la relación anterior. En esta los tensores pasan a ser pseudovectores formados por las componentes significativas de cada uno. Por lo tanto, definimos el pseudovector de tensiones de la siguiente forma: 11 12 13 1 6 5 12 22 23 6 2 4 13 23 33 5 4 3 ij                            ==             → 1 2 3 4 5 6 i           =      (2–2) En el caso del tensor de deformaciones se tiene un cambio en la definición de las componentes tangenciales. Se definen las deformaciones tangenciales de la siguiente forma: 1 En esta teoría haremos las suposiciones de pequeñas deformaciones y desplazamientos lo que nos permitirá, sin pérdida de generalidad, obtener las constantes elásticas de materiales compuestos. S 24 2 (i j) ij ij  = (2–3) Definida así, la deformada tangencial representa el ángulo que se alejan o se aproximan las caras de un cubo sometido a este tipo de deformación. Además, permite una relación directa con la tensión tangencial en materiales ortótropos a través del módulo de cizalladura, como se verá más adelante. Con este cambio de notación, definimos el pseudovector de deformaciones de la siguiente forma: 11 12 13 1 6 5 12 22 23 6 2 4 13 23 33 5 4 3 / 2 / 2 / 2 / 2 / 2 / 2 ij                            ==             → 1 2 3 4 5 6 i           =      (2–4) Por tanto, la relación (2–1) pasaría a definirse de forma general (esto es, para un material anisótropo) mediante 21 constantes: i ij j C  = (2–5) 1 11 12 13 14 15 16 1 2 22 23 24 25 26 2 3 33 34 35 36 3 4 44 45 46 4 5 55 56 5 6 66 6 C C C C C C C C C C C C C C C C C C Sim C C C                               =                                     (2–6) Puede hacerse un desarrollo análogo para el tensor de flexibilidad ij S , lo que permite establecer la relación inversa: -1 S =[C] (2–7) i ij j S  = (2–8) 1 11 12 13 14 15 16 1 2 22 23 24 25 26 2 3 33 34 35 36 3 4 44 45 46 4 5 55 56 5 6 66 6 S S S S S S S S S S S S S S S S S S Sim S S S                               =                                     (2–9) La simetría del tensor ij C viene dada por una propiedad de la energía de deformación en el modelo considerado: al no tener en cuenta las pérdidas de energía en forma de calor es posible definir una función potencial densidad de energía de deformación en la que la energía potencial almacenada no depende del camino recorrido si no 25 de las deformaciones finales alcanzadas. Además, cabe destacar que el tensor ij C viene dado en un sistema de referencia concreto, luego este depende del sistema de referencia elegido. Si se considera un material con simetría elástica, esto es, simetría en las constantes elásticas en dos direcciones perpendiculares al plano de simetría, es posible reducir el número de constantes que caracterizan su ley de comportamiento nada más que comparando los tensores expresados en dos sistemas de referencia equivalentes. La demostración de estas propiedades y la reducción del número de constantes puede encontrarse en [14] y [15]. Atendiendo a las propiedades de simetría elástica que presenten, podemos clasificar los materiales en diversos tipos: materiales monoclínicos, materiales ortótropos y materiales isótropos. 2.1.1 Materiales monoclínicos Se entienden por materiales monoclínicos aquellos que presentan un único plano de simetría elástica. Suponiendo el plano 2-3 de simetría, las leyes de comportamiento en tensiones y deformaciones quedan: 11 12 13 14 11 12 22 23 24 22 13 23 33 34 33 14 24 34 44 23 23 55 56 13 13 56 66 12 12 00 00 00 00 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C C C C CC CC                           =                              (2–10) 11 12 13 14 11 12 22 23 24 22 13 23 33 34 33 14 24 34 44 23 23 55 56 13 13 56 66 12 12 00 00 00 00 0 0 0 0 0 0 0 0 S S S S S S S S S S S S S S S S SS SS                           =                              (2–11) Por lo tanto, la ley de comportamiento de los materiales monoclínicos queda completamente definida mediante 13 constantes. 2.1.2 Materiales ortótropos Se entienden por materiales ortótropos aquellos que presentan tres planos de simetría elástica. Puede demostrarse que, teniendo dos planos de simetría, el tercer plano perpendicular a los anteriores también es de simetría (ver [15]). Con esta configuración las leyes de comportamiento quedan: 11 12 13 11 12 22 23 22 13 23 33 33 44 23 23 55 13 13 66 12 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C                           =                              (2–12) 32 Lo que se conoce como Ley de mezclas directa o simplemente Ley de mezclas. Un procedimiento similar puede seguirse para determinar el coeficiente de Poisson longitudinal 𝜈12 del compuesto a partir del de sus constituyentes, suponiendo nuevamente la misma deformación longitudinal en matriz y fibra: 12c f f m m vv    =+ (2–25) Considérese ahora el volumen sometido a un estado de carga en dirección transversal a la fibra. Relacionando el incremento de longitud del compuesto con los incrementos en la fibra y matriz por separado se llega a la relación: c c f f m m t t t    =+ (2–26) Donde 𝑡𝑐,𝑓,𝑚 y 𝜖𝑐,𝑓,𝑚 son los valores del espesor y de la deformación en el compuesto, fibra y matriz, respectivamente. Figura 2-4. Estado de carga para determinar la ley de mezclas inversa mediante un estado de isotensión. Figura extraída de [16]. Usando nuevamente la definición del módulo de Young y suponiendo un estado de isotensión en el compuesto puede encontrarse la siguiente ecuación para determinar el módulo elástico transversal del compuesto: 22 1 1 1 fm c f m vv E E E =+ (2–27) Lo que se conoce como Ley de mezclas inversa. Expresiones similares pueden hallarse para los módulos de cizalladura longitudinal y transversal, 𝐺12 y 𝐺23: 12 12 1 1 1 fm c f m vv G G G =+ (2–28) 23 23 1 1 1 fm c f m vv G G G =+ (2–29) Puede comprobarse que los resultados experimentales para 𝐸1 y 𝜎12 en materiales de fibra continua se ajustan con suficiente precisión a las estimaciones dadas por la ley de mezclas (2–24) y (2–25). En cambio, las estimaciones para 𝐸2, 𝐺12 y 𝐺23 no son lo suficientemente precisas, por lo que se han propuesto distintos modelos para predecir estas propiedades. 33 Figura 2-5. Comparación entre los resultados experimentales y la predicción de la ley de mezclas para el módulo de Young longitudinal 𝐸1. Figura extraída de [16]. Figura 2-6. Comparación entre los resultados experimentales y la predicción de la ley de mezclas inversa para el módulo de Young transversal 𝐸2. Figura extraída de [16]. 34 Tabla 2–2. Resultados de la ley de mezclas para un compuesto de fibra con 𝐸𝑚= 3.5𝐺𝑃𝑎, 𝜈𝑚= 0.35 y 𝐸𝑓=73𝐺𝑃𝑎, 𝜈𝑓= 0.22. 𝑣𝑓 𝐸1[𝐺𝑃𝑎] 𝜈12 𝐸2[𝐺𝑃𝑎] 𝐺12[𝐺𝑃𝑎] 𝐺23[𝐺𝑃𝑎] 0.10 10.45 0.3370 3.868 1.433 1.433 0.20 17.40 0.3240 4.323 1.603 1.603 0.30 24.35 0.3110 4.899 1.818 1.818 0.40 31.30 0.2980 5.653 2.100 2.100 0.50 38.25 0.2850 6.680 2.485 2.485 0.60 45.20 0.2720 8.163 3.043 3.043 0.70 52.15 0.2525 10.49 3.924 3.924 2.3.2 Ley de mezclas modificada: Modelo de Hopkins y Chamis El Modelo de Hopkins y Chamis, también conocido como la ley de mezclas modificada, surge por la necesidad de mejorar las predicciones obtenidas por la ley de mezclas inversa en compuestos de fibra continua. El modelo se basa en el presentado para la ley de mezclas (ver Figura 2-2), aunque realizado en dos etapas. En una primera etapa se considera un volumen representativo compuesto por un bloque de material ficticio 𝑏 rodeado de otro par de bloques de matriz. El espesor de este material se determina al pasar de una fibra circular a otra cuadrada, manteniendo la misma fracción volumétrica de fibra 𝜈𝑓 en el compuesto. Conocidas las propiedades 𝐸2𝑏, 𝐺12𝑏 y 𝐺23𝑏 del material ficticio es posible aplicar la ley de mezclas inversa usando dichas propiedades y 𝜈𝑓 ′=√𝜈𝑓. Las constantes del material ficticio pueden calcularse a su vez aplicando la ley de mezclas a otro modelo equivalente en el que un bloque de fibra está rodeado por dos bloques de matriz. Un esquema del modelo puede apreciarse en la Figura 2-7. Figura 2-7. Esquema del modelo de Hopkins y Chamis. Figura extraída de [17]. 35 Las ecuaciones de Hopkins y Chamis pueden resumirse en la siguiente formulación: 1 1 ff bm PPP  −  −  =+   (2–30) ( ) 1 b f f f m P P P  = + − (2–31) Siendo 𝑃 la propiedad 𝐸2, 𝐺12 o 𝐺23 que se desee estimar. En cuanto a 𝐸1 y 𝜎12 sigue pudiendo usarse la ley de mezclas. Tabla 2–3. Resultados de la ley de mezclas modificada para un compuesto de fibra con 𝐸𝑚= 3.5𝐺𝑃𝑎, 𝜈𝑚= 0.35 y 𝐸𝑓=73𝐺𝑃𝑎, 𝜈𝑓= 0.22. 𝑣𝑓 𝐸1[𝐺𝑃𝑎] 𝜈12 𝐸2[𝐺𝑃𝑎] 𝐺12[𝐺𝑃𝑎] 𝐺23[𝐺𝑃𝑎] 0.10 10.45 0.3370 4.791 1.785 1.785 0.20 17.40 0.3240 5.837 2.178 2.178 0.30 24.35 0.3110 7.012 2.620 2.620 0.40 31.30 0.2980 8.446 3.162 3.162 0.50 38.25 0.2850 10.30 3.866 3.866 0.60 45.20 0.2720 12.86 4.840 4.840 0.70 52.15 0.2525 16.65 6.300 6.300 2.3.3 Modelo de Hashin y Shtrikman El modelo de Hashin y Shtrikman [4] se basa en el empleo de principios variacionales para establecer límites en las propiedades de materiales macroscópicamente homogéneos e isótropos. Se trata de un método general que no tiene en cuenta la microestructura del compuesto. Por lo tanto, es útil para obtener los valores máximo y mínimo de las propiedades de materiales compuestos por inclusiones de cualquier tipo. Las ecuaciones de Hashin y Shtrikman establecen los siguientes límites en el módulo de rigidez volumétrica 𝐾 y en el módulo de cizalladura 𝐺: 1 max 3 1 (1 ) 34 f ff m f f f KK K K K G   −  = + − +  −+   (2–32) 1 min 3(1 ) 1 34 f mf f m m m KK K K K G   −  − = + +  −+   (2–33) 36 1 max 6 (K 2G ) 1 (1 ) 5G (3 4 ) f f f ff m f f f f GG G G K G   −  + = + − +  −+   (2–34) 1 min 6(1 )(K 2G ) 1 5G (3 4 ) f m m mf f m m m m GG G G K G   −  −+ = + +  −+   (2–35) Donde se ha supuesto que 𝐾𝑓> 𝐾𝑚 y 𝐺𝑓> 𝐺𝑚. Puede observarse que existe una simetría en las ecuaciones para ambas fases. Esto explica que no se haya tenido en cuenta una microestructura concreta para la obtención de las ecuaciones, y que la única distinción que hagamos entre las fases 𝑓 y 𝑚 sea a través de sus propiedades. Los límites que marcan las correlaciones de Hashin y Shtrikman son más estrechos que los dados por la ley de mezclas y la ley de mezclas inversas. Además, hay que tener en cuenta que los límites se acercan cuando hay una menor diferencia entre las propiedades de ambas fases. 2.3.4 Ecuaciones de Halpin-Tsai Las ecuaciones de Halpin-Tsai tratan de condensar los resultados obtenidos mediante modelos auto consistentes, modelos numéricos y resultados experimentales en una única formulación de la forma: 1 1 f mf v P Pv   + =− (2–36) 1 f m f m P P P P   − = + (2–37) Donde 𝑃 representa la propiedad que de desea calcular y 𝜉 es un parámetro adimensional de refuerzo, cuyo valor viene dado en función del modelo considerado. El parámetro 𝜉 puede tomar como límites los valores 0 e infinito, en cuyo caso llegamos a la ley de mezclas inversa y a la ley de mezclas, respectivamente. En [9] se obtienen las expresiones de 𝜉 para distintos modelos auto consistentes de fibra continua, como el modelo de Hill (en base a los resultados de Hermans [6]) y de inclusiones esféricas, como el modelo de Kerner. En la tabla siguiente se recogen los valores de 𝜉 de estos modelos, que serán los que se usarán posteriormente para la comparación de resultados. Las ecuaciones de Halpin-Tsai consideran válida la ley de mezclas para la estimación de las propiedades 𝐸1 y 𝜈12, en cuyo caso el valor de 𝜉 sería infinito. 37 Tabla 2–4. Valores del parámetro adimensional de refuerzo 𝜉 para distintos modelos micromecánicos. Datos extraídos de [9]. Modelo Propiedad elástica 𝜉 Hill (fibra continua) 𝑘 3 1 − 𝜈𝑚− 2𝜈𝑚 2 1 + 𝜈𝑚 𝐺23 1 + 𝜈𝑚 3 − 𝜈𝑚− 4𝜈𝑚 2 𝐺13 1 Kerner (inclusiones esféricas) 𝐾 2(1 − 2𝜈𝑚) 1 + 𝜈𝑚 𝐺 7 − 5𝜈𝑚 8 − 10𝜈𝑚 Tabla 2–5. Resultados del modelo de Hill para un compuesto de fibra con 𝐸𝑚= 3.5𝐺𝑃𝑎, 𝜈𝑚= 0.35 y 𝐸𝑓=73𝐺𝑃𝑎, 𝜈𝑓= 0.22. 𝑣𝑓 𝐸1[𝐺𝑃𝑎] 𝜈12 𝐸2[𝐺𝑃𝑎] 𝐺12[𝐺𝑃𝑎] 𝐺23[𝐺𝑃𝑎] 0.10 10.45 0.3370 4.288 1.558 1.513 0.20 17.40 0.3240 5.081 1.879 1.779 0.30 24.35 0.3110 6.025 2.280 2.113 0.40 31.30 0.2980 7.221 2.798 2.547 0.50 38.25 0.2850 8.809 3.491 3.133 0.60 45.20 0.2720 11.03 4.467 3.965 0.70 52.15 0.2525 14.38 5.943 5.243 3 El módulo de rigidez volumétrica en deformación plana k viene definido en [9]. 38 Tabla 2–6. Resultados del modelo de Kerner para un compuesto de inclusiones esféricas con 𝐸𝑚= 207𝐺𝑃𝑎, 𝜈𝑚= 0.30 y 𝐸𝑓=703𝐺𝑃𝑎, 𝜈𝑓= 0.22. 𝑣𝑓 𝐸[𝐺𝑃𝑎] 𝐺[𝐺𝑃𝑎] 0.10 231.3 89.44 0.20 258.5 100.5 0.30 289.1 113.0 0.40 323.9 127.3 0.50 363.8 143.8 0.60 410.1 163.1 39 3 RVE PARA LA IMPLEMENTACIÓN NUMÉRICA DE UN MODELO MICROMECÁNICO os capítulos nuevos pueden comenzarse de una forma cómoda copiando y pegando esta página (desde el principio hasta el final) en una página en blanco. Hay cuatro “profundidades” de Títulos que representan los capítulos, las secciones, las subsecciones y los apartados. Para asegurarse de que el capítulo comienza en página impar puede introducir un salto de sección impar. 3.1 Concepto de volumen representativo En el ámbito de la micromecánica de los materiales compuestos, un volumen representativo (Representative Volume Element, RVE) consiste en un extracto o modelo del material capaz de caracterizar su microestructura de forma global. El concepto de volumen representativo aparece por primera vez con el fin de generar modelos micromecánicos que permitan obtener las propiedades del compuesto [3]. (a) (b) Figura 3-1. Ejemplos de volúmenes representativos para compuestos de fibra corta. (a) Fibras alineadas distribuidas aleatoriamente. (b) Celda unitaria para modelo periódico. Modelos propuestos en [13] y [18]. Tradicionalmente [3], los volúmenes representativos se modelaban de tal forma que presentaban unas características comunes: • El RVE debe ser representativo, en media, del compuesto, por lo que este debe ser estadística o globalmente homogéneo. • Debe contener el suficiente número de inclusiones y estas deben ser lo suficientemente pequeñas para que las propiedades obtenidas sean independientes de los valores superficiales de tensión y desplazamiento (Principio de Saint Venant). De esta forma conseguimos la independencia con las condiciones de contorno aplicadas. L 40 En nuestro análisis, y con el fin de reducir el coste computacional, supondremos un compuesto periódico, lo que nos permitirá caracterizarlo mediante una celda unitaria como la mostrada en la Figura 3-1 (b). Este modelo está justificado cuando imponemos unas condiciones de contorno periódicas que representan de la forma más fiable posible el estado del volumen representativo embebido en el compuesto. Para la obtención de resultados se utilizarán en este documento RVE hexaédricos, cuyos ejes y dimensiones característicos se muestran en la Figura 3-2. Figura 3-2. Tipología de los RVE que se ensayarán en este documento durante la obtención de resultados. En el caso compuestos de fibra, la dirección de las fibras se corresponde con el eje 1. Las propiedades del compuesto se estiman mediante técnicas de homogenización [19] que nos permiten distinguir entre dos tipos de análisis del material: un análisis local, en lo referente a la microestructura y la solución elástica exacta en las fases; y un análisis global, relacionado con a la determinación de las propiedades a partir de los valores medios de tensión y deformación. 3.2 Análisis local-global A continuación, se va a desarrollar conceptualmente el procedimiento para la estimación de las propiedades de los materiales compuestos basándonos en un volumen representativo. Estos materiales permiten ser visualizados mediante dos escalas: una escala microscópica, en la que distinguimos las fases por las que está constituido el compuesto; y una escala macroscópica, en la que el material puede considerarse homogéneo. Figura 3-3. Esquema representativo de las distintas escalas en el análisis de un material compuesto. Para un punto 𝒙 en la escala macroscópica del material podemos considerar un RVE heterogéneo en la escala 41 microscópica con distintas propiedades en cada punto 𝐲. Figura extraída de [19]. La distinción entre ambas escalas nos permite plantear dos problemas completamente independientes. Considérese en la escala macroscópica un punto 𝒙 del compuesto. Mediante un análisis local podríamos modelar, para dicho punto, un RVE representativo de la microestructura. A partir de este modelo y mediante las condiciones de contorno adecuadas, podríamos obtener la solución al problema elástico para cada punto 𝐲 en la escala microscópica. De esta forma podríamos establecer una relación entre sus variables medias del tipo (2–21) en la escala macroscópica. Conocidas las propiedades efectivas del compuesto en cada punto podríamos realizar un análisis global, que nos permitiría ensayar el compuesto como un cuerpo homogéneo. Figura 3-4. Esquema sobre la obtención de las propiedades efectivas del compuesto a partir del RVE mediante un análisis local. Figura extraída de [15]. 3.3 Condiciones de contorno periódicas Puesto que estamos considerando un compuesto periódico, será necesario especificar las condiciones de contorno periódicas que hay que imponer a la celda unitaria. Considérense las siguientes definiciones de tensión y deformación media en el RVE: 1 RVE ii RVE V dV V  = (3–1) 1 RVE ii RVE V dV V  = (3–2) A partir del concepto de análisis global-local y la teoría de homogenización podemos desglosar los campos de tensión y deformación en el RVE en dos términos: un valor medio en todo el volumen y una perturbación periódica respecto a dicho valor medio. ( ) ( ) p i i i    =+xy (3–3) ( ) ( ) p i i i    =+xy (3–4) Debido a la periodicidad del compuesto, los valores medios de las variables elásticas pueden considerarse 48 (a) (b) Figura 4-4. Mallado fino del modelo para una fracción volumétrica de fibra pequeña 𝜈𝑓= 0.10. También se muestran los mallados para el valor máximo de 𝜈𝑓 que se considerará en este análisis. Puede apreciarse que el mallado grueso debe ser refinado en aquellas zonas de proximidad entre fibras. (a) (b) Figura 4-5. Mallado grueso del modelo para una fracción volumétrica de fibra elevada 𝜈𝑓= 0.75. 49 (a) (b) Figura 4-6. Mallado medio del modelo para una fracción volumétrica de fibra elevada 𝜈𝑓= 0.75. (a) (b) Figura 4-7. Mallado fino del modelo para una fracción volumétrica de fibra elevada 𝜈𝑓= 0.75. 50 4.1.2 Obtención de la matriz de rigidez Para una fracción volumétrica de fibra de 𝜈𝑓= 0.5 se resolverán cada uno de los 6 problemas con deformación media unidad 𝜖𝑖  = 1. Se calcularán cada una de las columnas de 𝑪 𝑭𝑬𝑴 obtenidas por separado para luego comprobar que el tensor se ajusta a la ley de comportamiento transversalmente isótropa, como se deduce de la disposición hexagonal. 4.1.2.1 Primera columna: deformación longitudinal 𝝐𝟏    40.9086 4.9409 4.9390 0.0000 0.0000 0.0000 iMPa      =     1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 i      =     Puede observarse en la Figura 4-8 que cada fase está sometida a una tensión constante e igual a su módulo de Young. Esto se debe al estado de deformación longitudinal (Figura 4-14) que es aproximadamente constante en todo el volumen e igual a la unidad. Vemos por tanto que el estado de deformación se aproxima con la suficiente exactitud a un estado de isodeformación, lo que es compatible con la ley de mezclas. Debido al estado de carga las tensiones tangenciales son prácticamente nulas en todos los elementos, aunque podemos ver en la Figura 4-11 valores no nulos de la tensión tangencial 𝜎23 en el plano perpendicular a la fibra. Aunque estos valores no son esperables en un material homogéneo, aquí aparecen debido a las diferentes propiedades que presentan matriz y fibra, particularmente debido al diferente módulo de Poisson. No obstante, todas las tensiones tangenciales son nulas en media, como puede apreciarse de los resultados de 𝜎𝑖 . Figura 4-8. Tensión 𝜎1 en cada elemento para una deformación media 𝜖1  = 1. Modelo de fibra continua. 51 Figura 4-9. Tensión 𝜎2 en cada elemento para una deformación media 𝜖1  = 1. Modelo de fibra continua. Figura 4-10. Tensión 𝜎3 en cada elemento para una deformación media 𝜖1  = 1. Modelo de fibra continua. 52 Figura 4-11. Tensión 𝜎23 en cada elemento para una deformación media 𝜖1  = 1. Modelo de fibra continua. Figura 4-12. Tensión 𝜎13 en cada elemento para una deformación media 𝜖1  = 1. Modelo de fibra continua. 53 Figura 4-13. Tensión 𝜎12 en cada elemento para una deformación media 𝜖1  = 1. Modelo de fibra continua. Figura 4-14. Deformación 𝜖1 en cada elemento para una deformación media 𝜖1  = 1. Modelo de fibra continua. 54 4.1.2.2 Segunda columna: deformación transversal 𝝐𝟐 4.9409 12.2836 5.5462 0.0000 0.0000 0.0000 i      =     0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 i      =     A diferencia del estado de carga anterior, de las siguientes figuras (Figura 4-16 y Figura 4-21) no podemos deducir ni un estado de isodeformación ni de isotensión. Por lo tanto, las hipótesis realizadas en la ley de mezclas y la ley de mezclas inversa no son válidas para aproximar las propiedades del material 𝐸2 y 𝐺23 con la suficiente exactitud. Puede verse nuevamente que las tensiones tangenciales son prácticamente nulas, lo que se ajusta al estado de carga considerado. Los valores no nulos de 𝜎23 se deben, al igual que en el caso anterior, a diferencias en el coeficiente de Poisson de las fases. Por su parte, en la Figura 4-21 se aprecia que la deformación que sufre la fibra es sustancialmente menor que la que sufre la matriz. Esto se debe a que la fibra es mucho más rígida que la matriz. Figura 4-15. Tensión 𝜎1 en cada elemento para una deformación media 𝜖2  = 1. Modelo de fibra continua. 55 Figura 4-16. Tensión 𝜎2 en cada elemento para una deformación media 𝜖2  = 1. Modelo de fibra continua. Figura 4-17. Tensión 𝜎3 en cada elemento para una deformación media 𝜖2  = 1. Modelo de fibra continua. 56 Figura 4-18. Tensión 𝜎23 en cada elemento para una deformación media 𝜖2  = 1. Modelo de fibra continua. Figura 4-19. Tensión 𝜎13 en cada elemento para una deformación media 𝜖2  = 1. Modelo de fibra continua. 57 Figura 4-20. Tensión 𝜎12 en cada elemento para una deformación media 𝜖2  = 1. Modelo de fibra continua. Figura 4-21. Deformación 𝜖2 en cada elemento para una deformación media 𝜖2  = 1. Modelo de fibra continua. 4.1.2.3 Tercera columna: deformación transversal 𝝐𝟑 4.9390 5.5462 12.2746 0.0000 0.0000 0.0000 i      =     0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 i      =     64 Figura 4-32. Tensión 𝜎23 en cada elemento para una deformación media 𝛾12     = 1. Modelo de fibra continua. Figura 4-33. Tensión 𝜎13 en cada elemento para una deformación media 𝛾12     = 1. Modelo de fibra continua. 65 Figura 4-34. Tensión 𝜎12 en cada elemento para una deformación media 𝛾12     = 1. Modelo de fibra continua. Figura 4-35. Deformación 𝜖12 en cada elemento para una deformación media 𝛾12     = 1. Modelo de fibra continua. Para formar el tensor 𝑪 𝑻𝑬𝑶 se han utilizado las siguientes componentes de 𝑪 𝑭𝑬𝑴: 𝐶11, 𝐶22, 𝐶21, 𝐶32 y 𝐶55. Puede observarse que el tensor obtenido se ajusta con gran exactitud a una ley de comportamiento transversalmente isótropa. 66 40.91 4.941 4.939 0 0 0 4.941 12.284 5.546 0 0 0 4.939 5.546 12.275 0 0 0 0 0 0 3.371 0 0 0 0 0 0 3.512 0 0 0 0 0 0 3.522 FEM     =     C 40.91 4.941 4.941 0 0 0 4.941 12.284 5.546 0 0 0 4.941 5.546 12.284 0 0 0 0 0 0 3.369 0 0 0 0 0 0 3.512 0 0 0 0 0 0 3.522 TEO     =     C 4.1.3 Obtención de las constantes del material En este apartado se representarán gráficamente las propiedades efectivas del material frente a la fracción volumétrica de fibra 𝜈𝑓. Al tratarse de un modelo transversalmente isótropo el comportamiento elástico del modelo vendrá dado por 5 constantes: el módulo de Young longitudinal 𝐸1, el coeficiente de Poisson longitudinal 𝜈12, el módulo de Young transversal 𝐸2, el módulo de cizalladura longitudinal 𝐺12 y el módulo de cizalladura transversal 𝐺13. Los resultados se comparan con los modelos micromecánicos de fibra continua expuestos en la sección 2.3. Los modelos elegidos para realizar las comparaciones son la ley de mezclas, la ley de mezclas modificada o modelo de Hopkins-Chamis y la ecuación de Halpin-Tsai basada en la solución de Hermans. Puesto que se han obtenido los resultados para 3 tamaños distintos del mallado también se podrá hacer un análisis de sensibilidad. 67 4.1.3.1 Módulo de Young longitudinal 𝑬𝟏 Figura 4-36. 𝐸1 frente a 𝜈𝑓 para los resultados del FEM en fibra continua con distintos refinamientos del mallado y los modelos micromecánicos analíticos especificados. Los resultados obtenidos para el módulo de Young longitudinal se ajustan con gran exactitud a los de la ley de mezclas, tal y como había predicho la literatura. Esto es coherente con las soluciones obtenidas para la tensión y la deformación en la dirección de la fibra en el problema de la primera columna, observando que es asumible un estado de isodeformación a través de la Figura 4-14. Vemos además, tal y como se comentó en el Capítulo 2, que las predicciones de los modelos de Hopkins-Chamis y Halpin-Tsai coinciden con los de la ley de mezclas. Apenas se puede hacer distinción entre los resultados para mallado grueso y fino, aunque este último es algo más preciso. 68 4.1.3.2 Coeficiente de Poisson longitudinal 𝝂𝟏𝟐 Figura 4-37. 𝜈12 frente a 𝜈𝑓 para los resultados del FEM en fibra continua con distintos refinamientos del mallado y los modelos micromecánicos analíticos especificados. En este caso vemos que el módulo de Poisson 𝜈12 sigue la tendencia lineal que predice la ley de mezclas, aunque hay una diferencia más notable en torno al 10%. Nuevamente, los resultados de Hopkins-Chamis y Halpin-Tsai coinciden con los de la ley de mezclas. 69 4.1.3.3 Módulo de Young transversal 𝑬𝟐 Figura 4-38. 𝐸2 frente a 𝜈𝑓 para los resultados del FEM en fibra continua con distintos refinamientos del mallado y los modelos micromecánicos analíticos especificados. Como se comentó previamente, los resultados del problema de la segunda columna no permitían considerar como válido un estado de isotensión. Esto se traduce en la discrepancia entre los resultados obtenidos y la ley de mezclas inversa, especialmente para fracciones volumétricas elevadas. En cambio, se observa una gran correspondencia con las ecuaciones de Halpin-Tsai basadas en la solución de Hermans para 𝜈𝑓 pequeños. Para 𝜈𝑓 mayores vemos que los resultados se acercan al modelo de Hopkins-Chamis. Vemos que el tamaño del mallado apenas influye en la solución salvo para fracciones volumétricas elevadas en las que el mallado grueso se vuelve demasiado inexacto. 70 4.1.3.4 Módulo de cizalladura longitudinal 𝑮𝟏𝟐 Figura 4-39. 𝐺12 frente a 𝜈𝑓 para los resultados del FEM en fibra continua con distintos refinamientos del mallado y los modelos micromecánicos analíticos especificados. Se pueden hacer comentarios similares a los hechos para la comparativa de 𝐸2. Se observa que, aunque la ley de mezclas inversa no predice correctamente los resultados, estos se ajustan a las ecuaciones de Halpin-Tsai con gran exactitud. Respecto al refinamiento del mallado, las diferencias son mínimas. 71 4.1.3.5 Módulo de cizalladura transversal 𝑮𝟐𝟑 Figura 4-40. 𝐺23 frente a 𝜈𝑓 para los resultados del FEM en fibra continua con distintos refinamientos del mallado y los modelos micromecánicos analíticos especificados. La comparativa en este caso es similar a la vista para 𝐸2: para valores pequeños de 𝜈𝑓 los resultados se ajustan al modelo de Halpin-Tsai, y para valores mayores se ajustan a Hopkins-Chamis. 4.2 Implementación de un modelo de inclusiones esféricas La segunda de las aplicaciones consistirá en un modelo de inclusiones esféricas. Siguiendo el criterio de elegir distribuciones de máxima compacidad se ha elegido una disposición cúbica centrada en las caras (Face Centered Cubic, FCC) que constituye la distribución de máxima compacidad en 3D (máximo en torno al 75% del volumen de las inclusiones). El volumen representativo consistirá por tanto en un cubo en el que las inclusiones se disponen tal y como se observa en la Figura 4-41. 72 Figura 4-41. Disposición FCC (Cúbica centrada en las caras). Debido a la simetría del RVE, el modelo nos conducirá a una ley de comportamiento isótropa, quedando sus propiedades definidas mediante 2 constantes (2–16): 1 11 12 12 1 2 12 11 12 2 3 12 12 11 3 23 11 12 23 13 11 12 13 12 11 12 12 000 000 000 0 0 0 ( )/ 2 0 0 0 0 0 0 ( )/ 2 0 0 0 0 0 0 ( )/ 2 C C C C C C C C C CC CC CC                               =       −             −       −             4.2.1 Modelo y mallado A continuación, se muestran los mallados propuestos para la resolución mediante el método de los elementos finitos del modelo de inclusiones esféricas. Los modelos se han realizado de forma paramétrica, obteniendo el diámetro de las esferas 𝑑𝑒 a partir de una fracción volumétrica de inclusiones 𝜈𝑓 mediante la ecuación: 1 2 3 312 f e v a a a d  = (4–3) También se hará un análisis de la sensibilidad del mallado, considerando dos tamaños distintos: un mallado más grueso y otro ligeramente más fino. Se dejará que ANSYS realice el mallado automáticamente a través de la opción smartsizing. Para conseguir una correcta adaptación a la geometría del modelo será necesario el uso de elementos tetraédricos. Tal y como se recomienda para esta geometría, se usarán elementos cuadráticos, en concreto el tipo de elemento SOLID186 que incorpora ANSYS. El uso de elementos cuadráticos nos permitirá reducir el número de elementos en comparación con el mallado realizado para el modelo de fibra continua. A continuación, se muestran ambos mallados para distintas fracciones volumétricas de fibra. 73 (a) (b) Figura 4-42. Mallado grueso del modelo para una fracción volumétrica de fibra pequeña 𝜈𝑓= 0.10. (a) (b) Figura 4-43. Mallado fino del modelo para una fracción volumétrica de fibra pequeña 𝜈𝑓= 0.10. 80 Figura 4-53. Tensión 𝜎1 en cada elemento para una deformación media 𝛾23     = 1. Figura 4-54. Tensión 𝜎2 en cada elemento para una deformación media 𝛾23     = 1. 81 Figura 4-55. Tensión 𝜎3 en cada elemento para una deformación media 𝛾23     = 1. Figura 4-56. Tensión 𝜎23 en cada elemento para una deformación media 𝛾23     = 1. 82 Figura 4-57. Tensión 𝜎13 en cada elemento para una deformación media 𝛾23     = 1. Figura 4-58. Tensión 𝜎12 en cada elemento para una deformación media 𝛾23     = 1. 83 Figura 4-59. Deformación 𝜖23 en cada elemento para una deformación media 𝛾23     = 1. Para formar el tensor 𝑪 𝑻𝑬𝑶 se han utilizado las componentes de 𝑪 𝑭𝑬𝑴: 𝐶11 y 𝐶21. Puede observarse que el tensor obtenido se ajusta con buena exactitud a una ley de comportamiento isótropa. 442.1 161.9 161.9 0 0 0 161.9 442.5 162.0 0 0 0 161.9 162.0 442.2 0 0 0 0 0 0 146.4 0 0 0 0 0 0 146.2 0 0 0 0 0 0 146.4 FEM     =     C 442.1 161.9 161.9 0 0 0 161.9 442.1 161.9 0 0 0 161.9 161.9 442.1 0 0 0 0 0 0 140.1 0 0 0 0 0 0 140.1 0 0 0 0 0 0 140.1 TEO     =     C 4.2.3 Obtención de las constantes del material En este apartado se representarán gráficamente las propiedades efectivas del material frente a la fracción volumétrica de fibra 𝜈𝑓. Al tratarse de un modelo isótropo el comportamiento elástico del modelo vendrá dado por 2 constantes: el módulo de Young 𝐸 y el módulo de cizalladura 𝐺. Los resultados se comparan con los modelos micromecánicos para compuestos macroscópicamente isótropos y homogéneos, así como con modelos auto consistentes de inclusiones esféricas, todos ellos expuestos en la sección 2.3. Los modelos utilizados han sido el modelo de Hashin-Shtrikman, que nos da unos límites superior e inferior, y el modelo de 84 4.2.3.1 Módulo de Young 𝑬 Figura 4-60. 𝐸 frente a 𝜈𝑓 para los resultados del FEM en inclusiones esféricas con distintos refinamientos del mallado y los modelos micromecánicos analíticos especificados. En primer lugar, notamos que el mínimo establecido por el modelo de Hashin-Shtrikman y los resultados del modelo de Halpin-Tsai basado en la solución de Kerner coinciden, consideración de la que ya se percataron en la literatura. Podemos decir por tanto que las inclusiones de tipo esféricos son las menos eficientes, pues dan los valores mínimos de 𝐸 y 𝐺 para una fracción volumétrica de refuerzo dada. Podemos ver que las soluciones obtenidas se ajustan con bastante precisión al modelo auto consistente de Kerner. 85 4.2.3.2 Módulo de cizalladura 𝑮 Figura 4-61. 𝐺 frente a 𝜈𝑓 para los resultados del FEM en inclusiones esféricas con distintos refinamientos del mallado y los modelos micromecánicos analíticos especificados. El módulo de cizalladura también se ajusta considerablemente bien a las ecuaciones de Halpin-Tsai basadas en el modelo de Kerner, aunque la precisión empeora para fracciones volumétricas mayores. 87 5 RESUMEN Y CONCLUSIONES 5.1 Resumen Se han presentado y explicado los distintos tipos de modelos micromecánicos que han permitido predecir el comportamiento de los materiales compuestos durante el último siglo. Se le ha dado especial relevancia a los modelos de fibra continua e inclusiones esféricas, que son las configuraciones con las que se han trabajado en este documento. Antes de pasar a la implementación numérica de los modelos, se han desarrollado el concepto de volumen representativo y la técnica de homogenización que nos permite obtener las propiedades del compuesto a partir de estas celdas unitarias. Una vez expuesto y cerrado el problema, se ha utilizado un software de elementos finitos para su resolución. Para automatizar esta tarea se han elaborado varias rutinas de obtención y análisis de datos mediante los programas comerciales. ANSYS y MATLAB. En el documento se recogen los resultados numéricos y gráficos del problema elástico, así como las gráficas comparativas entre las propiedades obtenidas mediante estos datos y las obtenidas a través de modelos micromecánicos analíticos equivalentes. Todo esto nos ha permitido llegar a conclusiones no solo sobre la validez de los modelos obtenidos por otros autores, sino también sobre la eficacia del método de los elementos finitos para obtener las propiedades mecánicas de distintos tipos de materiales compuestos. 5.2 Conclusiones • Las condiciones de contorno (3–7) nos permiten conseguir con gran precisión un estado de deformación cuyos valores medios vienen dados. Esto se desprende de los valores obtenidos del vector 𝜖𝑖  durante la obtención de cada columna, tanto para el caso de fibra continua como el de inclusiones esféricas. • El criterio de tomar las disposiciones de máxima compacidad en ambos modelos es el más adecuado para recrear una distribución aleatoria tanto de las fibras como de las inclusiones. En el caso del modelo de fibra continua, la disposición hexagonal le otorga un comportamiento transversalmente isótropo, característico de las distribuciones aleatorias. Además, las soluciones de Hermans y Kerner parten de un modelo auto consistente en el que fibra e inclusiones esféricas vienen rodeadas de matriz, no teniendo en cuenta cómo se disponen estas en el medio; esto acerca más estos modelos a predecir las distribuciones aleatorias, y se observa su gran cercanía con los resultados aquí obtenidos para una distribución hexagonal. • La ley de mezclas es válida para estimar las propiedades 𝐸1 y 𝜈12 en compuestos de fibra continua con suficiente precisión. Esto viene justificado no solo por los resultados obtenidos en las gráficas comparativas, sino por el campo de tensiones y deformaciones de la solución del problema que surge al imponer una deformación media 𝜖1 . Aquí se observa un estado de isodeformación en el volumen, que es la hipótesis que empleaba la ley de mezclas para obtener estas dos constantes ingenieriles. • La ley de mezclas inversa no es válida para estimar las propiedades 𝐸2, 𝐺12 y 𝐺23 en compuestos de fibra continua con suficiente precisión. A diferencia de la estimación de 𝐸1 y 𝜈12, las hipótesis supuestas por la ley de mezclas inversa no se cumplen en este caso. Por lo tanto, son preferibles usar los resultados de otros modelos, como los resultados de Hermans para fibra continua. • A la hora de estimar las propiedades 𝐸2, 𝐺12 y 𝐺23, los resultados de Hermans ofrecen una buena aproximación para valores de fracción volumétrica de fibra no demasiado elevados. Para 𝜈𝑓 mayores, la ley de mezclas modificada es preferible. 88 • El modelo de Kerner se ajusta con gran precisión a los resultados obtenidos para inclusiones esféricas en disposición FCC (Cúbica Centrada en las Caras). Además, las predicciones de este modelo coinciden con el límite inferior dado por el modelo de Hashin-Shtrikman para compuestos globalmente homogéneos e isótropos, luego se puede concluir que las inclusiones esféricas es el tipo de inclusión menos eficiente en cuanto a la rigidez que puede ofrecer el material para una misma 𝜈𝑓. 5.3 Trabajos futuros Los métodos aquí empleados para estimar las propiedades mecánicas de compuestos usando el método de los elementos finitos han demostrado ser una potente y versátil herramienta. Siguiendo la idea detrás de la modelización micromecánica se proponen los siguientes trabajos para sacarle el máximo partido: • Proponer modelos micromecánicos para compuestos de fibra corta y comparar los resultados numéricos con los obtenidos analíticamente. En este caso habría que tener también en cuenta la dependencia con el factor de forma de la fibra, 𝑙 𝑑. En el límite, podría comprobarse que nos aproximamos a los resultados obtenidos para fibra continua. • Proponer modelos micromecánicos para inclusiones de tipo elipsoidal y comparar los resultados numéricos con los obtenidos a través el método de Mori-Tanaka. En este caso, los límites corresponden al caso de una inclusión esférica y al de una fibra continua. • Modelar la interfase entre matriz e inclusión para analizar su influencia en las propiedades mecánicas del material. Esto conduciría a modelos micromecánicos más realísticos. Figura 5-1. Caracterización de la interfase entre fibra y matriz. • Elaborar modelos para compuestos basados en trenzado de fibras, como los mostrados en la figura. Figura 5-2. Ejemplos de RVE para compuestos trenzados. 89 REFERENCIAS [1] W. Voigt, «Ueber die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper,» Ann. Phys., nº 274, pp. 573-587, 1889. [2] A. Reuss, «Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle,» Z. angew. Math. Mech., vol. 9, pp. 49-58, 1929. [3] R. Hill, «Elastic Properties of Reinforced Solids: Some Theoretical Principles,» J. Mech. Phys. Solids, vol. 11, pp. 357-372, 1963. [4] Z. Hashin y S. Shtrikman, «A Variational Approach to the Theory of the Elastic Behavour of Multiphase Materials,» J. Mechs. Phys. Solids, vol. 11, pp. 127-140, 1963. [5] R. Hill, J. Mech. Phys. Solids, vol. 12, p. 199, 1964. [6] J. J. Hermans, Koninkl. Nederl. Akademie Van Wetenschappen-Amsterdam, Proceedings, Series B, vol. 70, p. 1, 1967. [7] E. N. Kerner, Proc. Phys. Soc. London, vol. 69B, p. 808, 1956. [8] H. Cox, «The Elasticity and Strength of Paper and other Fibrous Materials,» British Journal of Applied Physics, vol. 3, pp. 72-79, 1952. [9] J. Halpin y J. Kardos, «The Halpin-Tsai Equations: A Review,» Materials Research Laboratory, Missouri, 1976. [10] J. Aboudi, «Micromechanical analysis of composites by the method of cells,» Applied Mechanics Reviews, vol. 42, nº 7, pp. 193-221, 1989. [11] C. L. Tucker y E. Liang, «Stiffness predictions for unidirectional short-fiber composites: Review and evaluation,» Composites Science and Technology, nº 59, pp. 655-671, 1999. [12] S. B. R. Devireddy y S. Biswas, «Effect on Fiber Geometry and Representatve Volume Element on Elastic and Thermal Properties of Unidirectional Fiber-Reinforced Composites,» Hindawi Publishing Corporation, Rourkela, 2014. [13] K. P. Babu, P. Mohite y C. S. Upadhyay, «Development of an RVE and its stiffness predictions based on mathematical homogenization theory for short fibre composites,» Indian Institute of Technology, Kanpur, 2017. 96 /SOLU ceRVE,a1,a2,a3,0,0,0,1,0,0 !Cuarta columna: deformación unitaria ejes XY SOLVE FINISH /POST1 SRECOVER EPELRECOVER *VWRITE,EPELzz0,EPELxx0,EPELyy0,EPELxy0,EPELyz0,EPELxz0,Szz0,Sxx0,Syy0,Sxy0,Syz0,Sxz0 (12F10.4,$) FINISH /SOLU ceRVE,a1,a2,a3,0,0,0,0,1,0 !Quinta columna: deformación unitaria ejes YZ SOLVE FINISH /POST1 SRECOVER EPELRECOVER *VWRITE,EPELzz0,EPELxx0,EPELyy0,EPELxy0,EPELyz0,EPELxz0,Szz0,Sxx0,Syy0,Sxy0,Syz0,Sxz0 (12F10.4,$) FINISH /SOLU ceRVE,a1,a2,a3,0,0,0,0,0,1 !Sexta columna: deformación unitaria ejes XZ SOLVE FINISH /POST1 SRECOVER EPELRECOVER *VWRITE,EPELzz0,EPELxx0,EPELyy0,EPELxy0,EPELyz0,EPELxz0,Szz0,Sxx0,Syy0,Sxy0,Syz0,Sxz0 (12F10.4) FINISH *ENDDO *CFCLOSE 97 B. Rutinas de MATLAB Aquí se muestran los códigos de MATLAB que permiten procesar los datos obtenidos mediante ANSYS, así como los que implementan los distintos modelos micromecánicos usados para la comparativa. Las rutinas ANSYS_Results_Fibra.m y ANSYS_Results_Inclusiones.m son las que leen los datos de un archivo de texto y procesan la información para obtener el tensor de rigidez, las deformaciones medias de cada problema y las propiedades efectivas del material. Por su parte, cada modelo micromecánico está implementado en una función cuyas variables de entrada son la fracción volumétrica de fibra y las propiedades de los constituyentes. ANSYS_Results_Fibra.m function [Results,vf,E1,E2,G12,nu12,nu23,G23] = ANSYS_Results_Fibra(file) A=load(file); [m,~]=size(A); Results=struct; vf=zeros(1,m); E1=zeros(1,m); E2=zeros(1,m); G12=zeros(1,m); nu12=zeros(1,m); nu23=zeros(1,m); G23=zeros(1,m); for i=1:m Results(i).vf=A(i,1); Results(i).Strain(1:6,1)=A(i,2:7); Results(i).C(1:6,1)=A(i,8:13); Results(i).Strain(1:6,2)=A(i,14:19); Results(i).C(1:6,2)=A(i,20:25); Results(i).Strain(1:6,3)=A(i,26:31); Results(i).C(1:6,3)=A(i,32:37); Results(i).Strain(1:6,4)=A(i,38:43); Results(i).C(1:6,4)=A(i,44:49); Results(i).Strain(1:6,5)=A(i,50:55); Results(i).C(1:6,5)=A(i,56:61); Results(i).Strain(1:6,6)=A(i,62:67); Results(i).C(1:6,6)=A(i,68:73); 98 Results(i).C_=zeros(6); Results(i).C_(1,1)=Results(i).C(1,1); Results(i).C_(1,2)=Results(i).C(2,1); Results(i).C_(1,3)=Results(i).C(2,1); Results(i).C_(2,1)=Results(i).C(2,1); Results(i).C_(3,1)=Results(i).C(2,1); Results(i).C_(2,2)=Results(i).C(2,2); Results(i).C_(3,3)=Results(i).C(2,2); Results(i).C_(3,2)=Results(i).C(3,2); Results(i).C_(2,3)=Results(i).C(3,2); Results(i).C_(4,4)=(Results(i).C(2,2)-Results(i).C(3,2))/2; Results(i).C_(5,5)=Results(i).C(5,5); Results(i).C_(6,6)=Results(i).C(5,5); %Results(i).C_=abs((Results(i).C_-Results(i).C)); Results(i).E1=Results(i).C(1,1)-2*Results(i).C(2,1)^2/(Results(i).C(2,2)+Results(i).C(3,2)); Results(i).nu12=Results(i).C(2,1)/(Results(i).C(2,2)+Results(i).C(3,2)); Results(i).E2=(Results(i).C(1,1)*(Results(i).C(2,2)+Results(i).C(3,2))- 2*Results(i).C(2,1)^2)*(Results(i).C(2,2)-Results(i).C(3,2))/(Results(i).C(1,1)*Results(i).C(2,2)- Results(i).C(2,1)^2); Results(i).nu23=(Results(i).C(1,1)*Results(i).C(3,2)- Results(i).C(2,1)^2)/(Results(i).C(1,1)*Results(i).C(2,2)-Results(i).C(2,1)^2); Results(i).G12=Results(i).C(5,5); vf(i)=Results(i).vf; E1(i)=Results(i).E1; E2(i)=Results(i).E2; G12(i)=Results(i).G12; nu12(i)=Results(i).nu12; nu23(i)=Results(i).nu23; end G23=E2./(2*(1+nu23));