A Computational Approach for Prediction of the Damage Evolution and Mechanical Characteristics of Random Fiber Composites

In this study we developed a computational approach by implementing the damage models proposed by authors for simulating the damage evolution and mechanical properties of random short fiber composites according to the respective characteristics of the matrix, the reinforcement and the volume fraction. Material damage induced by fiber de-bonding is considered. A comparison between the different existing models of homogenization was performed to determine the model that best reflects the response of our study material. And furthermore a range of simulations was carried out to study the influence of various parameters of the composite for predicting the response of the material and damage evolution.


INTRODUCTION
Composite materials are increasingly used in several fields like automotive industry and aeronautic.
In the automotive field, the race for energy savings and limitations of pollutant emission is a growing interest on the part of manufacturers and suppliers for the use of materials with low densities.Meanwhile, concerns to protect passengers but also pedestrians in accidents are responsible for safety standards more stringent.Made of materials with complementary characteristics, composites appear as good candidates to meet all these requirements.It is in this context that the study of the mechanical behavior of composite materials under dynamic loading is necessary.The analysis of behavior with the damage of a composite subjected to such stresses requires theoretical tools and experimental approaches (Averous et al., 1997).A better understanding of the physical mechanisms of failure of the material would then optimize the modeling of the impact resistance of structures.Composite materials are generally characterized by greater specific strength than metallic materials and have particularly high energy dissipation (Lemaitre and Dufailly, 1987).
Among the most commonly used composite materials, the Sheet Molding Compound (SMC) is widely used in mass production because it meets the desired requirements (Main Lee et al., 1999).
This study is centered on predicting the properties and behavior by a multi-scale approach of behavior laws in dynamic for the study of composites with randomly oriented discontinuous reinforcements.
The laws of non-linear behavior were made through a homogenization type Mori-Tanaka based on the theory of Eshelby (1959) equivalent inclusion.The damage is then introduced at the local level through local criteria reflecting the physical phenomena of degradation.

MATERIALS AND METHODS
Raw material: SMC is a polyester matrix composite loaded with chalk and glass fiber reinforcement in the form of randomly oriented strand of 200 fibers.The fibers have a diameter of 15 µm and a length of 25 mm.The fibers represent a volume fraction of 18%.

Homogenization process:
The implementation of a model from a technical homogenization is carried out in three steps.
Representation step, where the constitution of the VER, shown in Fig. 1, is mechanically defined by the laws of behavior of components and geometrically by their shape or distributions.
Location step, which allows formalizing the relationship between the mechanical response at the microscopic level and the macroscopic scale.Is thus established, following a model of the behavior of each phase, the laws of location or concentration where the tensor strain localization, denoted A is introduced and the stress concentration tensor, denoted B. By definition, the laws of location A i and concentration B i for the "i" phase are given by (Okoli, 2001;Allix and Corigliano, 2000): (2) Homogenization Step, where representation and localization are used to construct the "macro mechanics" behavior law of ERV (Okoli, 2001): where, ‫ܥ‬ = Effective stiffness matrix ܵ = Effective compliance matrix

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 (1990) for a load imposed displacement and Reuss (1929) for traction imposed.These models, as will be shown later 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 model of Tanaka andMori (1970, 1972) and the selfconsistent 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.
To choose the model best suited to our material we made a comparison between Young modulus obtained by different methods with experimental results.
We note in Fig. 2 that's the model of Mori-Tanaka gives the closest to experimental results.
Mori-Tanaka's model: This model is the most suitable for our material because it has several advantages.This model can take into account in the actual behavior of the homogenized material, the presence of a large number of heterogeneities and the interactions between the local phases (Benveniste, 1987) where in 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.
Damage and failure criteria: Generally, damage is defined as a set of micro-structural changes in the material which causes more or less irreversible deterioration characteristics.In our material the damage observed phenomena are the de-bonding fiber/matrix and matrix cracking (Voigt, 1990), 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: With ‫ܣ‬ ாௌு : localization effectif tensor of Eshelby for phase i, given by: ‫ܣ‬ ாௌு : Effective localization tensor of Eshelby given by: where, ‫ܧ‬ ாௌு is Eshelby Shape Tensor The equivalent stiffness matrix of homogenized material by the method of Mori-Tanaka is expressed as: 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 de-bonding:
The behavior of composite materials is generally a damageable elastic or elastoplastic behavior.For our study material, the main source of nonlinearity is the damage of the fiber-matrix interface.To integrate this 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 us to calculate 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. 3.The vector ݊ ሬԦ is the outward normal vector of the fiber at the point A and T is the stress vector.
For a macroscopic stress Σ, the average stress tensor inside a fiber oriented in the direction θ i is given by the following equation (Voigt, 1990): where, ܳ is the localization Tensor for each family of reinforcements indexed "i".
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 a point of interface with the normal vector ݊ ሬԦ are obtained by simple projection of the stress tensor in the fiber: The vector ݊ ሬԦ in the axis system (1, 2, 3) is: Probability of failure in fiber/matrix interface: The process of breaking the fiber/matrix interface results of de-bonding on one or more sites (points on the equator of the fiber) emanating from a coupling between normal and tangential stress which produces different macroscopic damage levels for the different states of strain.
Several forms of interfacial failure criteria can be adopted.They consist of various combinations of normal and tangential stresses.There are two forms of failure criteria, linear and quadratic.For example the linear Coulomb criterion is given by Main Lee et al. (1999): where, ߚ and ܴ ௗ are parameters to identify.This criterion is attained when the value of R ୢ reached R ୢ .Fitoussi and Jendli et al. (2004) proposed a quadratic criterion with an elliptic form like: where, local constraints, normal and shear (ߪ , τ), 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.And two intrinsic parameters of the interface (ߪ , ߬ ) will be identified depending on the material.
In short fiber composites, several studies have shown that there are two general types of dispersions of the microstructure (Bay, 1991;Yen, 2002;Xiao and Mai, 1999).The first type is due to the difference in fiber concentration, this therefore uses the concept of distribution of the local fiber volume fraction.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 interface of fibers.For this, it is necessary to introduce a probabilistic aspect to the failure criterion.
The coupling between the probabilistic aspects and the Mori-Tanaka model is done through the introduction of statistical functions that can translate the kinetics of damage type of de-bonding.The kinetics of damage mechanisms involved is governed by a Weibull probability: The function ݂ሺߪ , ߬, ߪ , ߬ ሻ may take the form of one of the criteria of interfacial failure (linear or quadratic), 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 multi-scale model by using a local statistical criterion.We choose the quadratic form of the failure criterion for fiber-matrix interface which is written as follows: ሺσ , τ , 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.First, to evaluate the criterion, at each point of the interface to determine a percentage of failure of the interface.The second method is to evaluate the criterion to find the point where the criterion is maximum (Main Lee et al., 1999).
In this study we use the second method because we are interested in evaluating the volume fraction of the de-bonded fibers.

