Modeling and simulation of damage in elastomer structures at high strains

A new model is proposed in order to simulate correctly the behaviour of elastomers for very large deformations. The hyperelastic Hart–Smith s model is revisited in order to take into account the damage of elastomers and the opening and the closing of micro-defects. Applications concern structures with elastomer used in space industry. Shear loadings are typical solicitations for these structures. Whenever the material is able to undergo very large deformations (600– 700%), the simulation of this type of experiment up to so high strains is not easy. The failure of industrial finite element code leads us to develop a prototype software using a non-incremental method (large time increment method) and a material formulation with ‘‘rotated’’ quantities. Some applications are shown for elastic or hyperelastic materials, especially for nearly incompressible materials.


Introduction
The aim of this study is to simulate the behaviour of metal/elastomer structures used in space industry. The experiments performed on bulk elastomer specimens and on structures show the capacity of this material to undergo very large deformations. For instance, in a shear test case, rupture generally occurs at strains in range 600-700%. Obviously, simulation up to these high strains is not easy: first of all, because of numerical difficulties such as distortion of elements. An alternative developed by Chen et al. [4] is to work with meshless methods. The near incompressibility of elastomers is an additional important difficulty because of a very bad behaviour of classical elements. The complexity of the problem is confirmed by the number of attempts to develop elements with a good accuracy even for coarse meshes ( [2,7,9,19] for instance). Moreover, for the shear experiment case, local instabilities seem to occur (premature stop of calculations with the finite element code ABAQUS). Another difficulty is due to the modeling of material behaviour: only few hyperelastic laws are able to give good adequacy between calculations and experiments for very large strains. Moreover due to the high strains levels reached by this type of elastomer, damage such as chain damage or micro-void formation occur. The damage level can be evaluated through the degradation of mechanical material properties in tension. But in case of compression, no degradation of mechanical properties is observed: even if micro-defects exist, when micro-voids are closed, the initial properties are recovered. If some hyperelastic models taking into account the damage phenomena have been recently proposed [1,18], they do not include the unilateral aspect of damage which is nevertheless a worthy aspect. To be able to simulate the behaviour of the laminated structure, we 1 work on both constitutive relation and numerical simulation.
Concerning the damage modeling in elastomer, a new model using the concept of damage internal variables coupled to Hart-Smith energy is proposed. This classical hyperelastic model represents easily and correctly the hyperelastic behaviour of elastomers, especially in the case of very large strains since it enables the modeling of the increasing stiffness of elastomer. The unilateral damage is introduced through the split of this strain energy into two parts: a ''tension'' energy affected by damage and a ''compression'' energy not concerned with damage. Details are exposed in the second section. This model was identified through shear test results. The model has been implemented in the prototype code. A simple numerical example is presented.
Concerning the numerical simulation, as the industrial code ABAQUS has some difficulties to perform simulations in the shear test case, a prototype software was developed based on the LATIN method (LArge Time INcrement method, [12,15]). This non-incremental method is well adapted to strongly non-linear problems and has already demonstrated its capacity to crosscritical points without any particular technique (such as continuation techniques) in case of instability such as buckling [3]. This method is an iterative method which takes into account the whole loading process in a single time increment. It is recalled in the part Section 3.2 of this paper. It is combined with a formulation using varying ''rotated'' quantities and rated constitutive law. The treatment of incompressibility based on the B-Bar method is also presented. Some examples such as tension or shear experiment cases are shown in Section 4.

Experimental results
The behaviour of the elastomers is commonly represented by hyperelastic models, characterized by the strain energy function W, which is often considered to be dependent on the invariants I 1 ; I 2 ; I 3 of the Cauchy-Green deformation tensor C (or another strain tensor). Among the usual hyperelasticity laws (Mooney-Rivlin, Ogden, Varga,. . .), the best one for the studied material is the model of Hart-Smith, as shown in Fig. 1.
The agreement between experiments and model is quite good on a monotonic shear test, but some phenomena are not taken into account, essentially a decrease of the shear modulus during cycles and a decrease of the stiffening for high strains, as Figs. 2 and 3 show. The experimental results in a shear test case (realised with some loads and unloads) are shown in the axes G ¼ f ðcÞ and ln G ¼ f ðc 4 Þ (G is the shear modulus and c ¼ U d =h is the apparent shear strain).   Actually, the occurrence of micro-defects leads to a degradation of the mechanical material properties. Nevertheless it is worthy to note that for damaged material the type of load is crucial. In case of tension, the occurrence of micro-defects induces a decreasing in the material properties; on the contrary, in case of compression, when the micro-defects are closed, the initial properties are recovered. This non-symmetry in tension and compression appears clearly in the force-displacement response of the material plotted in Fig. 4, in the case of hydrostatic tension-compression test. The model exposed in the next part takes into account these experimental observations.

