Micromechanical Modeling of Random Short Fiber Reinforced Polymer Composites with Progressive Debonding Damage

Fiber-reinforced polymeric composites are widely used as structural materials in industry so there is an evident need for the prediction of their mechanical properties. This kind of materials shows various damage phenomena such as matrix cracking, interfacial debonding, fiber pullout and fracture. The predominant damage mechanism may vary according to fiber volume fraction, aspect ratio, strengths of the interface, fiber orientation and distribution. Progressive damage accumulation in the composite is known to affect the global mechanical properties [1]. In order to estimate the material overall response, the accumulated damage must then be included in the constitutive relations. Furthermore, in order to achieve a rigorous description, the constitutive equations must be derived from micromechanical considerations. Among the most commonly used composite materials, the Sheet Molding Compound (SMC) is widely used in mass production [2]. This study is centered on predicting the properties and mechanical behavior by a multi-scale approach of behavior laws in dynamic, for the study of composites with randomly oriented discontinuous reinforcements. The aim of this study is a contribution to mechanical modeling of damage in random fiber composites by combining to different models and adding a statistical parameter. This category of composites presents 3 different types of damage mechanism, matrix cracking, matrix/fiber debonding and fiber failure. The first mechanism in not took in consideration in most models for fiber reinforced composites damage prediction, the damage is not considered until fiber/matrix debonding. The first novelty of this work is in the combination of two models to estimate the global damage with more accuracy. The first one is the DSGZ model, named after its authors (Duan, Saigal, Greif, Zimmerman) [3] which is a phenomenological model for the prediction of polymeric matrix behavior. The second novelty is the use of a statistical parameter in the homogenization process, because the randomness in the fibers orientation causes a disparity in the results, and taking in consideration all the fiber orientation lead to heavy calculation. The laws of non-linear behavior were made through a Mori-Tanaka homogenization based on the theory of Eshelby equivalent inclusion. The damage is then introduced at the local level through local criteria reflecting the physical phenomena of degradation. The objective of the mechanics of heterogeneous materials is to estimate the macroscopic properties of an equivalent homogeneous material. Several authors include the well-known Voigt and Reuss bounds that take only the volume fraction of the composite materials into account, also Hashin and Shtrikmans introduced the notion of isotropic distribution of phases [1–2]. Other authors estimate the effective properties based up on the Mori Tanaka model or Self Consistent Scheme, consider only small elastic or elastic–plastic deformations on microscopic and macroscopic levels [4], [5]. Moreover, homogenization methods provide another way for the prediction of mechanical response of heterogeneous specimen replaced by a homogeneous equivalent material (HEM), which represents the material in an averaged sense [6–7]. In the present paper, the damage behavior of glass fiber-reinforced composite with polymeric matrix is investigated using a combined approach of micromechanical modeling and experimental characterization. The resulting micromechanical model is numerically implemented and used to simulate uniaxial loading for the sake of comparison with experimental results.


Introduction
Fiber-reinforced polymeric composites are widely used as structural materials in industry so there is an evident need for the prediction of their mechanical properties.This kind of materials shows various damage phenomena such as matrix cracking, interfacial debonding, fiber pullout and fracture.The predominant damage mechanism may vary according to fiber volume fraction, aspect ratio, strengths of the interface, fiber orientation and distribution.Progressive damage accumulation in the composite is known to affect the global mechanical properties [1].In order to estimate the material overall response, the accumulated damage must then be included in the constitutive relations.Furthermore, in order to achieve a rigorous description, the constitutive equations must be derived from micromechanical considerations.
Among the most commonly used composite materials, the Sheet Molding Compound (SMC) is widely used in mass production [2].This study is centered on predicting the properties and mechanical behavior by a multi-scale approach of behavior laws in dynamic, for the study of composites with randomly oriented discontinuous reinforcements.
The aim of this study is a contribution to mechanical modeling of damage in random fiber composites by combining to different models and adding a statistical parameter.
This category of composites presents 3 different types of damage mechanism, matrix cracking, matrix/fiber debonding and fiber failure.The first mechanism in not took in consideration in most models for fiber reinforced composites damage prediction, the damage is not considered until fiber/matrix debonding.The first novelty of this work is in the combination of two models to estimate the global damage with more accuracy.The first one is the DSGZ model, named after its authors (Duan, Saigal, Greif, Zimmerman) [3] which is a phenomenological model for the prediction of polymeric matrix behavior.The second novelty is the use of a statistical parameter in the homogenization process, because the randomness in the fibers orientation causes a disparity in the results, and taking in consideration all the fiber orientation lead to heavy calculation.
The laws of non-linear behavior were made through a Mori-Tanaka homogenization based on the theory of Eshelby equivalent inclusion.The damage is then introduced at the local level through local criteria reflecting the physical phenomena of degradation.
The objective of the mechanics of heterogeneous materials is to estimate the macroscopic properties of an equivalent homogeneous material.Several authors include the well-known Voigt and Reuss bounds that take only the volume fraction of the composite materials into account, also Hashin and Shtrikmans introduced the notion of isotropic distribution of phases [1][2].Other authors estimate the effective properties based up on the Mori Tanaka model or Self Consistent Scheme, consider only small elastic or elastic-plastic deformations on microscopic and macroscopic levels [4], [5].Moreover, homogenization methods provide another way for the prediction of mechanical response of heterogeneous specimen replaced by a homogeneous equivalent material (HEM), which represents the material in an averaged sense [6][7].
In the present paper, the damage behavior of glass fiber-reinforced composite with polymeric matrix is investigated using a combined approach of micromechanical modeling and experimental characterization.The resulting micromechanical model is numerically implemented and used to simulate uniaxial loading for the sake of comparison with experimental results.

Material of study
The composite sheet molding compound random (SMC-R) consists of an unsaturated polyester resin loaded with chalk (CaCO3) and glass fiber reinforcement randomly oriented with a weight content of 26%.These bundles of fibers are discontinuous and have a constant length (L = 25 mm).We define a family of reinforcement as a collection of fibers that exhibit the same volume fraction and orientation.The reinforcement orientation randomly distributed confers to the material a microscopic heterogeneous aspect and an overall transverse isotropic mechanical behavior.The SMC-R25 plates were prepared of thickness 2.7 mm and were cured at 140°C with an applied pressure averaging between 7 and 8 MPa for 2 min [8][9].Microscopic observations, using scanning electronic microscope, have been performed to investigate the material microstructure and to assess the random fibers orientation.Specimen cartography is then carried out using image analysis to quantify and characterize, at the microscopic scale [10][11].

Global modeling technique
The global modeling is carried out in three steps like shown in Fig. 1.The first step is the calculation of the strain tensor in the matrix and the damage if the elastic limit is reached.The second step is the calculation of the stress fiber/matrix interfaces stress using the stress tensor of the matrix.The last step is the homogenization process using the interfacial stress distribution and the strain tensor in the matrix.These steps will be detailed next.

. Homogenization technique
In general, engineering and natural materials are heterogeneous at a certain scale and this heterogeneous nature has a significant influence on the macroscopic behavior of multi-phase materials.To predict the macroscopic behavior of heterogeneous materials various homogenization techniques are used.
There are four basic steps to a homogenization process for a multi-phase composite [12]: the first stage is the definition of the RVE, the second is to analyze the mechanical behavior of each phase; then the third step is the description of the boundary conditions and interphase boundaries; and finally an homogenization strategy to predict the macroscopic behavior based on the mechanical response of the REV.
According to the Hill theory of average [13][14], the macroscopic strain and stress tensor components are defined as the average on the REV volume: 1. Representation step: in this step the laws of behavior of components, shape and distributions geometrically define the REV (Fig. 2).For our material, glass fibers are assimilated to ellipsoids.The laws of location Ai and concentration Bi for the "i" phase are given by [9]. . .

Homogenization
Step: where representation and localization are used to construct the "macro mechanics" behavior law of REV [15].The effective stiffness and compliance matrix, Ccomp and Scomp, are:   1 : .

Micromechanical approach
In general, the micromechanical modeling methods permeate research or actual overall behavior from the properties of the various components, interfaces and interactions.The first simple models have been proposed by Voigt (1889) for a load-imposed displacement, and Reuss (1929) for traction imposed stress.These models, as will be shown later by Hill (1952), respectively provide upper and lower bounds of effective elastic constants.Then came the study of Eshelby (1959) that served as the basis for the model of Tanaka andMori (1970, 1972) and the self-consistent scheme of Kroner (1958).The Voigt and Reuss models do not take into account the shape and orientation of reinforcements unlike the model of Eshelby and Mori-Tanaka [16].
The model of Mori-Tanaka is the most suitable for our material.It takes into account the actual behavior of the homogenized material, the presence of a large number of heterogeneities and the interactions between the local phases, where the matrix immersed in the heterogeneous medium is already disturbed by the presence of other heterogeneities.In addition, this model is an improved model of the equivalent inclusion of Eshelby.For the matrix behavior, we used the DSGZ model, [3] with the structure shown in Fig. 3.It is a phenomenolog-ical model for predicting the constitutive behavior and estimating elastic-plastic responses to static or dynamic solicitations of semi crystalline and amorphous polymers [12].
, ., Where:      Generally, damage is defined as a set of microstructural changes that causes more or less irreversible deterioration of the mechanical characteristics.In our material, the damage observed phenomena are the matrix cracking and fiber/matrix debonding; it is called an effective field in which we are reduced to solving the problem of heterogeneity with loading at infinity.Le localization tensor of phase i is given by: A EST is the Eshelby tensor given by: And Ai ESH is the localization effectif tensor of Eshelby for the i phasegiven by:   1 : .
And Ei ESH is the Eshelby Shape tensor.The equivalent stiffness tensor obtained by Mori-Tanaka's method is: After numerical application, we obtain the matrix and reinforcement stiffness tensors noted Cm and Cr.
The Eshelby shape tensor corresponding to ellipsoidal reinforcements is: 3077 Note that the stiffness tensor of the matrix and the reinforcement are isotropic, but the Eshelby tensor is not because of the ellipsoidal geometry of reinforcements.

Consideration of fiber/matrix interface debonding
The behavior of composite materials is in general a damageable elastoplastic behavior.For our study material, the main source of the damage nonlinearity is the fiber-matrix interface.To integrate the damage, we are now interested in the calculation of stress and strain fields in the fiber/matrix interface.The elastic model of Mori-Tanaka allows the calculation of the stress fields and deformation in each phase through the localization and concentration tensors.Fields in the fiber/matrix interface are determined by the constraints prevailing in the fiber.We consider a single fiber representing a family of fibers oriented (θ, φ) in the axis system represented in Fig. 2. The vector  ⃗ is the outward normal vector of the fiber at the point A and  ⃗ is the stress vector.For a macroscopic stress Σ, the average stress tensor in a fiber/matrix interface is given by the following equation (Voigt, 1889): where: Qi is the localization tensor for each family of reinforcements indexed "i" given by: The stress tensor in the interface can be calculated using the continuity condition of the normal stress at the interface.Normal and tangential stresses at an interfacial point with the normal vector  ⃗ are obtained by simple projection of the stress tensor in the fiber: .
  3.1.6.Failure probability in fiber/matrix interface The process of breaking the fiber/matrix interface (points on the equator of the fiber) emanates from a coupling between normal and tangential stress that produces different macroscopic damage levels for the different states of strain.Several forms of interfacial failure criteria can be adopted.There are two forms of failure criteria, linear and quadratic.In our study, we used a quadratic criterion with an elliptic form like: Where, local normal and shear constraints (σn, τ), at the same point of the fiber/matrix interface are function of the macroscopic load Σ, the fiber orientation (θ, φ), the point considered on the fiber/matrix and the mechanical properties of each phase.In addition, two intrinsic material parameters (σ0, τ0) of the interface will be identified.
In short fiber composites, there are two general types of microstructure dispersions [10], [15].The first type is the difference in fiber concentration.The second type is related to the presence of some relatively large pores, which can sometimes be very close to the fibers, thus affecting the state of stress at the fibers interface.For this, it is necessary to introduce a probabilistic aspect to the failure criterion.
The introduction of this aspects in the Mori-Tanaka model is done through the introduction of statistical functions f(σn,τ,σ0,τ0) that can translate the kinetics of damage type of de-bonding.The kinetics of damage mechanisms involved is governed by a Weibull probability [13]: The parameter m is a coefficient reflecting the sensitivity of the dispersion intensity of the damage.The initiation of the interfacial failure is introduced into the multiscale model by using a local statistical criterion.We choose the quadratic form of the failure criterion.σ0, τ0, m are parameters reflecting the strength of the fiber/matrix interface as well as the kinetics of damage to the statistical point of view.
There are two methods for assessing the failure criterion interface.The first method is to evaluate the criterion at each point of the interface to determine a percentage of interface failure.The second method is to evaluate the criterion to find the point where the criterion is maximum 2 .In this study we use the second method because we are interested in evaluating the volume fraction of thede-bonded fibers.

