# A new isotropic hyperelastic strain energy function in terms of invariants and its derivation into a pseudo-elastic model for Mullins effect: application to finite element analysis

Laurent Gornet, R. Desmorat, Gilles Marckmann, Pierre Charrier

► To cite this version:

Laurent Gornet, R. Desmorat, Gilles Marckmann, Pierre Charrier. A new isotropic hyperelastic strain energy function in terms of invariants and its derivation into a pseudo-elastic model for Mullins effect: application to finite element analysis. Stephen Jerrams, Niall Murphy. Constitutive Models for Rubber VII, VII, CRC Press, pp.265-271, 2012, 978-0-415-68389-0. 10.1201/b11687-49. hal-01004687

HAL Id: hal-01004687

https://hal.science/hal-01004687

Submitted on 6 Oct 2016

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>
# A new isotropic hyperelastic strain energy function in terms of invariants and its derivation into a pseudo-elastic model for Mullins effect: application to finite element analysis

L. Gornet & G. Marckmann

GeM, UMR CNRS 6183, Ecole Centrale Nantes, 1 rue de la Noë, BP 92101-44321 Nantes, Cedex 3, France.

R. Desmorat

LMT- Cachan ENS Cachan/CNRS/UPMC/PRES Univers Sud Paris, 61 avenue du Président Wilson, F-94230 Cachan Cedex, France.

P. Charrier

Modyn Trelleborg, Zone ind. de Carquefou, BP 419, 44474 Carquefou Cedex - France

**ABSTRACT:** The present paper focuses on static stiffness modelling of rubber materials for multiaxial loadings with a minimal number of material parameters in order to ensure robustness of both identification and Finite Element analysis. A physically motivated isotropic constitutive model is proposed for multiaxial large stretch rubber deformation. This model is expressed in terms of classical independent strain invariants. Only three parameters are needed to successfully represent both Treloar and Kawabata experiments (Kawabata et al. 1981) with a response equivalent to the one obtained by Ogden six parameters model (Ogden 1972). Mullins effect modelling is finally derived according to a damage mechanics approach (Chagnon et al. 2004).

## 1 INTRODUCTION

The aim of the present paper is to develop new constitutive models for the multiaxial loadings of elastomers which are both physically motivated, well adapted for numerical problems, and accurate for rubber materials. Hyperelastic models dedicated to rubber materials can be classified into two types of strain energy formulations. The first kind of models is issued from mathematical developments such as the well-known Rivlin series or as the Ogden model (Ogden 1972). The second kind of models is the one developed from physical motivations. Such models are based on both physics of polymer chains networks and statistical methods. In this work, we propose a strain energy function expressed in terms of independent strain invariants. The proposed model (GD) is dedicated to multiaxial loadings. This model is successfully identified on both Treloar and Kawabata experiments (Treloar 1944, Kawabata et al. 1981). We present a bridge between the phenomenological strain energy formulation and the physical motivation of the proposed (GD) model. In the Finite Element context, the proposed model can easily be implemented because of its strain invariants formulation. This has been done in computer codes like Cast3M-CEA (French Atomic Agency) and AB-AQUS. Finite Element predictions of a roll restrictor developed by Trelleborg automotive have been carried out with ABAQUS.

The continuous damage approach is applied to soften by damage the material parameters of this initial model in order to describe Mullins effect. A new model for Mullins effect (GDM model) is therefore derived. The identification of the GDM model is

made thanks to experimental data from (Marckmann & G., E. Verron 2006).

## 2 STRAIN ENERGY FUNCTIONS

In this section, the strain energy functions that define the proposed model (GD) are briefly recalled. Links between this model and the eight-chains model are highlighted. A phenomenological energy function expressed in term of the second invariant $I_2$ is also introduced based on physical motivations.

### 2.1 Incompressible GD strain energy density

Assuming that rubber materials are both isotropic and incompressible, the proposed strain energy function $W$ only depends on the two first invariants of the left Cauchy-Green stretch tensor $B$:

$$ 
\begin{aligned}
W_{GD}(I_1, I_2) &= W_1(I_1) + W_2(I_2) \\
&= h_1 \int e^{h_3(I_1-3)^2} dI_1 + 3h_2 \int \frac{1}{\sqrt{I_2}} dI_2
\end{aligned}
\quad (1) $$

where $h_1, h_2, h_3$ are the material parameters. In this strain energy, the $I_1$ part of $W$ describes the global response of the material and is equivalent to the Hart-Smith model (Hart-Smith 1966). The second term that involves $I_2$ improves the accuracy of the model for multiaxial loading conditions (a general power $I_2$-term has been introduced by Lambert-Diani & Rey, 1999). The true stress tensor is defined by the differentiation of the proposed strain energy with respect to $B$:
$$ \sigma = -pI + 2 \left( \frac{\partial W}{\partial I_1} + I_1 \frac{\partial W}{\partial I_2} \right) B - 2 \frac{\partial W}{\partial I_2} B^2 \quad (2) $$

$$ I_1 = \mathrm{Tr}(B), \quad I_2 = \frac{1}{2} \left[ \mathrm{Tr}(B)^2 - \mathrm{Tr}(B^2) \right] \quad (3) $$

## 2.2 Links between GD and eight-chains models

The present part is devoted to the comparison of the GD model to the eight-chains one (Arruda and Boyce 1993) in order to highlight the physical motivations of this proposed phenomenological model. As these constitutive equations are supposed to be qualitatively efficient for the entire range of strains, both small and large strain responses are compared. Comparison of these models is established using their polynomial expansion in terms of $I_1$ respectively (equations 4 and 5) for small strain. High order terms are then eliminated:

$$ W_1(I_1) = C_1 \sum_{i=0}^{\infty} \frac{C_3^i}{(2i+1)i!} (I_1 - 3)^{2i+1} \quad (4) $$

The first five terms of the development of the eight-chains model are:

$$ W_{8ch} = C_R \left[ \frac{1}{2}(I_1 - 3) + \frac{1}{20N}(I_1^2 - 3) \right] + C_R \left[ \frac{1}{1050N^2}(I_1^3 - 27) + \frac{19}{7000N^3}(I_1^4 - 81) \right] + C_R \left[ \frac{519}{673750N^4}(I_1^5 - 243) + \dots \right] \quad (5) $$

The small strain stiffness of the models are respectively defined by first terms of Equations (4) and (5). Then, both the strain energy functions reduce to the neo-Hookean expression. The both models are characterized by their ability to describe the strain-hardening of the material that takes place under large strains. This strain-hardening phenomenon is mainly due to the extensibility limit of polymer chains. Chagnon et al. (2004) established that the first part of the Hart-Smith strain energy $W_1(I_1)$, which is identical to the first part of the GD model, is equivalent to the eight-chains and Gent models for the entire range of strains.

$$ \frac{\partial W_1(I_1)}{\partial I_1} = h_1 e^{h_3(I_1-3)^2} \quad (6) $$

We proposed now to highlight the physical motivation of the second part of strain energy $W_2(I_2)$ of the GD model. As mentioned by Treloar (Treloar 1975), term as function of $I_2$ in constitutive equations can be seen as corrections of the phantom network theory (terms as function $I_1$). We propose to constrain the eight-chains model by a new network of chains

on the surface of the cube (figure 1). The confinement of the eight-chains model is governed by a strain energy potential. This potential constrains the eight-chains cube surface. Let us recall that the surface of the eight-chains model (a cube) is $I_2^{1/2}$ and that its increase under deformation is $I_2^{1/2} - 3^{1/2}$. We therefore define a pressure constrain of the eight-chains rubber network. This phenomenon is modeled by the second invariant energy part:

$$ W_2(I_2) = 3h_2 \int \frac{1}{\sqrt{I_2}} dI_2 \quad \text{or} \quad \frac{\partial W_2(I_2)}{\partial I_2} = \frac{3h_2}{\sqrt{I_2}} \quad (7) $$

where $3h_2$ stands for the pressure constrain of entanglement of the eight-chains rubber network.

<figure><img src="image_2.png" /><figcaption>Figure 1. Entanglement of the eight-chains rubber network. Eight-chains model completed by surrounding chains (in bold) is the physical motivations of the GD model.</figcaption></figure>

## 2.3 Compressible GD strain energy density

In this part we introduce the proposed strain energy function for compressible isotropic hyperelastic materials in terms of strain invariants. Based on kinematic assumption, we use a decoupled representation of the strain energy function originally proposed by (Flory 1961, Ogden 1984).

$$ W(\bar{I}_1, \bar{I}_2, J) = W_{iso}(\bar{I}_1, \bar{I}_2) + W_{vol}(J) \quad (8) $$

$$ J = \det(F), \quad \bar{F} = J^{-\frac{1}{3}}F, \quad \bar{B} = \bar{F}\bar{F}^T \quad (9) $$

