# Modeling and simulation of damage in elastomer structures at high strains

Xavier Aubard, Pierre-Alain Boucard, Pierre Ladevèze, Sylvain Michel

► To cite this version:

Xavier Aubard, Pierre-Alain Boucard, Pierre Ladevèze, Sylvain Michel. Modeling and simulation of damage in elastomer structures at high strains. Computers & Structures, 2002, 80 (27-30), pp.2289-2298. 10.1016/S0045-7949(02)00241-9. hal-01698052

HAL Id: hal-01698052

https://hal.science/hal-01698052

Submitted on 13 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling and simulation of damage in elastomer structures at high strains

X. Aubard<sup>a</sup>, P.-A. Boucard<sup>b</sup>, P. Ladevèze<sup>b,*</sup>, S. Michel<sup>b</sup>

<sup>a</sup> SNECMA-Moteurs, Le Haillan BP37, 33165 St-Medard en Jalles, France

<sup>b</sup> ENS Cachan/CNRS/Université Paris 6, Laboratoire de Mécanique et Technologie (LMT-Cachan), 61 Avenue du Pree’sident Wilson, 94235 Cachan Cedex, France

## Abstract

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.

*Keywords: Elastomers; Damage; Finite transformations; LATIN method*

## 1. 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

<footnote><sup>*</sup> Corresponding author. Tel.: +33-1-47-40-22-38; fax: +33-1-47-40-27-85.</footnote>

<footnote>E-mail addresses: ladeveze@lmt.ens-cachan.fr, pierre.ladeveze@lmt.ens-cachan.fr (P. Ladevèze).</footnote>
work on both constitutive relation and numerical sim-
ulation.

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 indus-
trial 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 cross-
critical 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.

## 2. Behaviour of the elastomer

2.1. *Experimental results*

The behaviour of the elastomers is commonly rep-
resented by hyperelastic models, characterized by the
strain energy function W, which is often considered to
be dependent on the invariants I₁, I₂, I₃ 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 phe-
nomena are not taken into account, essentially a de-
crease 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(γ) and ln G = f(γ⁴) (G is the shear modulus and
γ = U_d/h is the apparent shear strain).

<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison of several behaviour relations.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Experimental results in shear test case with successive loads and unloads in the plane G versus γ².</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Experimental results in shear test case with successive loads and unloads in the plane ln G versus γ⁴.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Force versus displacement for a hydrostatic tension-compression test.</figcaption></figure>

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.

## 2.2. Modeling

### 2.2.1. 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:

$$ \rho\psi = C_1 \int_0^{I_1^{-3}} \exp[C_3(I_1^2)] dI_1' + 3C_2 \ln\left(\frac{I_2}{3}\right) \quad (1) $$

