Micromechanics for Polymer Matrix Laminated Composites: A Literature Review of Fundamental Models, Advances and a Trend Analysis for Future Researches

1. Introduction

With the marked technical-scientific development in composite materials, mainly from the 1970s, several industrial sectors began to incorporate this class of materials in their components, in particular, because they present satisfactory specific mechanical properties, given their low density in comparison to steel, for example, until nowadays typically used in structural applications [1].

Currently, the great technological challenge of the energy transition needs to be faced by engineers and scientists and at the same time, there is also the process of greenhouse gas emission reduction. In this complex context, the application and research focused on laminated composites gains even more encouragement in the industrial and academic areas, as this class of multiphase materials is seen as an excellent alternative for structural parts and components in the renewable energy sectors, including, with regard to minimizing CO2eq emissions throughout the life cycle of their equipments [2].

Given this favorable scenario for structural laminated composites, it is of utmost importance to further develop the knowledge about the several thermomechanical behaviors that they can present. Therefore, as they are multiphase materials, this heterogeneity needs to be better accounted for the prediction of effective properties necessary for the macro analysis of a structure made up of these kind of composites. The purpose of micromechanical modeling focused on these laminates is to make a theoretical determination as to what would effectively be the homogenized characteristic tensors that should be used into the constitutive relations of these materials [3].

Aiming to provide a broad overview of micromechanical models for the laminated composite analyzes, this research is designed to review the fundamental modeling, that is, the most elementary and precursors models in the development of this branch of solid mechanics, as well as exploring the state of the art models whose proposals have been published in recent years and, finally, this study presents trends for future researches involving the topic.

In this sense, the research shows highlighted elementary models highlighted in the literature, such as: Voigt and Reuss models (mixture rule), the Chamis model, Bridging, the modeling by Eshelby equivalent inclusions, the Mori-Tanaka method, the model for predicting elastic constants in fabric-reinforced composites, as well as citations made to conventional numerical models.

Regarding the advances achieved in the past years and more recent proposals, the following stand out: reformulations for the high-fidelity generalized cell model, including the use of continuum damage mechanics to analyze the nucleation of discontinuities in the material, the multilayer concentric cylinder method, approaches including the cohesive zone model for evaluating the initiation of interlaminar failures, modeling that study the micromechanical mechanisms associated with delamination, models for inelastic behavior prior to failure, dynamic models, approaches for progressive microdamages analysis based on non-local continuum mechanics and stochastic models.

2. Literature Review

2.1. Fundamentals and Previous Studies

All material, at a certain scale of analysis, can be classified as heterogeneous. The major contribution of micromechanics is to consider the heterogeneity of the material while allowing it to be represented as a pseudo-homogeneous continuum in analyzes referring to larger length scales (structural sizes) [4].

The concepts of representative volume elemento (RVE) and repeated unit cell (RUC) are important in microanalysis. RVE concerns a volume with sufficient dimensions to consider the constituent microstructure of the material and, at the same time, that is representative of the material as a whole in terms of effective mechanical behavior, that is, when applying a tension or a displacement as boundary conditions to this RVE, its mechanical behavior would be the same as that of the material on a global scale (for a perfectly homogeneous material, this volume could be much smaller). The RUC concept establishes the existence of a standard microstructure that is repeated throughout the volume of the material, i. e, it considers the microstructure as periodic [4].

The properties of a lamina or assembled laminate can be determined experimentally or designed and predicted mathematically based on the properties of the constituents of the continuous and dispersed phases by a micromechanical modeling. As will be seen, the purpose of all smaller-scale analysis is to determine how the constituent elements of the material interact mechanically and contribute to the final properties of the composite, for example, the moduli of elasticity, compliance or elastic stiffness in terms of the properties of reinforcement and matrix elements, equation (1) [3]:

$$C_{ij}=C_{ij}(E_f,{\nu }_f,V_f,E_m,{\nu }_m,V_m)$$

(1)

where, $E_f$ is the modulus of elasticity of an isotropic fiber, ${\nu }_f$is Poisson’s ratio and $V_f$ is the ratio between the volume of fibers and the composite volume, being $E_m,{\nu }_m,V_m$ the respective analogues of the material that makes up the matrix [3].

By considering simplifying hypotheses that enable an effective modeling to the mechanical behavior of the composite, the mechanics of material approach to determining stiffness is common and practical. The most prominent hypothesis among those that set up this model is the one that considers the same strains for fiber and matrix when there is load in the direction of the reinforcing fiber in a unidirectional laminated composite. This hypothesis is reasonable, since if there were no this motion condition, there would be a detachment or rupture between the two components [3].

Based on the aforementioned hypotheses is the Voigt model, which leads to the rule of mixture for the composite stiffness tensor components. For a two-phase material, could be written [4]:

$$C=c_m{C}^m+c_f{C}^f$$

(2)

where, $c_m$ is the volume fraction of the matrix phase, $c_f$ is the volume fraction of the reinforcement phase, $C^m$ is the stiffness tensor of the matrix phase, $C^f$ constitutes the stiffness tensor of the reinforcement phase and C is the effective stiffness tensor of the composite [4].

The Reuss approximation assumes that the stresses acting on the different phases of the composite are uniform and equal to the average stresses, this implies the rule of mixture to the compliance tensors [4]:

$$S=s_m{S}^m+s_f{S}^f$$

(3)

where, $s_m$ is the volume fraction of the matrix phase, $s_f$is the volume fraction of reinforcement phase, $S^m$ is the compliance tensor of the reinforcement phase, $S^f$is the compliance tensor of the reinforcement phase and S is the effective compliance tensor of the composite [4].

The limitation for the Voigt model is found in the stresses arising from the same strains for fiber and matrix do not bring the boundaries of these phases to an equilibrium condition. Regarding the Reuss model, it appears that strains arising from the hypothesis of uniform stresses that makes the cohesion of phases physically impossible [4].

In the specific literature, there are some authors who seek to combine the Voigt and Reuss models to reasonably predict effective properties of certain composites in engineering, even though the two models do not satisfy the basic equations of the classical theory of elasticity [4]. It is demonstrated by the so-called Hill’s Theorem that the Voigt and Reuss models provide upper and lower limits on the effective properties of the composite material.

The bounding theorems propose the establishment of narrower limits than Voigt and Reuss models predict for the material constant values, within which it is known that the effective elastic properties of the composites are adjusted. Considering the interaction of elastic properties of the constituent phases and not just the volume fractions, it is clear that for a two-phase composite, in which the axis of symmetry for the transversely isotropic effective material coincides with the longitudinal direction, the relationship obtained for the elastic modulus was [4]:

$$\frac{c_m{c}_f}{\frac{c_m}{k_f}+\frac{c_f}{k_m}+\frac{1}{G_m}}\le \frac{E_L-c_m{E}_{Lm}-c_f{E}_{Lf}}{4({\nu }_m-{\nu }_f)²}\le \frac{c_m{c}_f}{\frac{c_m}{k_f}+\frac{c_f}{k_m}+\frac{1}{G_f}}$$

(4)

where, in addition to the already described volume fractions of the matrix and reinforcement constituents, there are:$k_m$>$k_f$ as the bulk modulus of matrix and fibers, respectively, $G_f>G_m$ shear moduli of the dispersed and cohesive phases, ${\nu }_m$ and ${\nu }_f$ are the respective Poisson’s ratios of matrix and fibers, E_L is the effective longitudinal elastic modulus of the composite and $E_{Lm}$ and $E_{Lf}$ are also the respective longitudinal elasticity moduli of matrix and reinforcement [4].

The Halpin-Tsai approach consists of a semi-empirical model that provides better fits between theoretical prediction and experimental results for both the effective elastic modulus transverse to fibers (E_T), as for effective shear modulus (G_LT) in the plane that contains the longitudinal and transverse directions to the fibers as indicated in Figure 1, that represents in-plane shear modulus [5].

According to the Halpin-Tsai model, could be written [5]:

$$E_T=E_m\left(\frac{1+\xi \eta c_f}{1-\eta c_f}\right)$$

(5)

With

$$\eta =\frac{\left(\frac{E_f}{E_m}\right)-1}{\left(\frac{E_f}{E_m}\right)+\xi }$$

(6)

$$G_{LT}=G_m\left(\frac{1+\xi \eta c_f}{1-\eta c_f}\right)$$

(7)

$$\eta =\frac{\left(\frac{E_f}{E_m}\right)-1}{\left(\frac{E_f}{E_m}\right)+\xi }$$

(8)

$\xi$ is the so-called reinforcing factor, which is a function of fiber geometry, fiber packing geometry and loading conditions [5].

Using the classical elasticity theory and the so-called average strain and average stress theorems for a RVE, it is shown using homogeneous boundary conditions and considering perfect cohesion at the interface and that the effective elasticity and compliance tensors for a two-phase material can be approximated by [4]:

$$C_{ijkl}{\epsilon }_{kl}^0=C_{ijkl}^m{\epsilon }_{kl}^0+c_f(C_{ijkl}^f-C_{ijkl}^m){\epsilon }_{kl}^{*f}$$

(9)

$$S_{ijkl}{\sigma }_{kl}^0=S_{ijkl}^m{\sigma }_{kl}^0+c_f\left(S_{ijkl}^f-S_{ijkl}^m\right){\sigma }_{kl}^{*f}$$

(10)

In which $C_{ijkl}$ and $S_{ijkl}$ are the effective stiffness and compliance tensors for the two-phase composite, ${\epsilon }_{kl}^0$ and ${\sigma }_{kl}^0$ are the strain and stress tensors induced by the homogeneous displacement boundary conditions and loading, respectively. ${\epsilon }_{kl}^{*f}$ and ${\sigma }_{kl}^{*m}$ are the average strain and stress tensors over the volume occupied by the reinforcement phase [4].

The so-called dilute approximation is a modeling that disregards the interaction between inclusions, then, the multiphase material can be assumed as a single inclusion embedded in continuous phase. Based on this hypothesis, for example, the effective bulk modulus, equation (11), for spherical inclusions can be demonstrated [4].

$$K=K_m+\frac{c_i(K_i-K_m)(3K_m+4{\mu }_m)}{(3K_i+4{\mu }_m)}$$

(11)

In which, $c_i$ is the volumetric fraction of concentration of the diluted inclusions, $K_m,K_i,{\mu }_m$ are the bulk modulus of the matrix, of the inclusion and one of the Lamé constants for the matrix, respectively.

The Mori-Tanaka Theory (MT) has proven to be relevant in many studies that search to determine the effective behavior of composite materials. This theory is based on the calculation of average internal stresses and strains in the matrix of a material containing precipitates using the concept of eigenstrains [4,70]. Then, the effective stiffness tensor for a two-phase composite can be described by:

$$C=C^m+c_f(C^f-C^m)A^f$$

(12)

where, $A^f$ is a concentration tensor associated with the linear transformation between the strain tensor induced by homogeneous displacement boundary conditions and the average strain tensor for the dispersed phase [4,70].

$${\epsilon }^{*f}=A^f{\epsilon }^0$$

(13)

Hence, the Mori-Tanaka theory originally considers that the interactions within the composite caused by the inclusions of the dispersed phase can be approximated by a single inclusion in the infinite matrix of the continous phase, which is subjected to the average strain tensor ${\epsilon }^{*m}$, i.e., the volumetric concentration of inclusions needs to be small for the model to have a satisfactory approximation [4,70].

The text written by Aboudi et al [4] even highlights an elementary micromechanical model called differential schemes (DS). In this differential model, a two-phase composite of initial volume $V_0$, first of all, has an increment of volume $\Delta V$ removed. This fact implies a corresponding removal of inclusion/reinforcement $\Delta V{c}_i$. Consecutively, instead of this volume of composite removed, the same amount of inclusion/reinforcement is added $\Delta V$. Consequently, the current volume of the dispersed phase is:

$$V_0{c}_{f(i+1)}=V_0{c}_{fi}+\Delta V-\Delta V{c}_{fi}$$

(14)

Which, in differential form, is [4]:

$$d{c}_f=\frac{dV}{V_0}(1-c_f)$$

(15)

And, also based on equation (11), what the differential model proposes is that the effective stiffness tensor of the composite in a new volume fraction of fibers $C\left(c_f+\Delta c_f\right)$ is approximated to [4]:

$$C\left(c_f+\Delta c_f\right)=C\left(c_f\right)-\frac{\Delta V}{V_0}\left[C_f-C\left(c_f\right)\right]A^f$$

(16)

where a small volume $\Delta V$ of inclusions/reinforcements, whose stiffness tensor is $C_f$, are homogeneously embedded in a continous phase of stiffness given by $C\left(c_f\right)$ , therefore, in differential form can be written [4]:

$$\frac{dC}{dc}=\frac{1}{1-c}(C_f-C)A^f$$

(17)

With $C=C^m$ when $c_f=0$ [4].

The Eshelby’s equivalent inclusion method works with inclusion in the ellipsoidal form embedded in an infinite matrix phase, i.e, the homogeneous boundary conditions act so far from the inclusion. The solution for the effective elastic properties based on the principle of equivalent ellipsoidal inclusion can be derived for different other geometries, such as: spheres, cylinders with an elliptical base, possible cracks, etc. [4,59,67].

The Eshelby approximation considers an infinite homogeneous matrix whose stiffness tensor is $C^m$ and which contains an ellipsoidal domain Ω with an eigenstrain ${\epsilon }^t$ (it could be a thermal strain, for example). As the natural change in the geometry of this domain Ω is restricted to matrix contours, a stress field develops that can be expressed by [4,59,67]:

$$\sigma =C^m(\epsilon -{\epsilon }^t)$$

(18)

where, $\epsilon$ is the total strain given by the sum of the eigenstrain and the elastic strain induced by the stresses. Eshelby showed that if it ${\epsilon }^t$ is uniform in this domain, then, $\epsilon$ it is also uniform in it and both strains are related through the so-called Eshelby tensor (P), whose components are functions of elastic constants and geometry of the inclusions [4,59,67]:

$$\epsilon =P{\epsilon }^t$$

(19)

When an infinite matrix containing this ellipsoidal domain Ω is subjected to a homogeneous boundary condition of displacements, the total stress tensor is given by [4,59,67]:

$$\sigma ={C^m{\epsilon }^0+C}^m(\epsilon -{\epsilon }^t)$$

(20)

Then, the Eshelby model replaces the ellipsoidal domain whose stiffness is given by the tensor$C^m$ by an equivalent stiffness inclusion $C^f$ such that [4,59,67]:

$$C^m\left({\epsilon }^0+\epsilon -{\epsilon }^t\right)=C^f\left({\epsilon }^0+\epsilon \right)\mathrm{in}\Omega .$$

(21)

Using the concept of equivalent inclusions of Eshelby in the Mori-Tanaka method, a three-phase model for this theory assumes the presence of an intermediate phase between the fiber and the matrix. Using multiple inclusions, the explicit formulation of this model is given by [5]:

$$C=\left[c_m{C}^{m}I+c_f{C}^f{A}^f+c_i{C}^i{A}^i\right]{\left[c_{m}I+c_f{A}^f+c_i{A}^i\right]}^{-1}$$

(22)

where the terms $c_i,C^i,A_i$ are the volume fraction, the stiffness tensor and the concentration tensor of the intermediate phase, respectively.$c_i,C^i,A_i$

$$A^i={\left[I+P^i\left({\left(C^m\right)}^{-1}\left(C^f-C^m\right)\right)\right]}^{-1}$$

(23)

$$A^f={\left[I+P^f\left({\left(C^m\right)}^{-1}\left(C^f-C^m\right)\right)\right]}^{-1}$$

(24)

In which $P^i$ and $P^f$ are the Eshelby tensors for the domains of the intermediate and reinforcement phases, respectively. For a cylindrical inclusion of reinforcement-intermediate phase in an isotropic matrix, it is expressed by [5]:

$$P=\left[\begin{array}{cccccc}{P}_{1111}& P_{1122}& P_{1133}& 0& 0& 0\\ P_{2211}& P_{2222}& P_{2233}& 0& 0& 0\\ P_{3311}& P_{3322}& P_{3333}& 0& 0& 0\\ 0& 0& 0& P_{2323}& 0& 0\\ 0& 0& 0& 0& P_{1313}& 0\\ 0& 0& 0& 0& 0& P_{1212}\end{array}\right]$$

(25)

For a cylindrical inclusion of reinforcement-intermediate phase in an isotropic matrix, the elements of $P^i$ are expressed by [5]:

$P_{1111}=0;P_{2222}=P_{3333}=\frac{5-4c_m}{8(1-c_m)};P_{2211}=P_{3311}=\frac{c_m}{2(1-c_m)};P_{2233}=P_{3322}=\frac{4c_{m-1}}{8(1-c_m)};P_{1122}=P_{1133}=0;P_{1313}=P_{1212}=0;P_{2222}=P_{2323}=\frac{3-4c_m}{8(1-c_m)}$

And for a cylindrical reinforcement/fiber inclusion in an intermediate phase, the elements of $P^f$ tensor are given by [5]:

$P_{1111}=0;P_{2222}=P_{3333}=\frac{5-4c_i}{8(1-c_i)};P_{2211}=P_{3311}=\frac{c_i}{2(1-c_i)};P_{2233}=P_{3322}=\frac{4c_{i-1}}{8(1-c_i)};P_{1122}=P_{1133}=0;P_{1313}=P_{1212}=0;P_{2222}=P_{2323}=\frac{3-4c_i}{8(1-c_i)}$

The cylindrical composite assemblage model aims to determine the effective elastic modulus of composites with hollow circular fiber alignments. For a hexagonal arrangement, the following prediction was obtained for the effective longitudinal modulus of elasticity of the composite [5]:

$$E_L=E_m\left(c_f\frac{E_f}{E_m}+c_m\right)\frac{E_m\left(A_1-A_3{B}_1\right)+E_f\left(A_2-A_4{B}_2\right)}{E_m\left(A_1-A_3\right)+E_f\left(A_2-A_4\right)}$$

(26)

Where:

$A_1=\frac{1+a^2}{1-a^2}-c_f;A_2=\frac{1+c_g}{c_m}+c_m;A_3=\frac{2c_f^2}{1-a^2};A_4=\frac{2c_m^2{c}_g}{c_m};$$B_1=\frac{c_m{c}_f{E}_f+c_f{c}_m{E}_m}{c_f{c}_f{E}_f+c_m{c}_m{E}_m};B_2=\frac{c_f{B}_1}{c_m};a=\frac{R_1}{R_2};c_g=\frac{R_2^2}{R^2};c_f=\left(1-a^2\right)c_g;=1-c_g$

Figure 2. Cylindrical hollow fiber as matrix reinforcement and respective internal and external diameters of the fiber and matrix [5].

Figure 2. Cylindrical hollow fiber as matrix reinforcement and respective internal and external diameters of the fiber and matrix [5].

Kumar et al [6] develop a model also considering three phases (matrix, interface and unidirectional fibers). For physical modeling, researchers assume that the interface is an isotropic and homogeneous material, thus, elastic constants such as $E_{11},E_{22},G_{12},G_{23}$ and damping coefficients ${\eta }_{11},{\eta }_{22},{\eta }_{12},{\eta }_{23}$ can be obtained from the related mathematical model. From the calculations performed, the study shows that the volume fraction of the interface, as well as, its determined properties, generate a significant effect on the loss factor of the composite and, in the end, conclude that the results of the proposed model present acceptable congruence with the modeling numerical performed by finite elements [6].

The method of cells approach assumes uniform and aligned spacing for the fibers in a matrix, meaning that the composite can be viewed as a periodic matrix arrangement in transverse directions to the fibers. In this model, each cell corresponds to the representative volume element and is composed of subcells as indicated in Figure 3 [5].

The volume of each cell can be expressed by [5]:

$$V_{bc}=l{B}_b{H}_c$$

(27)

Where $l{,B}_b,H_c$ are the length, width and height of each subcell identified by subindexes b and c. Therefore, the cell volume is given as follows [5]:

$$V=lBH$$

(28)

Where, $B=\sum _{b=1}^M{B}_{b}eH=\sum _{c=1}^N{H}_c$

Applying Hooke’s law to each subcell, can be written [5]:

$${\sigma }^{bc}=C^{bc}{\epsilon }^{bc}$$

(29)

As each cell represents a representative volume element for the composite, the average strains and stresses for each of them can be written as follows [5]:

$${\epsilon }^{cell}=\frac{1}{V}\sum _{b=1}^M\sum _{c=1}^N{V}_{bc}{\epsilon }^{bc}$$

(30)

$${\sigma }^{cell}=\frac{1}{V}\sum _{b=1}^M\sum _{c=1}^N{V}_{bc}{\sigma }^{bc}$$

(31)

Then, applying now Hooke’s law for a cell of the composite material, it can be expressed by equation (32) [5]:

$${\sigma }^{cell}=C{\epsilon }^{cell}$$

(32)

Where, $C=\frac{1}{V}\sum _{b=1}^M\sum _{c=1}^N{V}_{bc}{C}^{bc}{A}^{bc}{A}^{bc}$ is the effective stiffness tensor, with representing the strain concentration matrix of each subcell. An extension to global coordinates of a laminate can also be seen in [5].

In research aimed at predicting progressive failures in multidirectional laminated composites, Basu et al [7] (noting that local rotations combined with the continuous degradation of the matrix are the mechanisms responsible for fiber microbuckling) develop a modeling based on the so-called Schapery theory. Thermodynamically based, this theory assumes internal state variables in order to analyze the evolution of damage in laminated composites.

$$W_T=W+W_S$$

(33)

According to this theory, the total work $W_T$is composed of inelastic work $W_S$ and deformation work $W$. The portion of the total energy that is not returned ($W_S$) can be obtained from the material’s own response according to the stress-strain relationship displayed [7]. To satisfy the path independence condition for total work, these internal state variables ($S_{i’s})$ must be compliant to equation (35) [7].

By definition,

$$f_i\equiv -\frac{\partial W_{}}{\partial S_i}$$

(34)

$$f_i=\frac{\partial W_T}{\partial S_i}$$

(35)

Where, $f_i$ is the so-called thermodynamic force of the ith internal state variable $S_i$, whose value represents the magnitude of force available through a load. When this force exceeds $\frac{\partial W_T}{\partial S_i}$ (necessary intensity), the material undergoes the structural change associated with the variable $S_i$, then, there is $\frac{\partial S_i}{\partial t}>0$ (damage progress related to $S_i$ at a time), otherwise this temporal rate is null [7].

