Finite strain formulation of viscoelastic damage model for simulation of fabric reinforced polymers under dynamic loading

The use of fabric reinforced polymers in the automotive industry is growing significantly. The high specific stiffness and strength, the ease of shaping as well as the great impact performance of these materials widely encourage their diffusion. The present model increases the predictability of explicit finite element analysis and push the boundaries of the ongoing phenomenological model. Carbon fibre composites made up various preforms were tested by applying different mechanical load up to dynamic loading. This experimental campaign highlighted the physical mechanisms affecting the initial mechanical properties, namely intraand interlaminar matrix damage, viscoelasticty and fibre failure. The intralaminar behaviour model is based on the explicit formulation of the matrix damage model developed by the ONERA as the given damage formulation correlates with the experimental observation. Coupling with a Maxwell-Wiechert model, the viscoelasticity is included without losing the direct explicit formulation. Additionally, the model is formulated under a total Lagrangian scheme in order to maintain consistency for finite strain. Thus, the material frame-indifference as well as anisotropy are ensured. This allows reorientation of fibres to be taken into account particularly for in-plane shear loading. Moreover, fall within the framework of the total Lagrangian scheme greatly makes the parameter identification easier, as based on the initial configuration. This intralaminar model thus relies upon a physical description of the behaviour of fabric composites and the numerical simulations show a good correlation with the experimental results.


Introduction
In these last decades, the use of carbon fabric reinforced polymers (CFRP) in the automotive industry increased very significantly.Their high specific stiffness and strength, their great energy absorption as well as the reduced manufacturing cost widely encourage their diffusion.
Previously limited to small runs (premium vehicles, racing), last advances in highly productive manufacturing process lead to the use of CFRP for high volume automotive production.Therefore, the behaviour understanding and modelling of these materials become essential for their implementation into the design loop, needed for the deployment on mass-produced vehicles.In order to ensure the protection of pedestrians and drivers/passengers in case of collision with a CFRP panel, a model dedicated to the finite element analysis (FEA) of impacts is needed.
The impact properties of layered composite materials are well reviewed by Richardson [1] and Cantwell [2].The nonlinear material behaviour which leads to differences in the impact response of composites is attributed to fibre failure, intra-and interlaminar matrix cracking, fibrematrix debonding and strain rate sensitivity of the matrix.Additionally, the textile composite materials are capable of large shearing prior to the ultimate failure due to the sliding a Corresponding author: sylvain.treutenaere@univvalenciennes.fr and reorientation of yarns.The modelling of all these phenomena is essential to describe the impact behaviour of layered fabric composites.
To develop a material model at the ply scale achieves a balance between CPU-time and physical meaning needed for structural analysis.Indeed, the properties of the above mentioned phenomena (except for the interlaminar matrix cracking) are common to a ply but can vary from one ply to another.Among these, this study is focused on the finite strain modelling of the intralaminar matrix damage and the strain-rate dependency of fabric plies.
The use of the Continuum Damage Mechanic (CDM) for anisotropic and composite materials was introduced by Chaboche [3] and Ladevèze [4].Further works of Chaboche [5,6] introduced the damage deactivation after the closing of cracks.Marcin [7] adapted and extended the micromechanics based CDM model proposed by Chaboche [8] to the fabric reinforced materials for implicit simulations.
In addition, Marcin extends the matrix damage model by taking into account the resin-induced viscoelasticity.To do this, a bi-spectral viscoelastic model is added.However, this formulation is not adapted to an explicit analysis because of the Newton-Raphson iterative loop to achieve the equilibrium.
Accordingly, the spectral viscoelastic model is replaced by a generalized Maxwell model.It has been successfully extended and used in finite strain for anisotropic [9,10] materials.

EPJ Web of Conferences
In this paper, the matrix damage model proposed by Marcin is coupled to a generalised Maxwell model.This present damageable viscoelastic model is then expressed in the finite strain framework and implemented in the commercial finite element software LS-DYNA.A comparison between experimental and numerical results is done in the final part.

