TFM-2302 Alegre Sabatel

Vista previa del material en texto

0
75571601H_solicitud_
FDO.pdf
Equation Chapter 1 Section 1 
Trabajo Fin de Máster 
Máster en Ingeniería Industrial 
 
Análisis Micromecánico de Materiales Compuestos 
con Interfase Inclusión-Matriz Imperfecta mediante 
Elemento de Volumen Representativo 
Autor: Beatriz Alegre Sabatel 
Tutor: Luis Rodríguez de Tembleque Solano 
Dpto. Mecánica de Medios Continuos y 
Teoría de Estructuras 
Escuela Técnica Superior de Ingeniería 
Universidad de Sevilla 
 
Sevilla, 2022 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
iii 
 
Trabajo Fin de Máster 
Máster en Ingeniería Industrial 
 
 
 
 
 
 
Análisis Micromecánico de Materiales Compuestos 
con Interfase Inclusión-Matriz Imperfecta mediante 
Elemento de Volumen Representativo 
 
 
Autor: 
Beatriz Alegre Sabatel 
 
 
Tutor: 
Luis Rodríguez de Tembleque Solano 
Profesor Titular de Universidad 
 
 
 
Dpto. Mecánica de Medios Continuos y 
Teoría de Estructuras 
Escuela Técnica Superior de Ingeniería 
Universidad de Sevilla 
Sevilla, 2022 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
v 
 
Trabajo Fin de Máster: Análisis Micromecánico de Materiales Compuestos con Interfase Inclusión-Matriz 
Imperfecta mediante Elemento de Volumen Representativo 
 
 
Autor: Beatriz Alegre Sabatel 
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, 2022 
 
 
 
El Secretario del Tribunal 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
vii 
 
 
 
 
 
A mi familia 
A mis amigos 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
ix 
 
Agradecimientos 
Este proyecto va dedicado a todas aquellas personas que me han acompañado en este camino 
 
A mis padres, a mis hermanos, a mis abuelos, por quererme tanto. Por confiar en mí siempre y ayudarme 
 
A la familia que una elige, por los momentos que se recuerdan, por mostrarme lo que de verdad es 
importante. A Chávez, por ser esa amiga que todo el mundo quiere tener 
 
Al señor Abad, por su apoyo, por su confianza, por su cariño. Gracias por estar ahí siempre, en las buenas y 
en las no tan buenas. Gracias por hacerme tan feliz. Aquí va un pasito más de nuestro camino 
 
A mi tutor, Luis, por su enorme paciencia, dedicación y ayuda 
 
 Y por último, gracias a mí misma, por demostrarme una vez más que con esfuerzo todo se consigue 
 
Beatriz Alegre Sabatel 
Sevilla, 2022 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
xi 
 
Resumen 
n los últimos años, el interés por el estudio del comportamiento de los materiales compuestos se ha 
incrementado debido al auge en cuanto a sus posibles aplicaciones en infinidad de ámbitos de la 
ingeniería.. Entre ellos destacan los sectores automovilístico, naval y aeronáutico, donde cada vez se 
requieren estructuras más ligeras y eficientes, que permitan reducir así el consumo de combustible y las 
emisiones contaminantes. En todo ello, la alta resistencia específica de los materiales compuestos (como 
aquellos compuestos de carbono de matriz polimérica) permite introducir de manera progresiva elementos 
estructurales caracterizados por su ligereza y rigidez. 
Las propiedades de un material compuesto son el resultado de diversos factores: las propiedades de sus 
constituyentes, destacando el parámetro que define la fracción volumétrica del refuerzo, la distribución de las 
inclusiones o fibras, que determina la homogeneidad del material, y la geometría y orientación de éstas, que 
interfieren en la anisotropía del material. Sin embargo, un parámetro fundamental en los materiales compuestos 
es la interfase entre las fibras (o las inclusiones) y la matriz, ya que tiene un papel fundamental a la hora de 
controlar los mecanismos de fallo y su propagación, la tenacidad a fractura, así como el comportamiento y 
prestaciones de éstos. 
La caracterización y el diseño de materiales compuestos de altas prestaciones viene empleando modelos 
micromecánicos basados en esquemas analíticos que permiten determinar las propiedades de estos compuestos 
en función de la fracción en volumen y distribución de las fibras de refuerzo, además de las propiedades de fibra 
y matriz. Recientemente, algunos de estos modelos han empezado a incorporar las propiedades de la interfase, 
i.e., su rigidez – o flexibilidad –, a la hora de caracterizar o determinar las propiedades del compuesto. Sin 
embargo, dichos modelos analíticos están sujetos a restricciones en cuanto a la geometría de las inclusiones o 
refuerzos. 
En los últimos años, se vienen usando técnicas de simulación computacional basadas en modelos de celda 
unitaria o elemento de volumen representativo (RVE, por sus siglas en inglés) que permiten caracterizar y/o 
diseñar materiales compuestos en función de parámetros micromecánicos (i.e., fracción en volumen y 
distribución de las fibras de refuerzo o las propiedades de fibra y matriz). 
En este contexto se desarrolla este proyecto, en el que presenta un procedimiento computacional basado en la 
técnica de homogeneización numérica de RVE y el método de los elementos finitos (MEF) para determinar las 
propiedades mecánicas de materiales compuestos en función de parámetros micromecánicos, incluyendo la 
modelización de la interfase inclusión-matriz. De esta manera, podemos considerar las propiedades del material 
en condiciones de daño en la interfase. Esta metodología ha sido validada por comparación con modelos teóricos 
de referencia. Además, se han realizado estudios paramétricos que determinan su rango de aplicabilidad. 
Finalmente, la metodología propuesta se ha empleado en la caracterización de un compuesto de fibra de vidrio 
y matriz de epoxy, donde se puede ver la versatilidad de la metodología presentada. 
 
 
 
 
 
 
 
 
E 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
xiii 
 
Abstract 
n recent years, composite materials behaviour has become a study point of interest due to the amount of 
applications these type of materials can be used for in a several number of engineering sectors such as 
automobile, naval and aeronautical fields. The need to reduce fuel usage and pollutant emissions increases 
the requirement of lighter and more efficient structures. Composite materials are distinguished by their high 
resistance (carbon fiber polymer matrix) which allows the creation of light and stiff structures. 
Composite materials properties are determined by different factors: constituents’ properties such as the 
reinforcement volume fraction, the distribution of reinforcement that defines the material uniformity, and the 
reinforcement orientation and geometry, which influence in the anisotropy of the system. However, an essential 
parameter is the interphase that can be found between the reinforcement and the matrix, having an important 
role controlling the failure mechanisms and propagation, stress-strain behaviour to failure of the material and 
fracture toughness. 
The definition and design of composite materials for high performances has applied micromechanical models 
based on analytical procedures that allow to get their properties according to their fibervolume ratio and 
distribution apart from the fiber and matrix properties. Some of these models have introduced the interphase 
properties (i.e stiffness or flexibility) that also affects the composite properties, although their application is 
restricted to specific reinforcement geometries. 
Computational techniques based on unit cells or representative volume elements (RVE) have been used over the 
past few years to characterize composite materials according to micromechanics parameters (i.e. fiber volume 
ratio and fiber distribution, fiber and matrix properties). 
In this project, a computational procedure is developed based on the RVE numerical homogenization technique 
and the finite element method (MEF) to determine the composite materials mechanical properties according to 
micromechanics parameters, including the fiber-matrix interphase modelling. Therefore, we can consider these 
properties include the interphase damage. This methodology has been validated by the comparison with theorical 
models. In addition, parametric studies have been performed to determine the applicability range. Finally, the 
procedure proposed has been used to characterize the glass fiber reinforced epoxy, which shows us its versatility. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
I 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
xv 
 
Índice 
Agradecimientos ix 
Resumen xi 
Abstract xiii 
Índice xv 
Índice de Tablas xvii 
Índice de Figuras xix 
Notación xxiii 
1 Introducción 11 
1.1 Introducción 11 
1.2 Objetivos 15 
1.3 Organización del documento 16 
2 Ecuaciones constitutivas y modelos micromecánicos 17 
2.1 Leyes de comportamiento 17 
2.1.1 Materiales ortótropos 18 
2.1.2 Materiales transversalmente isótropos 19 
2.1.3 Materiales isótropos 19 
2.1.4 Relaciones constitutivas 21 
2.2 Modelos teóricos micromecánicos: fibra matriz e inclusión-matriz 22 
2.2.1 Modelo de Halpin-Tsai 22 
2.2.2 Modelo de Mori-Tanaka de dos fases 24 
2.2.3 Modelo de Mori-Tanaka de tres fases 26 
2.3 Modelo de Hopkins y Chamis 28 
2.4 Comparativa entre Halpin-Tsai, Mori-Tanaka y Hopkins y Chamis 31 
2.4.1 Fibra continua en matriz isótropa 31 
2.4.2 Inclusiones esféricas en matriz isótropa 35 
3 Modelo de homogeneización numérica mediante Elemento de Volumen Representativo 39 
3.1 Introducción 39 
3.2 Condiciones de contorno 41 
3.3 Implementación numérica 42 
3.4 Comparativa modelo de homogeneización con modelos analíticos 44 
3.4.1 Fibra continua en matriz isótropa 45 
3.4.2 Inclusión esférica en matriz cuadrada isótropa 48 
4 Modelización de interfase inclusión-matriz imperfecta 51 
4.1 Introducción 51 
4.2 Modelo de Mori-Tanaka para modelar interfase inclusión-matriz imperfecta 53 
4.2.1 Modelo de Mori-Tanaka con daño en la región interfacial 53 
4.2.2 Implementación 54 
4.3 Modelo basado en RVE para modelar interfase inclusión-matriz imperfecta 55 
4.3.1 Modelo de Hashin 56 
4.3.2 Implementación 57 
 
 
 
5 Resultados 59 
5.1 Validación y estudio paramétrico 59 
5.1.1 Descripción del modelo 59 
5.1.2 Análisis de sensibilidad del mallado 60 
5.1.3 Influencia del parámetro η 61 
5.1.4 Influencia de las propiedades de los materiales inclusión-matriz 70 
5.2 Estudio de la influencia del daño de la interfase fibra/matriz en las propiedades de materiales 
compuestos transversalmente isótropos 71 
6 Resumen y Conclusiones 79 
6.1 Resumen 79 
6.2 Conclusiones 79 
6.3 Líneas futuras 80 
Referencias 81 
Anexo I 83 
Rutinas de Matlab 83 
Hopkins_Chamis_ESF.m 83 
Two_Phase_MT_ESF.m 83 
Halpin_Tsai_Kerner.m 84 
Hopkins_Chamis_FIB.m 84 
Two_Phase_MT_FIB.m 85 
Halpin_Tsai_FIB.m 85 
Anexo II 87 
Rutinas de ANSYS 87 
Fibra_continua.mac 87 
Inclusiones_esfericas.mac 101 
EPELRECOVER.mac 115 
ceRVE.mac 116 
Rutinas de Matlab 129 
ANSYS_Results_fibra.m 129 
ANSYS_Results_Inclusiones.m 130 
Comparación modelos analíticos y modelos homogeneización 130 
Anexo III 137 
Rutinas de ANSYS 137 
Fastcompleta_vf0p3_eta0p005_IE 137 
FCCI_2completa_eta0p01_vf0p1_CC2.mac 151 
Rutinas de Matlab 180 
Two_Phase_intD_MT_Spherinc2.m 180 
Modelo_Hashin_interfase_ vf0p1_etaVariable.m 181 
Influencia_propiedades_materiales.m 186 
 
 
 
 
 
 
 
 
xvii 
 