Damage evolution by interface de-bonding:
The debonding occurs as the criterion is reached on one of the points of the interface (Lataillade et al., 1996).This debonding decreases the participation of the fiber to the reinforcement of the composite.This effect is translated by the replacement of the fiber de-bonded by an "equivalent reinforcement".
One way of taking into account the fiber/matrix debonding is to consider de-bonded fiber as anellipsoidal cavity of stiffness zero (Zairi et al., 2008).The debonded fibers, shown schematically in Fig. 4, no longer contribute to the strengthening of the composite (Zouhaier et al., 2009).
After this, the compositeis consisted of three phases: matrix, reinforcing fibers and de-bonded fibers.The interface de-bonding occurs after the threshold for nonlinearity when the local de-bonding criterion for fiber-matrix interface is reached.It is therefore necessary to consider the probabilistic aspects to estimate the rate of damaged fibers (Zaïri, 2008;Lachaud and Michel, 1997).
For the i-th family of fibers with the orientation (θ, φ), the volume fraction of damaged fibers ܸ݂ ௗ is determined by the percentage of damage relative to the total volume fraction of fibers Vf i .For each family of reinforcement, the maximum value of the failure probability maxPr (α) is calculated on achieving interface with each increment charge.(Depending on the angle α around the fiber): Later it is this proportion of damaged fibers that will be replaced with ellipsoidal cavities like Fig. 5 shows.
The Fig. 6 shows the homogenization of the composite material which is then carried out in two successive stages.This schema takes into account the redistribution of stresses related damage.The model takes into account 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.

Damage model algorithm for fiber/matrix interface:
All the steps of the algorithm are summarized in Fig. 7.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 de-bonding of fiber/matrix, then is made the two-step homogenizing procedure.
The composite has three phases: matrix (V m ), reinforcing fibers (Vf r ) and damaged fiber (Vf d ).For the increment n we have: For n = 0: ܸ݂ ௗ = 0, ܸ݂ = ܸ݂, initial fiber fraction of the formulation.At each load increment n+1 we have: 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 drop in Young's modulus as a function of loading.

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 in Fig. 8 shows that from a certain value of Nf = 3, the calculated modules tend to a limit.So we can say that the fibers divide into 3 families is sufficient to represent the material.
The equivalent stiffness tensors obtained by combining 3 and 10 families of inclusions are:  These two tensors verify the symmetry properties of the tensor rigidities.The first tensor is monoclinic (13 independent coefficients).In the second case, the number of families has increased; the equivalent tensor is transverse isotropic (5 independent coefficients).
If we increase the number of families of reinforcements and we multiply their combinations (θ, φ), 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 begin the study with a wide range of form ratio (1-1000).The results are shown in Fig. 9.
Note that the calculated modules tend to a limit value from a value of 120.The diameter of the fibers being constant, it is in fact the length of the fibers increases and the values 120 of shear and tensile module tend to those of a long fiber composite.This implies that R must be in the range (1-120), we then made a more refined study of changes in elastic properties in this interval.The results are given in Fig. 10.More generally, Fig. 10 shows that the tensile modulus is more sensitive to the value of the form ratio than the shear modulus for values of R<40.