Initially, Schapery considered two internal state variables, one of them for the energy associated with microcracks in the matrix (S) and another for intra-layer transverse cracks (S_c), so, the inelastic work is written as follows [7]:

$$W_s=S+S_c$$

(36)

The study of Basu et al [7] only considers the irreversibility related to S, i. e, the energy dissipated due to the accumulation of microcracks in the matrix. The effects on the elastic modulus E₂₂ and shear modulus G₁₂ are given by:

$$E_{22}=E_{220}{e}_S\left(S\right)$$

(37)

$$G_{12}=G_{120}{g}_S\left(S\right)$$

(38)

In which $E_{220}$ and $G_{120}$ are the respective elastic moduli for the intact material, $e_S\left(S\right)$ and $g_S\left(S\right)$ are polynomial functions of the internal state variable S associated with microcracks [7].

Regarding the strain energy density or strain work, Basu et al [7] write:

$$W=\frac{1}{2}\left(C_{11}{\epsilon }_{11}^2+C_{22}{\epsilon }_{22}^2\right)+C_{12}{\epsilon }_{11}{\epsilon }_{22}+\frac{1}{2}{C}_{66}{\gamma }_{12}^2$$

(39)

To consider geometric nonlinearities, the Green strain tensor and the second Piola-Kirchoff tensor must be used in equation (38). It can be demonstrated that [7]:

$$\frac{{\epsilon }_{11}^2}{2}\frac{\partial C_{11}}{\partial S}+\frac{{\epsilon }_{22}^2}{2}\frac{\partial C_{22}}{\partial S}+\frac{{\gamma }_{12}^2}{2}\frac{\partial G_{12}}{\partial S}+{\epsilon }_{11}{\epsilon }_{22}\frac{\partial C_{12}}{\partial S}=-1$$

(40)

In [7], the physical hypotheses of the model consider${\nu }_{12}{\nu }_{21}\ll 1,$so, $C_{ij}=E_{ij}$. It is assumed that the stiffness in the direction of the fibers is not influenced by S, then, $\frac{\partial E_{11}}{\partial S}\cong 0$ and also that ${\nu }_{12}$ varies while remaining $C_{12}={\nu }_{12}{E}_{22}$ constant, thus, the equation (39) for damage evolution reduces to [7]:

$$\frac{{\epsilon }_{22}^2}{2}\frac{\partial E_{22}}{\partial S}+\frac{{\gamma }_{12}^2}{2}\frac{\partial G_{12}}{\partial S}=-1$$

(41)

Obviously, in inelastic processes the rate of entropy generation is non-negative ($\frac{\partial S_{}}{\partial t}\ge 0$), given that physically a reversible damage is not allowed by the assumed mechanism [7].

Once established the theoretical basis for the model, Basu et al [7] suggest that the microbuckling phenomenon in fibers under compression load is accounted for by the rotation of the fiber at a given point in the material.

$${\gamma }_{12}={\gamma }_{12}^{\infty }+\phi -\beta$$

(42)

Where, ${\gamma }_{12}^{\infty }$ is the contribution of applied external shear, if any. The angular variables $\phi$ and $\beta$can change as a function of the far-field stresses [7].

After implementing of subroutine in ABAQUS and the good approximation of the simulation results with experimental data, the authors concluded that this model dispels the notion of a fixed compressive strength, specifically, it uses the mechanism of the failure process to provide the local compressive strength in given point of a lamina [7].

When studying specific models for determining effective elastic constants in fabric-reinforced composites, the modeling developed by Vasiliev and Morozov stands out. [8], which describe the typical structuring of knitted fabrics as a fabric in which warp yarns and weft or fill yarns are woven and to demonstrate the influence of the angle (α) on the modulus of elasticity of the material, the researchers consider a tension in the longitudinal direction and the rule of mixture for the fabric itself. It is worth mentioning that the magnitude of this elastic modulus of the fabric is, as a general rule, tens of times greater than that of polymer matrices, therefore, the authors assume it to be the apparent elastic modulus of the composite itself in that direction [8].

$$E_l{A}_l=E_{warp}{A}_{warp}+E_{weft}{A}_{weft}$$

(43)

Where, $A_l$ is the apparent cross-sectional area, oriented with normal in the longitudinal direction. $A_{weft}$ is the area of the section containing the transverse yarns (weft) and $A_{warp}$ is the area of the section with longitudinal yarns (warp). After defining the geometric measurements associated with each area ($t_1$ e $t_2$), the longitudinal elastic modulus of the composite is expressed as follows [8]:

$$E_l=\frac{E_2}{2}+\frac{E_1({4t}_1+t_2)}{4\left({2t}_1+t_2({cos}^4\alpha +\frac{E_1}{E_2}{sen}^4\alpha +\left(\frac{E_1}{G_{12}}-2{\nu }_{21}\right){sen}^2\alpha {cos}^2\alpha )\right)}$$

(44)

It is noted that the fibers that make up the transverse direction yarns are orthogonal to the loading, so, it is assumed that the transverse modulus of elasticity in the fabric is equal to the transverse Young’s modulus in a unidirectional composite, i.e, $E_{weft}=E_2$ and $E_1$ is the longitudinal elastic modulus of a unidirectional composite reinforced by the same fibers [8].

In according to Vasiliev and Morozov [8], if the same fibers in a generic arrangement of woven fabric form a cross-ply [0/90] composite, for example, they will give it an effective longitudinal elastic modulus greater than that of the woven fabric, just because the fibers, in cross-ply, are not under the pre-bending as in the woven fabric.

Despite this general decrease in the stiffness and resistance of composites reinforced by woven fabrics in relation to those of the [0/90] arrangement, in two particular types of weaves, practically the same properties are observed as in that cross-ply composites. These are plain and twill weaves. In plain weave, very common in the formation of woven fabrics, the longitudinal fibers are intertwined with the transverse fibers, passing repeatedly above and below them, so that, globally, there is a tendency to cancel out these pre-bendings of the fibers caused by the weaving [8].

Figure 4. Photo taken by the first author. Illustrates the plain weave profile in fabric made with ramie fibers.

Figure 4. Photo taken by the first author. Illustrates the plain weave profile in fabric made with ramie fibers.

The Bridging model works with a tensor $A_{ij}$(bridging tensor), which represents how able the continous phase is in sharing the external loading to the composite in relation to the dispersed phase [9], in other words:

$${\sigma }_i^m=A_{ij}{\sigma }_i^f$$

(45)

Therefore, the effective elastic compliance tensor of the composite is determined by [9]:

$$S_{ij}=\left(c_f{S}_{ij}^f+c_m{S}_{ij}^m{A}_{ij}\right){(c_{f}I+c_m{A}_{ij})}^{-1}$$

(46)

$$A_{ij}=\left[\begin{array}{ccc}\begin{array}{cc}{a}_{11}& a_{12}\\ 0& a_{22}\end{array}& \begin{array}{cc}{a}_{12}& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}{a}_{33}& 0\\ 0& a_{44}\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}{a}_{55}& 0\\ 0& a_{66}\end{array}\end{array}\right],\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mathrm{i}\mathrm{j}=1,2,...6$$

(47)

In which, $a_{ij}\ne 0$ and they can be easily consulted, as well as specifications for their parameters, in the Wang and Huang paper [9].

Chamis is a model modified from the mixture rule, in which the longitudinal elastic modulus and in-plane Poisson’s ratio are identical to those determined by the aforementioned rule, the other constants are mathematically expressed as follows [9]:

$$E_{22}=\frac{E^m}{1-\sqrt{c_f}\left(1-\frac{E^m}{E_{22}^f}\right)}$$

(48)

$$G_{12}=\frac{G^m}{1-\sqrt{c_f}\left(1-\frac{G^m}{G_{22}^f}\right)}$$

(49)

$$G_{23}=\frac{G^m}{1-\sqrt{c_f}\left(1-\frac{G^m}{G_{23}^f}\right)}$$

(50)

Concluding the initial approach of this research upon fundamental micromechanical models for laminated composites, it is necessary to comment about some numerical models that are used in the calculations of the homogenizations proposed by micromechanical modeling. Well, a micromechanical finite element analysis (FEM) has acceptable computational costs in linear elasticity problems, however, when the matrix phase material presents non-linear behaviors, some other specific methods become more advantageous, such as: Voronoi cellular finite elements (VCFEM), the generalized method of cells (GMC), finite volume direct averaging micromechanics (FVDAM) and the asymptotic method for unit cell homogenization (VAMUCH) [9].

Confronting each othersome of the micromechanical models described previously, Table 1 below, in which the total average error involved in predicting five elastic constants ($E_{11},E_{22},G_{12},G_{21},{\nu }_{12}$) for each of these models is listed in ascending order in relation to their respective lack of accuracy [9].

It is worth noting that for the longitudinal elastic modulus ($E_{11}$) all models provide practically the same values, neverthless, for the other constants, it is observed that some approaches are sometimes superior, sometimes inferior compared to each other [9]. Therefore, the choice of one model or another to predict values of effective elastic constants in a composite must involve consideration of the degree of relevance for those constants that are more and less important in view of the final application of the material to be developed.

Table 1. Total average error of some micromechanical models in predicting five elastic constants for composite materials [9].

model	Average Error	Rank
Bridging	10.38%	1
FVDAM	12.83%	2
FEM	13.08%	3
Chamis	14.09%	4
GMC	15.07%	5
Halpin-Tsai	19.24%	6
Mori-Tanaka	19.59%	7
Rule of Mixture	28.40%	8
Eshelby	30.72%	9

Table 1. Total average error of some micromechanical models in predicting five elastic constants for composite materials [9].

model	Average Error	Rank
Bridging	10.38%	1
FVDAM	12.83%	2
FEM	13.08%	3
Chamis	14.09%	4
GMC	15.07%	5
Halpin-Tsai	19.24%	6
Mori-Tanaka	19.59%	7
Rule of Mixture	28.40%	8
Eshelby	30.72%	9

2.2. Advances and Recent Researches

In order to include the failure mechanisms of fiber-reinforced structural composites in processes of numerical simulation, a reformulated high-fidelity generalized method of cells (HFGMC) is used with a finite element code. This approach, at the microstructural level, uses the model of cells and subcells to calculate the stress field within a repeating unit cell (RUC) of the microstructure of the composite [10]. The displacement field in the subcell is approximated by the following Legendre polynomial [10]:

$$u_i^{b,c}={\epsilon ’}_{ij}{x}_j+W_{i\left(00\right)}^{b,c}+{y^{*}}_2^b{W}_{i\left(10\right)}^{b,c}+{y^{*}}_3^c{W}_{i\left(01\right)}^{b,c}+\frac{1}{2}\left(3{{y^{*}}_2^b}^2-\frac{h_b²}{4}\right)W_{i\left(20\right)}^{b,c}+\frac{1}{2}\left(3{{y^{*}}_3^c}^2-\frac{l_c²}{4}\right)W_{i\left(02\right)}^{b,c}$$

(51)

Where l and h are dimensional terms of the subcells, $y_2^{*}$ and $y_3^{*}$ are the coordinates of the subcell in the local coordinate system. The term ${\epsilon ’}_{ij}{x}_j$ represents the average strain contribution (homogenized) that is identical for all subcells in the unit cell. The remaining terms represent the fluctuations of the displacement fields that are distributed throughout the unit cell [10].