ÍNDICE DE TABLAS 
 
 
Tabla 2-1. Resumen tipos de materiales y leyes de comportamiento [6] 20 
Tabla 2-2. Resumen de constantes ingenieriles para tipos de materiales y leyes de comportamiento 21 
Tabla 2-3. Valores factor de refuerzo adimensional ξ para modelos micromecánicos de Hermans y 
Kerner [13] 23 
Tabla 2-4. Propiedades materiales del modelo fibra circular alineada 31 
Tabla 2-5. Propiedades materiales del modelo inclusiones esféricas en matriz isótropa 35 
Tabla 5-1. Resultados de las propiedades efectivas para distintos mallados 60 
Tabla 5-2. Comparativa valores obtenidos a partir de la solución teórica y la solución numérica para el 
Módulo de Young transversal 𝐸2 en el modelo de fibras continuas en matriz isótropa (para RVE con 
interfase se aplica η =0.01) 74 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
xix 
 
ÍNDICE DE FIGURAS 
 
 
Figura 1-1. Esquema de material compuesto, en el que se observan las fases principales: matriz, 
refuerzo e interfase. 11 
Figura 1-2. Fotografía del B787 durante su primer vuelo. 12 
Figura 1-3. Uso de materiales compuestos en el modelo B787. 13 
Figura 1-4. Modelo BMW i3 realizado en fibra de vidrio. 13 
Figura 1-5. Raqueta de tenis actual realizada en materiales compuestos. 13 
Figura 1-6. Chaqueta de motociclismo tejida con kevlar. 14 
Figura 1-7. Fibras naturales para refuerzo de plástico. 14 
Figura 2-1. Material con tres planos de simetría [7] 18 
Figura 2-2. Tensiones normales en materiales ortótropos y deformaciones [7] 19 
Figura 2-3. Material transversalmente isótropo [7] 19 
Figura 2-4. Material isótropo [7] 21 
Figura 2-5. Modelo fibras alineadas (a) y modelo inclusiones esféricas [14] (b) 23 
Figura 2-6. Material de dos fases con fibra continua en matriz isótropa [10] 25 
Figura 2-7. Material de dos fases con inclusiones esféricas en matriz isótropa 25 
Figura 2-8. Material de tres fases con fibra continua en matriz isótropa [10] 27 
Figura 2-9. Material de tres fases con inclusiones esféricas en matriz isótropa 27 
Figura 2-10. Ejemplos de elementos y su fracción volumétrica [19] 28 
Figura 2-11. Volumen sometido a una fuerza en dirección longitudinal de la fibra [7] 29 
Figura 2-12. Volumen sometido a una fuerza en dirección transversal de la fibra [7] 29 
Figura 2-13. Esquema del modelo de Hopkins y Chamis [7] 30 
Figura 2-14. Módulo de elasticidad 𝐸1 normalizado respecto al módulo de elasticidad de la matriz 
(𝐸𝑚) en función del porcentaje de fibra. 31 
Figura 2-15. Coeficiente de Poisson 𝜐12 normalizado respecto al coeficiente de Poisson de la matriz 
(𝜐𝑚) en función del porcentaje de fibra. 32 
Figura 2-16. Módulo de elasticidad 𝐸2 normalizado respecto al módulo de elasticidad de la matriz 
(𝐸𝑚) en función del porcentaje de fibra. 33 
Figura 2-17. Módulo de cizalladura 𝐺12 normalizado respecto al módulo de cizalladura de la matriz 
(𝐺𝑚) en función del porcentaje de fibra. 34 
Figura 2-18. Módulo de elasticidad 𝐸 normalizado respecto al módulo de elasticidad de la matriz 
(𝐸𝑚) en función del porcentaje de fibra. 35 
Figura 2-19. Coeficiente de Poisson 𝜐 normalizado respecto al Coeficiente de Poisson de la matriz 
(𝜐𝑚) en función del porcentaje de fibra. 36 
Figura 2-20. Módulo de cizalladura 𝐺 normalizado respecto al módulo de cizalladura de la matriz 
(𝐺𝑚) en función del porcentaje de fibra. 37 
Figura 3-1. RVE para compuesto formado por fibras continuas [5] 39 
Figura 3-2. Selección de un RVE representativo del material [13] 39 
 
 
 
Figura 3-3. Simplificación del cálculo de las propiedades efectivas del material [21] 40 
Figura 3-4. Esquema de descomposición del análisis del cálculode las propiedades efectivas 
mediante RVE [3] 40 
Figura 3-5. Modelo de inclusión esférica con tres fases para vf = 0.1 y vf = 0.5 42 
Figura 3-6. Ejemplo de modelo de inclusiones esféricas discretizado en distinto número de elementos
 43 
Figura 3-7. Esquema implementación modelo de homogeneización numérica 43 
Figura 3-8. Ejemplo de modelo de elementos finitos de fibra continua en matriz isótropa discretizado 
en elementos 44 
Figura 3-9. Comparativa modelos analíticos y modelos mediante RVE para 𝐸1 45 
Figura 3-10. Comparativa modelos analíticos y modelos mediante RVE para 𝐸2 46 
Figura 3-11. Comparativa modelos analíticos y modelos mediante RVE para 𝜈12 47 
Figura 3-12. Comparativa modelos analíticos y modelos mediante RVE para 𝐺12 48 
Figura 3-13. Comparativa modelos analíticos y modelos mediante RVE para 𝐸 49 
Figura 3-14. Comparativa modelos analíticos y modelos mediante RVE para 𝐺 50 
Figura 4-1. Esquema interfase 51 
Figura 4-2. Representación de discontinuidad de desplazamiento en interfase 2D [25] 52 
Figura 4-3. Direcciones interfase [26] 52 
Figura 4-4. Esquema del modelo de interfase con resorte [28] 53 
Figura 4-5. Esquema de implementación del modelo Mori Tanaka-Lee con interfase imperfecta 
propuesto por Lee [27] [28] 54 
Figura 4-6. Descomposición del RVE de inclusiones esféricas en matriz cuadrada isótropa con 
interfase imperfecta definida por el parámetro 𝜂 55 
Figura 4-7. Esquema de implementación del modelo de daño de la interfase en un RVE usando el 
MEF 57 
Figura 5-1. Mallado resultante con 4 divisiones (izquierda) y 24 divisiones (derecha) 60 
Figura 5-2. Influencia del parámetro η sobre el módulo de elasticidad efectivo del compuesto para 
una fracción de volumen 𝑣𝑓 = 0.1 61 
Figura 5-3. Influencia del parámetro η sobre el módulo de elasticidad efectivo del compuesto para 
una fracción de volumen 𝑣𝑓 = 0.3 62 
Figura 5-4. Influencia del parámetro η sobre el módulo de elasticidad efectivo del compuesto para 
una fracción de volumen 𝑣𝑓 = 0.5 63 
Figura 5-5. Influencia del parámetro η en el caso de una interfase perfecta. 64 
Figura 5-6. Mallado del RVE para una fracción de volumen de fibra 𝑣𝑓 = 0.3 y parámetro η=
0.005 (a); Mallado del RVE para una fracción de volumen de fibra RVE para 𝑣𝑓 = 0.3 y η= 0.1 
(b) 65 
Figura 5-7. RVE para 𝑣𝑓 = 0.1, 0.3 𝑦 0.5 65 
Figura 5-8. Valores de módulo de Young efectivo del compuesto inclusión matriz para distintos 
valores de daño de la interfase para modelo teórico y de RVE de inclusiones esféricas en matriz 
isótropa (η= 0.005) 66 
Figura 5-9. Resultados de RVE de inclusión esférica en matriz cuadrada con interfase perfecta 
(𝐺𝑚𝛾𝑠𝑅=10 − 4) y η =0.005. (a)Desplazamiento para deformación unitaria impuesta en el eje 𝑋1; 
(b) Tensión obtenida para deformación unitaria impuesta en el eje 𝑋1. 67 
Figura 5-10. Resultados de RVE de inclusión esférica en matriz cuadrada con interfase imperfecta 
xxi 
 
(𝐺𝑚𝛾𝑠𝑅=100) y η =0.005. (a)Desplazamiento para deformación unitaria impuesta en el eje 𝑋1; (b) 
Tensión obtenida para deformación unitaria impuesta en el eje 𝑋1. 68 
Figura 5-11. Resultados de RVE de inclusión esférica en matriz cuadrada con interfase totalmente 
dañada (𝐺𝑚𝛾𝑠𝑅=104) y η =0.005. (a)Desplazamiento para deformación unitaria impuesta en el eje 
𝑋1; (b) Tensión obtenida para deformación unitaria impuesta en el eje 𝑋1. 69 
Figura 5-12. Influencia de las propiedades de los materiales inclusión-matriz 70 
Figura 5-13. RVE fibra circular continua con volumen de fibra 𝑣𝑓 = 0.3 y espesor definido por η 
=0.01 71 
Figura 5-14. Valores de módulo de Young longitudinal en RVE de fibra circular continua en matriz 
isótropa 72 
Figura 5-15. Valores de módulo de Young transversal E2 en RVE de fibra circular continua en matriz 
isótropa con un valor de parámetro η =0.01 para las fracciones volumétricas de fibra 𝑣𝑓 =
0.1, 0.3 𝑦 0.5. 73 
Figura 5-16. Resultados de RVE de fibra continua en matriz isótropa con interfase perfecta 
(𝐺𝑚𝛾𝑠𝑅=10 − 4) y η =0.01. (a)Desplazamiento para deformación unitaria impuesta en el eje 𝑋2 
Ec.(82); (b) Tensión obtenida para deformación unitaria impuesta en el eje 𝑋2 Ec.(82). 75 
Figura 5-17. Resultados de RVE de fibra continua en matriz isótropa con interfase perfecta 
(𝐺𝑚𝛾𝑠𝑅=10 − 4) y η =0.01. (a) Desplazamiento para deformación unitaria impuesta en el eje 𝑋1 
Ec.(81); (b) Tensión para deformación unitaria impuesta en el eje 𝑋1 Ec.(81). 76 
Figura 5-18. Resultados de RVE de fibra continua en matriz isótropa con interfase totalmente dañada 
(𝐺𝑚𝛾𝑠𝑅=104) y η = 0.01 (a) Desplazamiento para deformación unitaria impuesta en el eje 𝑋2 
Ec.(82); (b) Tensión para deformación unitaria impuesta en el eje 𝑋2 Ec.(82). 77 
Figura 5-19. Resultados de RVE de fibra continua en matriz isótropa con interfase totalmente dañada 
(𝐺𝑚𝛾𝑠𝑅=104) y η = 0.01 (a) Desplazamiento para deformación unitaria impuesta en el eje 𝑋1 
Ec.(81); (b) Tensión para deformación unitaria impuesta en el eje 𝑋1 Ec.(81). 78 
Figura 6-1. Compuesto piezoeléctrico [22] 80 
 
 
 
 
 
 
xxiii 
 
Notación 
𝜎𝑖𝑗 Tensor de tensiones 
𝜀𝑘𝑙 Tensor de deformaciones 
𝐶𝑖𝑗𝑘𝑙 Tensor de rigidez 
𝑆𝑖𝑗𝑘𝑙 Tensor de flexibilidad 
𝛾𝑖𝑗 Deformaciones tangenciales 
{𝜎} Pseudovector de tensiones 
{𝜀} Pseudovector de deformaciones 
[𝐶] Matriz de rigidez 
[𝑆] Matriz de flexibilidad 
𝐸𝑖 Módulo de Young en la dirección 𝑖 
𝜐𝑖𝑗 Coeficiente de Poisson 
𝐺𝑖𝑗 Módulo de cizalladura 
𝑃 Propiedad del compuesto 
𝑃𝑚 Propiedad del refuerzo 
𝑃𝑓 Propiedad de la matriz 
𝜉 Factor de refuerzo adimensional 
𝑣𝑓 Fracción volumétrica de la fibra/inclusión 
𝑣𝑚 Fracción volumétrica de la matriz 
𝑣𝑖 Fracción volumétrica de la interfase 
[𝑆𝑒] Expresión matricial del tensor de Eshelby 
[𝑆𝑒𝑓] Expresión matricial del tensor de Eshelby de la fibra 
[𝑆𝑒𝑖] Expresión matricial del tensor Eshelby de la interfase 
𝑉 Volumen 
𝑉𝑓 Volumen de fibra 
𝑉𝑚 Volumen de matriz 
𝑡 Espesor 
𝜀𝑖 Deformación 
𝜎𝑖 Tensión 
𝑣𝑏 Fracción volumétrica de material ficticio b 
𝑃𝑏 Propiedad material ficticio b 
𝐸𝑚 Módulo de Young de la matriz 
𝐸𝑓 Módulo de Young de la fibra/inclusión 
𝐸𝑖 Módulo de Young de la interfase 
𝜐𝑚 Coeficiente de Poisson de la matriz 
𝜐𝑓 Coeficiente de Poisson de la fibra/inclusión 
𝜐𝑖 Coeficiente de Poisson de la interfase 
𝐺𝑚 Módulo transversal de la matriz 
 
 
 