Continuum matrix damage model
Experimentally the matrix cracks are observed to be oriented in the directions of reinforcement or transverse to them.As the Onera Damage Model MicroStructure (ODM MS) [7] is based on this assumption, the present model relies on a close version of the ODM MS.

Constitutive relation
The model is formulated in strain-space to ensure a good efficiency for FEA.The constitutive relation, by using the Voigt notation, is defined as follows where σ and ε are respectively the infinitesimal stress and strain tensor.C 0 is the initial stiffness tensor and Cm is the effective stiffness tensor.The stored strains ε s are needed in order to avoid the discontinuity of the (σ, ε) response for bi-axial loading, and to ensure the recovery of the initial stiffness after the damage deactivation.
The evolution of the residual strains ε r is linearly dependent on the damage evolution.

Effective stiffness tensor
The damage is introduced by adding additional compliance to the initial compliance tensor S 0 .The effective stiffness tensor is thus defined with ( Cm ) −1 = S 0 + S m .The additional compliance tensor due to the matrix damage is given by where H m i is the compliance tensor associated with the damage variable d m i .η m i represents the crack closure index which varies from 0 (closed crack) to 1 (opened crack).

Damage evolution
In order to describe separately the different fracture modes, two driving forces affect the evolution of each damage variable: the normal y nm i and the tangential y tm i to the damage directions.The thermodynamic forces so defined are dependent on the effective strain tensor, which corresponds to the positive part of the spectral decomposition of the strain tensor.
Therefore, the damage criterion is defined as follows with f the cumulative distribution function of Weibull.This formulation -compared to a standard frame where the thermodynamic forces are based on the derivative of the Helmholtz free energy -makes easier the parameter identification, retains the explicitness of the formulation and ensures the thermodynamic coherence.

Viscoelastic damageable continuum model
The generalised Maxwell model is well adapted to the explicit finite element method because of the direct dependence of the total stress with the strain.It is symbolically represented by a spring of stiffness C ∞ paralleled with N Maxwell elements (Fig. 1).Hence, the total stress, provided by the coupling between the above mentioned matrix damage and the generalised Maxwell viscoelastic models, is given by where σ ∞ is the long time stress, defined by means of the ODM MS with the equation (1), and h j is the effective viscoelastic stress of the j-th branch.This tensor is determined thanks to the Boltzman superposition principle and is expressed in an iterative form [9,10] by and Cm j the effective viscoelastic stiffness tensor.This tensor is defined in a same manner as the effective stiffness tensor presented in the subsection 2.2.Consequently, where d m i is the same damage variable as that used during the computation of the long time stress in the Eq. ( 1), S 0 j is the initial viscoelastic compliance tensor and H ve,m i, j is the viscoelastic compliance tensor associated with the damage variable d m i .To use the same damage variables is not unreasonable as the matrix damage and the viscoelasticity are both inherent in the resin.

Material model formulation for finite strain
Given the shear behaviour of the fabric composite ply with the possibility of large rotations of the yarns, the base (or undeformed) and the current (deformed) configurations cannot be assumed identical.It is therefore essential to extend the model in finite strain.
With few exceptions, Lagrangian (or Material) coordinates are used to describe the deformations of a solid.This approach facilitates the formulation of material constitutive models since the position and the physical properties of solid particles are described according to a reference position of these material particles and time.

Expression of tensors
As the objective is the formulation of a constitutive model in finite strain for structural analysis, it is appropriate to recall the definition of strain and stress tensors.
Let d X be a position vector which describes material points in the undeformed configuration.The material points in the deformed configurations are now described by dx (Fig. 2).The change of the material points is defined by the function x i = m i (X I , t).Hence, the differential is which leads to the definition of the deformation gradient tensor F by The transformation of a volume element dv 0 to a volume element dv is given by the relation The elongation of a position vector in the N direction is defined as Lastly, the shear angle is determined by From now, a distinction should be made between the tensors expressed according to the base configuration or the current configuration.First of all, suppose that the reference is the base configuration.The scalar product of the material vectors in the current configuration, evaluated according to these material vectors in the base configuration, is given by where C = F T F is the right Cauchy-Green deformation tensor.Similarly, the variation of the scalar product according to the material vectors in the base configuration is defined with with δ J K = 1 when J = K , otherwise zero, and where E = 1 2 (C − 1) is called the Green-Lagrange strain tensor with 1 the identity tensor.
In a consistent manner, but by considering the current configuration as reference, the strain can be described by the Almansi tensor A = 1 2 (1 − B −1 ) where B = FF T is the left Cauchy-Green deformation tensor.Now that the different strain tensors are defined, the stresses have to be evaluated according to both configurations too.d f is a force which acts on the current body and is the only one measurable from the experiments.By expressing the force vector according to the base or the current configuration, this leads to with T the Cauchy stress tensor and the first Piola-Kirchoff stress tensor (PK1).Note that whereas T is symmetric, because of its mixed nature is not.
Let d f 0 be a virtual force, seen as the equivalent of d f which may act on the reference configuration.d f 0 has no physical existence and is the transposition of d f in the base configuration: Hence, a new stress tensor integrally based on the reference configuration can be defined through the relation where S is called the second Piola-Kirchoff stress tensor (PK2) and S is symmetric.As a result, the behaviour law in finite strain can be expressed either according to the current configuration by finding the relation T = t(A) or according to the reference configuration by finding the relation S = s(E).

Iterative formulations
By means of the Finite Element Method (FEM), the final solution is obtained from the base configuration 04011-p.3C 0 by stepping computations in an incremental solution process.At step n, a current configuration C n is computed based on the reference configuration C R .The choice of this reference configuration leads to different Lagrangian formulations.

Updated Lagrangian formulation
In the updated Lagrangian formulation, the reference is assumed to be the previous current configuration C n−1 and is usually updated after each incremental step.Therefore, the new current configuration C n is directly getting from the previous current configuration C n−1 (Fig. 3).
The reference configuration being at a deformed state, the behaviour law is formulated in terms of the conjugate stress-strain pair T : A.
But this scheme leads to additional computations to achieve the incremental objectivity.In addition, objective stress rate as the commonly employed Green-Naghdi or Jaumann, do not correctly follow the directions of anisotropy of the composite material [11,12].LS-DYNA, for example, uses the updated Lagrangian framework as stress-return algorithm.

Total Lagrangian formulation
Otherwise, in the total Lagrangian formulation, the reference is assumed to be the base configuration at each step.Therefore, the current configuration C n is always getting from the base configuration C 0 , regardless the current configuration computed at the previous step (Fig. 4).
The reference configuration being at the undeformed state, the behaviour law is formulated in terms of the conjugate stress-strain pair (S, E).
With this formulation the objectivity is ensured by the use of total Lagrangian tensors and the direction of fibres is well followed.This formulation is used to formulate the present model in finite strain.

Coupling
The formulation used by the finite element program differs from the chosen scheme for the material model.The links between the total and the updated Lagrangian frameworks need to be established.This is done by converting the second Piola-Kirchoff stress tensor to the Cauchy stress tensor through the pushforward operation by using the formula

Constitutive relation
The matrix damage and the viscoelastic models share the same damage variables.Therefore, it is essential that both models are expressed by using the same formulation.Due to the anisotropy and as explained in the Sect.4.2.2, the model is formulated in the total Lagrangian framework.
The second Piola Kirchoff stress tensor S and the Green-Lagrange strain tensor E are symmetric and can be represented through the Voigt notation, respectively s and ε G .By following the methods given by Kaliske [9,10] and according to the damageable vicoelastic model formulated in the Sect.3, the constitutive relation becomes with and The parameter identification have to be adapted to this finite strain framework.

Shear locking
In the particular case of textile fabric, which is when the yarns are interlaced, a locking phenomenon appears after the shear reaches a critical angle.Accordingly, a significant rise of the elastic shear stiffness is observed.To take into account this phenomenon the constitutive relation becomes for the in-plane shear stress component with x = x if x > 0 and 0 if x < 0 and where α sl and ε sl G are material parameters.