Damage evolution by interface de-bonding
The debonding occurs as the criterion is reached on one of the points of the interface [16].This debonding decreases the participation of the fiber to the reinforcement of the composite.This effect is translated by the replacement of the de-bonded fiber, that no longer contribute to the strengthening of the composite, by an "equivalent reinforcement" as an ellipsoidal cavity of stiffness zero [15], [17].After this, the composite is consisted of three phases: matrix, reinforcing fibers and debonded fibers.The interface debonding occurs after the threshold for nonlinearity when the local de-bonding criterion is reached.Later it is this proportion of damaged fibers that will be replaced with ellipsoidal cavities like shown in Fig. 4.
It is then necessary to consider the probabilistic aspects to estimate the rate of damaged fibers [6], [18].For the i-th family of fibers with the orientation (θ, φ), the volume fraction of damaged fibers Vfi is determined by the percentage of damage relative to the total volume fraction of fibers Vfi.For each family of reinforcement, the maximum value of the failure probability maxPr(α) is calculated on achieving interface with each increment.(Depending on the angle α around the fiber): The Fig. 5 shows the homogenization process which is carried out in two successive stages.

Fig. 6 Structure of the Mori-Tanaka damage model
All the steps of the algorithm are summarized in Fig. 6.This schema, takes into account the redistribution of stresses and the evolving anisotropy of the material related to the evolution of the damage.It also allows predicting the evolution of the loss of elastic modulus for all the coefficients of the stiffness tensor of the composite.
Initially, the mechanical properties of the matrix and the reinforcement are introduced.Then a first homogenization is done before load application.If the applied load is less than the damage limit, it is within the elastic field.If we go into the plastic range then we calculate interfacial constraints and the failure probability, after, if the damage criterion is satisfied, that is to say that there was fiber/matrix interfacial de-bonding, then is made the two-step homogenizing procedure.
The composite has three phases: matrix (Vm), reinforcing fibers (Vfr) and damaged fiber (Vfd).For the increment n we have: for n=0: 00 0, dr Vf Vf Vf  initial fiber fraction of the formulation. .
At each load increment n+1 we have: 1 Pr .
We define a new material: At each load increment, the introduction of damaged fibers of characteristic zero produces a gradual decrease in the stiffness of the composite.The damage is characterized by a loss in Young modulus.

Results and analysis
The model built on the basis of the homogenization method of Mori-Tanaka is best suited for estimating the elastic properties of the material studied.This homogenizing method allows to take into account the shape and size of the reinforcement and to establish a coherent behavior with a law of distribution of fiber orientation.