The standard Hart-Smith model
As the Hart-Smith law gives good results, this will be the basis of the model elaborated and which takes into account the material damage. In the case of incompressibility, the strain energy function is: where I 1 and I 2 are the invariants of tensor C (see Section 3.1.1). C 1 , C 2 and C 3 are three parameters, which must be identified. Thus, the behaviour relation is

Introduction of damage
If we examine the behaviour of the Hart-Smith law for small and for large strain, we find that the constants materials can be defined by Moreover, since the damage affects both initial shear modulus and stiffening at large strains, it is easy to take it into account by introducing three scalar damage variables, d 1 ; d 2 and d 3 , which affect respectively C 1 ; C 2 and C 3 . Thus, The damage is now introduced but we do not yet take into account its unilateral aspect. The solution is based on the strain energy function splitted into two parts: one for tensile for which damage variables act, the other one for compression for which initial characteristics are used, also: But this split is not trivial. The solution proposed in [11] is to replace classical invariants by new quantities obtained with the positive and negative parts of the Green-Lagrangian strain E. They are defined by writing the incompressibility constraint, that is In terms of the new quantities, I 0 1 and I 0 2 , this equation becomes where G is the Almansi-Lagrange tensor such as Taking into account that for small and mean strains, the terms det E and det G can be neglected, the quantities used to separate tension and compression (I þ 1 , I À 1 , I þ 2 and I À 2 ) are defined by the following equations: where hi À is the negative part (associated to negative eigenvalues) of and hi þ the positive part (associated to positive eigenvalues). And the new strain energy function for hyperelastic model with unilateral damage, can be written The state laws associated with this model are classically obtained: Y i are the thermodynamical forces associated with d i . By the addition of evolution law for the damage variables , the model would be complete.

Identification
Before using the proposed model, the three elastic parameters C 1 , C 2 and C 3 and the parameters describing the evolution laws for damage have to be identified. It was performed through shear test results (monotonic for elastic parameters and non-monotonic for damage parameters) specially by exploiting shear modulus properties (Eq. (4)). The different stages of the procedure are summarized below. Let us precise that as a first approximation, the term with I 2 has been neglected since the shear test case do not enable to split the contribution of the two first invariants and moreover simulation with Hart-Smith model with only the first invariant and experimental results are in good agreement, also here C 2 ¼ 0.
• As the gradient deformation in a shear test case can be written the new quantities defined in (8) and (9) verify Thus, stress tensor r and the thermodynamical forces Y 1 and Y 3 can be written in terms of the shear rate c through Eqs. (13), (14) and (16).
• Let s be the measure of the global force. The shear modulus can be expressed in terms of c, since: As in the case of the Hart-Smith model, the study of the asymptotic behaviour of this modulus enable us to estimate the level of damage for a given shear strain rate since: Also, d 1 and d 3 can be obtained through shear experiment performed with successive loadings and unloadings by using the decrease of the initial shear modulus (i.e. at c ! 0) and the decrease of the slope of the curve ln G ¼ f ðc 4 Þ for c ! 1. • For every level of damage calculated, the corresponding thermodynamical force can be evaluated through the integration of Eqs. (14) and (16). Knowing the couples (d i ; Y i ) for different levels of strain, an interpolation enables to propose the evolution laws As in the case of composites, the evolution laws can be written where Y 0 and Y c are material parameters (Y 0 is the threshold under which d i may not grow).
The identification was performed with Mathematica. A verification procedure was also performed with the simulation of a shear experiment with loads and unloads. Because of calculationÕs complexity, only the first damage parameter was used. The equations of the problems are • The state laws • The evolution law with the material constants: The results of simulation are plotted in Fig. 5 and are compared with experimental results. Symbols are used for experimental data, whereas full lines are used for simulation. This curve show the ability of the model to estimate the degradation of the material properties. The stiffening effect for high strain is here a little bit overestimated. For a more precise modeling, the solution is either to let d 3 free or to identify C 3 on a large range of strain to obtain a ''mean value''.