Identification
The large rotations of the yarns appear under inplane shearing and other loading do not lead to large deformations.Therefore, the identification procedure presented in this article only concerns the parameters determined from cyclic in-plane shear test which is derived from the standardized NF EN ISO 14129 monotonic inplane shear test.
All measurements are carried out in the global coordinate system and provide: ε X = ε( X ) the elongation in the centre of the coupon according to the X direction, -ε Y = ε( Y ) the elongation in the centre of the coupon according to the Y direction, -f = f • X the measured force applied to the coupon.
As a result, the deformation gradient tensor F and the first Piola-Kirchoff stress tensor are in the global frame After a change of basis to be placed in the material frame and by using the Eqs.( 9) and (11), the Lagrangian tensors are determined and given by Hence, it becomes possible to identify the parameters of the total Lagrangian behaviour law S = s(E).
The damage evolution law is determined through the evolution of the shear stiffness on quasi-static tests and without considering any viscoelastic effects.However in the case of textile preforms and when the angle of locking is reached, it becomes impossible to distinguish the effects of the damage and of the shear locking.So after this critical shear angle, the damage evolution law is approximated to best fit the evolution of the residual strains.Indeed, the model is formulated in such a way that the residual strains are directly dependent on the damage variables.The gap in the shear stiffness is eventually considered as be due to the shear locking.The initial viscoelastic parameters are identifies through Dynamic Mechanical Analysis (DMA) tests [13] and the viscoelastic compliance tensors associated to the damage variables are determined by an iterating process to best fit the shear behaviour.

Numerical examples and discussions
The formulation proposed is illustrated by the simulation of in-plane shear tests carried out on two different fabric preforms.
Both are made up carbon fibres but the process to get the fabric preforms differ.The first one is a Non-Crimp Fabric (NCF) where the carbon fibre are arranged in two uni-directional orthogonal plies and then stitched together by a glass fibre tow (Fig. 6a).The second one is a plain-weave preform made up interlaced yarn including 3000 fibres each (Fig. 6b).
The experimental tests are carried out by means of a high-speed jack facility on coupons defined by the standard NF EN ISO 14129.Three speed loadings were set: 1.7 mm/s, 41 mm/s and 1000 mm/s.The shear strain is measured by an optical extensometer.The reference is taken as the mean value of the experimental results of a given configuration.
The Fig. 7 shows the comparison between the experimental results and the finite element analysis results of in-plane shear tests on both fabric materials at various strain-rate.The values are not stated because of confidentiality reasons.The model correlates well the material behaviours, even for completely different fabric preforms.

Concluding remarks
In this paper, a finite strain damageable viscoelastic model for fabric reinforced polymers has been presented.
The model is a combination between the Onera Damage MicroStructure model and a generalised Maxwell viscoelastic model.Because of the large rotation of the yarns for shear loading, the model is extended in finite strain.The total Lagrangian formulation is used for both matrix damage and viscoelastic models in order to well track the fibre orientation and ensure the objectivity.Moreover, the damage variables can be shared between the both element of the finite strain model.The phenomenon of shear locking in the particular case of textile preforms is also taking into account by a adding a non-linearity to the elastic shear stiffness.From the numerical point of view, the model is implemented in the commercial finite element software LS-DYNA.It is validated through experimental tests carried out at various strain-rate and reaching the loading rate of 1 m/s.The simulation results show the good efficiency of the proposed model.

04011-p.5 EPJ Web of Conferences
However, in order to fully simulate the behaviour of layered fabric composites, additional physical phenomenon have to be taking into account.Such is the case, for instance, of the fibre failure or the intralaminar damage.

Figure 2 .
Figure 2. Deformation of a body with representation of undeformed and current configurations.

Figure 3 .
Figure 3.The major configurations used in Updated Lagrangian framework.

Figure 4 .
Figure 4.The major configurations used in Total Lagrangian framework.

Figure 6 .
Figure 6.View of the upper face of the preforms.

Figure 7 .
Figure 7. Resultant stress-strain curves of in-plane shear test simulation compared to experimental results at various strainrates.