Influence of the number of reinforcement family
A family is a group of reinforcing fibers that share the same geometric properties.For our study, we varied the number of families from 1 to 20 to see the influence of this parameter on the material properties.The evolution of tensile and shear modulus shows that from a value of Nf = 3, the calculated modules tend to a limit.Therefore, we can say that the fibers divide into 3 families is sufficient to represent the material.
The equivalent stiffness tensor obtained by combining three families of inclusions is: This tensor is monoclinic and verifies the symmetry properties of the tensor rigidities.The first tensor (  family's number; the equivalent tensor is transverse isotropic (5 independent coefficients).We find the elastic behavior of the equivalent homogeneous material.Thus, the random orientation of non-spherical reinforcements led to a transverse isotropic composite.

Influence of the form ratio of reinforcement
The evolution of Young's and shear modulus based on the form ratio R (fiber length/fiber diameter) is directly related to the length of the fiber as the fiber diameter is constant.We note, in Fig. 7, that the calculated modules tend to a limit from a value of form ratio R=120.The diameter of the fibers being constant, it is in fact the length of the fibers.
More generally, we note that the tensile modulus is more sensitive to the value of the form ratio than the shear modulus.

Influence of the volume fraction
We study the influence of the volume fraction of reinforcement, taking a ratio of R = 25000/15 (fiber length Lf = 25 mm, fiber diameter df = 15 mm).We see from Fig. 8 that the extreme values of E1 are logically those properties of the matrix or fibers, respectively, for Vf = 0 we have Em = 3000 MPa and for Vf= 1 we obtain Ef = 73000 MPa.Note that the modules are highly sensitive to variations in fiber volume fraction.We find that for a volume fraction of 18%, corresponding to our material, tensile modulus E1=14,64 MPa.
In addition, this choice also involves restricting the range of variations in study modules depending on the volume fraction Vf.Indeed, the fiber volume fraction, which can go up to 60-70% in the long-fiber reinforced composites, cannot practically exceed 40% in the case of short fibers and polymer matrix.More, it is important to note that beyond a value greater than 40%, the method of Mori-Tanka no longer reflects correctly the behavior of the composite.The restriction comes from the method itself, it is not capable of taking into account the interactions "medium distance" between heterogeneities.

Parameter identification of interfacial de-bonding criterion
To take into account the de-bonding fiber/matrix in the model, we identify the parameters (σ0, τ0, m) which are the characteristics of the interface and the statistics of interfacial failure.The local probability of inter-facial failure is given by:     (  ,, ).
The first step is the identification of the local failure criterion parameters; it is performed by inverse method.Indeed, for a set of given parameters (σ0, τ0, m), the model predicts the evolution of density dmicro of created cracks in fiber-matrix interface in accordance with the imposed deformation.The results are shown in Table 1.From Table 2, we see that the maximum value of the error is 10.55% for the yield stress.For the young modulus and εr the maximum value is 3.29% which is very acceptable.

Predicting the behavior of SMC-R25 at different strain rates
Fig. 9 shows a comparison experiment-simulation for the tensile and compression response to different deformation rates.We can see a good concordance regardless of the strain rate.It is easy to see that the phase of the elastic behavior (ε≤0.50%)remains insensitive to the effect of strain rate.On a macroscopic scale, increasing the strain rate creates a rise of the elastic limit.This then leads to an increased level of non-linearity as well as the stress and strain at rupture.Once matched the criteria according to the strain rate, the model predicts the stress-strain regardless of the strain rate.

Predicting the macroscopic damage
At the macroscopic scale, the modulus loss is translated through a damage variable Dmacro: We see from the Fig. 10 that for a strain rate of 20s -1 , the degradation occurs at a macroscopic level deformation of 0.2%, while for a strain rate of 160 s -1 , the first stiffness reduction occurs at a deformation of 0.35%.
In addition, we can see that the damage is relatively minor when the strain rate increases.This phenomenon can be explained by the fact that the local damage growth is changed in terms of deformation and exhibits a reduced kinetic due to the effect of strain rate.
Both aspects are closer to the effect of viscosity produced by the delay reaching the dissipation interfacial areas.Therefore, we can notice a delay in the macroscopic damage thresholds.Microscopically, damage resulted through the variable dmicro representing the amount of micro-cracks d(θ), introduced at each increment of stress, and is directly functions of interface failure probabilities Pr(θ), calculated for each family orientation at each increment of macroscopic stress (ΔΣREV).
We have defined 3 families of orientation: θ1(0°-30°), θ2(31°-60°), θ3(61°-90°).From the results shown in Fig. 11 we can say that for a tensile test, the fiber-matrix interfaces of the θ3 family are mainly subjected to normal stress interfaces of θ1 and θ2 families are subject to a tensile shear coupling.It was noted that the level of deformation increases when the orientation of the fiber decreases.We can conclude that the damage coming from the fiber/matrix interface de-bonding are mainly caused by pure normal interfacial stress.This leads to specify that the fiber-matrix interface damage is anisotropic in nature.Therefore, compositesSMC-R25 are characterized by a damageable visco-elastic behavior.

Influence of form ratio
We conduct our numerical simulations with a volume fraction of reinforcement Vf=18%.Fig. 12 illustrates the effect of the reinforcement's geometry on the equivalent behavior.It is found that the elasto-plastic behavior depends on the reinforcement's elongation.More reinforcement is long more modules and the yield increase.Fig. 12 Effect of the form ratio on the SMC-R25 behavior

Conclusion
The main objective of this study is to establish the behavior of composite materials with short fiber in order to predict the mechanical properties and damage evolution.The modeling technique adopted in consists on the combination of a phenomenological model for the matrix behavior and homogenization technique with a statistical media.It was shown that this approach can be successfully applied to predict the macroscopic behavior and damage evolution of this range of materials with complex material behavior and microstructure.
The comparison of the simulation results and experimental response provided an extra assessment of the proposed modeling technique.The material parameters obtained for the matrix material have been used with satisfactory agreement.
The modeling of the macroscopic behavior using the viscoplastic model and homogenization method demonstrates an accurate description of the deformation behavior.Finally, additional developments to improve the modeling are still required.First, a cavity nucleation criterion is necessary to describe the damage initiation and the uses of a computational homogenization method with incremental scheme.
The proposed modeling strategy seem to be an efficient tool to establish the micro-macro relationship.


Fig.3Structure of the DSGZ model .4.Damage and failure criteria

Fig. 9
Fig. 9 Comparison experimental-simulation for different strain rates

Fig. 10
Fig. 10 Evolution of the macroscopic as function of the strain for different strain rates 5.2.Predicting the evolution of microscopic damage