𝐺𝑖 Módulo de cizalladura de la interfase 
�̅� Tensiones promedio 
𝜀 ̅ Deformaciones promedio 
𝜎𝑖
𝑝
 Perturbación teórica de tensión 
𝜀𝑖
𝑝
 Perturbación teórica de deformación 
𝑢𝑖 Campo de desplazamiento 
𝑢𝑖
𝑃(𝑦) Desplazamientos asociados a las perturbaciones 
𝑎𝑖 Parámetro geométrico del RVE 
𝑅 Radio de la inclusión 
𝑡𝑖 Espesor de la interfase 
𝜂 Parámetro de la interfase 
𝑡𝑛 Tracción normal de la interfase 
𝑡𝑡 Tracción tangencial de la interfase 
𝑢𝑛 Discontinuidad en el desplazamiento normal de la interfase 
𝑢𝑡 Discontinuidad en el desplazamiento tangencial de la interfase 
𝑘𝑛 Parámetro de rigidez normal de la interfase 
𝑘𝑡 Parámetro de rigidez tangencial de la interfase 
α Flexibilidad de la interfase en dirección tangencial de la interfase 
β Flexibilidad de la interfase en dirección normal de la interfase 
𝛾𝑠 Flexibilidad de la interfase 
[𝐼] Matriz identidad 
𝜆𝑖 Parámetro de Lamé 
𝑑𝑒 Diámetro de la inclusión 
𝑋1 Eje de dirección del material 1 
𝑋2 Eje de dirección de material 2 
𝑋3 Eje de dirección de material 3 
 
 
 
 
 
 
 
 
 
 
 
1 INTRODUCCIÓN 
1.1 Introducción 
 saac M. Daniel and O. Ishai [1] define un material compuesto como aquel material constituido por dos o 
más fases en escala macroscópica que se caracteriza por poseer propiedades superiores a las que poseen sus 
constituyentes de forma aislada. De manera generalizada, una de las fases se corresponde con una fase 
discontinua llamada "refuerzo",que aporta rigidez, y una segunda fase menos rígida, que envuelve a la anterior 
y denominamos "matriz" (Figura 1-1). En algunas ocasiones, para controlar la interacción entre ambas, existe 
un tercer constituyente conocido como “interfase”. 
Por tanto, las propiedades de un material compuesto son el resultado de diversos factores: las propiedades de sus 
constituyentes, destacando el parámetro que define la fracción volumétrica del refuerzo, la distribución de éste 
que define la homogeneidad del material (a menor uniformidad, más heterogéneo se considera el material), y la 
geometría y orientación que interfieren en la anisotropía del material. 
Los constituyentes estarán determinados por el tipo de aplicación del material compuesto; para materiales con 
poco rendimiento, el refuerzo consistirá en fibras cortas que aportarán cierto endurecimiento, envueltas en una 
matriz como constituyente principal; en el caso inverso de materiales de alto rendimiento, el refuerzo presentado 
en fibras continuas conformará la estructura principal del material, acompañado de una matriz como medio de 
transición entre fibras. La interfase, que siempre será de pequeñas dimensiones, tiene un papel importante a la 
hora de controlar los mecanismos de fallo y su propagación, la tenacidad a fractura y el comportamiento tensión-
deformación del material [1]. 
 
Figura 1-1. Esquema de material compuesto, en el que se observan las fases principales: matriz, refuerzo e interfase. 
I 
 
 Introducción 
12 
 
12 
Los materiales compuestos se encuentran presentes en distintos ámbitos de la ingeniería al caracterizarse por 
tener una gran capacidad de adaptación y por tanto, mucha versatilidad. Destaca su amplio uso en la industria 
aeroespacial, automovilística, construcción naval, sector energético o deportivo. 
 En la industria aeronáutica, centrada en la reducción de costes de mantenimiento y de operación, y encaminada 
desde sus inicios a la reducción de peso y aumento de prestaciones y alcance/tiempo de vuelo, el uso de 
materiales compuestos (i.e., plásticos reforzados con fibra de carbono) se ha incrementado de forma exponencial 
comenzando con la introducción de este tipo de materiales en estructuras secundarias hasta llegar al caso del 
B787 Dreamliner (Figura 1-2) en el que representan el 50% del total de componentes del avión (Figura 1-3). 
Con el paso del tiempo se ha comprobado, además de su ligereza y consecuente reducción del consumo, la baja 
capacidad de corrosión y alta resistencia a la fatiga. 
En el sector automovilístico, el uso de la fibra de carbono o la fibra de vidrio no solo ha aportado la reducción 
de masa de los vehículos, sino también una eficacia en cuanto a la simplificación de los procesos productivos y 
los consumos energéticos requeridos (Figura 1-4). Debido a las propiedades de los materiales compuestos, como 
la alta resistencia a la oxidación, se reducen los tratamientos posteriores a la fabricación necesarios para alcanzar 
la funcionalidad de la pieza. Se minimiza así la inversión en maquinaria, alcanzando un abaratamiento 
considerable. 
Si se analizan aplicaciones recreativas, muchos deportes se han visto beneficiados por el avance de los materiales 
compuestos. Elementos como la raqueta de tenis han evolucionado, no tanto en el aspecto geométrico, como en 
el material que la conforma, siendo inicialmente de madera y evolucionando hasta las actuales raquetas de fibra 
de carbono, aportando un gran salto de calidad en cuanto a precisión, además de la ligereza ya mencionada 
(Figura 1-5). El klevar, otro material compuesto muy resistente al impacto y al calor, se utiliza en prendas de 
motociclismo (Figura 1-6). 
Este campo se encuentra en constante evolución. En los últimos años ha habido un desarrollo en el área de la 
investigación y aplicación de los materiales compuestos reforzados con fibras naturales, como juta, sisal o 
henequén, coco, lino, piña, etc (Figura 1-7). Estos materiales de refuerzo están siendo estudiados para sustituir 
parcialmente y/o en su totalidad a algunas fibras sintéticas, para algunas aplicaciones de matriz termoestable, 
principalmente en prestaciones menos severas. La importancia del empleo de las fibras naturales, para la 
sustitución de materiales sintéticos, se basa principalmente en que no son tóxicas ni perjudiciales a la salud, 
poseen baja densidad, se aprovechan como materia prima biodegradable. 
La abundante disponibilidad de fibras de materiales naturales muestra un gran abanico de posibilidades para ser 
explotadas de forma inteligente en varios sectores industriales. La utilización de estas fuentes naturales en 
aplicaciones de ingeniería apoya al desarrollo global, en los aspectos socioeconómicos, de medio ambiente y de 
tecnologías verdes [2]. 
 
 
Figura 1-2. Fotografía del B787 durante su primer vuelo. 
 
13 
 
 
Figura 1-3. Uso de materiales compuestos en el modelo B787. 
 
 
Figura 1-4. Modelo BMW i3 realizado en fibra de vidrio. 
 
Figura 1-5. Raqueta de tenis actual realizada en materiales compuestos. 
 
 Introducción 
14 
 
14 
 
 
 
Figura 1-6. Chaqueta de motociclismo tejida con kevlar. 
 
 
 
Figura 1-7. Fibras naturales para refuerzo de plástico. 
 
15 
 