Implementation of the model
The behaviour relation is implemented in the prototype code we are developing. It requires its transformation in a rated law form and some modifications of the algorithm because of the strong non-linearity of this behaviour relation. Details are not exposed here (the interested reader could refer to [11,12]). But a result of simulation for a tension test is exposed. Thus Fig. 6 shows the comparison of the response in terms of force versus stretch (k ¼ L=L 0 ) for the classical model of Hart-Smith and with introduction of the first variables of damage d 1 . Fig. 7 shows the evolution of this variable versus stretch. It is important to note that the chosen evolution law is a law with delay effect (see [10,13]), which is an extension of dynamical model and which enables to regularize the damage model. The general form is where d ¼ f ðY Þ is the ''static'' evolution law.

General formulation for large strain problems
The key of the method is the introduction of a formulation which gives a common framework to the constitutive laws and the equilibrium equations. Rated constitutive laws and varying ''rotated'' quantities, especially the displacement rate V are used. More details can be found in [12]. where subscript S and AS denotes respectively the symmetric and antisymmetric part of a tensor, and I d is the identity tensor.
The chosen framework is the corotational space frame defined by We can define a strain rate _ R R and its conjugate stress r, such that the virtual power of internal forces can be written; These new quantities defined by are X 0 -material since they are invariant under a change in reference frame. A motivation for such a choice is given in [17].

Formulation in terms of rotated quantities
One of the advantages of the corotational frame is that the problem can be completely formulated in terms of rotated quantities. The second is that the resulting formulation is quite similar to the same problem in small strains. We split the strain rate into the sum of an elastic part R e and an inelastic part R p . The inelastic strain is to be separated from the other internal variables (denoted collectively by X). We denote the conjugate of X by Y so that the mechanical dissipation is given by: tr½r _ R R p À Y _ X X. Let quðr; XÞ be the volumic energy, where q is the density and uðr; XÞ the free energy. Using classical hypotheses, the following state equations are obtained: Since the elasticity is linear, _ R R e ¼ K À1 _ r r, where K is HookeÕs tensor. Moreover, under the conventional hypothesis: where the operator G is a given material characteristic. Models with internal variables are of particular interest because of their simplicity. Their evolution law is a differential equation: where B is a given operator which depends on the material and which must be positive (BðX Þ; X > 0) in order for the second law of thermodynamics to hold. The formulation in terms of rotated quantity can then be written: 2. The state and constitutive equations

The kinetic equations defining the configuration
where and the operator curl is defined by 4. The compatibility equation Once the reference problem is solved, it is possible to find the displacement U by solving the two differential equations: 3.2. Resolution of the problem by means of the large time increment method

Presentation of the algorithm
The LATIN method is an iterative method, which takes into account the whole loading process in a single time increment. It rests on three principles.
The first principle. The idea is to separate the difficulties by splitting the equations into two groups. The first one ðA d Þ contains equations which are linear and possibly global in space, the second one ðCÞ contains equations which are local and possibly non-linear. Typically, we put in ðA d Þ the compatibility equation, the equilibrium equation and the state laws and we put in ðCÞ the constitutive equation and the equations defining the configuration.
The second principle. We use an iterative two stage scheme to solve successively the first and the second group over the whole loading process. It is represented in Fig. 8.
It can be detailed as follow: At a given iteration, the data are s n , defined over X 0 Â ½0; T . In the first stage or the local stage, we determine the quantities s nþ 1 2 which satisfy the constitutive relations exactly. Then, we can find A and Q which define the new configuration. In this stage, we only solve local differential equations.
In the second stage or the global stage, we ''freeze'' the configuration parameters and then we find s nþ1 which verify the equilibrium equation. In this stage, we only solve linear problems. For both stage, we need search directions, which are parameters of the method.
The third principle. The key is to introduce an appropriate space-time representation for the unknowns. In finite strains, we seek correction of the unknowns (for instance V nþ1 À V n ) as a product of time functions by space functions, thus The problem is solved iteratively: we alternatively compute the best space functions when time functions are fixed and the best time functions when space functions are fixed. The initial solution is obtained through a calculation in small strain.