$$ \bar{I}_1 = \mathrm{Tr}(\bar{B}), \quad \bar{I}_2 = \frac{1}{2} \left[ \mathrm{Tr}(\bar{B})^2 - \mathrm{Tr}(\bar{B}^2) \right] \quad (10) $$

Where $W_{iso}(I_1, I_2)$ and $W_{vol}(J)$ stand for the isochoric and volumetric elastic strain energy of the material, respectively. The isochoric part is equivalent to the equation (1) with modified invariants. The volumetric part is chosen to enforce a nearly quasi incompressible behaviour.

$$ W_{iso}(\bar{I}_1, \bar{I}_2) = h_1 \int e^{h_3(\bar{I}_1-3)^2} d\bar{I}_1 + h_2 \int \frac{1}{\sqrt{\bar{I}_2}} d\bar{I}_2 \quad (11) $$

$$ W_{vol} = \frac{1}{D}(J-1)^2 \quad (12) $$

The volumetric properties part is presented in (Doll and Schweizerhof 2000). The true stress tensor
is defined by the differentiation of the proposed strain energy with respect to B:

$$ \sigma_{iso} = 2J^{-\frac{5}{3}} \left( \frac{\partial W_{iso}}{\partial I_1} + \bar{I}_1 \frac{\partial W_{iso}}{\partial I_2} \right) B - 2J^{-\frac{7}{3}} \frac{\partial W_{iso}}{\partial I_2} B^2 \quad (13) $$

$$ \sigma_{vol} = \frac{\partial W_{vol}}{\partial J} I, \quad \sigma = \sigma_{vol} + \sigma_{iso} \quad (14) $$

### 3 EXPERIMENTAL DATA AND IDENTIFICATION OF PARAMETERS

In order to compare the efficiency of the models, we choose two complementary data sets issued from classical references (Marckmann & Verron 2006). The first set is due to Treloar (1944). In the current study, data from Treloar for unfilled natural rubber (cross-linked with 8 parts of S phr) were used. This material exhibits highly reversible elastic response and no stretch-induced crystallization up to 400%. Thus it is well-modeled by hyperelastic constitutive equations. Experimental measures were performed for four different loading conditions: equibiaxial extension of a sheet (EQE), uniaxial tensile extension (UE), pure shear (PS) and biaxial extension (BE). The second data set is due to Kawabata et al. (1981). It was obtained using an experimental apparatus for general biaxial extension testing. In terms of stretch ratios, unfilled polyisoprene specimens were stretched from 1.04 to 3.7 in the first direction ($\lambda_1$) and from 0.52 to 3.1 in the perpendicular direction ($\lambda_2$). These values correspond to moderate strain but lead to deformation conditions from uniaxial extension to equibiaxial extension. Here, both experimental data sets are simultaneously considered to compare the models because the two materials are quite similar. Thus, for a given model, a unique set of material parameters must be able to reproduce these data with a good agreement. The parameter identification is performed using genetic algorithms as presented in (Marckmann and Verron 2006). The GD model includes entanglement of the eight-chains rubber network and identified responses are depicted on figures 2-4. The GD model responses on figures 2 to 4 are almost equivalent to Ogden six parameters results for bi-axial extensions presented on figures 5 to 7. Harts-Smith identified responses are depicted on figures 8-10.

<figure><img src="image_3.png" /><figcaption>Figure 2. Experimental data (- -) and GD model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible GD model response.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Experimental data (- -) and GD model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible GD model response.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Experimental data (●) and GD model identification for biaxial extension (BE), equibiaxial extension (EQE), pure shear (PS), uniaxial tensile extension (UE). Incompressible GD model response.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Experimental data (- -) and Ogden model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible Ogden six model response.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. Experimental data (- -) and Ogden model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible Ogden six model response.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Experimental data (•) and Ogden model identification for biaxial extension (BE), equibiaxial extension (EQE), pure shear (PS), uniaxial tensile extension (UE). Incompressible Ogden six model response.</figcaption></figure>

The GD parameters are: $h_1=0.142$ MPa, $h_2=1.585 \ 10^{-2}$ MPa and $h_3=3.495 \ 10^{-4}$. The proposed model GD is able to accurately reproduce the whole "S" shaped response of the material. The model behaves satisfactory under all the presented loadings (figures 2-4).

The Ogden six parameter strain energy density is classically