1.2 Objetivos 
Como se ha comentado en el apartado anterior, las propiedades de un material compuesto son el resultado de 
diversos factores: las propiedades de sus constituyentes, destacando el parámetro que define la fracción 
volumétrica del refuerzo, la distribución de las inclusiones o fibras, que define la homogeneidad del material, y 
la geometría y orientación de éstas, que interfieren en la anisotropía del material. Sin embargo, un parámetro 
fundamental en los materiales compuestos es la interfase entre las fibras (o las inclusiones) y la matriz, ya que 
tiene un papel fundamental a la hora de controlar los mecanismos de fallo y su propagación, la tenacidad a 
fractura, así como el comportamiento y prestaciones de éstos. 
La caracterización y el diseño de materiales compuestos de altas prestaciones viene empleando modelos 
micromecánicos basados en esquemas analíticos que permiten determinar las propiedades de éstos en función 
de la fracción en volumen y distribución de las fibras de refuerzo o las propiedades de fibra y matriz. Ente los 
mencionados modelos analíticos destacan los modelos de Mori-Tanaka [3], Halpin-Tsai [4] y Hopkins y Chamis 
[5]. Recientemente, algunos modelos de estos modelos han empezado a incorporar las propiedades de la 
interfase, i.e., su rigidez – o flexibilidad –, a la hora de caracterizar o determinar las propiedades del compuesto. 
Sin embargo, dichos modelos analíticos están sujetos a restricciones en cuanto a la geometría de las inclusiones 
o refuerzos (ver modelo basado en Mori-Tanaka desarrollado por Lee et al. [7,8]. 
En los últimos años, se vienen usando técnicas de simulación computacional basadas la técnica numérica de 
celda unitaria [3,4,5] o elemento de volumen representativo (RVE, por sus siglas en inglés), que permiten 
caracterizar y/o diseñar materiales compuestos en fusión de parámetros micromecánicos (i.e., fracción en 
volumen y distribución de las fibras de refuerzo o las propiedades de fibra y matriz). 
En este contexto se desarrolla este proyecto, que tiene como objetivo presentar un procedimiento computacional 
basado en el RVE y el método de los elementos finitos (MEF), que permita determinar las propiedades 
mecánicas de materiales compuestos en función de parámetros micromecánicos, incluyendo la modelización de 
la interfase inclusión-matriz. De esta manera, podremos considerar las propiedades del material en condiciones 
de daño en la interfase refuerzo matriz. 
Para conseguir este propósito se plantearán los siguientes objetivos: 
1.- Se describirán e implementarán en MATLAB distintos modelos micromecánicos analíticos que 
permitan validar, más adelante, el modelo RVE-MEF desarrollado.Para ello se compararán las 
propiedades efectivas de varios materiales compuestos. 
2.- Desarrollo de un modelo de homogeneización numérica mediante elemento de volumen 
representativo, usando el MEF, que además incorpore la modelización de interfase inclusión-matriz 
imperfecta. Dicho modelo consiste en una interfase de espesor finito (y muy reducido), que tiene en 
cuanta la posible existencia de discontinuidad en el campo de desplazamientos en la interfase fibra-
matriz, debido a un posible daño en la misma. Por tanto, en función de la rigidez – o flexibilidad – de 
dicha interfase podremos caracterizar las propiedades del material sin dañar en condiciones de daño 
total o parcial. 
3.- Validación de esta metodología mediante comparación con modelos teóricos de referencia. 
4.- Realizar estudios paramétricos que determinan su rango de aplicabilidad. Finalmente, la metodología 
propuesta se ha empleado en la caracterización de un compuesto de fibra de vidrio y matriz de epoxy, 
donde se puede ver la versatilidad de la metodología presentada. 
 
 
 Introducción 
16 
 
16 
1.3 Organización del documento 
El presente proyecto está dividido en diferentes capítulos, cuyos contenidos se describen a continuación: 
▪ Capítulo 2 - Ecuaciones constitutivas y modelos micromecánicos: en este capítulo se introducirán 
conceptos como el tensor de rigidez y tensor de flexibilidad a partir de los cuales se obtienen las 
ecuaciones constitutivas de las leyes de comportamiento que definen los tipos de materiales atendiendo 
a las propiedades de simetría elástica. Se describirán distintos modelos teóricos micromecánicos y 
compararán las propiedades efectivas obtenidas a partir de ellos para dos tipos de materiales 
compuestos. 
▪ Capítulo 3 - Modelo de homogeneización numérica mediante elemento de volumen representativo 
(RVE): en este tercer apartado se define el Elemento de Volumen Representativo (RVE) y las técnicas 
de homogeneización utilizadas para el cálculo de las propiedades efectivas mediante el Método de 
Elementos Finitos (MEF). Posteriormente se procede a comparar y validar los resultados obtenidos a 
partir del RVE con los modelos analíticos previos para el mismo caso de materiales compuestos. 
▪ Capítulo 4 – Modelización de interfase inclusión-matriz imperfecta: este cuarto capítulo presenta la 
interfase imperfecta entre la inclusión y la matriz y se detalla la implementación de los modelos 
analíticos y numéricos que incluyen la modelización del daño de este tercer constituyente. 
▪ Capítulo 5 – Resultados: en este apartado se desarrolla el modelo de homogeneización numérica que 
incorpora la interfase inclusión-matriz imperfecta, y el modelo teórico de referencia descrito en el 
capítulo previo para su validación y estudio paramétrico para determinar su rango de aplicabilidad. De 
forma adicional, se implementa esta metodología en un compuesto de fibra de vidrio y matrix epoxy 
(material transversalmente isótropo) comprobando así su versatilidad. 
▪ Capítulo 6 – Resumen y conclusiones: en este último capítulo, se obtienen una serie de conclusiones 
además de proponer líneas futuras de estudio a partir de este trabajo. 
 
: 
 
17 
 
2 ECUACIONES CONSTITUTIVAS Y MODELOS 
MICROMECÁNICOS 
 n este capítulo se van a presentar las ecuaciones constitutivas de las leyes de comportamiento, 
introduciendo conceptos como el tensor de rigidez y el tensor de flexibilidad. Se definirán los tipos de 
materiales atendiendo a las propiedades de simetría elástica y las respectivas leyes de comportamiento. 
Por último, se describirán modelos teóricos micromecánicos por distintos autores particularizando para dos 
geometrías diferentes: fibra-matriz e inclusión-matriz. 
2.1 Leyes de comportamiento 
Las leyes de comportamiento de los materiales compuestos se definen a partir del tensor de rigidez y el tensor 
de flexibilidad. Al estar dentro del régimen elástico-lineal, se puede hablar de la relación entre las tensiones 𝜎𝑖𝑗 y 
las deformaciones 𝜀𝑘𝑙 a través del tensor de rigidez 𝐶𝑖𝑗𝑘𝑙: 
 𝜎𝑖𝑗 = 𝐶𝑖𝑗𝑘𝑙 𝜀𝑘𝑙 (1) 
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: 
 𝛾𝑖𝑗 = 2𝜀𝑖𝑗(𝑖 ≠ 𝑗) (2) 
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. Con este cambio de notación, definimos el pseudovector de deformaciones 
de la siguiente forma [6]: 
 
𝜀𝑖𝑗 = [
𝜀11 𝜀12 𝜀13
𝜀12 𝜀22 𝜀23
𝜀13 𝜀23 𝜀33
] = [
𝜀1 𝜀6/2 𝜀5/2
𝜀6/2 𝜀2 𝜀4/2
𝜀5/2 𝜀4/2 𝜀3
] → 𝜀𝑖 = 
[
 
 
 
 
 
𝑙
𝜀1𝑙
𝜀2𝑙
𝜀3𝑙
𝜀4𝑙
𝜀5𝑙
𝜀6𝑙]
 
 
 
 
 
 (3) 
Como consecuencia de la simetría de los tensores de deformación y tensión, se usará la formulación compacta 
quedando la relación definida como: 
 {𝜎} = [𝐶]{𝜀} (4) 
donde [𝐶] es la matriz de rigidez (ó matriz de comportamiento) del material y, {𝜎} y {𝜀} son los pseudovectores 
tensión y deformación respectivamente, que incluyen las componentes 𝜎𝑖𝑗 y 𝜀𝑘𝑙. 
Desarrollando la anterior expresión, redefinimos la relación entre las tensiones y deformaciones con la siguiente 
ecuación: 
 
[
 
 
 
 
 
𝑙𝜎1𝑙
𝑙𝜎2𝑙
𝑙𝜎3𝑙
𝑙𝜎4𝑙
𝑙𝜎5𝑙
𝑙𝜎6𝑙]
 
 
 
 
 
= 
[
 
 
 
 
 
𝐶11 𝐶12 𝐶13
𝐶11 𝐶22 𝐶23
𝐶11 𝐶11 𝐶33
 
𝐶14 𝐶15 𝐶16
𝐶24 𝐶25 𝐶26
𝐶34 𝐶35 𝐶36
𝐶11 𝐶11 𝐶11
𝐶11 𝑆𝑖𝑚 𝐶11
𝐶11 𝐶11 𝐶11
 
𝐶44 𝐶45 𝐶46
𝐶11 𝐶55 𝐶56
𝐶11 𝐶11 𝐶66]
 
 
 
 
 
 
[
 
 
 
 
 
𝑙
𝜀1𝑙
𝜀2𝑙
𝜀3𝑙
𝜀4𝑙
𝜀5𝑙
𝜀6𝑙]
 
 
 
 
 
 (5) 
E 
 
 Ecuaciones constitutivas y modelos micromecánicos 
18 
 
18 
De forma análoga, puede realizarse el mismo razonamiento para el tensor de flexibilidad [𝑆], lo que permite 
establecer la siguiente relación inversa (ver Ref. [7]): 
 [𝑆] = [𝐶]−1 (6) 
 {𝜀} = [𝑆]{𝜎} (7) 
Por tanto, las matrices [𝐶] y [𝑆] constituyen, respectivamente, las matrices de rigidez y flexibilidad del material. 
Un material con simetría elástica será aquel cuyas constantes elásticas presenten simetría en dos direcciones 
perpendiculares al plano de simetría [6][8]. En ese caso, la ley de comportamiento que caracteriza dicho material 
podrá definirse con un menor número de constantes [9]. Así pues, conforme a las propiedades de simetría 
elástica, pueden definirse diversos tipos de materiales que estarán caracterizados por una ley de comportamiento 
concreta. 
2.1.1 Materiales ortótropos 
Un material ortótropo es aquel que posee tres planos de simetría coincidentes con los planos de coordenadas con 
respecto a la alineación de las fibras. Se puede demostrar que si se tienen dos planos ortogonales simétricos, 
existe un tercer plano ortogonal de simetría [9]. 
 
Figura 2-1. Material con tres planos de simetría [7] 
 
Debido a esta particularidad, la ley de comportamiento de los materiales ortótropos queda definida únicamente 
a partir de 9 constantes. 
 
 
(8) 
Como se puede observar en la ecuación anterior, las tensiones normales no producen deformaciones tangenciales 
siempre y cuando sean aplicadas en las direcciones ortotrópicas [7]. 
19 
 
 
Figura 2-2. Tensiones normales en materiales ortótropos y deformaciones [7] 
 
2.1.2 Materiales transversalmente isótropos 
Se definen los materiales transversalmente isótropos como aquellos que contienen tres planos de simetría (al 
igual que los materiales ortótropos), pero además en uno de los planos se considera un material isótropo. En la 
siguiente imagen el plano de isotropía coincide con el plano 𝑥2-𝑥3: 
 
Figura 2-3. Material transversalmente isótropo [7] 
 
Como consecuencia, la ley de comportamiento de este tipo de material queda definida con 5 constantes [6]: 
 
 
(9) 
2.1.3 Materiales isótropos 
En los materiales isótropos cada uno de los planos es un plano de simetría, por lo que las propiedades son las 
mismasen todas las direcciones (ver Figura 2-4). 
 
 
 
 Ecuaciones constitutivas y modelos micromecánicos 
20 
 
20 
La ley de comportamiento de este tipo de material queda completamente definida con 2 constantes como se 
muestra en la siguiente ecuación: 
 
 
(10) 
A continuación, se muestra una tabla resumen que recoge lo explicado en apartados anteriores: 
 Ortótropo Transversalmente isótropo Isótropo 
Constantes 
ingenieriles 
independientes 
 
Constantes 
ingenieriles 
dependientes 
 
 
Tensor de 
flexibilidad 
 
 
Tabla 2-1. Resumen tipos de materiales y leyes de comportamiento [6] 
 
 
21 
 
 
Figura 2-4. Material isótropo [7] 
 
2.1.4 Relaciones constitutivas 
A modo de síntesis, se muestra la Tabla 2-1, donde se incluyen las expresiones de la matriz de flexibilidad para 
cada una de las leyes de comportamiento analizadas y la relación de sus términos con las constantes ingenieriles. 
De forma adicional, las expresiones que permiten obtener las constantes ingenieriles a partir de los términos de 
la matriz de rigidez [𝐶] se muestran en la Tabla 2-2. 
 Ortótropo Transversalmente isótropo Isótropo 
Constantes 
ingenieriles a 
partir de 
matriz de 
rigidez 
 
 
Tabla 2-2. Resumen de constantes ingenieriles para tipos de materiales y leyes de comportamiento 
 
 
 Ecuaciones constitutivas y modelos micromecánicos 
22 
 
22 
2.2 Modelos teóricos micromecánicos: fibra matriz e inclusión-matriz 
Los modelos micromecánicos permiten obtener las propiedades homogéneas de los materiales compuestos 
constituidos por microestructuras heterogéneas que poseen distintas propiedades y leyes de comportamiento. 
Estas propiedades homogéneas del compuesto son llamadas propiedades efectivas [6]. Los modelos 
micromecánicos difieren unos de otros en las técnicas de caracterización de las propiedades efectivas 
considerando diferentes parámetros influyentes como la estructura de la fibra, orientación, ondulación, interfase 
fibra-matriz, etc [10]. A lo largo de los siguientes apartados, se presentan varios modelos micromecánicos 
particularizados para los casos de fibra continua-matriz e inclusión-matriz. 
2.2.1 Modelo de Halpin-Tsai 
El modelo de Halpin-Tsai aparece como un modelo semiempírico, siendo un punto intermedio entre técnicas 
más complejas, como los métodos numéricos, y aproximaciones de resistencia de materiales [10]. Está basado 
en los trabajos desarrollados por Hill y posteriormente Hermans, que se fundamentan en un modelo 
autoconsistente generalizado (“self -consistent micromechanics method”) [11]. Las ecuaciones de Halpin-Tsai 
se caracterizan por su versatilidad pues permiten su utilización para una gran variedad de geometrías de fibra, 
porcentajes volumétricos y propiedades elásticas [10][11]. El modelo de Halpin-Tsai se expresa a partir de la 
siguiente formulación compacta: 
 𝑃
𝑃𝑚
=
1 + 𝜉𝜂𝑣𝑓
1 − 𝜂𝑣𝑓
 (11) 
 
𝜂 =
(
𝑃𝑓
𝑃𝑚
) − 1
(
𝑃𝑓
𝑃𝑚
) + 𝜉
 (12) 
donde: 
▪ 𝑃: propiedad del compuesto. 
▪ 𝑃𝑚: propiedad del refuerzo. 
▪ 𝑃𝑓: propiedad de la matriz. 
▪ 𝜉: factor de refuerzo adimensional 
▪ 𝑣𝑓: fracción volumétrica de la fibra. 
Aplicando las ecuaciones para las propiedades 𝐸1, 𝐸2, 𝐺12, 𝐺23 y 𝑣23, se obtienen las expresiones que se 
muestran a continuación: 
 𝐸1 = 𝐸𝑓v𝑓 + 𝐸𝑚v𝑚 (13) 
 𝐸2
𝐸𝑚
=
1 + 𝜉𝜂v𝑓
1 − 𝜂v𝑓
 (15) 
 𝜐12 = 𝜐𝑓v𝑓 + 𝜐𝑚v𝑚 (14) 
 𝐺12
𝐺𝑚
=
1 + 𝜉𝜂v𝑓
1 − 𝜂v𝑓
 (16) 
 
𝜂 =
(
𝐺𝑓
𝐺𝑚
) − 1
(
𝐺𝑓
𝐺𝑚
) + 𝜉
 (17) 
23 
 
En las ecuaciones anteriores, aparece el parámetro adimensional 𝜉 cuyo valor depende de la geometría de la 
fibra, empaquetamiento, condiciones de carga y la propiedad a calcular. En la siguiente tabla se muestran las 
estimaciones atendiendo a los resultados obtenidos por los modelos de Hermans para fibra continua y Kerner 
para inclusiones esféricas y las distintas propiedades [12]: 
 𝑷 𝑷𝒇 𝑷𝒎 𝝃 
Hermans 
 
𝑘23 𝑘𝑓 𝑘𝑚 
1 − 𝑣𝑚 − 𝑣𝑚
2
1 + 𝑣𝑚
 
𝐺23 𝐺𝑓 𝐺𝑚 
1 + 𝑣𝑚
3 − 𝑣𝑚 − 4𝑣𝑚
2 
𝐺12 𝐺𝑓 𝐺𝑚 1 
Kerner 
𝐾 𝐾𝑓 𝐾𝑚 
2(1 − 2𝑣𝑚)
1 + 𝑣𝑚
 
𝐺 𝐺𝑓 𝐺𝑚 
7 − 5𝑣𝑚
8 − 10𝑣𝑚
 
Tabla 2-3. Valores factor de refuerzo adimensional ξ para modelos micromecánicos de Hermans y Kerner [13] 
 
El modelo desarrollado por Hermans aplica a fibras continuas (ver Figura 2-5a), mientras que el modelo de 
Kerner se estudia para aquellos materiales compuestos con inclusiones esféricas (ver Figura 2-5b). 
 
Figura 2-5. Modelo fibras alineadas (a) y modelo inclusiones esféricas [14] (b) 
(a) 
(b) 
 
 Ecuaciones constitutivas y modelos micromecánicos 
24 
 
24 
2.2.2 Modelo de Mori-Tanaka de dos fases 
En este apartado se presenta el modelo de Mori-Tanaka, definido como un modelo tipo Eshelby [15] en el que 
se contempla la interacción con el resto de fibras [10][16]. Este método es uno de los más efectivos a la hora de 
calcular las propiedades efectivas ortótropas del material [10]. 
Para aquellos casos en los que el material esté constituido por dos fases (fibra y matriz), el cálculo de la matriz 
de comportamiento efectiva se realiza a partir de las siguientes expresiones [10]: 
 [𝐶] = (𝑣𝑚[𝐶𝑚] + 𝑣𝑓[𝐶𝑓][𝐴])[𝐴1] (18) 
donde 
 [𝐴1] = (𝑣𝑚[𝐼] + 𝑣𝑓[𝐴])
−1
 (19) 
 
[𝐴] = [[𝐼] + [𝑆𝑒][𝐶𝑚]
−1([𝐶𝑓] − [𝐶𝑚])]
−1
 (20) 
siendo: 
▪ 𝑣𝑓 tanto por ciento de volumen de fibra 
▪ 𝑣𝑚 tanto por ciento de volumen de la matriz 
▪ [𝑆𝑒] es el tensor de Eshelby expresado en forma matricial y cuyas componentes se determinarán a partir 
de las propiedades de la matriz y de la forma de la fibra. 
 
El tensor de Eshelby [𝑆𝑒] queda expresado en forma matricial de la siguiente manera: 
 
 
(21) 
donde sus términos toman diversos valores en función de las propiedades de la matriz y la forma de las fibras. 
 
Si se tiene una fibra continua en una matriz isótropa (ver Figura 2-6) los valores son [17]: 
 
𝑆1111 = 0, 𝑆2222 = 𝑆3333 =
5 − 4𝑣𝑚
8(1 − 𝑣𝑚)
, 𝑆2211 = 𝑆3311 =
𝑣𝑚
2(1 − 𝑣𝑚)
, 
 𝑆2233 = 𝑆3322 =
4𝑣𝑚 − 1
8(1 − 𝑣𝑚)
, 𝑆1122 = 𝑆1133 = 0, 𝑆1313 = 𝑆1212 = 0.5, 
 𝑆2323 = 2
3 − 4𝑣𝑚
8(1 − 𝑣𝑚)
 
(22) 
 
25 
 
 
Figura 2-6. Material de dos fases con fibra continua en matriz isótropa [10] 
 
Para el caso concreto de un material formado por dos fases (matriz+inclusión), en el que las inclusiones tienen 
forma esférica y la matriz es isótropa, las componentes del tensor de Eshelby [𝑆𝑒] tienen los siguientes valores 
[17]: 
 
𝑆1111 = 𝑆2222 = 𝑆3333 =
7 − 5𝑣𝑚
15(1 − 𝑣𝑚)
, 
𝑆1122 = 𝑆2233 = 𝑆3311 =
5𝑣𝑚 − 1
15(1 − 𝑣𝑚)
, 
𝑆1212 = 𝑆2323 = 𝑆3131 = 2
4 − 5𝑣𝑚
15(1 − 𝑣𝑚)
, 
(23) 
Una vez obtenidas las componentes de la matriz de comportamiento efectiva, se pueden calcular las propiedades 
del compuesto a partir de las relaciones constitutivas indicadas en la Tabla 2-2. 
 
Figura 2-7. Material de dos fases con inclusiones esféricas en matriz isótropa 
 
 
 Ecuaciones constitutivas y modelos micromecánicos 
26 
 
26 
2.2.3 Modelo de Mori-Tanaka de tres fases 
El modelo de Mori-Tanaka, además de aplicarse para el caso de materiales formados por dos fases, puede 
utilizarse para calcular las propiedades efectivas de materiales constituidos por tres fases (fibra/inclusión, matriz 
e interfase). En esos casos, el cálculo de la matriz efectiva se efectúa a partir de las siguientes expresiones: 
 [𝐶] = (𝑣𝑚[𝐶𝑚][𝐼] + 𝑣𝑓[𝐶𝑓][𝐴𝑓] + 𝑣𝑖[𝐶𝑖][𝐴𝑖])[𝐴1] (24) 
donde 
 [𝐴1] = (𝑣𝑚[𝐼] + 𝑣𝑓[𝐴𝑓] + 𝑣𝑖[𝐴𝑖])
−1
 (25) 
 
[𝐴𝑓] = [[𝐼] + [𝑆𝑓][𝐶𝑚]
−1([𝐶𝑓] − [𝐶𝑚])]
−1
 (26) 
 [𝐴𝑖] = [[𝐼] + [𝑆𝑖][𝐶𝑚]
−1([𝐶𝑖] − [𝐶𝑚])]
−1
 (27) 
siendo 
▪ 𝑣𝑓 fracción volumétrica de la fibra 
▪ 𝑣𝑚 fracción volumétrica de la matriz 
▪ 𝑣𝑖 fracción volumétrica de la interfase 
▪ [𝑆𝑒𝑓] tensor de Eshelby de la fibra 
▪ [𝑆𝑒𝑖] tensor de Eshelby de la interfase 
Al igual que en apartados anteriores,el valor de las componentes de los tensores de Eshelby [𝑆𝑒𝑓] y [𝑆𝑒𝑖] está 
supeditado a la geometría de la fibra e interfase y propiedades de la matriz. Para el caso concreto de una fibra 
continua en una matriz isótropa (ver Figura 2-8), los elementos del tensor [𝑆𝑒𝑓] se calcularán a partir de las 
expresiones [18]: 
 
 
(28) 
 
𝑆1111 = 0, 𝑆2222 = 𝑆3333 =
5 − 4𝑣𝑖
8(1 − 𝑣𝑖)
, 𝑆2211 = 𝑆3311 =
𝑣𝑖
2(1 − 𝑣𝑖)
, 
 𝑆2233 = 𝑆3322 =
4𝑣𝑖 − 1
8(1 − 𝑣𝑖)
, 𝑆1122 = 𝑆1133 = 0, 𝑆1313 = 𝑆1212 = 0.50, 
 𝑆2323 = 2
3 − 4𝑣𝑖
8(1 − 𝑣𝑖)
 
(29) 
 
 
27 
 
De la misma manera, para calcular las componentes del tensor [𝑆𝑒𝑖] se utilizan las siguientes ecuaciones: 
 
𝑆1111 = 0, 𝑆2222 = 𝑆3333 =
5 − 4𝑣𝑚
8(1 − 𝑣𝑚)
, 𝑆2211 = 𝑆3311 =
𝑣𝑚
2(1 − 𝑣𝑚)
, 
 𝑆2233 = 𝑆3322 =
4𝑣𝑚 − 1
8(1 − 𝑣𝑚)
, 𝑆1122 = 𝑆1133 = 0, 𝑆1313 = 𝑆1212 = 0.5, 
 𝑆2323 = 2
3 − 4𝑣𝑚
8(1 − 𝑣𝑚)
 
(30) 
 
Figura 2-8. Material de tres fases con fibra continua en matriz isótropa [10] 
 
Si el material está constituido por tres fases y una disposición geométrica de inclusiones esféricas en una matriz 
isótropa (ver Figura 2-9), los elementos de los tensores de Eshelby [𝑆𝑒𝑓] y [𝑆𝑒𝑖] se obtendrán aplicando las 
expresiones que se muestran a continuación: 
[𝑺𝒆𝒇] [𝑺𝒆𝒊] 
𝑆1111 = 𝑆2222 = 𝑆3333 =
7 − 5𝑣𝑖
15(1 − 𝑣𝑖)
, 
𝑆1122 = 𝑆2233 = 𝑆3311 =
5𝑣𝑖 − 1
15(1 − 𝑣𝑖)
, 
𝑆1212 = 𝑆2323 = 𝑆3131 = 2
4 − 5𝑣𝑖
15(1 − 𝑣𝑖)
, 
𝑆1111 = 𝑆2222 = 𝑆3333 =
7 − 5𝑣𝑚
15(1 − 𝑣𝑚)
, 
𝑆1122 = 𝑆2233 = 𝑆3311 =
5𝑣𝑚 − 1
15(1 − 𝑣𝑚)
, 
𝑆1212 = 𝑆2323 = 𝑆3131 = 2
4 − 5𝑣𝑚
15(1 − 𝑣𝑚)
, 
(31) 
 
 
Figura 2-9. Material de tres fases con inclusiones esféricas en matriz isótropa 
 
 Ecuaciones constitutivas y modelos micromecánicos 
28 
 
28 
2.3 Modelo de Hopkins y Chamis 
En este último apartado se va a describir el modelo micromecánico de Hopkins y Chamis, utilizado para obtener 
las propiedades mecánicas de materiales compuestos de fibra continua y también conocido como Ley de Mezclas 
Modificada. 
Para poder comprender el modelo en cuestión, se presenta de forma breve la Ley de Mezclas, uno de los modelos 
analíticos más sencillos para el cálculo de las propiedades del material compuesto. Este método supone un estado 
de isodeformación o isotensión [6] no siendo muy preciso para ciertas propiedades. 
Considerando un volumen 𝑉, se cumplen las siguientes relaciones [7]: 
 𝑉 = 𝑉𝑓 + 𝑉𝑚 (32) 
 
1 =
𝑉𝑓
𝑉
+ 
𝑉𝑚
𝑉
 (33) 
 𝑣𝑓 + 𝑣𝑚 = 1 (34) 
donde 𝑣𝑓 y 𝑣𝑚 son las fracciones volumétricas de la fibra y la matriz respectivamente, que seguirán el modelo 
representado en la Figura 2-10: 
 
Figura 2-10. Ejemplos de elementos y su fracción volumétrica [19] 
 
En el caso en el que el volumen 𝑉 esté sometido a una carga en dirección de la fibra, se aplica la siguiente 
ecuación: 
 𝜎1𝐴 = 𝜎𝑓𝐴𝑓 + 𝜎𝑚𝐴𝑚 (35) 
 𝐹1 = 𝜎1𝐴 (36) 
29 
 
 
Figura 2-11. Volumen sometido a una fuerza en dirección longitudinal de la fibra [7] 
 
Siguiendo el razonamiento planteado en la ecuación (34), se puede calcular el módulo de Young en la dirección 
longitudinal 𝐸1 y el coeficiente de Poisson longitudinal 𝜐12 (suponiendo un estado de isodeformación) [6] 
también conocido como Ley de Mezclas Directa: 
 𝐸1 = 𝑣𝑓𝐸𝑓1 + 𝑣𝑚𝐸𝑚1 (37) 
 𝜐12 = 𝑣𝑓 𝜐𝑓+ 𝑣𝑚 𝜐𝑚 (38) 
Si el volumen está sometido a una carga transversal, se plantea la Ley de Mezclas Inversa que parte de la relación 
entre los espesores 𝑡 y deformaciones ε del compuesto, fibra y matriz: 
 𝑡 = 𝑡𝑓𝜀𝑓 + 𝑡𝑚𝜀𝑚 (39) 
 
Figura 2-12. Volumen sometido a una fuerza en dirección transversal de la fibra [7] 
 
Se estiman propiedades como el módulo de Young en dirección transversal 𝐸2 y módulos de cizalladura 
longitudinal y transversal 𝐺12 y 𝐺23 [7]: 
 1
𝐸2
= 𝑣𝑓
1
𝐸2𝑓
+ 𝑣𝑚
1
𝐸𝑚
 (40) 
 1
𝐺12
= 𝑣𝑓
1
𝐺12𝑓
+ 𝑣𝑚
1
𝐺𝑚
 (41) 
 1
𝐺23
= 𝑣𝑓
1
𝐺23𝑓
+ 𝑣𝑚
1
𝐺𝑚
 (42) 
 
 Ecuaciones constitutivas y modelos micromecánicos 
30 
 
30 
 
El modelo de Hopkins y Chamis surge como una mejora de la Ley de Mezclas, aportando más precisión a la 
hora de calcular las propiedades en dirección transversal. Se obtiene a partir del razonamiento anterior, aunque 
se añade al modelo una fase intermedia donde la fibra pasa a tener una sección rectangular y a ser de un material 
ficticio 𝑏 [7]: 
 
Figura 2-13. Esquema del modelo de Hopkins y Chamis [7] 
 
La relación entre la fracción de volumen del material ficticio 𝑏 y de la fibra 𝑓 es [7]: 
 𝑣𝑏 = √𝑣𝑓 (43) 
Y por tanto, para estimar las propiedades del material compuesto 𝐸2, 𝐺12 o 𝐺23, Hopkins y Chamis plantean las 
siguientes ecuaciones (siendo 𝑃 la propiedad a calcular) [6]: 
 
𝑃 = (
√𝑣𝑓
𝑃𝑏
+
1 − √𝑣𝑓
𝑃𝑚
)
−1
 (44) 
 𝜐12 = √𝑣𝑓 𝑃𝑓 + (1 − √𝑣𝑓)𝑃𝑚 (45) 
 
 
 
 
 
 
31 
 
2.4 Comparativa entre Halpin-Tsai, Mori-Tanaka y Hopkins y Chamis 
En este apartado se aplican los modelos teóricos micromecánicos previamente descritos para la obtención de las 
constantes ingenieriles de materiales compuestos de dos fases con disposiciones geométricas diferentes: fibra 
continua en matriz isótropa e inclusiones esféricas en matriz isótropa. Se compararán los resultados obtenidos 
en el cálculo de las propiedades 𝐸1, 𝐸2, 𝐺12 y 𝜐12 a partir de los modelos de Halpin-Tsai, Mori-Tanaka y Hopkins 
y Chamis. 
2.4.1 Fibra continua en matriz isótropa 
Para el caso del compuesto de fibra continua alineada en una matriz isótropa (ver Figura 2-5a), los materiales 
empleados consisten en una fibra de vidrio y una matriz de resina epoxy. Las propiedades de ambas fases se 
recogen en la siguiente tabla: 
 Módulo de Young 
Coeficiente 
de Poisson 
Matriz 𝐸𝑚 = 3.5𝐺𝑃𝑎 𝜈𝑚 = 0.35 
Fibra 𝐸𝑓 = 73𝐺𝑃𝑎 𝜈𝑓 = 0.22 
 
Tabla 2-4. Propiedades materiales del modelo fibra circular alineada 
Se obtienen los siguientes resultados de constantes ingenieriles para un amplio rango de fracciones de volumen 
de fibra 𝑣𝑓: 
 
Figura 2-14. Módulo de elasticidad 𝐸1 normalizado respecto al módulo de elasticidad de la matriz (𝐸𝑚) en función del 
porcentaje de fibra. 
 
 Ecuaciones constitutivas y modelos micromecánicos 
32 
 
32 
Los modelos de Halpin-Tsai y Hopkins y Chamis siguen la ley de mezclas para las propiedades 𝐸1 y 𝜐12, de ahí 
que sean totalmente coincidentes. Con respecto al modelo de Mori-Tanaka, los valores obtenidos para el módulo 
de Young longitudinal 𝐸1 coinciden exactamente con los dos anteriores. 
Si se representan los valores del coeficiente de Poisson, de nuevo coinciden los modelos de Halpin-Tsai y 
Hopkins y Chamis. Los resultados obtenidos a partir de Mori-Tanaka mantienen la linealidad y se desmarcan 
conforme el porcentaje de volumen de fibra es mayor. 
 
 
Figura 2-15. Coeficiente de Poisson 𝜐12 normalizado respecto al coeficiente de Poisson de la matriz (𝜐𝑚) en función del 
porcentaje de fibra. 
33 
 
 
Figura 2-16. Módulo de elasticidad 𝐸2 normalizado respecto al módulo de elasticidad de la matriz (𝐸𝑚) en función del 
porcentaje de fibra. 
 
Para el caso de la constante ingenieril 𝐸2, los modelos de Mori-Tanaka y Halpin-Tsai tienen valores muy 
similares. Sin embargo, el modelo de Hopkins y Chamis difieren de los anteriores en cuanto aumenta el 
porcentaje de volumen de fibra. 
 
 
 Ecuaciones constitutivas y modelos micromecánicos 
34 
 
34 
 
Figura 2-17. Módulo de cizalladura 𝐺12 normalizado respecto al módulo de cizalladura de la matriz (𝐺𝑚) en función del 
porcentaje de fibra. 
 
Por último, comparando los resultados obtenidos para el módulo de cizalladura 𝐺12, los modelos de Halpin-Tsai 
y Mori-Tanaka vuelven a ser iguales, y no muy distintos a los obtenidos en el modelo de Hopkins y Chamis. 
Todos los modelos son coincidentes para un volumen de fibra nulo. Dependiendo de la propiedad calculada, las 
similitudesentre modelos varían, aunque todos se mantienen en la misma tendencia. Se puede concluir que los 
tres modelos micromecánicos son válidos para el caso de un compuesto de fibras circulares continuas en una 
matriz isótropa. 
 
35 
 
2.4.2 Inclusiones esféricas en matriz isótropa 
Se estudia un segundo caso, donde el compuesto contiene inclusiones de carburo de tungsteno en una matriz de 
cobalto (ver Figura 2-7). Las propiedades de ambas fases se recogen en la siguiente tabla: 
 Módulo de Young 
Coeficiente 
de Poisson 
Densidad 
Matriz 𝐸𝑚 = 207𝐺𝑃𝑎 𝜈𝑚 = 0.3 𝜌𝑚 = 1150 𝑘𝑔/𝑚
3 
Inclusión 𝐸𝑓 = 703𝐺𝑃𝑎 𝜈𝑓 = 0.22 𝜌𝑓 = 1780 𝑘𝑔/𝑚
3 
 
Tabla 2-5. Propiedades materiales del modelo inclusiones esféricas en matriz isótropa 
 
Se trata de un material compuesto isótropo, por lo que las propiedades a calcular se simplifican a 𝐸, 𝐺 y 𝜐. A 
continuación, se muestran los resultados obtenidos en la comparativa de los modelos de Mori-Tanaka y Halpin-
Tsai (en concreto el modelo de Kerner) mencionados para distintos porcentajes de volumen de fibra: 
 
Figura 2-18. Módulo de elasticidad 𝐸 normalizado respecto al módulo de elasticidad de la matriz (𝐸𝑚) en función del 
porcentaje de fibra. 
 
 Ecuaciones constitutivas y modelos micromecánicos 
36 
 
36 
 
Los modelos de Halpin-Tsai y Mori-Tanaka para el caso del material con inclusiones esféricas en la matriz 
isótropa son totalmente coincidentes a lo largo de los distintos volúmenes de fibra. 
Con respecto a los valores obtenidos al calcular el coeficiente de Poisson, los modelos de Halpin-Tsai y Mori-
Tanaka empiezan obteniendo resultados similares para pequeños porcentajes de volumen de fibra, pero se 
desmarcan conforme aumenta el porcentaje de volumen de fibra. 
 
 
Figura 2-19. Coeficiente de Poisson 𝜐 normalizado respecto al Coeficiente de Poisson de la matriz (𝜐𝑚) en función del 
porcentaje de fibra. 
 
 
. 
37 
 
 
Figura 2-20. Módulo de cizalladura 𝐺 normalizado respecto al módulo de cizalladura de la matriz (𝐺𝑚) en función del 
porcentaje de fibra. 
 
Por último, en la obtención de los valores de la constante ingenieril 𝐺12, los modelos de Halpin-Tsai y Mori-
Tanaka vuelven a ser idénticos para el material en cuestión. 
De nuevo, se puede concluir que los modelos micromecánicos de Halpin-Tsai, Mori-Tanaka obtienen resultados 
parecidos para el cálculo de las constantes ingenieriles mencionadas para el caso de un material con inclusiones 
esféricas en una matriz isótropa y cuadrada. 
En el Anexo I se encuentran las rutinas de MATLAB utilizadas para la generación de resultados presentados en 
las figuras: Fig. 2-14 a Fig. 2-20. 
 
 Ecuaciones constitutivas y modelos micromecánicos 
38 
 
38 
 
39 
 
3 MODELO DE HOMOGENEIZACIÓN NUMÉRICA 
MEDIANTE ELEMENTO DE VOLUMEN 
REPRESENTATIVO 
ras explicar las leyes de comportamiento y algunos de los modelos micromecánicos teóricos en los 
apartados anteriores, en este capítulo se analiza la homogeneización numérica mediante una celda unitaria 
ó Elemento de Volumen Representativo (en adelante RVE, por sus siglas en inglés), modelizada mediante 
el Método de los Elementos Finitos (MEF). Para ello se introduce el RVE y las técnicas de homogeneización 
que permiten calcular las propiedades efectivas del material compuesto. Se analizan las condiciones de contorno 
a aplicar y finalmente se detalla el proceso de implementación numérica a seguir para resolver el problema 
elástico. 
3.1 Introducción 
Como se ha comentado en el capítulo anterior, el comportamiento de los materiales compuestos puede ser 
conocido a partir del estudio de sus constituyentes, siendo de gran importancia la disposición de las fibras en la 
matriz. Ejemplificando para el caso de la fibra continua a escala microscópica, su representación más común es 
la disposición cuadrada de fibras en la matriz. Debido a su simetría y su configuración periódica, seleccionando 
únicamente una porción rectangular se puede analizar el material a escala microscópica. Esta selección es 
conocida como un Elemento de Volumen Representativo (RVE) [5]. 
 
Figura 3-1. RVE para compuesto formado por fibras continuas [5] 
Por consiguiente, un RVE es un elemento diferencial del material que se caracteriza por [20]: 
▪ Ser representativo del compuesto. 
▪ Contener un número suficiente de inclusiones de manera que las propiedades obtenidas presenten 
independencia con respecto a las tensiones y desplazamientos aplicados en la superficie. 
 
Figura 3-2. Selección de un RVE representativo del material [13] 
T 
 
 Modelo de homogeneización numérica mediante Elemento de Volumen Representativo 
40 
 
40 
Para definir las propiedades del material compuesto, denominadas propiedades efectivas, se hace uso de las 
técnicas de homogeneización. Estas técnicas descomponen el análisis del material en dos escalas: una primera 
escala microscópica donde se estudian los detalles microestructurales teniendo en cuenta las propiedades de cada 
una de las fases y su heterogeneidad; y un segundo análisis a escala macroscópica que obtiene las tensiones y 
desplazamientos promedios al considerar el material homogéneo. 
 
Figura 3-3. Simplificación del cálculo de las propiedades efectivas del material [21] 
Existe una gran variedad de técnicas de homogeneización tanto numéricas como analíticas. En este proyecto se 
va a hacer uso de la homogeneización numérica, fundamentada en un análisis local-global del material 
empleando un RVE [22]. En ese primer análisis local, se selecciona un RVE de la microestructura y se resuelve 
el problema elástico tras aplicar unas condiciones de contorno apropiadas. Posteriormente, se procede a realizar 
un análisis global estableciendo una relación entre las propiedades calculadas y las tensiones y deformaciones 
promedio �̅� y 𝜀,̅ alcanzando las propiedades efectivas del material. 
 [𝑆] = [𝐶]−1 (46) 
 {𝜀} = [𝑆]{𝜎} (47) 
siendo [𝐶] , [𝑆] los tensores de rigidez y flexibilidad calculadas respectivamente promediados. 
 
Figura 3-4. Esquema de descomposición del análisis del cálculo de las propiedades efectivas mediante RVE [3] 
41 
 
3.2 Condiciones de contorno 
La periodicidad en la disposición de las distintas fases que conforman los materiales compuestos es una 
característica clave. Asumir esa periodicidad permite, a la hora de analizar las propiedades de los constituyentes, 
poder simplificar su estructura hasta reducirla al RVE. De forma general, los campos de tensión y deformación 
de un RVE se rigen por las siguientes expresiones: 
 𝜎𝑖 = 𝜎�̅�(𝑥) + 𝜎𝑖
𝑝
(y) (48) 
 𝜀𝑖 = 𝜀�̅�(𝑥) + 𝜀𝑖
𝑝
(y) (49) 
donde: 
▪ 𝜎�̅� y 𝜀�̅� son los valores promedios de tensión y deformación en el volumen 
▪ 𝜎𝑖
𝑝
 y 𝜀𝑖
𝑝
 corresponden con las perturbaciones periódicas con respecto a los términos de tensión y 
desplazamiento promedios [23]. 
Profundizando en su definición y particularizándolos al RVE, los términos de tensión y deformación media se 
expresan como: 
 
 
(50) 
 
 
(51) 
Con respecto a los términos de perturbaciones periódicas, al tratarse el RVE de un volumen periódico, adquieren 
un valor nulo [23]. 
 
 
(52) 
 
 
(53) 
Considerando los razonamientos anteriores, las condiciones de contorno adecuadas para poder solventar el 
problema elástico en el RVE deben de ser condiciones de contorno periódicas. 
El procedimiento utilizado para resolver el problema planteado en la escala microscópica consiste en imponer 
una deformación media conocida 𝜀�̅� y, aplicando la Ec. (47), obtener el valor de la tensión media resultante 𝜎�̅� 
[6]. Este método es conocido como el método de la imposición de la deformación media [9]. Conforme a lo 
desarrollado por Suquet (1987), para el caso anterior, el campo de desplazamientos 𝑢𝑖 puede expresarse como 
[23]: 
 𝑢𝑖 = 𝜀𝑖𝑗̅̅ ̅ · 𝑦𝑗 + 𝑢𝑖
𝑃(𝑦) (54) 
siendo 𝑢𝑖
𝑃(𝑦) los desplazamientos asociadosa las perturbaciones. Teniendo en cuenta la periodicidad y 
continuidad de los desplazamientos, se plantean las siguientes condiciones de contorno: 
 
 
 
 
 
 Modelo de homogeneización numérica mediante Elemento de Volumen Representativo 
42 
 
42 
 
𝑢𝑖
𝑃(𝑎1, 𝑦2, 𝑦3) = 𝑢𝑖
𝑃(−𝑎1, 𝑦2, 𝑦3) 
∀𝑦2 ∈ [−𝑎2, 𝑎2] 
(55) 
 ∀𝑦3 ∈ [−𝑎3, 𝑎3] 
 
𝑢𝑖
𝑃(𝑦1, 𝑎2, 𝑦3) = 𝑢𝑖
𝑃(𝑦1, −𝑎2, 𝑦3) 
∀𝑦1 ∈ [−𝑎1, 𝑎1] 
 ∀𝑦3 ∈ [−𝑎3, 𝑎3] 
 
𝑢𝑖
𝑃(𝑦1, 𝑦2, 𝑎3) = 𝑢𝑖
𝑃(𝑦1, 𝑦2, −𝑎3) 
∀𝑦1 ∈ [−𝑎1, 𝑎1] 
 ∀𝑦2 ∈ [−𝑎2, 𝑎2] 
siendo 𝑎1, 𝑎2 y 𝑎3 los parámetros geométricos del RVE mostrados en la Figura 3-5. 
 
Figura 3-5. Modelo de inclusión esférica con tres fases para vf = 0.1 y vf = 0.5 
Tal y como desarrolla Luciano [23], las deformaciones medias impuestas y, por tanto, las condiciones de 
contorno del problema se expresan con los siguientes términos [9]: 
 
𝑢𝑖(𝑎1, 𝑦2, 𝑦3) − 𝑢𝑖(−𝑎1, 𝑦2, 𝑦3) = 2𝜀𝑖1̅̅̅̅ 𝑎1 
∀𝑦2 ∈ [−𝑎2, 𝑎2] 
(56) 
 ∀𝑦3 ∈ [−𝑎3, 𝑎3] 
 
𝑢𝑖(𝑦1, 𝑎2, 𝑦3) − 𝑢𝑖(𝑦1, −𝑎2, 𝑦3) = 2𝜀𝑖2̅̅̅̅ 𝑎2 
∀𝑦1 ∈ [−𝑎1, 𝑎1] 
 ∀𝑦3 ∈ [−𝑎3, 𝑎3] 
 
𝑢𝑖(𝑦1, 𝑦2, 𝑎3) − 𝑢𝑖(𝑦1, 𝑦2, −𝑎3) = 2𝜀𝑖3̅̅̅̅ 𝑎3 
∀𝑦1 ∈ [−𝑎1, 𝑎1] 
 ∀𝑦2 ∈ [−𝑎2, 𝑎2] 
Con las ecuaciones anteriores y una condición de contorno adicional (restricción completa de uno de los puntos 
del RVE para evitar el movimiento de sólido rígido), el problema elástico queda completamente definido. 
3.3 Implementación numérica 
Una vez definido el RVE y las condiciones de contorno periódicas a aplicar, las ecuaciones de la elasticidad se 
resuelven a partir del método de elementos finitos (MEF) con ayuda del software Ansys APDL. Este método 
consiste en modelar un volumen y discretizarlo en un número de elementos finito que están conectados entre sí 
a través de nodos. El conjunto de elementos y nodos se denomina malla. Conociendo la geometría, las 
propiedades físicas de los materiales del modelo y las condiciones de contorno el problema computacional queda 
resuelto. 
43 
 
 
Figura 3-6. Ejemplo de modelo de inclusiones esféricas discretizado en distinto número de elementos 
Para poder obtener las propiedades efectivas del material compuesto, es necesario resolver el sistema de 
ecuaciones planteado en la ecuación (46). 
Por tanto, el procedimiento a seguir es el siguiente: 
 
 
Figura 3-7. Esquema implementación modelo de homogeneización numérica 
 
 Modelo de homogeneización numérica mediante Elemento de Volumen Representativo 
44 
 
44 
Para facilitar el cálculo del tensor de rigidez, se imponen deformaciones unitarias. El proceso se repite 6 veces 
generando 6 problemas independientes, uno por cada una de las deformaciones medias a aplicar. Una vez 
obtenido el tensor de rigidez por completo, se procede a estimar las propiedades efectivas conforme a la Tabla 
2-1. 
3.4 Comparativa modelo de homogeneización con modelos analíticos 
Una vez definidos los pasos a seguir para la implementación del método de elementos finitos y las respectivas 
condiciones de contorno para la obtención del tensor de rigidez, se procede a la aplicación de éste a los dos 
ejemplos de materiales compuestos analizados en apartados anteriores. Por un lado, se tiene un compuesto de 
fibra de vidrio en una matriz isótropa de resina epoxy representado con una geometría de fibras continuas en 
una matriz isótropa, como se muestra en la Figura 3-6, y por otro lado se modela un RVE de inclusiones esféricas 
de carburo de tungsteno en una matriz de cobalto cuadrada (ver Figura 3-8). 
El primero de los compuestos se trata de un material transversalmente isótropo, por lo que tras obtener el tensor 
de rigidez, se calculan las propiedades 𝐸1, 𝐸2, 𝜈12 y 𝐺12 que definen el material. Para el caso de las inclusiones 
esféricas, se hallan las constantes ingenieriles 𝐸 y 𝐺. 
A continuación, se comparan los valores obtenidos para cada una de las propiedades efectivas calculadas 
mediante el método de homogeneización numérica con aquellos valores generados por los modelos analíticos 
de Halpin-Tsai, Mori-Tanaka y Hopkins y Chamis definidos en el capítulo 2 de este proyecto. En caso de 
resultados similares quedaría validado el método de homogeneización numérica implementado. 
 
 
Figura 3-8. Ejemplo de modelo de elementos finitos de fibra continua en matriz isótropa discretizado en elementos 
 
45 
 
3.4.1 Fibra continua en matriz isótropa 
Los resultados obtenidos para el material compuesto de fibras continuas en una matriz isótropa para las 
constantes ingenieriles 𝐸1, 𝐸2, 𝜈12 y 𝐺12 se representan en las siguientes gráficas: 
 
 
Figura 3-9. Comparativa modelos analíticos y modelos mediante RVE para 𝐸1 
 
 Modelo de homogeneización numérica mediante Elemento de Volumen Representativo 
46 
 
46 
 
Figura 3-10. Comparativa modelos analíticos y modelos mediante RVE para 𝐸2 
 
Los resultados obtenidos en el modelo de homogeneización numérica mediante el RVE por método de 
elementos finitos presentan valores muy similares con respecto a los modelos teóricos, en concreto con el 
modelo de Mori-Tanaka para porcentajes bajos de volumen de fibra (Figura 3-9 y Figura 3-10). Esta semejanza 
disminuye conforme aumenta la fracción volumétrica de fibra puesto que se mantiene el mismo mallado o 
discretización de elementos, para todo el rango de porcentajes de volumen de fibra, lo cual desvirtúa el resultado. 
47 
 
 
Figura 3-11. Comparativa modelos analíticos y modelos mediante RVE para 𝜈12 
 
Para el cálculo de 𝜈12 (Figura 3-11), los resultados obtenidos mediante el método de elementos finitos son 
prácticamente iguales a los alcanzados por el método analítico de Mori-Tanaka, para todo el intervalo analizado 
de porcentajes de volumen de fibra. 
Por último, analizando los resultados hallados para la constante ingenieril 𝐺12 en la Figura 3-12, coincide el 
modelo de homogeneización numérica mediante el RVE con los modelos teóricos micromecánicos de Halpin-
Tsai y Mori-Tanaka. Se mantiene de forma generalizada una mayor similitud para fracciones volumétricas de 
fibra bajas, debido al mallado común de la fibra para los distintos elementos de volumen representativos 
estudiados (a mayor tamaño de elementos, menor precisión). 
Se puede concluir que los resultados obtenidos por el método de elementos finitos mediante el elemento de 
volumen representativo son similares a los alcanzados mediante los modelos analíticos, y por lo tanto, se trata 
de un método válido. 
 
 
 
 Modelo de homogeneización numérica mediante Elemento de Volumen Representativo 
48 
 
48 
3.4.2 Inclusión esférica en matriz cuadrada isótropa 
A continuación, se implementa el modelo de homogeneización numérica para el material compuesto constituido 
por inclusiones esféricas en la matriz isótropa (ver Figura 2-5b). Se comparan los resultados obtenidos para las 
propiedades 𝐸 y 𝐺, obteniendo poca diferencia entre los modelos teóricos micromecánicos y el modelo de 
elementos finitos. Los valores no son tan dispares en los casos de volumen de fibra elevados, debiéndose a una 
geometría más homogénea del modelo y por tanto, un mallado más regular. 
 
 
Figura 3-12. Comparativa modelos analíticos y modelos mediante RVE para 𝐺12 
 
49 
 
 
Figura 3-13. Comparativa modelos analíticos y modelos mediante RVE para 𝐸 
 
 
 
 Modelo de homogeneización numérica mediante Elemento de Volumen Representativo 
50 
 
50 
 
Figura 3-14. Comparativa modelos analíticos y modelos mediante RVE para 𝐺 
 
En el Anexo II se encuentra el listado de comandos de ANSYS y las rutinas de MATLAB usadas para la 
generación de resultados presentados en las figuras: Fig. 3-9 a Fig. 3-14. 
51 
 
4 MODELIZACIÓN DE INTERFASE INCLUSIÓN-
MATRIZ IMPERFECTA 
n este capítulo se introduce el concepto de la interfase imperfecta en el material y cómo puede definirse 
su comportamiento a partir de los parámetrosde rigidez. Además, se presentan varios modelos 
micromecánicos que incluyen la interfase imperfecta y la caracterizan, tanto modelos teóricos como 
modelos de homogenización numéricos basados en el elemento de volumen representativo (RVE). Por último, 
para cada uno de los modelos mostrados, se detalla cómo se ha realizado su implementación en este proyecto. 
4.1 Introducción 
Dada una matriz de material (2) y una inclusión de material (1) delimitadas por una interfase i, se tiene el 
siguiente esquema: 
 
Figura 4-1. Esquema interfase 
 
Se define el parámetro η como el cociente entre el espesor de la interfase y el radio de la inclusión: 
 𝜂 =
𝑡𝑖
𝑅⁄ → parámetro de la interfase 
(57) 
Se considera una interfase perfecta como la capa intermedia entre la matriz y la inclusión por la que se transmiten 
las tensiones y desplazamientos de forma continua. Sin embargo, existen una serie de modelos que indican que 
las propiedades de la interfase deben ser ligeramente diferentes a las propiedades del resto de constituyentes, 
pudiendo ser más flexibles. 
 
E 
R 
𝑡𝑖 
Matriz (2) 
Interfase (i) 
Inclusión (1) 
Siendo 
R: radio de la inclusión 
𝑡𝑖: espesor de la interfase 
 
 
 Modelización de interfase inclusión-matriz imperfecta 
52 
 
52 
La imperfección de la interfase conlleva que las tensiones de la matriz y la inclusión se mantengan continuas 
pero se produzca un salto en los desplazamientos [24]. Para facilitar el modelaje de la interfase imperfecta (o 
también conocida como interfase dañada), se aplica una metodología que convierte un modelo 3D, con una 
transición de propiedades entre ambos materiales, en una interfase 2D para así reducir el número de incógnitas 
que definen dicha interfase [25]. 
 
Figura 4-2. Representación de discontinuidad de desplazamiento en interfase 2D [25] 
Estas discontinuidades de los desplazamientos en las direcciones normales y tangenciales están relacionadas con 
las respectivas tracciones de la interfase mediante los parámetros de rigidez de la interfase: 
 𝑡𝑛 = 𝑘𝑛𝑢𝑛 (58) 
 
Figura 4-3. Direcciones interfase [26] 
 𝑡𝑡 = 𝑘𝑡𝑢𝑡 (59) 
siendo: 
𝑡𝑖 tracciones 
𝑢𝑖 discontinuidades en los desplazamientos normales y tangenciales 
𝑘𝑖 parámetros de rigidez de la interfase 
Se diferencian dos casos límites en la interfase imperfecta, dependiendo del valor de los parámetros de rigidez: 
A. Interfase completamente imperfecta: los parámetros de rigidez normal 𝑘𝑛 y tangencial 𝑘𝑡 tienen un 
valor nulo. Las discontinuidades son tal que la matriz (2) y la inclusión (1) se encuentran despegadas. 
 𝑘𝑛 = 𝑘𝑡 = 0 (60) 
B. Interfase perfecta: en este caso, los parámetros de la rigidez tienden a infinito, obteniéndose unas 
discontinuidades de desplazamientos nulas y por tanto, una unión rígida entre la inclusión (1) y la matriz 
(2). 
 𝑘𝑛 = 𝑘𝑡 = ∞ (61) 
Por tanto, la interfase queda definida de manera imperfecta para cualquier valor finito de los parámetros de 
rigidez. 
En los siguientes apartados, se presentan distintos modelos micromecánicos que incluyen la imperfección de la 
interfase, o también conocida como daño en la interfase, y relacionan los parámetros de rigidez de la interfase 
con las propiedades mecánicas de ésta. 
53 
 
4.2 Modelo de Mori-Tanaka para modelar interfase inclusión-matriz imperfecta 
4.2.1 Modelo de Mori-Tanaka con daño en la región interfacial 
El modelo de Mori-Tanaka adaptado por Lee [7] [8] asemeja el comportamiento de la interfase imperfecta al 
comportamiento de un resorte. En concreto, Lee considera la imperfección de la interfase como un muelle lineal 
sin espesor [27]. La flexibilidad del muelle a través de la interfase está definida en un rango entre cero e infinito, 
y se corresponde con una interfase perfecta o completamente dañada respectivamente (aquella que no transmite 
cargas entre la matriz y la inclusión) [28]. 
 
Figura 4-4. Esquema del modelo de interfase con resorte [28] 
 
Como se muestra en el esquema anterior, la flexibilidad del muelle a través de la interfase puede expresarse 
como la combinación de la flexibilidad del muelle en dirección tangencial (α) y la flexibilidad del muelle en 
dirección normal (β) a la interfase. 
Según el modelo presentado por Zhong y Meguid [29], el campo de desplazamientos es uniforme siempre y 
cuando se cumplan dos condiciones; la primera, que las flexibilidades del muelle a través de la interfase sean 
iguales en ambas direcciones y en segundo lugar, que la inclusión tenga una forma geométrica perfectamente 
esférica [28]. 
 𝛾𝑠 = 𝛼 = 𝛽 (62) 
Asumiendo las hipótesis previas, Lee define la matriz efectiva teniendo en cuenta el daño de la interfase a partir 
de las siguientes ecuaciones: 
 [𝐶] = (𝑣𝑚[𝐶𝑚][𝐼] + 𝑣𝑓[𝐶𝑓][𝐴𝑓] + 𝑣𝑖[𝐶𝑖][𝐴𝑖])[𝐴1] (63) 
 
[𝐴1] = (𝑣𝑚[𝐼] + 𝑣𝑓[Â] + 𝒗𝒇
𝜸𝒔
𝑹
[𝑪𝒇][Â])
−1
 (64) 
 
[Â] = [[𝐼] + [𝑆𝑚][𝐶𝑚]
−1([𝐶𝑓] − [𝐶𝑚]) +
𝜸𝒔
𝑹
([𝑰] − [𝑺𝒎])[𝑪𝒇]]
−1
, (65) 
donde se añaden una serie de términos a las ecuaciones convencionales del modelo de Mori-Tanaka que incluyen 
el parámetro que indica la flexibilidad de la interfase 𝛾𝑠 y por tanto, el daño en la interfase . 
 
 
 Modelización de interfase inclusión-matriz imperfecta 
54 
 
54 
El parámetro del daño de la interfase queda definido con la siguiente expresión: 
 𝐺𝑚
𝛾𝑠
𝑅
 (66) 
siendo 
R el radio de la inclusión 
𝐺𝑚 el módulo transversal de la matriz 
El valor de este parámetro queda comprendido en el intervalo [10−3, 103], por lo que se determinan dos casos 
extremos de la interfase a partir del valor del parámetro de daño: 
▪ Condición de interfase imperfecta: 𝛾𝑠 = 10
3 𝑅
𝐺𝑚
⁄ 
▪ Condición de interfase perfecta: 𝛾𝑠 = 10
−3 𝑅
𝐺𝑚
⁄ 
4.2.2 Implementación 
Para poder implementar el modelo de Mori Tanaka adaptado por Lee [27] [28] y modelar la interfase inclusión-
matriz imperfecta, se ha generado un código en Matlab. A continuación, se muestra un esquema de los pasos 
seguidos a lo largo del código de implementación que se encuentra en el Anexo IV. 
 
Figura 4-5. Esquema de implementación del modelo Mori Tanaka-Lee con interfase imperfecta propuesto por Lee [27] [28] 
55 
 
4.3 Modelo basado en RVE para modelar interfase inclusión-matriz imperfecta 
En el apartado anterior, se ha analizado el modelo de Mori Tanaka presentado por Lee como un modelo teórico 
que caracteriza la interfase imperfecta y su comportamiento a partir del parámetro de daño tal y como se muestra 
en la ecuación 65. Sin embargo, también existen otros modelos que definen el comportamiento de la interfase 
imperfecta basados en un RVE por el modelo de elementos finitos. 
En los apartados posteriores se presenta el modelo de Hashin que indica cómo a partir de los valores aplicados 
a las propiedades de la interfase en el RVE-MEF, se obtiene el comportamiento de interfase imperfecta. 
 
 
Figura 4-6. Descomposición del RVE de inclusiones esféricas en matriz cuadrada isótropa con interfase imperfecta 
definida por el parámetro 𝜂 
 
 
 
 Modelización de interfase inclusión-matriz imperfecta 
56 
 
56 
4.3.1 Modelo de Hashin 
El modelo de Hashin explica que las condiciones de la interfase imperfecta son definidas como una relación 
lineal entre la discontinuidad en los desplazamientos de la interfase y las respectivas componentes de las 
tensiones a tracción [24]. Además, añade que estas condiciones representan el efecto de una interfase muy 
delgada y elástica. Finalmente concluye que los parámetros de rigidez en dirección normal (𝑘𝑛) y tangencial 
(𝑘𝑡) se calculan a partir del espesor de la interfase y sus propiedades elásticas [30]. 
 
𝑘𝑛 =
𝜆𝑖 + 2𝐺𝑖
𝑡𝑖
 (67) 
 
𝑘𝑡 =
𝐺𝑖
𝑡𝑖
 (68) 
Estas expresiones son válidas tanto para interfases muy flexibles como para interfases más rígidas y de cualquier 
forma geométrica [24]. Si se reordenan las ecuaciones