First remark
Concerning the space aspects, the problem is solved using the finite element method. The discretized quantity is here the displacement rate V , and not the displacement U as it is usually the case in finite elements.
As a consequence, in the global stage: • if the third principle is not applied, we solve for every time t 2 ½O; T a linear problem: KðtÞV ðtÞ ¼ F ðtÞ where K, is a matrix similar to the one obtained in the classical finite element method, • if P3 is applied, we solve: -one differential equation to obtain aðtÞ, -one linear problem where K int contains integrals on both time and space.

Treatment of incompressibility
Simulation are done for plane strain with four-node quadrilateral elements. As in the standard case, the use of such elements to discretize the rotated displacement rate V in the case of quasi-incompressibility can lead, to ''volumetric locking'' or in the context of LATIN method to non-convergence of the algorithm. Many techniques have been developed to overcome this difficulty: reduced integration with stabilization matrix [2,14], assumed strain method or B-Bar method which can derive from mixed method [8,19], assumed stress method [16], enhanced strain element [5,6]. The chosen method is the B-Bar method, developed by Hughes [8] for its simplicity. The method is applied to the matrix derivation B, which in our case links the strain rate _ R R to the rotated displacement rate V through Eq. (38): The matrix is split into deviatoric and dilatational parts. The deviatoric part is classically computed contrarily to the dilatational one for which we take only the mean contribution on the element. The new matrix used to solve the finite element problem is

Tensile test
First of all, in order to make the reader understand the specificity of the LATIN method and especially the non-incremental aspect, results in the case of a tensile test are shown. The material is supposed to be linear isotropic elastic with a YoungÕs modulus E ¼ 7:2 Â 10 10 Pa and m ¼ 0:3. The conditions of experiment and the solution of the simulation are shown in Figs. 9 and 10. Fig. 4 enables to follow the behaviour of the algorithm.
In Fig. 11, the force versus displacement curve is plotted for several iterations. It enables to see the initial solution obtained with a calculation with the small strain assumptions and the successive corrections until convergence is attained. At every iteration, we have an approximated solution on the whole loading process. Fig. 12 shows the evolution of the energy error indicator which provides a measure of convergence of the method.

Shear test
In this part, a simulation of a shear experiment is performed, as described in Fig. 13. The used behaviour is the same than above (elastic) but nearly incompressible with m ¼ 0:499. Two cases are shown: ratio h=L ¼ 1 and ratio h=L ¼ 0:4.
Figs. 14 and 15 show respectively the most deformed solution we can obtain with the industrial code ABA-QUS and a more deformed solution obtained with the LATIN method.   8 QUS and a more deformed solution obtained with the LATIN method.
A premature stop of calculations (dependent on the mesh) combined with the occurrence of some ''cusps'' make presume the occurrence of local instabilities. Similar results have been obtained with other models such as Hart-Smith one. Thus such an instability is not essentially due to the material behaviour. The development of a criterion in order to detect this type of instabilities is in progress.

Conclusion
The studied structures with elastomer are able to undergo very large deformation up to 600% or 700%, which involves damage of the elastomer. The proposed model, which rests on the continuum damage mechanic, also a phenomenological approach, is both simple with few parameters to identify and efficient to take into account the degradation of material properties. Moreover, it enables to take into account the non-symmetry of the response of the damaged material in tension or in compression which is therefore a worthy aspect never taken into account to our knowledge. Concerning the simulation, even if the LATIN method enables to perform shear experiment simulation up to higher strains than a classical incremental code ABAQUS, all the difficulties are not overcome. Hence current work is done on a criterium to detect the local instability which seems to occur. Of course, one more difficulty is always due to element distorsion which could appear when deformation becomes large, even if the sensibility is diminished through the use of a material formulation.