The solution of the global system of equations allows the calculation of the displacement field within the unit cell. After this obtained solution, the microvariables W can be calculated allowing to determine the components of the strain tensor in each subcell [10]. From this, the strain concentration tensor, the average stress field in the subcell, as well as the effective or equivalent stiffness tensor can be determined [10].

Meshi et al. [11] also use HGFMC to predict the behavior of a cross-ply laminated composite. In this research, the authors work with experimental data from specific literature and their own experiments in order to obtain mechanical properties of the continuos and dispersed phases. Therefore, applying linear and non-linear constitutive laws, solutions can be determined for the global and local systems of the studied laminate and conclude that modeling using the high-fidelity generalized method of cells is able to simulate the non-linear behavior of soft composites (low stiffness to bending ~10⁰ MPa) [12].

The micromechanical theory of the generalized cell method, using the continuum damage mechanics model of mixed-mode and multiaxial for the composite matrix search to predict the evolution of micro damage in a simple unidirectional polymer matrix laminated composite [13]. This modeling considers that the nucleation of damage in each subcell is determined by quadratic functions called damage strain (${\epsilon }_i^D$), in which the subindex i indicates the direction of the material, with i = 1 being the direction of alignment of the fiber with the lamina, so [13]:

$${\epsilon }_1^D={\left({\left(\frac{{\epsilon }_{11}}{X_{\epsilon }}\right)}^2+{\left(\frac{{\gamma }_{12}}{S_{\epsilon }}\right)}^2+{\left(\frac{{\gamma }_{13}}{S_{\epsilon }}\right)}^2\right)}^{1/2}$$

(52)

Where ${\epsilon }_{11}$ is the normal engineering strain component, ${\gamma }_{12}$ and ${\gamma }_{13}$ are the corresponding shear strains. The two other damage strains are defined by analogy to equation (52), $X_{\epsilon }$ and $S_{\epsilon }$ are defined as the allowable normal and shear strains, respectively. The model assumes that damage begins when any of the damage strains exceed unity, that is [13]:

$${\epsilon }_i^D\ge 1,\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{i}\in \{1,2,3\}$$

(53)

The researchers in [13] still work with Schapery’s theory for modeling microdamages at the lamina level. As already discussed in this research, Schapery’s theory assumes the non-occurrence of irreversible strains in the matrix, therefore, the variable S is associated with the work potential that contributes to structural changes that generate microdamages in the material.

Two studies carried out by Ivancevic and Smojver [14], as well as by Massarwa et al [15], in which researchers also use continuum damage mechanics and nonlinear multiscale damage modeling, respectively, to modify the HGFCM and account for, in properties predicted for study composites, the progressive damage due to mechanical stress.

The method of N-layer concentric cylinders (NCYL) can be used to model the microstructure of unidirectional laminae in laminated composites [16]. The assemblage of fibers within the ply can be assumed to be a hexagonal arrangement in the transverse plane, in which all fibers have identical cross sections, or a random arrangement, where the fibers have different diameters, however, the volume fraction is maintained. In both cases, the fibers are considered infinitely long and the effective response of the composite is macroscopically homogeneous and transversely isotropic, requiring only five independent constants to form the composite stiffness tensor [16].

$$C^c=\left[\begin{array}{cccccc}{E}_1^c+4{\nu }_{12}^c²K_{23}^c& 2{\nu }_{12}^c{K}_{23}^c& 2{\nu }_{12}^c{K}_{23}^c& 0& 0& 0\\ 2{\nu }_{12}^c{K}_{23}^c& K_{23}^c+G_{23}^c& K_{23}^c-G_{23}^c& 0& 0& 0\\ 2{\nu }_{12}^c{K}_{23}^c& K_{23}^c-G_{23}^c& K_{23}^c+G_{23}^c& 0& 0& 0\\ 0& 0& 0& G_{23}^c& 0& 0\\ 0& 0& 0& 0& G_{12}^c& 0\\ 0& 0& 0& 0& 0& G_{12}^c\end{array}\right]$$

(54)

Although the choice of these constants is not unique, the researchers worked in [16] with $E_1^c$ (axial modulus of elasticity), ${\nu }_{12}^c$ (axial Poisson’s ratio), $G_{12}^c$ (in plane shear modulus), $K_{23}^c$ (plane strain bulk modulus), $G_{23}^c$ (out-of-plane shear modulus) [16].

Basically, this multilayer model considers a core of linear and transversely isotropic elastic fibers wrapped in N-1 concentric cylindrical layers of matrix [16]. The stiffness tensor ($C^f)$of the fibers can be written with fibers properties and has the following form:

$$C^f=\left[\begin{array}{cccccc}{C}_{11}^f& C_{12}^f& C_{12}^f& 0& 0& 0\\ C_{12}^f& C_{22}^f& C_{23}^f& 0& 0& 0\\ C_{12}^f& C_{23}^f& C_{22}^f& 0& 0& 0\\ 0& 0& 0& C_{44}^f& 0& 0\\ 0& 0& 0& 0& C_{55}^f& 0\\ 0& 0& 0& 0& 0& C_{55}^f\end{array}\right]$$

(55)

To analysis of microdamage in the matrix, the authors in [16] consider the global non-linear behavior of the stress-strain relationship exhibited by polymer composites. Although this nonlinearity is due to the coalescence of microdamages that evolve into macroscopic cracks in the continous phase of the composite, the researchers mathematically model the problem according to the J₂ theory of plastic strain through a secant modulus approximation. In order to use tensile test results, for example, to determine the damage response of the matrix phase material when subjected to combined loading, the equivalent stress and strain in the matrix (${\sigma }_{eq}^m\mathrm{a}\mathrm{n}\mathrm{d}{\epsilon }_{eq}^m$) can be determined by [16]:

$${\sigma }_{eq}^m=\sqrt{\frac{1}{2}\left[{\left({\sigma }_{11}^m-{\sigma }_{22}^m\right)}^2+{\left({\sigma }_{22}^m-{\sigma }_{33}^m\right)}^2+{\left({\sigma }_{33}^m-{\sigma }_{11}^m\right)}^2\right]+3\left[{\left({\sigma }_{12}^m\right)}^2+{\left({\sigma }_{13}^m\right)}^2+{\left({\sigma }_{23}^m\right)}^2\right]}$$

(56)

$${\epsilon }_{eq}^m=\frac{1}{1+{\nu }_s^m}\sqrt{\frac{1}{2}\left[{\left({\epsilon }_{11}^m-{\epsilon }_{22}^m\right)}^2+{\left({\epsilon }_{22}^m-{\epsilon }_{33}^m\right)}^2+{\left({\epsilon }_{33}^m-{\epsilon }_{11}^m\right)}^2\right]+\frac{3}{4}\left[{\left({\epsilon }_{12}^m\right)}^2+{\left({\epsilon }_{13}^m\right)}^2+{\left({\epsilon }_{23}^m\right)}^2\right]}$$

(57)

$${\nu }_s^m=\frac{1}{2}+\frac{E_s^m}{E_e^m}\left({\nu }_e^m-\frac{1}{2}\right)$$

(58)

Where, ${\nu }_s^m,E_s^m=\frac{{\sigma }_{eq}^m}{{\epsilon }_{eq}^m},E_e^m,{\nu }_e^m$are the secant Poisson’s ratio of the matrix, the secant elastic modulus, elastic modulus and the Poisson’s coefficient of the matrix, respectively [16]. Finally, there is the following non-linear relationship between stress and strain [16]:

$${\sigma }_{eq}^m={\sigma }_y^m-\frac{k_1}{k_2}\left(e^{{-k}_2{\epsilon }_{eq}^m}-e^{{-k}_2\frac{{\sigma }_y^m}{E_e^m}}\right)$$

(59)

As ${\sigma }_y^m$ the yield stress of the matrix (from tensile testing), $k_1\mathrm{a}\mathrm{n}\mathrm{d}{k}_2$ are two parameters that govern the progress of microdamages in the matrix, which can be determined from the experimental response to the stress-strain relationship [16].

A micromechanical analytical model for predicting effective elastic properties in multidirectional laminated composites with wavy defect layers is described and analyzed in relation to previously published numerical and experimental data [17]. The analytical model proposed by this research begins to be constituted as follows:

$${\psi }_{uw}\left(x\right)=A\cos \left(\frac{2\pi }{L}x\right);\left|x\right|\le \frac{L}{2}$$

(60)

With A being the amplitude and L the wavelength associated with the extension of the representative volume elemento (RVE) in the longitudinal x direction. The equation (60) describes wavy and uniform shapes throughout the thickness of the laminate [17].