where $I_1$ and $I_2$ are the invariants of tensor $\mathbb{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

$$ S = \frac{\partial W}{\partial \mathbb{E}} = 2C_1 \exp(C_3 I_1^2) + \frac{C_2}{I_2} (I_1 \mathbb{I} - \mathbb{C}) \quad (2) $$

### 2.2.2. 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

$$ \lim_{\gamma \to 0} G = 2C_1 + C_2 \quad (3) $$

$$ \lim_{\gamma \to \infty} \left( \frac{\ln G}{\gamma^4} \right) = C_3 \quad (4) $$

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, $C_i$ becomes $C_i^0(1 - d_i)$: if the material is safe $d_i = 0$ so $C_i = C_i^0$; if the material is completely damaged $d_i = 0$ and so $C_i = 0$.

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 split into two parts: one for tensile for which damage variables act, the other one for compression for which initial characteristics are used, also:

$$ W(\mathbb{E}, d_i) = W^+(\mathbb{E}, d_i) + W^-(\mathbb{E}) $$

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 $\mathbb{E}$. They are defined by writing the incompressibility constraint, that is

$$ \det(\mathbb{F}^\dagger \mathbb{F}) = 1 \quad \text{or} \quad \det((\mathbb{F}^\dagger \mathbb{F})^{-1}) = 1 \quad (5) $$

In terms of the new quantities, $I'_1$ and $I'_2$, this equation becomes

$$ (I'_1)^2 + 2I'_1 - 4\operatorname{tr}[\mathbb{E}\mathbb{E}] + 16\det\mathbb{E} = 0 \quad (6) $$

and

$$ (I'_2)^2 + 2I'_2 - 4\operatorname{tr}[\mathbb{G}\mathbb{G}] + 16\det\mathbb{G} = 0 \quad (7) $$

where $\mathbb{G}$ is the Almansi-Lagrange tensor such as

$$ \mathbb{G} = \frac{1}{2}(\mathbb{I}_d - (\mathbb{A}^\dagger \mathbb{A}^{-1})) $$

Taking into account that for small and mean strains, the terms $\det \mathbb{E}$ and $\det \mathbb{G}$ can be neglected, the quantities used to separate tension and compression ($I_1^\dagger, I_1^-, I_2^\dagger$ and $I_2^-$) are defined by the following equations:

$$ (I_1^+)^2 + 2I_1^+ - 4\operatorname{tr}[\langle\mathbb{E}\rangle_+ \langle\mathbb{E}\rangle_+] = 0 \quad (8) $$

$$ (I_1^-)^2 + 2I_1^- - 4\operatorname{tr}[\langle\mathbb{E}\rangle_- \langle\mathbb{E}\rangle_-] = 0 \quad (9) $$

$$ (I_2^+)^2 + 2I_2^+ - 4\operatorname{tr}[\langle\mathbb{G}\rangle_+ \langle\mathbb{G}\rangle_+] = 0 \quad (10) $$

$$ (I_2^-)^2 + 2I_2^- - 4\operatorname{tr}[\langle\mathbb{G}\rangle_- \langle\mathbb{G}\rangle_-] = 0 \quad (11) $$

where $\langle\bullet\rangle_-$ is the negative part (associated to negative eigenvalues) of $\bullet$ and $\langle\bullet\rangle_+$ the positive part (associated to positive eigenvalues). And the new strain energy
function for hyperelastic model with unilateral damage, can be written

$$
\rho\psi(\mathbb{E}, d_i) = \int_0^{I_1^+} (1 - d_1)C_1 \exp[C_3(1 - d_3)u^2] du \\
+ 3(1 - d_2)C_2 \log \left(1 + \frac{I_2^+}{3}\right) + \int_0^{I_1^-} C_1 \\
\times \exp[C_3 u^2] du + 3C_2 \log \left(1 + \frac{I_2^-}{3}\right) \quad (12)
$$

The state laws associated with this model are classically
obtained:

$$
\begin{align}
\mathcal{S} &= \frac{\partial W}{\partial \mathbb{E}} \nonumber \\
&= \frac{4C_1(1-d_1)}{I_1^+ + 1} \exp \left[ C_3(1-d_3)(I_1^+)^2 \right] \langle \mathbb{E} \rangle_+ \nonumber \\
&\quad + \frac{4C_1}{I_1^- + 1} \exp \left[ C_3(I_1^-)^2 \right] \langle \mathbb{E} \rangle_- \nonumber \\
&\quad + \frac{12C_2(1-d_2)}{(3+I_2^+)(I_2^+ + 1)} \langle \mathbb{G} \rangle_+ (1+2\mathbb{E})^{-2} \nonumber \\
&\quad + \frac{12C_2}{(3+I_2^-)(I_2^- + 1)} \langle \mathbb{G} \rangle_- (1+2\mathbb{E})^{-2} \tag{13}
\end{align}
$$

and with $Y_i = -\partial W / \partial d_i$

$$
Y_1 = \int_0^{I_1^+} C_1 \exp[C_3(1 - d_3)u^2] du \quad (14)
$$

$$
Y_2 = 3C_2 \log \left(1 + \frac{I_2^+}{3}\right) \qquad (15)
$$

$$
Y_3 = \int_0^{I_1^+} (1 - d_1) C_1 C_3 u^2 \exp [C_3 (1 - d_3) u^2] du \quad (16)
$$

$Y_i$ are the thermodynamical forces associated with $d_i$. By the addition of evolution law for the damage variables ($\dot{d}_i = f(Y_i)$), the model would be complete.

2.3. Identification

Before using the proposed model, the three elastic
parameters C₁, C₂ and C₃ 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 prop-
erties (Eq. (4)). The different stages of the procedure are
summarized below. Let us precise that as a first approxi-
mation, the term with I₂ 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 ex-
perimental results are in good agreement, also here
C₂ = 0.

• As the gradient deformation in a shear test case can be written

$$
\mathbb{F} = \begin{bmatrix} 1 & \gamma & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix} \tag{17}
$$

the new quantities defined in (8) and (9) verify

$$
\begin{align}
I_1^+ &= -1 + \sqrt{1 + 4 \left( \frac{\gamma^2}{4} + \frac{\gamma}{2} \sqrt{1 + \gamma^2} \right)^2} \nonumber \\
I_1^- &= -1 + \sqrt{1 + 4 \left( \frac{\gamma^2}{4} - \frac{\gamma}{2} \sqrt{1 + \gamma^2} \right)^2} \tag{18}
\end{align}
$$

Thus, stress tensor σ and the thermodynamical forces
Y₁ and Y₃ can be written in terms of the shear rate γ
through Eqs. (13), (14) and (16).

• Let $\tau$ be the measure of the global force. The shear modulus can be expressed in terms of $\gamma$, since:

$$
G_{\text{sec}} = \frac{\tau}{\gamma} = S_{22} + \frac{1}{\gamma} S_{12} \quad (19)
$$

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:

$$
\lim_{\gamma \to 0} G = C_1(2 - d_1) \quad (20)
$$

$$
\lim_{\gamma \to \infty} \left( \frac{\ln G}{\gamma^4} \right) = C_3(1 - d_3) \quad (21)
$$

Also, $d_1$ and $d_3$ can be obtained through shear ex-
periment performed with successive loadings and
unloadings by using the decrease of the initial shear
modulus (i.e. at $\gamma \to 0$) and the decrease of the slope
of the curve $\ln G = f(\gamma^4)$ for $\gamma \to \infty$.

* 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 $d_i = f(Y_i)$. As in the case of composites, the evolution laws can be written

$$
d = \min \left( \left\langle \frac{\sqrt{Y} - \sqrt{Y_0}}{\sqrt{Y_c}} \right\rangle_+ , 1 \right) \quad (22)
$$

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 un-
loads. Because of calculation's complexity, only the first
damage parameter was used. The equations of the problems are

* The state laws

$$ \sigma = \Delta \left\{ \frac{4C_1}{I_1^+ + 1} \exp [C_3(I_1^+)^2] \langle \mathbb{E} \rangle_+ \right\} \quad (23) $$

$$ + \frac{4C_1}{I_1^- + 1} \exp [C_3(I_1^-)^2] \langle \mathbb{E} \rangle_- \right\} \Delta^t \quad (24) $$

$$ Y_1 = \int_0^{I_1^+} \exp[C_3 u^2] du \quad (25) $$

* The evolution law

$$ d_1 = \min \left[ \left\langle \frac{\sqrt{Y_1} - \sqrt{Y_0}}{\sqrt{Y_c}} \right\rangle_+ , 1 \right] \quad (26) $$

with the material constants: $C_1 = 0.185$ MPa, $C_3 = 7.75 \times 10^{-4}$, $Y_c = 2.164$ MPa, $Y_0 = 0.856$ MPa. 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 shows 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”.

## 2.4. 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

<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison experiment results/calculation: load/unload in shear test case.</figcaption></figure>

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 ($\lambda = 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

$$ \dot{d} = \frac{k}{a} \exp [-a \langle f(Y) - d \rangle_+] \quad (27) $$

where $d = f(Y)$ is the “static” evolution law.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Force versus displacement in traction.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Evolution of $d_1$.</figcaption></figure>
### 3. Numerical approach

#### 3.1. 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 $\underline{\mathbf{V}}$ are used. More details can be found in [12].

##### 3.1.1. Material rotated quantities

The used notations are as follows: $M_0$ and $M$ are respectively the initial and the final positions of a point, $\underline{\mathbf{U}}$ is the displacement. The classical variables used in the context of finite deformation are:

* the deformation gradient: $\mathbf{F} = \mathbf{I}_d + (d\underline{\mathbf{U}}/dM_0)$,

* the Cauchy-Green stress tensor: $\mathcal{C} = \mathbf{F}^t\mathbf{F}$,

* the Green-Lagrangian strain: $\mathbf{E} = \frac{1}{2}(\mathbf{F}^t\mathbf{F} - \mathbf{I}_d)$,

* the usual strain rate (Eulerian): $\mathcal{D} = (\mathbf{F}\mathbf{F}^{-1})_S$,

* the second Piola-Kirchoff stress (Lagrangian): $\mathcal{S}$,

* the Cauchy stress tensor (Eulerian): $\boldsymbol{\sigma}_c$,

where subscript S and AS denotes respectively the symmetric and antisymmetric part of a tensor, and $\mathbf{I}_d$ is the identity tensor.

The chosen framework is the corotational space frame defined by

$$ \dot{\mathbf{R}}\mathbf{R}^t = (\mathbf{F}\mathbf{F}^{-1})_{\text{AS}} \quad \text{with} \quad \mathbf{R}|_{t=0} = \mathbf{I}_d \qquad (28) $$

We can define a strain rate $\dot{\boldsymbol{\Sigma}}$ and its conjugate stress $\boldsymbol{\sigma}$, such that the virtual power of internal forces can be written;

$$ P_i^* = - \int_{\Omega_0} \mathrm{tr}[\boldsymbol{\sigma} \dot{\boldsymbol{\Sigma}}^*] d\Omega_0 \qquad (29) $$

These new quantities defined by

$$ \dot{\boldsymbol{\Sigma}} = \mathbf{R}^t \mathcal{D} \mathbf{R} \quad \text{and} \quad \boldsymbol{\sigma} = \det(\mathbf{F}) \mathbf{R}^t \boldsymbol{\sigma}_c \mathcal{R} \qquad (30) $$

are $\Omega_0$-material since they are invariant under a change in reference frame. A motivation for such a choice is given in [17].

##### 3.1.2. 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 $\boldsymbol{\Sigma}_e$ and an inelastic part $\boldsymbol{\Sigma}_p$. The inelastic strain is to be separated from the other internal variables (denoted collectively by $\mathcal{X}$). We denote the conjugate of $\mathcal{X}$ by $\mathcal{Y}$ so

that the mechanical dissipation is given by: tr[$\boldsymbol{\sigma}\dot{\boldsymbol{\Sigma}}_p$] = $\mathcal{Y}\mathcal{X}$. Let $\rho\varphi(\boldsymbol{\sigma}, \mathcal{X})$ be the volumic energy, where $\rho$ is the density and $\varphi(\boldsymbol{\sigma}, \mathcal{X})$ the free energy. Using classical hypotheses, the following state equations are obtained:

$$ \boldsymbol{\Sigma}_e = \rho \frac{\partial \varphi(\boldsymbol{\sigma}, \mathcal{X})}{\partial \boldsymbol{\sigma}} \quad \text{and} \quad \mathcal{Y} = \rho \frac{\partial \varphi(\boldsymbol{\sigma}, \mathcal{X})}{\partial \mathcal{X}} \qquad (31) $$

Since the elasticity is linear, $\dot{\boldsymbol{\Sigma}}_e = \mathbf{K}^{-1}\dot{\boldsymbol{\sigma}}$, where $\mathbf{K}$ is Hooke's tensor. Moreover, under the conventional hypothesis:

$$ \mathcal{Y} = \rho \frac{\partial \varphi(\boldsymbol{\sigma}, \mathcal{X})}{\partial \mathcal{X}} = \mathbf{G}(\mathcal{X}) \qquad (32) $$

where the operator $\mathbf{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:

$$ \begin{bmatrix} \dot{\boldsymbol{\Sigma}}_p \\ -\mathcal{X} \end{bmatrix} = \mathbf{B} \begin{pmatrix} \boldsymbol{\sigma} \\ \mathcal{Y} \end{pmatrix} \qquad (33) $$

where $\mathbf{B}$ is a given operator which depends on the material and which must be positive ($\mathbf{B}(\mathcal{X})$, $X > 0$) in order for the second law of thermodynamics to hold.

The formulation in terms of rotated quantity can then be written:

Find $(\boldsymbol{\Sigma}_p, \mathcal{X}, \mathcal{Y}, \boldsymbol{\sigma}, \mathcal{A}, \mathcal{Q})$ on $\Omega_0 \times [0, T]$ satisfying

1. Equilibrium equations

$$ \forall \mathcal{V}^* \in \mathcal{V}_0^* = \{ \mathcal{V}^*, \mathcal{V}^*|_{\mathcal{I}_0, \mathcal{I}_0} = 0 + \text{"regularity"} \} \times \int_{\Omega_0} \mathrm{tr}[\boldsymbol{\sigma}\dot{\boldsymbol{\Sigma}}(\mathcal{V}^*)] d\Omega_0 + P_e^* \qquad (34) $$

2. The state and constitutive equations

$$ \dot{\boldsymbol{\Sigma}}_e = \mathbf{K}^{-1} \dot{\boldsymbol{\sigma}} \quad \mathcal{Y} = \mathbf{G}(\mathcal{X}) \quad \dot{\boldsymbol{\Sigma}} = \dot{\boldsymbol{\Sigma}}_e + \dot{\boldsymbol{\Sigma}}_p \qquad (35) $$

$$ \begin{bmatrix} \dot{\boldsymbol{\Sigma}}_p \\ -\mathcal{X} \end{bmatrix} = \mathbf{B} \begin{pmatrix} \boldsymbol{\sigma} \\ \mathcal{Y} \end{pmatrix} $$

with $\boldsymbol{\Sigma}_p|_{t=0} = \mathcal{O}$ and $\mathcal{X}|_{t=0} = \mathcal{O}$.

3. The kinetic equations defining the configuration $\mathcal{A}, \mathcal{Q}$

$$ \dot{\mathcal{A}} = \dot{\boldsymbol{\Sigma}} \mathcal{A} \quad \mathcal{A}|_{t=0} = \mathbf{I}_d \qquad (36) $$

$$ \dot{\mathcal{Q}} = \mathcal{L}(\dot{\boldsymbol{\Sigma}}, \mathcal{Q}) \quad \mathcal{Q}|_{t=0} = \mathcal{O} \qquad (37) $$

where

$$ \mathcal{L}(\dot{\boldsymbol{\Sigma}}, \mathcal{Q}) = \overline{\mathrm{curl}} \dot{\boldsymbol{\Sigma}} - \{[\mathrm{tr}[\dot{\boldsymbol{\Sigma}}\mathcal{Q}] - \mathrm{tr}[\dot{\boldsymbol{\Sigma}}]]\mathrm{tr}[\mathcal{Q}]\}|_{\mathcal{I}_d} + [2\mathrm{tr}[\mathcal{Q}]|_{\mathcal{I}_d} - \mathcal{Q}]\dot{\boldsymbol{\Sigma}} + [\mathrm{tr}[\dot{\boldsymbol{\Sigma}}]|_{\mathcal{I}_d} - \dot{\boldsymbol{\Sigma}}]\mathcal{Q} \} $$

and the operator $\overline{\mathrm{curl}}$ is defined by

$$ [\overline{\mathrm{curl}} \boldsymbol{\bullet}]_{\wedge} = 2 \left[ \frac{d \boldsymbol{\bullet}}{d \mathcal{X}} \mathcal{A}^{-1} \right]_{\text{AS}} $$
4. The compatibility equation

$$
\begin{gathered}
\exists \underline{V} \in \mathcal{V} = \{ \underline{V}, \underline{V}|_{\partial_1 \Omega} = \mathbb{R}^t \dot{\underline{U}}_d + \text{"regularity"} \}, \\
\dot{\underline{\Sigma}} = \dot{\underline{\Sigma}}(\underline{V}) = \left[ \frac{\partial \underline{V}}{\partial \underline{M}_0} \mathbb{A}^{-1} - [\underline{V} \wedge \mathbb{Q}] \right]_S
\end{gathered}
\qquad (38)
$$

Once the reference problem is solved, it is possible to
find the displacement $\underline{U}$ by solving the two differential
equations:

$$
\begin{align}
\dot{\underline{R}} &= \mathbb{R}[\operatorname{curl} \underline{V} + \frac{1}{2}(\mathbb{Q}\underline{V} - \operatorname{det}(\mathbb{Q})\underline{V})]_{\wedge}, \quad \underline{R}|_{t=0} = \underline{l}_d \nonumber \\
\dot{\underline{U}} &= \underline{R}\underline{V} \tag{39}
\end{align}
$$

3.2. Resolution of the problem by means of the large time increment method

3.2.1. 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 diffi-
culties 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 (Γ) 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
(Γ) 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 $\Omega_0 \times [0, T]$. In the first stage or
the local stage, we determine the quantities $s_{n+\frac{1}{2}}$ which
satisfy the constitutive relations exactly. Then, we can

<figure><img src="image_9.png" /><figcaption>Fig. 8. Scheme of the method.</figcaption></figure>

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 $\underline{V}_{n+1} - \underline{V}_n$) as a product of time functions by space functions, thus

$$
\underline{V}_{n+1} - \underline{V}_n = a(t) \underline{W}(\underline{M}_0)
$$

The problem is solved iteratively: we alternatively
compute the best space functions when time functions
are fixed and the best time functions when space func-
tions are fixed. The initial solution is obtained through a
calculation in small strain.

3.3. Remarks on the numerical resolution

3.3.1. First remark

Concerning the space aspects, the problem is solved
using the finite element method. The discretized quantity
is here the displacement rate $\underline{V}$, and not the displace-
ment $\underline{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 \in [O, T]$ a linear problem:

$$
\mathcal{H}(t) \underline{V}(t) = \underline{F}(t)
$$

where $\mathcal{H}$, 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

$$
\mathcal{H}_{\text{int}} \underline{V} = \underline{F}
$$

where $\mathcal{H}_{\text{int}}$ contains integrals on both time and space.

3.3.2. 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 $\underline{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 diffi-
culty: 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 $\mathcal{B}$, which in our case links the strain rate $\dot{\boldsymbol{\Sigma}}$ to the rotated displacement rate $\underline{\boldsymbol{V}}$ through Eq. (38):

$$ \begin{bmatrix} \dot{\boldsymbol{\Sigma}}_{11} \\ \dot{\boldsymbol{\Sigma}}_{22} \\ \dot{\boldsymbol{\Sigma}}_{33} \\ 2\dot{\boldsymbol{\Sigma}}_{12} \end{bmatrix} = \mathcal{B} \begin{bmatrix} V_1 \\ V_2 \end{bmatrix} = \begin{bmatrix} \mathcal{B}11 & \mathcal{B}12 \\ \mathcal{B}21 & \mathcal{B}22 \\ 0 & 0 \\ \mathcal{B}22 - \mathcal{B}12 & \mathcal{B}11 - \mathcal{B}21 \end{bmatrix} \begin{bmatrix} V_1 \\ V_2 \end{bmatrix} \quad (40) $$

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

$$ \bar{\mathcal{B}} = \mathcal{B}_{\text{dev}} + \bar{\mathcal{B}}_{\text{dil}} $$

# 4. Numerical examples

## 4.1. 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 \times 10^{10}$ Pa and $v = 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.

## 4.2. 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 $v = 0.499$. Two cases are shown: ratio $h/L = 1$ and ratio $h/L = 0.4$.

*First case:* ratio $h/L = 1$.

Figs. 14 and 15 show respectively the most deformed solution we can obtain with the industrial code ABAQUS and a more deformed solution obtained with the LATIN method.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Description of the tension experiment.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Solution for $U_d = 8$ cm.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Evolution of force versus displacement during iterations.</figcaption></figure>

*Second case:* ratio $h/L = 0.4$.

As below, Figs. 16 and 17 show the most deformed solution we can obtain with the industrial code ABA-
<figure><img src="image_13.png" /><figcaption>Fig. 12. Evolution of the error indicator versus iterations.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Description of the shear experiment.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Solution at maximal strain ($U_{max}$ = 0.99) with ABAQUS.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Solution with Latin $U_d$ = 1.36.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Solution at maximal strain ($U_{max}$ = 0.46) with ABAQUS.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Solution with Latin $U_d$ = 0.5.</figcaption></figure>

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.

## 5. 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.

## Acknowledgement

The authors are grateful to CNES for its financial support.
References

[1] Andrieux F, Saanouni K, Sidoroff F. Sur les solides hyperélastiques à compressibilité induite par l'endommagement. CRAcad Sci Paris 1997;24(Série IIb):281-8.

[2] Belytschko T, Bindeman LP. Assumed strain stabilization of the 4-node quadrilateral with 1-point quadrature for non-linear problems. Comput Meth Appl Mech Eng 1991;88:311-40

[3] Boucard P-A, Ladevèze P, Poss M, Rougée P. A non incremental approach for large displacement problems. Comput Struct 1997;64(1-4):499-508.

[4] Chen J-S, Pan C, Wu C-T. Large deformation analysis of rubber based on a reproducing kernel particle method. Comput Mech 1997;(19):211-27.

[5] Crisfield MA, Moita GF, Gelenic G, Lyons LPR. Enhanced lower-order elements formulations for large strains. In: Owen DRJ, Onate E, Hinton E, editors. Proceedings of the Fourth International Conference on Computational Plasticity: Fundamentals and applications, Barcelona 3-6 April 95. Swansea: Pineridge Presse; 1995.

[6] De-Souza-Neto EA, Peric D, Dutko M, Owen DRJ. Design of simple low order finite elements for large strain analysis of nearly incompressible solids. Int J Solids Struct 1996;33(20-22):3277-96.

[7] Glaser S, Armero F. On the formulation of enhanced strain finite elements in finite deformations. Eng Comput 1997;14(7):759-91.

[8] Hugues TJR. Generalization of selective integration procedures to anisotropic and non-linear media. Int J Numer Meth Eng 1980;15:1413-8.

[9] Koh BC, Kikuchi N. Improved hourglass control in anisotropic linear elasticity. Comput Meth Appl Mech Eng 1987;65:1-46.

[10] Ladevèze P. A damage computational approach for composites: basic aspects and micromechanical relations. Comput Mech 1995;(17):142-50.

[11] Ladevèze P. On the theory of large (visco)-plastic deformations: a new class of models with and without damage and a computational strategy. Internal report 220 (in French), LMT, 1998.

[12] Ladevèze P. Nonlinear computational structural mechanics. New York: 1999.

[13] Ladevèze P, Moes N, Douchin B. Constitutive relation error estimators for (visco) plastic finite element analysis with softening. Comput Meth Appl Mech Eng 1999;176: 247-64.

[14] Liu WG, Ong JS-J, Uras RA. Finite element stabilization matrices—a unification approach. Comput Meth Appl Mech Eng 1985;53:13-46.

[15] Michel S. Modélisation et simulation de structures élastomères endommageables en grandes deformations. AUM-AFM, 14th ed. Congrès français de mécanique, Toulouse 30 août-3 septembre 1999.

[16] Pian THH, Sumihara K. Rational approach for assumed stress finite elements. Int J Numer Meth Eng 1984;20:1685-95.

[17] Rougée P. A new Lagrangian intrinsic approach to large deformations in continuous media. Eur J Mech A/Solids 1991;10(1):15-39.

[18] Simo JC. On a fully three-dimensional finite-strain visco-elastic damage model: formulation and computational aspects. Comput Meth Appl Mech Eng 1987;60:153-73.

[19] Simo JC, Rifai MS. A class of mixed assumed strain methods and the method of incompatible modes. Int J Numer Meth Eng 1990;29:1595-638.