$$W = \sum_{n=1}^{3} \frac{\mu_n}{\alpha_n} (\lambda_1^{\alpha_n} + \lambda_2^{\alpha_n} + \lambda_3^{\alpha_n} - 3) \quad (15)$$

Parameters are: $\mu_1=0.63$ MPa, $\mu_2=1.2 \ 10^{-3}$ MPa, $\mu_3=-1 \ 10^{-2}$ MPa, $\alpha_1=1.3$, $\alpha_2=5$ and $\alpha_3=-2$. The conditions $\mu_i\alpha_i>0$ ensures the positive definite character of the strain energy. In order to recall the influence of $W_2(I_2)$, let us recall the classical strain energy of Hart-Smith (equation 16).

$$ 
\begin{aligned}
W_{HS}(I_1, I_2) &= W_1(I_1) + W_2(I_2) \\
&= a_1 \int e^{a_3(I_1-3)^2} dI_1 + a_4 a_3 \ln(I_2/3)
\end{aligned} 
\quad (16) 
$$

The Hart-Smith parameters are: $a_1=0.140$ MPa, $a_2=5.25410^{-4}$ and $a_3=1.290$.

<figure><img src="image_9.png" /><figcaption>Figure 8. Experimental data (- -) and Hart-Smith model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible Hart-Smith model response.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Experimental data (- -) and Hart-Smith model identification for biaxial tensile tests. Piola-Kirchhoff function of extension for several transverse extensions. Incompressible Hart-Smith model response.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. Experimental data (•) and Hart-Smith model identification for biaxial extension (BE), equibiaxial extension (EQE), pure shear (PS), uniaxial tensile extension (UE). Incompressible Hart-Smith model response.</figcaption></figure>

## 4 MULLINS EFFECT

Elastomers present a loss of stiffness after the first loading cycle of a fatigue experiment (Mullins, 1969). It has been observed that this phenomenon is only dependent on the maximum deformation previously reached in the history of the material. It is quite important to model it because the mechanical behaviour of rubber products is highly modified by this softening phenomenon. Moreover, as the Mullins effect depends on the maximum deformation endured previously, material REV of the product are not identically affected. As a consequence, it is not acceptable to determine experimentally an accommodated hyperelastic constitutive equation for the material. The stress-softening should be explicitly included in the model. The continuum damage me-
chanics has often been used to model the Mullins effect even if phenomenon undergoing Mullins effect is not a strictly speaking an irreversible damage phenomenon and can even be modelled without damage (Cantournet et al, 2009). For example, it can be recovered with time and annealing accelerates this recovery. A thermodynamic variable *d* is introduced to represent the loss of stiffness and the corresponding stress-softening. The general theory of Continuum Damage Mechanics is detailed in Lemaitre and Chaboche (1990) book. When applied to GD hyperelastic density, it yields to GDM model. This new strain energy function for hyperelastic model with damage variables can be written (eq 17) :

$$W_{\text{GDM}}(I_1, I_2) = \tilde{h}_1 \int e^{\tilde{h}_3(I_1-3)^2} dI_1 + 3\tilde{h}_2 \int \frac{1}{\sqrt{I_2}} dI_2 \quad (17)$$

$$\tilde{h}_1 = h_1(1-d_1), \tilde{h}_2 = h_2(1-d_2), \tilde{h}_3 = h_3(1-d_3)$$

The incompressible state laws (eq. 2) associated with this model are classically obtained by equation 18.

$$\sigma = \left. \frac{\partial W_{\text{GDM}}}{\partial B} \right|_{I_i=1} \qquad (18)$$

We define the thermodynamic forces $Y_{d_i}$ associated with damage internals variables $d_i$ by:

$$-Y_{d_i} = \frac{\partial W_{\text{GDM}}}{\partial d_i} \qquad (19)$$

A non standard damage model is build here as the damage thermodynamics forces are not used to describe damage evolution. We consider instead that the quantity governing the damage evolution laws is the maximum of the first invariant. This assumption is achieved according to the physical motivations of maximum strain state endured during the history of the deformation (Marckmann et al. 2002). In order to do this we can introduce damage criterion functions

$$f_i = I_1 - k_i(d_i) \qquad (20)$$

such as $f_i<0$ implies no damage evolution and as damage evolves at $f_i=0$ with then $d_i - k_i^{-1}(I_1)$. Choosing particular expressions for $k_i$-functions allows to derive the damage evolution laws (equations 19-20) and make the model complete. The proposed GDM model is based on an improved method already proposed in (Chagnon et al. 2004). The evolution equation of the damage variable is expressed thanks to the first strain invariant and presents an exponential form (equations 21). This model is able to represented unloadings Mullins effect responses (figures 11 and 12).

$$d_1 = d_{1\infty} \left( 1 - \exp\left( - \frac{I_1^{\text{max}}}{\eta_1} \right) \right) \qquad (21)$$

$$d_2 = d_{2\infty} \left( 1 - \exp\left( - \frac{I_2^{\text{max}}}{\eta_2} \right) \right)$$

where $d_{1\infty}$, $\eta_1$, $d_{2\infty}$, $\eta_2$ and $b$ are material parameters. $I_1^{\text{max}}$ represents the maximum value of the first strain invariant obtained during loadings. The coupling with damages $d_1$ and $d_2$ is similar to the discontinuous damage part of the constitutive equations proposed by Miehe (1995). Considering the mass conservation of the polymer network, which implies that the number of monomer segments per unit volume $N$ must remain constant (Marckmann et al. 2002), where $N$ is the mean number of monomer per chain and $n$ the number of chain per unit of volume. According to relationship between parameters derived in (Chagnon et al. 2004), this leads to link the mechanical properties $\tilde{h}_1$ and $\tilde{h}_3$ and shows that $d_3$ is not an independent thermodynamics damage variable,

$$d_3 = 1 - F(d_1), \quad F(d_1) = \frac{1}{\left(3(b \tilde{h}_1 - 1)\right)^2} \qquad (22)$$

The GDM is identified on idealized experimental data as presented in (Chagnon et al. 2004). Behaviour is considered time-independent. Figure 11 represents the response of the model for a tensile cyclic test. Pure Shear response is presented by Figure 12. These figures shows how the GDM model is able to accurately reproduce the response of successive loadings for different maximal stretches which characterizes the Mullins effect.

The GDM parameters are: $h_1=0.0157$ MPa, $h_2=0.0098$ MPa, $h_3=0.000561$, $d_{1\infty}=0.5525$, $\eta_1^{-1}=0.0119$, $d_{2\infty}=1.0 \eta_2^{-1}=0.3645$ and $b=62.69$
<figure><img src="image_12.png" /><figcaption>Figure 11. Tensile behaviour of a hyperelastic incompressible GDM model with stress-softening. Identification is performed for a material used on a roll restrictor.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12. Pure Shear behaviour of a hyperelastic incompressible GDM model with stress-softening. Identification is performed for a material used on a roll restrictor.</figcaption></figure>

## 5 FINITE ELEMENT PREDICTIONS

In the finite element context, the proposed GD and GDM models can easily be implemented because of their strain invariants formulations. This has been done in computer codes Cast3M-CEA (French Atomic Agency) and ABAQUS. The two previous constitutive models were implemented in the finite element context, thanks to the UMAT and UHYPER facilities.

<figure><img src="image_14.png" /><figcaption>Figure 13. Engine roll restrictor made of steel and rubber parts.</figcaption></figure>

Finite element predictions of a roll restrictor (figure 13) developed by Trelleborg automotive have been performed with ABAQUS and the GDM model. This behavior law is able to describe the local loss of stiffness of the material and the non-homogeneity of the structure after the first “first gear full torque” engine acceleration. The material properties used for finite element simulations are not the ones used for real structures. Predictions are performed with the GDM model discussed in previous sections. The 3D model is composed of C3D8H elements. It corresponds to an incompressible finite element hybrid formulation. It appears that the stress softening level is heterogeneous in the roll restrictor rubber parts (figures 14-16). Damage is located on the surface for most part of the structure excepted in the both snubbers.

<figure><img src="image_15.png" /><figcaption>Figure 14. White parts represent the damage level $d_1$ (Mullins effect) pattern in the roll restrictor during the first “first gear full torque” engine's acceleration.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. White parts represent the damage level $d_2$ (Mullins effect) pattern in the roll restrictor during the first “first gear full torque” engine's acceleration.</figcaption></figure>

Figure 17 illustrates the above Mullins observations in respect of a GDM model with isotropic damage evolution law (equation 23). Mullins effect is located on the surface for most part of the structure excepted in the both snubbers.

$$d = d_1 = d_2, d_2 = 0 \\ d = d_\infty \left( 1 - \exp \left( - \frac{I_1^{\max}}{\eta} \right) \right) \quad (23)$$
<figure><img src="image_17.png" /><figcaption>Figure 16. White parts represent the damage level $d_3$ (Mullins effect) pattern in the roll restrictor during the first “first gear full torque” engine's acceleration.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17. Isotropic damage model ($d = d_1 = d_2, d_3 = 0$). White parts represent the damage area "d" (Mullins effect) during the first “first gear full torque” engine's acceleration.</figcaption></figure>

## 6 CONCLUSIONS

We propose here a simple isotropic hyperelastic model (GD) expressed in terms of classical independent strain invariants of the symmetric Cauchy-Green tensor. The strain-energy part as function of $I_1$ is taken identical to the Hart-Smith one. This part is equivalent to the Eight chain model. Concerning the function of the second invariant $I_2$ a squareroot part is proposed. The corresponding energy density contribution is connected to the non-affine deformation of the entanglement Eight-chains network. The proposed model is successfully identified on both Treloar and Kawabata experiments. Only 3 parameters are needed to describe the experimental results. The response quality is equivalent to the one of the Ogden six parameter model.

The Mullins effect is finally taken into account by coupling with damage the GD model. The GDM model is thus derived stating that the loss of stiffness depends on the maximum value of the first invariant $I_1$. In the finite element context the proposed models GD and GDM are implemented in the finite element codes Cast3M and ABAQUS. Damage predictions observed on the FE simulations of the Engine roll are in good agreement with experimental data.

## 7 REFERENCES

* Arruda E. & Boyce M.C. 1993. A three dimensional constitutive model for the large stretch behavior of rubber elastic materials. *J. Mech. Phys. Solids*, 41, 2, 389-412

* Chagnon G. Marckmann G. and Verron E. 2004. A comparison of the hart-smith model with arruda-boyce and gent formulations for rubber elasticity, *Rubber Chemistry and Technology*, 77, 724-735

* Chagnon, G. Verron, E. Gornet, L. Marckmann, G. & Charrier, P. 2004. On the relevance of continuum damage mechanics as applied to the Mullins effect in elastomers. *J. Mech. Phys. Solids*, 52, 1627-1650.

* Cantournet S., Desmorat R., Besson J., Müllins effect and cyclic stress softening of filled elastomers by internal sliding and friction thermodynamics model, *International Journal of Solids and Structures*, 46, pp. 2255-2264, 2009.

* Doll S. and Schweizerhof K. 2000. On the development of volumetric strain energy function, *J. Appl Mech T ASME* 67, 17-21

* Flory, P.J. 1961. Thermodynamic relations for highly elastic materials, Transactions of the Faraday Society 57 (6,7), 829-838.

* Hart-Smith, L. J. 1966. Elasticity parameters for finite deformations of rubber-like materials. *Z. angew. Math. Phys.* 17, 608-626

* Kawabata, S., Matsuda, M., Tei, K., and Kawai, H. 1981, *Macromolecules*, 154-162.

* Lambert-Diani J., Rey C., New phenomenological behavior laws for rubbers and thermoplastic elastomers, *Eur. J. Mech. A/Solids*, 18, 1027-1043, 1999.

* Lemaitre, J. & Chaboche, J. L. 1990. *Mechanics of solid materials*, Cambridge University Press.

* Marckmann G., E. Verron 2006. Comparison of hyperelastic models for rubberlike materials, *Rubber Chemistry and Technology*. 5, 835-858

* Marckmann, E. Verron, L. Gornet, G. Chagnon, P. Charrier and P. Fort, A theory of network alteration for the Mullins effect, *J. Mech. Phys. Solids* **50** (2002), pp. 2011–2028

* Miehe, C. 1995. Discontinuous and continuous damage evolution in Ogden type large strain elastic materials. *Eur. J. Mech., A/Solids* 14, 697-720.

* Mullins, L. 1969. Softening of rubber by deformation. *Rubber Chem. Technol.* 42, 339-362.

* Ogden, R. W., 1984. Recent advances in the phenomenological theory of rubber elasticity. Rubber Chemistry and Technology. 59, 361-383.

* Ogden, R. W., 1972. Large deformation isotropic elasticity - on the correlation of theory and experiment for incompressible rubberlike solids. Proc. R. Soc. Lon. A. 326, 565-584.

* Treloar, L. R. G. 1975, The Physics of Rubber Elasticity, Oxford Classic Texts.

* Treloar, L. R. G. 1944, *Trans. Faraday Soc.*, 59-70.

* Cast3M, www-cast3m.cea.fr

* ABAQUS, www.simulia.com