$${\psi }_{gw}\left(x,z\right)=\left\{\begin{array}{c}A\left(\frac{\frac{H_w}{2}-\left|z\right|}{\frac{H_w}{2}}\right)\cos \left(\frac{2\pi }{L}x\right);\left|x\right|\le \frac{L}{2},\left|z\right|\le \frac{H_w}{2}\\ 0;\left|x\right|\le \frac{L}{2},\left|z\right|>\frac{H_w}{2}\end{array}\right.$$

(61)

While equation (61) mathematically models the gradual, non-uniform ripples in the laminate thickness. It becomes evident that the functions in (60) and (61) change the z coordinate in a ground laminate to [17]:

$$z_w=\left\{\begin{array}{c}{z}_0+{\psi }_{uw}\left(x\right);\mathrm{for}\mathrm{uniform}\mathrm{ripples}\\ z_0+{\psi }_{gw}\left(x,z\right);\mathrm{for}\mathrm{gradual}\mathrm{ripples}\end{array}\right.$$

(62)

Where $z_w$ is the z coordinate of the lamina points in a wavy layer configuration for the composite. Since the mid-plane of each rectified lamina is parallel to the xy plane, the angular orientations of the point of the corrugations in relation to the xz plane, for uniform and gradual ripples, are respectively expressed by equations (63) and (64) [17]:

$$\phi ={\tan }^{-1}\left(\frac{d{\psi }_{uw}\left(x\right)}{dx}\right);\left|x\right|\le \frac{L}{2}$$

(63)

$$\phi =\left\{\begin{array}{c}{\tan }^{-1}\left(\frac{d{\psi }_{gw}\left(x,z_0\right)}{dx}\right);\left|x\right|\le \frac{L}{2},\left|z\right|\le \frac{H_w}{2}\\ 0;\left|x\right|\le \frac{L}{2},\left|z\right|>\frac{H_w}{2}\end{array}\right.$$

(64)

Takeda [17] performs a homogenization process of the representative volume element (RVE) in two steps. As a first step, the researcher considers that the out-of-plane xy $\left({\sigma }_{zz},{\sigma }_{zy},{\sigma }_{zx}\right)$ and the strains in-plane xy $\left({\epsilon }_{xx},{\epsilon }_{yy},{\epsilon }_{xy}\right)$ are uniform in discrete strips in the RVE, see Figure 5. From this discretization, the effective stiffness tensor for each strip is represented by contracted equation below:

$$C^{*}=C_{ij}^{*}\left(\text{6x6}\mathrm{matrix}\right)$$

(65)

In which the upper index (*) indicates that it is an effective or average quantity in each strip.

The second homogenization step assumes that stresses $\left({\sigma }_{yy}^{*},{\sigma }_{zz}^{*},{\sigma }_{yz}^{*}\right)$ and strains $\left({\epsilon }_{yy}^{*},{\epsilon }_{zz}^{*},{\epsilon }_{yz}^{*}\right)$ are continuous across the strip interfaces. In this way, the effective or average stiffness tensor for the homogenized representative volume element (Figure 6) is represented by a 6x6 matrix and indicated in a contracted form by equation (66) [17]:

$$C^{**}=C_{ij}^{**}$$

(66)

Where the upper index (**) indicates an effective or average quantity for the RVE.

Figure 5. N strips of the first homogenization step of the representative volume element of the composite with ripples in the constituent layers [17].

Figure 5. N strips of the first homogenization step of the representative volume element of the composite with ripples in the constituent layers [17].

Figure 6. Representative volume element homogenized according to the second step [17].

Figure 6. Representative volume element homogenized according to the second step [17].

Na analysis combining the so-called micromechanical theory of failure with a cohesive zone model, in which a criterion for the initiation of interlaminar failures is established is the research proposal [18]. According to the authors of this work, to better understand the theory used in this study it is necessary to know its basic hypotheses, which are [18]:

• The fibers are assumed to be linear elastic, transversely isotropic and with uniform distribution.
• The matrix is considered linear elastic, isotropic and without voids.
• Every micro strain/stress component at any fiber or matrix point is linearly dependent on macro strains/stresses and on the temperature increment.

  $${\sigma }_i=M_{\sigma }^i{\sigma }^{’}+A_{\sigma }^i\Delta T$$

  (67)

Where, ${\sigma }_i$ is the micro stress at point i (has six stress components), ${\sigma }^{’}$ is the macro stress in the periodic unit cell representative of the composite (has six stress components), $M_{\sigma }^i$ is the amplification factor of the mechanical stress at point i (6x6 matrix ), $A_{\sigma }^i$ is the thermal stress amplification factor at point i (6x1 matrix) [18].

Regarding the failure criteria used for fiber and matrix, the researchers in [18] emphasize that fibers are the main components of mechanical strength of laminated polymer composites and that damage to them is the primary indicator for the failure of the material when subjected to longitudinal tension and compression stresses. When longitudinal tension occurs, the stress componente ${\sigma }_{11}^f$ of each critical point in the fiber is adopted as a parameter to evaluate whether or not the reinforcing element failure, such as:

$$\frac{{\sigma }_{11}^f}{{\sigma }_{cr}^t}=1$$

(68)

And for the compression criterion, it has:

$$\frac{{\sigma }_{eq}^f}{{\sigma }_{cr}^c}=1$$

(69)

In which, ${\sigma }_{eq}^f$ it is analogous to equation (56) for fibers,${\sigma }_{cr}^t$ and ${\sigma }_{cr}^c$ are the critical values of ${\sigma }_{11}^f$ and ${\sigma }_{eq}^f$ for when fiber failures occur under tensile and compressive loads [18].

In continuity, the model for the progressive decrease in the effective compliance tensor as damage to the composite increases and a criterion for the strength to delamination of the material are stand out, they are also addressed in the researches [18,19]. In the first case, from the expression of the strain energy density for the damaged composite, equation (70), it can be written as:

$${En}_d=\frac{1}{2}\left[\frac{{\sigma }_{11}^2}{E_1(1-d_f)}+\frac{{\sigma }_{22}^2}{E_2(1-d_m)}+\frac{{\sigma }_{33}^2}{E_3(1-d_m)}-\frac{2{\sigma }_{11}{\sigma }_{22}{\nu }_{12}}{E_1}-\frac{2{\sigma }_{11}{\sigma }_{33}{\nu }_{13}}{E_1}-\frac{2{\sigma }_{22}{\sigma }_{33}{\nu }_{23}}{E_2}+\frac{{\sigma }_{12}^2}{G_{12}(1-d_{12})}+\frac{{\sigma }_{13}^2}{G_{13}(1-d_{13})}+\frac{{\sigma }_{23}^2}{G_{23}(1-d_{23})}\right]$$

(70)

Where, $d_{ij}$ are the shear damage parameters in the ij plane [18,19]. So:

$$S_{ijkl}=\frac{\partial {En}_d}{\partial {\sigma }_{ij}}=\left[\begin{array}{ccc}\begin{array}{cc}\frac{1}{E_1\left(1-d_f\right)}& -\frac{{\nu }_{12}}{E_1}\\ -\frac{{\nu }_{12}}{E_1}& \frac{1}{E_2\left(1-d_m\right)}\end{array}& \begin{array}{cc}-\frac{{\nu }_{13}}{E_1}& \begin{array}{c}0\\ \end{array}\\ -\frac{{\nu }_{23}}{E_2}& 0\end{array}& \begin{array}{cc}\begin{array}{c}0\\ \\ 0\end{array}& \begin{array}{c}0\\ \\ 0\end{array}\\ & \end{array}\\ \begin{array}{cc}& \\ & \end{array}& \begin{array}{cc}\frac{1}{E_3\left(1-d_m\right)}& 0\\ \begin{array}{c}\\ \end{array}& \frac{1}{G_{12}\left(1-d_{12}\right)}\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ symmetry& & \begin{array}{cc}\frac{1}{G_{13}\left(1-d_{13}\right)}& 0\\ & \frac{1}{G_{23}\left(1-d_{23}\right)}\end{array}\end{array}\right]$$

(71)

About the criterion for the initiation of delamination (cohesive zone model), it can be expressed as [19]:

$${\left(\frac{〈{\sigma }_n〉}{N}\right)}^2+{\left(\frac{{\sigma }_s}{S}\right)}^2+{\left(\frac{{\sigma }_t}{T}\right)}^2=1$$

(72)

In which, ${\sigma }_n,{\sigma }_s,{\sigma }_t$are the normal and shear stresses, respectively, associated with the respective delamination crack opening modes, while N, S and T are the respective strengths of the interfaces to each of the aforementioned modes. The parentheses “< >” indicate that under pure compression delamination does not begin [19].

For failure analysis in curved laminated composites subjected to bending, Andraju and Raju in [20] use micromechanical models to describe some intra-laminar (matrix failure) and inter-laminar (delamination) damages. According to the research, two models are used to analyze the initiation of failure in the composite matrix:

$$F_{mt}=\sqrt{{\left(\frac{{\epsilon }_{33}}{\frac{Z_T}{C_{33}}}\right)}^2+{\left(\frac{{\gamma }_{13}}{\frac{{\sigma }_{c13}}{C_{55}}}\right)}^2+{\left(\frac{{\gamma }_{23}}{\frac{{\sigma }_{c23}}{C_{66}}}\right)}^2}\ge 1$$

(73)

Where, $F_{mt}$ is the parameter for matrix failure under tension when ${\epsilon }_{33}>0$, $Z_T$ is the tensile strength in the thickness/depth direction, ${\sigma }_{c13}\mathrm{a}\mathrm{n}\mathrm{d}{\sigma }_{c23}$ are shear strengths of the laminate [20,21].

$$F_{mt}=\frac{\left({\sigma }_{22}+{\sigma }_{33}\right)²}{Y_T²}+\frac{\left({\sigma }_{23}^2-{{\sigma }_{22}\sigma }_{33}\right)}{{\sigma }_{c23}²}+\frac{\left({\sigma }_{12}^2+{\sigma }_{13}²\right)}{{\sigma }_{c12}²}\ge 1$$

(74)

In equation (74), $F_{mt}$ is the parameter for matrix failure also under tension, but when ${\sigma }_{22}+{\sigma }_{33}>0$, $Y_T$ it indicates the tensile strength in the direction transverse to the fibers and ${\sigma }_{c12}$ is the shear strength in-plane of the laminae [20,21].

Still in the study [20], the authors use an exponential model for the progressive evolution of damage expressed as follows:

$$d_{mt}=1-\frac{1}{F_{mt}}exp\left(\frac{C_{ij}{\epsilon }_i²(1-F_{mt})L_c}{G_{mt}}\right)$$

(75)

Where $L_c$is the characteristic length of the representative volume elemento (RVE) of the laminate and $G_{mt}$ is the fracture energy corresponding to the failure of the matrix subjected to tension [20], thus:

$$C^e=C_{ij}^e \left(\mathrm{symmetrical}\text{6x6}\right)$$

(76)

It is the elastic stiffness tensor accounting for progressive damage $d_{mt}$, in which [20]:.$C_{11}^e=C_{11};C_{12}^e={(1-d_{mt})C}_{12};C_{13}^e={(1-d_{mt})C}_{13};C_{22}^e={(1-d_{mt})²C}_{22};C_{23}^e={(1-d_{mt})²C}_{23};C_{33}^e={(1-d_{mt})²C}_{33};C_{44}^e=C_{44};C_{55}^e={(1-d_{mt})C}_{55};C_{66}^e={(1-d_{mt})²C}_{66}$; with other symmetric terms being null.

Therefore, the stress-strain relationship can be expressed conditionally as in equation (77) below [20]:

$$=C\epsilon \left\{\begin{array}{c}se(0<d_{mt}\le 1;C=C^e\\ se\left(d_{mt}=0\right);C=C^0\end{array}\right.$$

(77)

Where, $C^0$ is the stiffness tensor for the undamaged composite. Finally, the authors in [20] use the cohesive zone model in equation (72) as a criterion to evaluate the opening of delamination cracks considering the stresses for mode I (normal) and modes II and III (shear).

Micromechanical analyzes considering progressive damage in laminated composites have been explored in more recent research. As a final emphasis, Liu et al [22] develop a study analyzing progressive damage to compression of laminated composites with stress concentrations. The authors propose a subroutine element of the commercial software ABAQUS that is constituted using the third-order shear deformation theory in conjunction with a two-dimensional progressive damage model and micromechanical criteria for fiber and matrix failure analyses. In order to verify the validity of the modeling, compression tests on laminated composites with cavities were carried out and the failure load predicted by theoretical modeling differed only 6.1% from that presented in the experimental analysis [22].

An interesting study conducted by Hu et al [23] propose a peridynamic micromechanical model (PD models), thus, based on the non-local continuum mechanics to analyze damage mechanisms in laminated composites. According to these authors, mathematical modeling is efficient in comparison to experimental literature data in terms of predictions that take into account important micromechanical characteristics and the initiation of damage. In this way, it makes it possible to investigate both the effective elastic properties and progressive mechanisms of deterioration of the material with good accuracy.

More recently, Li et al [24] developed a work in which they propose a more computationally efficient peridynamic micromechanical model (PD models) to characterize the strain and progressive damage behaviors in a laminated composite modeled with a single ply of material points. At the end of the research, as described by the authors, the alternative peridynamic model actually appears to have lower computational cost than conventional modeling using non-local continuum mechanics.

A study which presents micromechanical modeling for the fiber-matrix interface is developed by Lei et al [25]. As described by these researchers, the phenomenon of interfacial stress transfer is a relevant mechanical problem presented by fiber-reinforced composites. This behavior involves four successive stages, these are: when the interface is still intact and unites fiber and matrix, the process of detachment and separation between the two phases, the interface with fibers and matrix completely separated and the removal of fibers from within the matrix by tension loading (pull out). A model for the transfer of elastic stress can be developed through the balance of forces on the fiber and expressed as follows [25]:

$$\mathit{\tau }=-\frac{\mathit{r}}{2}\left(\frac{\mathit{d}\mathit{\sigma }}{\mathit{d}\mathit{x}}\right)$$

(78)

Where, $\mathit{\tau }$ is the shear stress along the fiber immersed in the matrix given by the so-called Cox model, $\mathit{\sigma }$ is the axial stress in the fiber and r is the fiber radius [25].

According to [25] the so-called Piggott model is used for the distribution of elastic axial stress in the fiber before interface separation, which is given by:

$$\mathit{\sigma }={\mathit{\sigma }}_{\mathit{a}\mathit{p}\mathit{p}}\frac{\mathit{s}\mathit{i}\mathit{n}\mathit{h}\left[\mathit{n}(\mathit{L}-\mathit{x})/\mathit{r}\right]}{\mathit{s}\mathit{i}\mathit{n}\mathit{h}\left(\mathit{n}\mathit{s}\right)}$$

(79)

In (79), ${\mathit{\sigma }}_{\mathit{a}\mathit{p}\mathit{p}}$ it is the stress acting on the fiber section outside the matrix, x is the distance to the fiber entry into the matrix, L is the effective stress transfer length, that is, it is the fiber length at which the axial stress will decay to zero, s is the fiber aspect ratio (L/r) and n is a parameter related to the geometry and material of the fiber and matrix [25].

Deriving the equation (79) using equation (78), it can be written [25]:

$$\mathit{\tau }={\mathit{\sigma }}_{\mathit{a}\mathit{p}\mathit{p}}\frac{\mathit{n}\mathit{c}\mathit{o}\mathit{s}\mathit{h}\left[\mathit{n}(\mathit{L}-\mathit{x})/\mathit{r}\right]}{2\mathit{s}\mathit{i}\mathit{n}\mathit{h}\left(\mathit{n}\mathit{s}\right)}$$

(80)

Which $\mathit{\tau }$ is the interfacial shear stress in the fiber section immersed in the matrix [25].

If continuous strains are applied, failure due to fiber fracture and failure due to detachment between fiber and matrix occur, respectively, when [25]:

$$\mathit{\sigma }=\frac{2{\mathit{\tau }}_{\mathit{m}}}{\mathit{n}}\ge {\mathit{\sigma }}_{\mathit{b}}$$

(81)

$${\mathit{\tau }}_{\mathit{m}}=\frac{\mathit{n}\mathit{\sigma }}{2}\ge {\mathit{\tau }}_{\mathit{b}}$$

(82)

With ${\mathit{\tau }}_{\mathit{m}}$ being the maximum interfacial shear stress (for which debonding between fiber and matrix occurs), ${\mathit{\sigma }}_{\mathit{b}}$ and ${\mathit{\tau }}_{\mathit{b}}$ are the limit of fiber tensile strength and interfacial shear strength, respectively [25].

As the load increases in a pull out test, detachment initially occurs at the entry point of the fiber into the matrix (x = 0) and propagates along the interface, using the two previously mentioned stress models (Cox and Piggott), the authors of the research [25] also make demonstrations related to the transfer of shear stress due to friction in the detachment region between the continuos and dispersed phases.

Well, phenomena associated with delamination such as fiber bridging (Figure 7) can also be modeled by micromechanical analyses. According to Daneshjoo et al [26], the main micromechanisms involved in this phenomenon are three:

• Fragmentation of the matrix (Matrix spalling): Before fiber bridging occurs, the fibers are peeled-off from the matrix due to the opening of delamination cracks. After fiber bridging, the ends of the fiber section peeled from within the matrix continue to be inserted into the matrix and, under new loads, these ends also tend to leave the matrix, however, the continous phase begins to fragment due to the stress concentration generated by the phenomenon of fiber bridging prevail.
• Detachment of fibers (Fiber pull-out): Weak unions with detachments between fiber and matrix contribute to the loads being sufficient to remove the fibers from within the matrix. Then, after fiber bridging, the detachment conditions between the composite phases are favorable to fiber pull out.
• Fiber fracture: Under combined flexural and axial loads, fibers that have undergone the fiber bridging phenomenon tend to break.

The practical effects of fiber bridging end up being useful, given that they improve the fracture toughness of the laminate, in other words, the material becomes more resistant to the crack propagations. Micromechanical models for calculating the absorption energy in the fiber bridging zone, involving each of the three failure micromechanisms cited above, can be analyzed from [26], however, Mirsayar in [27] works on the micromechanical modeling of fiber bridging in a more detailed.

Figure 7. Representation of a fiber bridging phenomenon in delamination crack in the composite material [27].

Figure 7. Representation of a fiber bridging phenomenon in delamination crack in the composite material [27].

In research that uses a criterion based on the stress/energy combination for mixed modes of fracture in laminated composites, taking into account the micromechanical analysis of the fiber bridging [27], the author develops the following approach to this phenomenon:

$${\xi }_{IC,f}={\int }_0^{û}\frac{F\left(u\right)}{dA}du={\int }_0^{û}p\left(u\right)du$$

(83)

Where, ${\mathsf{\xi }}_{\mathrm{I}\mathrm{C},\mathrm{f}}$ it is defined as the fiber bridging tenacity, $\mathrm{F}\left(\mathrm{u}\right)$ it is the strength to the opening u of the crack that is developed along the so-called fiber bridging length â (see Figure 7), $\mathrm{p}\left(\mathrm{u}\right)$ it is defined as the crack closure pressure and û is the opening at the end of the fiber bridging zone. The number of fibers involved in the bridging mechanism is a function of u, therefore, considering this fact, the crack closing pressure can be expressed by [27]:

$$\mathrm{p}\left(\mathrm{u}\right)=\mathrm{f}\left(\mathrm{u}\right).\mathrm{n}\left(\mathrm{u}\right)$$

(84)

In which, $f\left(u\right)$ is the force per unit of fibers and $\mathrm{n}\left(\mathrm{u}\right)$ is the number of fibers per unit area. Using Weibull statistical analysis, $\mathrm{n}\left(\mathrm{u}\right)$ can be given by [27]:

$$n\left(u\right)=\left\{\begin{array}{c}{n}_0,u<u^b\\ n_{0}exp\left(-\frac{C_b}{l_r}{\left(\frac{{\sigma }_b^{*}}{{\sigma }_r}\right)}^m\frac{l}{u}(u-u^b)\right),u\ge u^b\end{array}\right.$$

(85)

Where m is the Weibull modulus, ${\mathrm{u}}^{\mathrm{b}}\cong \frac{{\mathrm{u}}_{\mathrm{p},\mathrm{m}\mathrm{á}\mathrm{x}}}{2}$ is considered the crack opening displacement, below which the number of fiber failures becomes negligible, ${\mathrm{u}}_{\mathrm{p},\mathrm{m}\mathrm{á}\mathrm{x}}$is the crack opening displacement corresponding to the maximum closure pressure, ${\mathsf{\sigma }}_{\mathrm{r}}$ is the stress level in which 63.2% of the fibers along the ${\mathrm{l}}_{\mathrm{r}}$ length fail, ${\mathsf{\sigma }}_{\mathrm{b}}^{\mathrm{*}}$ is the average maximum bending stress at the root of the fibers that were peeled, l is the bridging span and ${\mathrm{C}}_{\mathrm{b}}<1$ consists of a correction factor that is dimensionless [27].

The calculations of p(u) and f(u) must be done considering two steps. In an initial phase $(\mathrm{n}\left(\mathrm{u}\right)={\mathrm{n}}_0)$ and in a steady state, considering that the phenomenon governing micromechanism and the evolution of f(u) correspond to each phase [27]. In the initial phase, after the fibers are peeled from the matrix and before the tensile stress acts on them due to the opening of the crack, the beam behavior is predominant over these reinforcement elements that were released from the layers of the laminated composite [27]. Then, the strength to crack opening is determined from the beam theory, as expressed by equation below [27]:

$$f\left(u\right)=\frac{3\pi E_f{r}^4}{4l_0^2}\left(\frac{u}{l_0}\right)$$

(86)

According to the displacement u of the crack opening increases, a significant normal tensile stress σ(u) is developed in the fiber that has been peeled (see Figure 8) [27]. By equilibrium condition, it is shown that the tensile stress due to the resistive shear τ and the debonding toughness of the fiber/matrix interface $K_{IIC,i}$ in a fiber of radius r can be expressed by [27]:

$$\sigma \left(u\right)=\sqrt{\frac{\tau E_f{l}_0}{r}}\left(\frac{u}{l_0}\right)+\frac{2}{\sqrt{r}}{K}_{IIC,i}$$

(87)

For θ << 1 rad, the total closure pressure is given by [27]:

$$p\left(u\right)=\pi r^{2}n\left(u\right)\left[\left(\frac{3\pi E_f{r}^2}{4l_0^2}+\frac{2}{\sqrt{r}}{K}_{IIC,i}\right)\left(\frac{u}{l_0}\right)+\sqrt{\frac{\tau E_f{l}_0}{r}}{\left(\frac{u}{l_0}\right)}^2\right]$$

(88)

Therefore, the mathematical model for toughness to the fiber bridging phenomenon is determined as follows [27]:

$${\xi }_{IC,f}=\pi r^2{\int }_0^{û}n\left(u\right)\left[\left(\frac{3\pi E_f{r}^2}{4l_0^2}+\frac{2}{\sqrt{r}}{K}_{IIC,i}\right)\left(\frac{u}{l_0}\right)+\sqrt{\frac{\tau E_f{l}_0}{r}}{\left(\frac{u}{l_0}\right)}^2\right]du=\pi r^2{\int }_0^{u^b}{n}_0\left[\left(\frac{3\pi E_f{r}^2}{4l_0^2}+\frac{2}{\sqrt{r}}{K}_{IIC,i}\right)\left(\frac{u}{l_0}\right)+\sqrt{\frac{\tau E_f{l}_0}{r}}{\left(\frac{u}{l_0}\right)}^2\right]du+\pi r^2{\int }_{u^b}^{û}{n}_0\left[\left(\frac{3\pi E_f{r}^2}{4l_0^2}+\frac{2}{\sqrt{r}}{K}_{IIC,i}\right)\left(\frac{u}{l_0}\right)+\sqrt{\frac{\tau E_f{l}_0}{r}}{\left(\frac{u}{l_0}\right)}^2\right]exp\left(-\frac{C_b}{l_r}{\left(\frac{{\sigma }_b^{*}}{{\sigma }_r}\right)}^m\frac{l}{u}(u-u^b)\right)du$$

(89)

In which, during the initial phase (u <$u^b$) only the non-exponential term is considered. So, it is noted here that as u increase, the fiber bridging thoghness also increases in this initial phase.

Experimental studies indicate that there is a limit to the integrity of the matrix in relation to the increase of the crack opening due to delamination [27]. The matrix begins to spall in the steady state and a micromechanical model relating fracture toughness to the maximum normal component of the force exerted by the fiber is indicated in equation (90), where ${\mathsf{\beta }\mathsf{\alpha }}_0$ is the average length of the crack that is parallel to the fiber in the which acts ${\mathrm{f}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ [27].

$$\frac{f_{max}}{{\beta \alpha }_0}=K_{IC}^m\sqrt{\pi {\alpha }_0}$$

(90)

Figure 9. Action of the maximum normal component exerted by the fiber on the matrix and the average length of the crack that fragments the matrix [27].

Figure 9. Action of the maximum normal component exerted by the fiber on the matrix and the average length of the crack that fragments the matrix [27].

Therefore, associating equations (88) and (90), the maximum force generated on the matrix by flexion and traction of the fiber can be described as [27]:

$${f\left(u\right)}_{max}=\pi r^2{\left[\left(\frac{3\pi E_f{r}^2}{4l_0^2}+\frac{2}{\sqrt{r}}{K}_{IIC,i}\right)\left(\frac{u}{l_0}\right)+\sqrt{\frac{\tau E_f{l}_0}{r}}{\left(\frac{u}{l_0}\right)}^2\right]}_{max}={\beta \alpha }_0{K}_{IC}^m\sqrt{\pi {\alpha }_0}$$

(91)

From the equation (91) it is possible to determine the maximum displacement for opening of the crack before the matrix spalls, $u^{sp}$ . Therefore, the fiber bridging toughness for the steady-state phase is determined by equation (92) [27]:

$${\xi }_{IC,f}={\beta \alpha }_0{n}_0{K}_{IC}^m\sqrt{\pi {\alpha }_0}{\int }_{u^{sp}}^{û}exp\left(-\frac{C_b}{l_r}{\left(\frac{{\sigma }_b^{*}}{{\sigma }_r}\right)}^m\frac{l}{u}(u-u^b)\right)du$$

(92)

Can be observed from the equation (92) by this review that as u e l increases (l > u), considering that spalling of the matrix occurs for this stage, the toughness to the fiber bridging phenomenon during the steady state tends to zero.

Closing the analysis about the Mirsayar’s research [27], stands out the modeling of the stress field near the tip of an existing crack in an orthotropic laminated composite:

$$\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{12}\end{array}=\frac{1}{\sqrt{2\pi r}}Re\left[\frac{K_I}{{\mu }_1-{\mu }_2}\left(\begin{array}{c}{\mu }_1{\mu }_2({\mu }_2{M}_2-{\mu }_1{M}_1)\\ ({\mu }_1{M}_2-{\mu }_2{M}_1)\\ {\mu }_1{\mu }_2(M_1-M_2)\end{array}\right)+\frac{K_{II}}{{\mu }_1-{\mu }_2}\left(\begin{array}{c}({\mu }_2{²M}_2-{\mu }_1²M_1)\\ (M_2-M_1)\\ ({\mu }_1{M}_1-{{\mu }_{2}M}_2)\end{array}\right)\right]$$

(93)

Where, ${\sigma }_{ij}$ are the components of the stress tensor in a rectangular coordinate system $x_1{x}_2$ defined along the principal directions of the orthotropic material, $K_I\mathrm{a}\mathrm{n}\mathrm{d}{K}_{II}$ are the stress intensity factors for modes I and II, respectively, $M_1\mathrm{a}\mathrm{n}\mathrm{d}{M}_2$ are terms defined by [27]:

$$M_j=\frac{1}{\sqrt{cos\theta +{\mu }_{j}sen\theta }}\mathrm{c}\mathrm{o}\mathrm{m}j=\mathrm{1,2}$$

(94)

And the parameters ${\mu }_1\mathrm{a}\mathrm{n}\mathrm{d}{\mu }_2$are the pair of complex conjugate roots of the characteristic equation (95) [27]:

$$S_{11}{\mu }^4-2S_{16}{\mu }^3+\left(2S_{12}+S_{66}\right){\mu }^2-2S_{26}\mu +S_{22}=0$$

(95)

Where $S_{ij}$ are the components of the elastic compliance tensor for the orthotropic composite [27].

A dynamic constitutive micromechanical model is suggested by Seyedalikhani et al [28] to predict mechanical behaviors as functions of strain rate in a epoxy matrix laminated composite. According to these authors, the generalized constitutive model for the polymer matrix phase dependent on the strain rate must encompass three aspects. The first one refers to the elastic behavior of polymers when subjected to normal and shear stresses, secondly, the inelastic behavior of polymers, which is prescribed based on the Johnson-Cook constitutive model and finally, a model to predict the strength limit of the polymer [28].

Regarding the elastic behavior of the polymer matrix, the following modeling is presented in study [28]:

$$H\left({\epsilon }^{’}\right)={\alpha }_1{\left(ln{\epsilon }^{’*}\right)}^{{\beta }_1}+{\gamma }_1$$

(96)

In which, ${\alpha }_1,{\beta }_1$ and ${\gamma }_1$ are scaling constants of the material, ${\epsilon }^{’*}=\frac{{\epsilon }^{’}}{{\epsilon }_0^{’}}$, with ${\epsilon }^{’}$ and ${\epsilon }_0^{’}$ the applied and reference strain rates, respectively, and $H$ is the theoretical dynamic elastic modulus [28].

The modified Johnson-Cook ’s constitutive model, equation (97), is used to model the inelastic strains presented by the polymer matrix [28]:

$$\sigma =(A+a{\epsilon }^3+b{\epsilon }^2+c\epsilon )(1+C_1\left(ln{\epsilon }^{’*}\right)+C_2{\left(ln{\epsilon }^{’*}\right)}^2)$$

(97)

Where, $\sigma$ and $\epsilon$ are the plastic stress and strain equivalent to the inelastic phenomenon occurring in polymers, the constants involved are unknowns to be determined in mechanical tests with different strain rates [28].

To the strength of polymer matrix as a function of the strain rate in the polymer, the mathematical model can be described as the equation (98) [28]:

$$S_m\left({\epsilon }^{’m}\right)={\alpha }_2{\left(ln{\epsilon }^{’*}\right)}^{{\beta }_2}+{\gamma }_2$$

(98)

Where, $S_m\left({\epsilon }^{’m}\right)$ is the ultimate tensile strength and ${\alpha }_2{,\beta }_{2,}{\gamma }_2$ are constants of the material [28].

With regard to the tensile strength of the glass fibers also as a function of the strain rate, in [28], the model is determined as follows:

$$S_f\left({\epsilon }^{’f}\right)={\alpha }_3+{{\beta }_3\left({\epsilon }^{’*}\right)}^{{\gamma }_3}$$

(99)

In equation (99), $S_f\left({\epsilon }^{’f}\right)$ is the ultimate tensile strength of a unidirectional composite at the strain rate shared by the fibers and ${\alpha }_3{,\beta }_3$ and ${\gamma }_3$ are material constants to be determined experimentally [28].

With meaningful highlight, the work [28] describes the theory of plasticity at the micromechanical level. According to the authors, most polymer composites reinforced by fibers unidirectionally present relevant inelastic strains before failing under in-plane shear and, in some situations, when subjected to transverse tensile loads too. The following model shows the relationship between the stress increments in the fiber and in the matrix for a unidirectional composite [28]:

$${d\sigma }_i^m=A_{ij}{d\sigma }_j^f$$

(100)

Where, $A_{ij}$ is the bridging tensor whose physical meaning has already been described in the previous section (2.1) of this research.

Still within this non-linear modeling, the elastic and elasto-plastic compliance tensor of the matrix phase is indicated by [28]:

$$S_{ij}^m=\left\{\begin{array}{c}{S}_{ij}^{me},{\tau }_0^m\le \frac{\sqrt{2}}{3}{\sigma }_Y^m\\ S_{ij}^{me}+S_{ij}^{mp},{\tau }_0^m>\frac{\sqrt{2}}{3}{\sigma }_Y^m\end{array}\right.$$

(101)

Where, ${\tau }_0^m$ and ${\sigma }_Y^m$ the octahedral shear stress and the uniaxial yield strength of the matrix phase, respectively,$S_{ij}^{me}$ are the components of the elastic compliance tensor and $S_{ij}^{mp}$ are the components of the plastic compliance tensor given by the matrix equation below [28]:

$$S_{ij}^{mp}=\frac{1}{2M_T^m{\tau }_0^2}\left[\begin{array}{ccc}\begin{array}{cc}{\sigma }_{11}^d{\sigma }_{11}^d& {\sigma }_{22}^d{\sigma }_{11}^d\\ & {\sigma }_{22}^d{\sigma }_{22}^d\end{array}& \begin{array}{cc}{\sigma }_{33}^d{\sigma }_{11}^d& 2{\sigma }_{23}^d{\sigma }_{11}^d\\ {\sigma }_{33}^d{\sigma }_{22}^d& 2{\sigma }_{23}^d{\sigma }_{22}^d\end{array}& \begin{array}{cc}2{\sigma }_{13}^d{\sigma }_{11}^d& 2{\sigma }_{12}^d{\sigma }_{11}^d\\ 2{\sigma }_{13}^d{\sigma }_{22}^d& 2{\sigma }_{12}^d{\sigma }_{22}^d\end{array}\\ \begin{array}{cc}& \\ & \end{array}& \begin{array}{cc}{\sigma }_{33}^d{\sigma }_{33}^d& 2{\sigma }_{23}^d{\sigma }_{23}^d\\ & 4{\sigma }_{23}^d{\sigma }_{23}^d\end{array}& \begin{array}{cc}2{\sigma }_{13}^d{\sigma }_{33}^d& 2{\sigma }_{12}^d{\sigma }_{33}^d\\ 4{\sigma }_{13}^d{\sigma }_{23}^d& 4{\sigma }_{12}^d{\sigma }_{23}^d\end{array}\\ \begin{array}{cc}& \\ & \end{array}symmetry& \begin{array}{cc}& \\ & \end{array}& \begin{array}{cc}4{\sigma }_{13}^d{\sigma }_{13}^d& 4{\sigma }_{12}^d{\sigma }_{13}^d\\ & 4{\sigma }_{12}^d{\sigma }_{12}^d\end{array}\end{array}\right]$$

(102)

Where, ${\sigma }_{ij}^d={\sigma }_{ij}-\frac{1}{3}{\sigma }_{kk}{\delta }_{ij}$ are the components of the deviator stress tensor and ${\tau }_0$ is the octahedral shear stress, $M_T^m=\frac{E^m{E}_T^m}{E^m-E_T^m}$ ,where $E^m$ and $E_T^m$ are the respective moduli of elasticity and the equivalent of the strain hardening for the matrix [28].

Hence, stresses and strains in the nonlinear regime can be ccalculated iteratively by [28]:

$${\sigma }_{ij}^{k+1}={\sigma }_{ij}^k+{d\sigma }_{ij}$$

(103)

$${\epsilon }_{ij}^{k+1}=S_{ij}^{k+1}{\sigma }_{ij}^{k+1}$$

(104)

With k = 0,1,2... as the counter that indicates the respective loading step for which the stress and strain fields are calculated iteratively.

A synthesis about the study [28] ends by standing out two topics of the modeling proposed by the researchers. The first of them is about the sharing of the strain rate in matrix and fibers for a unidirectional composite, which is mathematically modeled based on the rule of mixture as follows:

$${\epsilon }_i^{’f}=S_{ij}^f{\left(V_f{S}_{ij}^f+V_m{S}_{ij}^f{A}_{ij}\right)}^{-1}{\epsilon }_j^{’}$$

(105)

$${\epsilon }_i^{’m}=S_{ij}^m{A}_{ij}{\left(V_f{S}_{ij}^f+V_m{S}_{ij}^m{A}_{ij}\right)}^{-1}{\epsilon }_j^{’}$$

(106)

In which ${\epsilon }_i^{’f}$ and ${\epsilon }_i^{’m}$ are the strain rates shared by the reinforcement and matrix phases, respectively, in relation to the total strain rate (${\epsilon }_j^{’})$ for the unidirectional composite [28].

The second topic would be the descriptions of dynamic progressive failure criteria in laminated composites, which can be expressed as follows [28]:

$${\left(\frac{{\sigma }_x}{X_t\left({\epsilon }^{’}\right)}\right)}^2+{\left(\frac{{\sigma }_{xy}}{S_{xy}\left({\epsilon }^{’}\right)}\right)}^2=\left\{\begin{array}{c}\ge 1,failure\\ <1,safety\end{array}\right.$$

(107)

The equation (107) is the failure criterion due to dynamic tensile loading on the fibers, where ${\sigma }_x,X_t\left({\epsilon }^{’}\right),{\sigma }_{xy},S_{xy}\left({\epsilon }^{’}\right)$ are: the normal tensile stress, the tensile strength, the in-plane shear stress and the in-plane shear strength, respectively [28].

For dynamic compression loading, can be written [28]:

$${\left(\frac{{\sigma }_x}{X_c\left({\epsilon }^{’}\right)}\right)}^2=\left\{\begin{array}{c}\ge 1,failure\\ <1,safety\end{array}\right.({\sigma }_x<0)$$

(108)

With $X_c\left({\epsilon }^{’}\right)$ being the compressive strength [28].

Similarly, for the matrix, it is used [28]:

$${\left(\frac{{\sigma }_y}{Y_t\left({\epsilon }^{’}\right)}\right)}^2+{\left(\frac{{\sigma }_{xy}}{S_{xy}\left({\epsilon }^{’}\right)}\right)}^2=\left\{\begin{array}{c}\ge 1,failure\\ <1,safety\end{array}\right.$$

(109)

$${\left(\frac{{\sigma }_y}{Y_c\left({\epsilon }^{’}\right)}\right)}^2+{\left(\frac{{\sigma }_{xy}}{S_{xy}\left({\epsilon }^{’}\right)}\right)}^2=\left\{\begin{array}{c}\ge 1,failure\\ <1,safety\end{array}\right.$$

(110)

In equation (109), ${\sigma }_y$ and $Y_t\left({\epsilon }^{’}\right)$ represent the tensile stress transverse to the fibers and their respective strength and in equation (110) $Y_c\left({\epsilon }^{’}\right)$ is the compressive strength limit. It is worth mentioning that tensile and compression failure modes in the fibers are catastrophic, in which the composite completely loses its ability to support loads, while the so-called non-catastrophic failures (material retains partially the capacity to support loads) caused by other modes are also discussed briefly in [28].

Jarali et al [29] developed micromechanical modeling aimed at three-phase composites with shape memory polymer matrices. In this research, the authors extend Eshelby’s method by considering two distinct inclusions immersed separately in a heterogeneous shape memory matrix. The effective average elastic modulus assuming dual inclusion modeling for the analyzed composite was given by:

$$\overline{E}=E^{SMP}\left[I-f^F{\left[{{(E}^{SMP}-E^F)}^{-1}{E}^{SMP}-S^F\right]}^{-1}-f^N{\left[{{(E}^{SMP}-E^N)}^{-1}{E}^{SMP}-S^N\right]}^{-1}\right]$$

(111)

Equation (111) is obtained by the authors in [29] as a new dilute distribution relationship based on Eshelby’s theory, which contemplates elastic properties of the composite when the representative volume element has two inclusions in a heterogeneous matrix. In that equation, $E^{SMP}$ is the elastic modulus of the shape memory polyme matrix, I is the identity tensor, $f^F,f^N$, $E^F,$ $E^N,S^F{\mathrm{a}\mathrm{n}\mathrm{d}S}^N$ are the volume fractions of the fiber and carbon nanotube inclusions, the respective elastic moduli and the fourth order Eshelby tensors also for fiber and nanotube. With the novel proposal summarized above, Jarali et al [29] obtain a satisfactory comparison in relation to fundamental models such as the rule of mixture, for example, showing the consistency of the suggested model.

A novel micromechanical modeling to estimate effective properties in unidirectional composites reinforced by circular cross-section fibers was developedby Vignoli et al [30]. Based on the rule of mixture, the following equations were obtained:

$$E_1=E_1^f{V}_f+E^m(1-V_f)$$

(112)

$${\nu }_{12}={\nu }_{12}^f{V}_f+{\nu }^m(1-V_f)$$

(113)

$$E_2=E^m\left[1+2{\left(\frac{V_f}{\pi }\right)}^{1/2}\left(\frac{\pi }{2a_{22}}-\frac{\ln \left(a_{22}+\sqrt{a_{22}^2-1}\right)}{a_{22}\sqrt{a_{22}^2-1}}-1\right)\right]$$

(114)

$$G_{12}=G^m\left[1+2{\left(\frac{V_f}{\pi }\right)}^{1/2}\left(\frac{\pi }{2a_{12}}-\frac{\ln \left(a_{12}+\sqrt{a_{12}^2-1}\right)}{a_{12}\sqrt{a_{12}^2-1}}-1\right)\right]$$

(115)

$$G_{23}=G^m\left[1+2{\left(\frac{V_f}{\pi }\right)}^{1/2}\left(\frac{\pi }{2a_{23}}-\frac{\ln \left(a_{23}+\sqrt{a_{23}^2-1}\right)}{a_{23}\sqrt{a_{23}^2-1}}-1\right)\right]$$

(116)

In which the parameters, $a_{12}=2\sqrt{\frac{V_f}{\pi }}\left[\left(\frac{G^m}{G_{12}^f}\right)-1\right],a_{23}=2\sqrt{\frac{V_f}{\pi }}\left[\left(\frac{G^m}{G_{23}^f}\right)-1\right],a_{22}=2\sqrt{\frac{V_f}{\pi }}\left[\left(\frac{E^m}{E_2^f}\right)-1\right]$, with $E_1^f,E_2^f{V}_f,E^m,{\nu }_{12}^f,{\nu }^m,G_{12},G_{23},G^m$ being the longitudinal and transverse modulus of elasticity of the fibers, the volumetric fraction of fibers, the modulus of elasticity of the matrix (considered homogeneous and isotropic), in-plane Poisson’s ratio of the fibers, Poisson’s ratio of the matrix, the in-plane and out-plane shear moduli of lamina and the matrix shear modulus, respectively [30].

The researchers in [30], through a comparative analysis with other fundamental models, were able to conclude that the model called VSPKc showed the best results based on experimental data from literature and finite element analysis.

The technique known as asymptotic homogenization and the unit cell concept for multiphase materials are used in [31] in order to develop a general micro and nanomechanical model for an anisotropic smart nanocomposite with quantum dots embedded and piezoelectrically active constituents. In summary, the concept of unit cell is used for the microstructure of the material and a power series expansion of the variables u_i(x,y) (displacement), σ_i(x,y) (stress), D_i(x,y ) (electrical displacement), φ_i(x,y).

The effective properties of the multiphase material are determined from the set of equations (117)-(122) and, according to the authors of the study [31], the analytical results present good agreement with those obtained numerically with the aid of the finite element software ABAQUS.

$${\overline{C}}_{ijkl}\frac{\partial ²u_k^0}{\partial x_j\partial x_l}+{\overline{e}}_{kij}\frac{\partial ²{\phi }^0}{\partial x_j\partial x_k}=f_i\left(x\right)$$

(117)

$${\overline{e}}_{ikl}\frac{\partial ²u_k^0}{\partial x_i\partial x_l}+{\overline{k}}_{ik}\frac{\partial ²{\phi }^0}{\partial x_i\partial x_k}=0$$

(118)

On what:

$${\overline{C}}_{ijkl}=\frac{1}{\left|Y\right|}{\int }_Y^{}\left(C_{ijkl}\left(y\right)+{\tau }_{ij}^{kl}\left(y\right)\right)dv$$

(119)

$${\overline{e}}_{kij}=\frac{1}{\left|Y\right|}{\int }_Y^{}\left(e_{kij}\left(y\right)+{\lambda }_{ij}^k\left(y\right)\right)dv$$

(120)

$${\overline{e}}_{ikl}=\frac{1}{\left|Y\right|}{\int }_Y^{}\left(e_{ikl}\left(y\right)+d_i^{kl}\left(y\right)\right)dv$$

(121)

$${\overline{k}}_{ik}=\frac{1}{\left|Y\right|}{\int }_Y^{}\left({-k}_{ik}\left(y\right)+g_{ik}^{}\left(y\right)\right)dv$$

(122)

In which, ${\overline{C}}_{ijkl},{\overline{e}}_{kij},{\overline{k}}_{ik}$ are the tensors of elastic, piezoelectric and dielectric permittivity coefficients, in addition, y is the variable that accounts for the periodicity of the composite’s microstructure and ${\tau }_{ij}^{kl},{\lambda }_{ij}^k,d_i^{kl},g_{ik}^{}$ are terms associated with the rate of variation of the aforementioned tensors, as can be seen in [31] .

Finally, this literature review must stands out an interesting probabilistic micromechanical modeling carried out by Naskar et al [32], where a stochastic concept about representative volume element (RVE) is proposed. The authors of that research use a structural domain divided into N micro subdomains in which the material properties are presented randomly in space, in other words, the fundamental volumetric unit is considered within a stochastic context instead of the homogenization proposed by a conventional representative volume element. Hence, effects of the random spatial distribution of mechanical properties for different regions of a structure can be accounted in the modeling [32].

3. Final Remarks

Notably, micromechanical analyses have had the opportunity to deal with increasingly sophisticated modeling from the point of view of mathematics and appropriated applied sciences. However, traditional or reformulated elementary methods and approaches are still widely used in studies on the mechanics of composites involving their constituent phases, in other words, these fundamental models have also served as the basis for advances made in the area and have not fallen into disuse with the emergence of more elaborate proposals.

In fact, for all the content analyzed in this systematic review of literature, the micromechanics of polymer laminated composites constitutes a branch of solid mechanics that is completely open and is undergoing frank and broad expansion with regard to the development and consolidation of new knowledges within engineering applied sciences, thus, offering the opportunity for several physical and mathematical modelings to be proposed, duly evaluated and improved with the aim of predicting effective properties to describing progressive damage mechanisms in these multiphase materials with greater accuracy.

Exposed all this context, a trend analysis is constructed for future work that may come in the line of increasingly developing numerical models in order to optimize computational costs and with approaches that are able of considering strain fields influenced by the neighborhood, in other words, those are based on non-local theories and methods, for example, to multiphase laminates that also feature nanoparticle reinforcement.

In addition to the perspectives described previously, novel stochastic micromechanical approaches are expected, i. e., modeling that assumes that random character of the composite’s properties in more realistic representative volumes and, finally, studies considering characteristics of fibers, matrices and interfaces that can model the geometric path described by the tip a delamination crack through the existing duality with the stress wave fronts that advance through the material after the dynamic loading, for example, it also is a part of projections for future analyses. In case of a satisfactory approximation in relation to experimental observations, this last type of modeling mentioned could contribute significantly to the design of polymer matrix laminates with greater interlaminar strengths.