Influence of the volume fraction of reinforcement:
we study the influence of the volume fraction of reinforcement, taking a ratio of R = 25000/15 (fiber length L f = 25 mm, fiber diameter d f = 15 mm).
We see from Fig. 11 that the extreme values of E 1 are logically those properties of the matrix or fibers, respectively, for V f = 0 we have E m = 3000 MPa and for V f = 1 we obtain E f = 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 E 1 is equals to 14,640 MPa.
In addition, this choice also involves restricting the range of variations in study modules depending on the volume fraction V f .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 for the fiber-matrix interface; it's performed by inverse method.Indeed, for  a set of parameters (σ 0 , τ 0 , m) given, the model predicts the evolution of densitie ݀ of created cracks in the fiber-matrix interface in accordance with the imposed deformation.
From Fig. 12, we obtain the values shown in Table 1.
Once matched the criteria according to the strain rate, the model predicts the stress-strain regardless of the strain rate curves.
Predicting the behavior of SMC at different strain rates: Figure 13 shows a comparison experimentsimulation for the tensile response to different rates of deformation.We can see a good concordance regardless of the strain rate.
From Fig. 13 it is easy to see that the phase of the elastic behavior (ε≤0.45%)remains insensitive to the effect of strain rate.The slope of this phase remains relatively unchanged.On a macroscopic scale, increasing the strain rate creates a delay in the onset of non-linearity behavior.This then leads to an increased level of non-linearity as well as the stress and strain at rupture.
From Table 2 we see that the maximum value of the error is 11.15% for the yield stress and for the young modulus and the ߝ ோ the maximum value is 3.59 which is very acceptable.
Predicting the macroscopic damage: At the macroscopic scale, the loss modulus is translated through a damage variable D macro : We see from the Fig. 14 that for a strain rate of 20 s 1 , the degradation occurs at a macroscopic level Fig. 14: Evolution of the macroscopic damage in the SMC-R depending on the strain for different strain rates 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 small when the strain rate increases.To explain this phenomenon, we can say that the local growth of damage 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.In addition, the evolution of slopes also changed considerably.

Predicting the evolution of microscopic damage:
Microscopically, damage resulted through the variable d micro representing the amount of micro-cracks d (θ), introduced at each increment of stress, is directly functions of interface failure probabilities Pr (θ), calculated for each family orientation at each increment of macroscopic stress (∆∑ ோா ) is given by: We have defined 3 families of orientation: θ 1 (0°-30°), θ 2 (31°-60°), θ 3 (61°-90°) and in Fig. 15 we see the evolution of microscopic damage depending on the strains for the strain rate of 20 and 150 s -1 .
From the results shown in Fig. 16 we can say that for a tensile test, the fiber-matrix interfaces of the θ 3 family are mainly subjected to normal stress while interfaces of θ 1 and θ 2 families are subject to a stress of tensile shear coupling.It was noted that the level of deformation increases when the orientation of the fiber decreases.The fiber-matrix interface failure is developed for θ 3 family than other families of orientation.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, composites SMC-R 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%.Figure 17 illustrates the effect of the reinforcement's geometry on the equivalent behavior.It is found that the elastic-plastic behavior depends on the reinforcement's elongation.More reinforcement is long more modules and the yield increase.

CONCLUSION
The main objective of this study is to establish the behavior of composite materials with short fiber SMC-R in order to predict the elastic properties and damage.
We have developed a general application that provides the mechanical properties of a heterogeneous material by introducing the respective characteristics of the matrix, the reinforcement and the volume fraction.And also to predict the damage of this material as a function of strain rate, orientation angles, the number of families and the volume fraction.
The laws of behavior used in this study were applied to the finite element analysis of the behavior under impact.The objective is to enable the prediction of damage of composite structures with short fibers.

Fig. 4 :
Fig. 4: Schematic of the composite system with de-bonding interface

Fig. 9 :
Fig. 9: Evolution of tensile and shear modules depending on the form ratio

Fig. 10 :
Fig. 10: Tensile and shear modulus depending on the form ratio R in the range (1, 120)

Fig. 15 :Fig. 16 :Fig. 17 :
Fig. 15: Evolution of cracks density as a function of fiber orientation for a strain rate of 20 s -1 Allix, O. and A. Corigliano, 2000.Some aspects of interlaminar degradation in composites.Comput.Method Appl.M., 185: 203 Calculation of the local stresses in the matrix around the fiber After numerical applications carried out for the studied material, the stiffness tensors of the matrix and the reinforcement are, respectively: