# Development of new constitutive equations for the Mullins effect in rubber using the network alteration theory

Grégory Chagnon, Erwan Verron, Gilles Marckmann, Laurent Gornet

► To cite this version:

Grégory Chagnon, Erwan Verron, Gilles Marckmann, Laurent Gornet. Development of new constitutive equations for the Mullins effect in rubber using the network alteration theory. International Journal of Solids and Structures, 2006, 43 (22-23), pp.6817-6831. 10.1016/j.ijsolstr.2006.02.011. hal-01006742

HAL Id: hal-01006742

https://hal.science/hal-01006742

Submitted on 24 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Development of new constitutive equations for the Mullins effect in rubber using the network alteration theory

G. Chagnon<sup>a,*</sup>, E. Verron<sup>b</sup>, G. Marckmann<sup>b</sup>, L. Gornet<sup>b</sup>

<sup>a</sup> Laboratoire Sols Solides Structures, UMR CNRS 5521, Université Joseph Fourier, Institut National Polytechnique de Grenoble, BP 53, 38041 Grenoble Cedex 9, France

<sup>b</sup> Institut de Recherche en Génie Civil et Mécanique, UMR CNRS 6183, École Centrale de Nantes, BP 92101, 44321 Nantes Cedex 3, France

## Abstract

During cyclic loading, both natural and synthetic elastomers exhibit a stress-softening phenomenon known as the Mullins effect. In the last few years, numerous constitutive equations have been proposed. The major difficulty lies in the development of models which are both physically motivated and sufficiently mathematically well defined to be used in finite element applications. An attempt to reconcile both physical and phenomenological approaches is proposed in this paper. The network alteration theory of Marckmann et al. [Marckmann, G., Verron, E., Gornet, L., Chagnon, G., Charrier, P., Fort, P., 2002. A theory of network alteration for the Mullins effect. J. Mech. Phys. Solids 50, 2011–2028] is considered and modified. The equivalence between three different strain energy functions is then used to develop two new constitutive equations. They are founded on phenomenological strain energy densities which ensure simple numerical use, but the evolution of their material parameters during stress-softening is based on physical considerations. Basic examples illustrate the efficiency of this approach.

**Keywords:** Rubber material; Stress-softening; Finite strain; Constitutive equation

## 1. Introduction

Both natural rubber and synthetic elastomers are widely used in industrial design. In order to reduce the time necessary to develop new parts, numerical simulations are nowadays revealed of fundamental importance. In this context, one of the major difficulties encountered by engineers consists in the choice of a well-adapted constitutive model which satisfactorily reproduces the large strain inelastic response of rubbers. Indeed, they exhibit a time-dependent behaviour (relaxation, creep, hysteresis and Payne effect) and a particular stress-softening

<footnote>* Corresponding author. Tel.: +33 4 76 82 70 85; fax: +33 4 76 82 70 43.</footnote>

<footnote>E-mail address: gregory.chagnon@hmg.inpg.fr (G. Chagnon).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Stress-stretch response of a pseudo-elastic material. The sample is first stretched to $\lambda_1$ and the stress follows path I. Then the unloading path from $\lambda_1$ to 0 follows path I'. The second loading path from 0 to $\lambda_{II'} > \lambda_1$ first follows path I' until $\lambda = \lambda_1$ then it follows path II. The second unloading path from stretch ratio $\lambda_{II'}$ to 0 follows path II' which is different than path I'. It is to note that for a given stretch, stress on II' is lower than stress on I'. Repeating this process, the loading path from 0 to $\lambda_{III'}$ is the path that joins paths II' and III. Finally, the corresponding unloading response follows path III'.</figcaption></figure>

phenomenon which takes place during the five first cycles of a fatigue experiment. This phenomenon is known as the Mullins effect (Mullins, 1969).

The present paper only focuses on the Mullins effect; all other phenomena are not considered. Thus, it is admitted that the response of rubber is pseudo-elastic as schematically described in Fig. 1. For 60 years, numerous models have been proposed to predict stress-softening of rubbers. In the 1950s and 1960s, first works were only qualitative and mono-dimensional models were developed; three-dimensional constitutive equations have been proposed some years later. In the following, only three-dimensional models are reviewed. As there is no unanimous explanation of the physical origins of the Mullins effect, different approaches have been considered by authors.

First, some authors follows the phenomenological two-network theory of Mullins and Tobin (1957) who consider that the material is composed of a soft and a hard phase and that a part of hard phase is broken and transformed into soft phase during deformation. This approach was adopted by Wineman and Huntley (1994) and Huntley et al. (1996, 1997) who consider that the stress should be corrected by a scalar reformation function which depends on a measure of the deformation classically defined as a scalar function of principal stretches. Similarly, Beatty and Krishnaswamy (2000) and Zuñiga and Beatty (2002) proposed a measure in terms of strain invariants. These phenomenological models successfully describe the stress-softening phenomenon and their formulation is well adapted to finite element simulations. Nevertheless, these constitutive equations are not well adapted to material which admit large stress-hardening. Recently, using a similar approach, Qi and Boyce (2004) took into account the quantity of fillers considering that stress-softening only occurs in rubber bulk. Finally, Dorfmann and Ogden (2003) and Kazakevičiute-Makovska and Kačianauskas (2004) proposed an original theory in which the stress-softening function evolves during the unloading path of cycles instead of during loading path as in previous models.

Second, the general theory of Continuum Damage Mechanics has been applied to the Mullins effect since the end of the 1980s. Simo (1987) introduced the thermodynamical framework and proposed a three dimensional constitutive equation for the stress-softening. Later, many authors proposed phenomenological equations to predict the stress-softening phenomenon with this theory (DeSouza Neto et al., 1994; Miehe, 1995; Lion, 1996; Bikard and Desoyer, 2001). It has been recently demonstrated that the Continuum Damage Mechanics and the two-network approaches are similar (Chagnon et al., 2004b). It should be mentioned that all these phenomenological models are often easily implementable in finite element applications but cannot describe very large stress-softening response of elastomers because of the coupling between the strain energy density and the evolution laws adopted for damage variables.

Third, some physical-based constitutive equations were developed considering the evolution of the polymer network under loading as proposed by Bueche (1960, 1961). Govindjee and Simo (1991) proposed a three dimensional formulation that mimics the evolution of the chain network. They consider that chains break
progressively with the maximum deformation imposed. More recently, Drozdov and Dorfmann (2001) considered that the macromolecular chains are composed of two kinds of bonds: the flexed and the extended ones. When the material is deformed, the macromolecular chains extends and the quantity of bonds evolves. This permits to represent stress-softening and permanent set but the analytical form of the model is not easily usable for numerical problems. Marckmann et al. (2002) consider the evolution of the cross-links in the polymer network during deformation. Authors studied the evolution of the network without distinguishing the different natures of links breakage (between chains, between fillers and chains, or internal chain fracture). Thus, this leads to a constitutive equation well adapted to both filled and unfilled materials. More recently, Horgan et al. (2004) used a similar approach based on chain extensibility limit to derive a constitutive model based on the Gent (1996) strain energy density. The model is revealed hyperelastic during the first loading and pseudo-elastic during the subsequent loading.

The aim of the present paper is to develop new constitutive models for the Mullins effect which are both physically-motivated and well adapted for numerical problems. In this way, the network alteration theory of Marckmann et al. (2002) is retained for its physical meaning and the strain energy densities of Gent (1996) and Hart-Smith (1966) are considered for their mathematical simplicity and their ability to predict large strain response of elastomers. In the next section, the Marckmann et al. (2002) model is first briefly described. Then, a minor modification based on physical considerations of the network alteration theory is proposed. Implementation difficulties of this approach are exhibited, different solutions to simplify it are analyzed and are revealed not efficient. Section 3 is devoted to the derivation of two new constitutive equations. They are based on the equivalence between different strain energy densities: the 8-chain (Arruda and Boyce, 1993), the Gent (1996) and the Hart-Smith (1966) models. Finally, the implementation of these models in the finite element software Abaqus is presented in Section 4. Some examples illustrates their ability to predict the Mullins effect in elastomers.

## 2. The network alteration theory for the Mullins effect

### 2.1. Basic formulation of the Marckmann et al. model

Marckmann et al. (2002) studied the evolution of the rubber network due to stress-softening. Rubber materials are composed of long macromolecular chains connected by strong and weak links; the number of these links is increased by the presence of fillers. When the material is stretched for the first time, some links are broken. Thus, the Mullins effect can be explained by the evolution of the macromolecular network under loading. First, the mean distance between cross-links, i.e. the mean length of chains involved in elasticity, increases. Second, the number of chains per unit of volume decreases. Then, considering that the driving parameter of the Mullins effect is the maximum deformation previously endured by the material, the number of chains per unit of volume and the mean number of monomers per chain are respectively a decreasing and an increasing function of the maximum previous deformation. Moreover, authors made an additional assumption: they consider that the number of active monomers per unit of volume, i.e. monomers of chains that lie between two cross-links and are involved in elasticity, remains unchanged under loading.

These physical considerations are introduced in the 8-chain model because the corresponding strain energy density only depends on two material parameters which are the mean number of monomers per chain $N$ and a scalar proportional to the chain density $C_r$:

$$W = C_r N \left[ \frac{\lambda_{\text{chain}}}{\sqrt{N}} \beta \ln \left( \frac{\beta}{\sinh \beta} \right) \right], \quad (1)$$

where

$$\beta = \mathcal{L}^{-1} \left( \frac{\lambda_{\text{chain}}}{\sqrt{N}} \right) \quad \text{and} \quad \lambda_{\text{chain}} = \sqrt{\lambda_1^2 + \lambda_2^2 + \lambda_3^2}. \quad (2)$$

In this equation, $\mathcal{L}$ is the Langevin function defined by $\mathcal{L}(\beta) = \coth(\beta) - 1/\beta$ and $(\lambda_i)_{i=1,3}$ are the principal stretch ratios. In the approach of Marckmann et al. (2002), the material parameters depend on the maximum deformation through exponential functions:
$$C_r = C_{r0} \exp(-C_{r1} \bar{x}) \quad \text{and} \quad N = N_0 \exp(N_1 \bar{x}), \qquad (3)$$

where $C_{r0}$, $C_{r1}$, $N_0$ and $N_1$ are the material parameters, and $\bar{x}$ is the maximum of a deformation measure $\alpha$:

$$\bar{x} = \bar{x}(t) = \max_{\tau \in [0, t]} \alpha(\tau). \qquad (4)$$

As explained in Chagnon et al. (2004b), the first strain invariant $I_1$, i.e. the trace of the Cauchy–Green strain, can be chosen because of its ability to simply represent three-dimensional deformation conditions. The additional assumption, i.e. the number of active monomer remains constant under loading, leads to

$$C_{r1} = N_1. \qquad (5)$$

Thus the model admits only three material parameters.

## 2.2. A minor modification of the model

The strong assumption which states that the number of active monomers remains unchanged under loading is discussed in this paragraph and a more relevant assumption is proposed. During the first loading, some links are broken in the bulk material. The nature of these broken links is not well established; nevertheless, one can think that many weak (van der Waals) and a few strong links (covalent) are broken: the chains can break or can be detached from fillers. Thus, new dangling chains take place. In the classical theory of rubber elasticity, only chains between cross-links are involved in the elastic response of the macromolecular network, dangling chains are not active in the elastic response of the material. Consequently, the number of active monomers does not remain constant during loading: it decreases with the maximum deformation endured by the material (Chagnon, 2003). The condition (5) which relies the material parameters in the approach of Marckmann et al. (2002) is no longer relevant. The number of active monomers in the network is defined as the product of the chain density by the mean number of monomers per chain and it is proportional to

$$C_r N(\bar{x}) = C_{r0} N_0 \exp[(N_1 - C_{r1})\bar{x}], \qquad (6)$$

and this function decreases if and only if $N_1 - C_{r1} \le 0$.

Finally, the new model, called M2-model through the rest of the paper, admits four material parameters and two of them should fulfill the previous inequality. The parameters are fitted using experimental data of Chagnon et al. (2004b): $C_{r0} = 0.0575$, $C_{r1} = 0.0282$, $N_0 = 4.85$ and $N_1 = 0.0278$. The comparison between experimental data and model predictions are presented in Fig. 2 for both uniaxial extension and pure shear data. Results obtained with this model are similar to those previously presented in Marckmann et al. (2002). Nevertheless, the physical assumptions adopted for the present model are more realistic than those proposed in the above-mentioned paper. It is to note that the mechanical response presented in Fig. 2 admits

<figure><img src="image_3.png" /><figcaption>Fig. 2. Comparison of the predictions of the M2-model (—) with experimental data (· · · ○ · · ·): (a) uniaxial extension, (b) pure shear.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Evolution of the material parameters of M2-model with the maximum deformation: (a) chain density (proportional), (b) chain length (proportional) and (c) number of active monomers in the network (proportional).</figcaption></figure>

an upturn for stretch ratios of about 5.8 which corresponds to the maximum stretch ratio of the experimental database considered here. In order to improve the fit of the model, experimental data for stretch ratios greater than 6.0 should have been available. Nevertheless, this upturn is quite interesting because it corresponds to the maximum stretch above which the previous material parameters should not be used for numerical applications. In fact, it defines the deformation range in which Drucker convexity conditions are fulfilled.

The evolution of the material parameters as a function of the first strain invariant is shown in Fig. 3. Obviously, the number of active monomers in the network is a decreasing function of the deformation level (see Fig. 3(c)). It is noticeable that the slope of this curve is small: only few monomers are transferred from active to dangling chains under loading.

## 2.3. Limitations

In the framework of hyperelasticity, numerous forms of the strain energy density have been proposed for many years. Nevertheless, only few of them are actually used for finite element simulations. As an example, in the finite element software Abaqus©, the following models are implemented: Rivlin series (Rivlin and Saunders, 1951), the van der Waals model (Kilian, 1981), the Ogden model (Ogden, 1972) and the 8-chain model (Arruda and Boyce, 1993). More precisely, the later is not implemented under its original form but only its fifth order series expansion is considered:

$$W = C_r \sum_i C_i (I_i^i - 3^i), \qquad (7)$$

where $C_r$ is a material parameter proportional to the chain density and $C_i$ are the material parameters which only depend on the mean number of monomers per chain $N$. Their values are given in Table 1 up to the 10th order series expansion. In fact, the original form of this model is not easy to use in finite element applications because it is written in the principal strain directions and it explicitly contains the inverse Langevin function.

<table><caption>Table 1
Material parameters of the $I_1$-series expansion for the 8-chain model (Arruda and Boyce, 1993)</caption><tr><td>$C_1$</td><td>1/2</td></tr><tr><td>$C_2$</td><td>1/20N</td></tr><tr><td>$C_3$</td><td>11/1150N<sup>2</sup></td></tr><tr><td>$C_4$</td><td>19/7000N<sup>3</sup></td></tr><tr><td>$C_5$</td><td>519/673750N<sup>4</sup></td></tr><tr><td>$C_6$</td><td>59991/262762500N<sup>5</sup></td></tr><tr><td>$C_7$</td><td>105771/1532781250N<sup>6</sup></td></tr><tr><td>$C_8$</td><td>3123763/148898750000N<sup>7</sup></td></tr><tr><td>$C_9$</td><td>54543778207/8577405555468750N<sup>8</sup></td></tr><tr><td>$C_{10}$</td><td>74301767/38899798437500N<sup>9</sup></td></tr></table>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Predictions of the M2-model: original formulation (—), 10th order $I_1$-series expansion (---), fifth order $I_1$-series expansion (···). (a) Uniaxial extension, (b) pure shear.</figcaption></figure>

Even if the inverse Langevin function can be replaced by approximated functions (Bergstrom, 1999; Perrin, 2000), difficulties remains unchanged. The discrepancy between analytical and series expansion formulations is not significant for small and medium strain but it cannot be neglected for large strain.

For the original formulation, the curve admits an asymptote which corresponds with the extensibility limit of chains; in the two other cases, the function is polynomial. So, the difference is not very significant for classical elastic computations, but for the prediction of the Mullins effect it is quite different: responses of the Marckmann et al. and M2-models highly depend on the position of the strain-hardening part of secondary loading curves and the evolution $N$ is directly relied to the position of this asymptotic part of the curve. As an example, predictions for the M2-model are performed with the original formulation and the fifth and 10th series expansions for both uniaxial tensile and pure shear loading conditions. Here, the previous values of the material are considered and the corresponding results are presented in Fig. 4. This example reveals that $I_1$-series expansions are not relevant to correctly describe the stress-softening phenomenon. Fitting the model under its $I_1$-series forms does not lead to results as good as those obtained with the original formulation (with the inverse Langevin function): strain-hardening parts of the response are not correctly reproduced. Thus, the original formulation is the only one which leads to good correlation between phenomena and modelling, but it is not easily usable in finite element softwares. In fact, the use of the concept of network alteration for the Mullins effect necessitates alternative formulations; this is the aim of the following section.

### 3. Development of new constitutive equations

In order to avoid the difficulties due to the mathematical complexity of the original 8-chain model and the limitations of its series expansions, different strain energy densities should be considered and associated with the concept of network alteration. Nevertheless, this theory being based on the evolution of network parameters, i.e. chain density and length, relationships between the material parameters of these strain energy densities and the network parameters have to be investigated. Here, the Gent and the Hart-Smith models are examined. The former is considered because it is explicitly expressed in terms of the first strain invariant. The latter is chosen because its strain-hardening response is approached by an exponential function contrary to both previous constitutive equations which large strain responses are defined by asymptotic functions (corresponding to the limit of extensibility). The Gent strain energy density is defined by

$$W = -\frac{E}{6} J_m \ln \left[ 1 - \frac{I_1 - 3}{J_m} \right], \qquad (8)$$

where $E$ and $J_m$ are the two material parameters. The response of this model presents a similar shape as the one of the 8-chain model: it is defined by an initial stiffness, through $E$, and the position of an asymptote, through $J_m$, which is related to the extensibility limit of chains and reproduces the strain-hardening phenomenon at large
strain. A few years ago, Boyce (1996) demonstrated the equivalence between the 8-chain and the Gent models by establishing the following relationships between their material parameters:

$$E = 3C_r \quad \text{and} \quad J_m = 3(N-1). \qquad (9)$$

So, very simple equations relate the initial stiffness E to the chain density and the maximum deformation Jm to the mean number of monomers per chain.

A similar study was recently conducted by Chagnon et al. (2004a) to compare the I₁-part of the Hart-Smith and the 8-chain constitutive equations. The original Hart-Smith strain energy density is defined as follows:

$$W = C_1 \int \exp\{C_3(I_1 - 3)^2\} dI_1 + C_2 \ln\left(\frac{I_2}{3}\right), \qquad (10)$$

where (Cᵢ)ᵢ=1,3 are the three material parameters. The parameter C₂ is related to the second strain invariant I₂ and permits to well reproduce the stiffness change observed for moderate strain, similarly to the classical Moo-ney (1940) theory. Here, this parameter is set to 0 in order to define I₁-part of the Hart-Smith model which is revealed comparable with the 8-chain and Gent models. The only difference between the I₁-part of the Hart-Smith model and the two other models is that its large strain response is driven by an exponential function instead of an asymptote. Finally, this model involves only two parameters: the first one C₁ describes the initial stiffness and the second one C₃ the large strain response. Chagnon et al. (2004b) derived the following relationships between the parameters:

$$C_r = 2C_1 \quad \text{and} \quad C_3 J_m^2 = k, \qquad (11)$$

with k ∈ [2.83; 2.98]. So, similarly to the Gent model, the material parameters of the Hart-Smith model can be related to network characteristics: C₁ to the chain density and C₃ to the mean number of monomers per chain.

Thanks to the physical meaning of material parameters involved in the phenomenological strain energy densities of Gent and Hart-Smith, two new constitutive equations for the Mullins effect can be derived; they are studied in the following.

### 3.1. A network alteration model based on the Gent strain energy density

Thanks to the constitutive equation (3) of the M2-model and to the relationships between material parameters of the 8-chain and Gent models (9), the alteration equations of the Gent material parameters can be easily established:

$$E = 3C_{r0} \exp(-C_{r1}\alpha) \quad \text{and} \quad J_m = 3[N_0 \exp(N_1\alpha) - 1]. \qquad (12)$$

Considering the material parameters of the M2-model fitted previously, Fig. 5 presents the comparison between responses of the M2- and the altered Gent models. Even if the new constitutive equation is revealed

<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison of the altered Gent (—) and the M2- (---) models using the same values of the material parameter: (a) uniaxial extension, (b) pure shear.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Comparison of the 8-chain (---) and Gent (—) responses for uniaxial extension.</figcaption></figure>

slightly stiffer than the M2-model for large strain, response curves of both models are quite similar. The discrepancy between them is due to the fact that the equivalence relationships are only based on the initial stiffness and the extensibility limit. As an example, Fig. 6 shows both elastic 8-chain (with $C_r = 1$ and $N = 10$) and Gent (with $E = 3$ and $J_m = 27$) responses under uniaxial tensile loading conditions. Both curves are superimposed for small strain and close to for large strain, in the neighborhood of the asymptote. But for moderate strain they slightly differs: the Gent model is stiffer. This explains why stresses predicted with the altered Gent model are greater than those corresponding to the M2-model in Fig. 5.

To improve the model, it is possible to directly fit the material parameters of the altered Gent model using experimental data. Then, the evolution of parameters is considered exponential:

$$E = E_0 \exp(-E_1 \alpha) \quad \text{and} \quad J_m = J_{m0} \exp(J_{m1} \alpha). \qquad (13)$$

The following values of the parameters are obtained: $E_0 = 0.196$, $E_1 = 0.0384$, $J_{m0} = 12.3$ and $J_{m1} = 0.0304$. It is noticeable that they are different than those predicted with the equivalence relationships. The corresponding stress-softening model, defined by Eqs. (8) and (13), is called the M3-model through the rest of the paper.

Predictions of this model for both uniaxial extension and pure shear are in good agreement with experiments as shown in Fig. 7. Results are similar to those previously obtained with the M2-model (see Fig. 2). Finally, the evolution of network parameters, i.e. $E/3$, $J_m/3 + 1$ and $E/3 \times (J_m/3 + 1)$ which respectively correspond to the chain density, the chain length and the number of active monomers in the network with deformation are presented in Fig. 8. Even if the values of the parameters differ from those of the M2-model, their evolutions are similar. Moreover, it is shown that the number of active monomers is a decreasing function of the maximum deformation; in fact, $E_1 \ge J_{m1}$ is a sufficient but not a necessary condition to ensure that the number of active monomers decreases with deformation.

## 3.2. A network alteration model based on the Hart-Smith strain energy density

The method employed in this section is similar to the previous one which conducted to the M3-model. The major difference rests in the fact that the equivalence relationships between 8-chain and Hart-Smith models are
<figure><img src="image_8.png" /><figcaption>Fig. 7. Comparison of the M3-model (—) with experimental data (· · · ○ · · ·): (a) uniaxial extension, (b) pure shear.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Evolution of the material parameters for the M3-model: (a) chain density (proportional), (b) chain length (proportional) and (c) number of active monomers in the network (proportional).</figcaption></figure>

more complicated than those examined previously. Thanks to Eq. (3) and relationships between material parameters (Eqs. (9) and (11)), the evolution of the material parameters of the Hart-Smith model is given by

$$
C_1 = \frac{C_{r0}}{2} \exp(-C_{r1}\bar{x}) \quad \text{and} \quad C_3^{-1} = \frac{9}{k} [N_0 \exp(N_1\bar{x}) - 1]^2. \qquad (14)
$$

The corresponding mechanical (stress–strain) response of the model is identical to the one presented in Fig. 6.
Stress for large strain are overestimated due to the discrepancy observed between the 8-chain and Hart-Smith
models for moderate strain (Chagnon et al., 2004a). Thus, in order to improve the model, material parameters
should be directly determined by fitting experimental data as performed above with the Gent model. The
material functions adopted are then

$$
C_1 = C_{10} \exp(-C_{11} \bar{x}) \quad \text{and} \quad C_3 = C_{30} \exp(-C_{31} \bar{x}), \qquad (15)
$$

where C₁₀, C₁₁, C₃₀ and C₃₁ are the material parameters. The corresponding model, defined by Eq. (10) with
C₂ = 0 and (15), is referred as to the M4-model in the following. Values of the parameters obtained with exper-
imental data are: C₁₀ = 0.0376, C₁₁ = 0.0434, C₃₀ = 0.0135 and C₃₁ = 0.0515.

The comparison of this new constitutive equation with experiments for both uniaxial extension and pure
shear is shown in Fig. 9. Results are in good accordance with experiments and predictions are similar to those
of the M2- and M3-models (see Figs. 2 and 7). Finally, evolution of the network parameters, i.e. 2C₁,
√k/9C₃ + 1 and 2C₁(√k/9C₃ + 1) which respectively correspond to the chain density, the chain length
and the number of active monomers in the network, are shown in Fig. 10. It is verified that the number of
active monomers is a decreasing function of the maximum deformation: a sufficient but not necessary condi-
tion to ensure this decrease is C₁₁ ≥ C₃₁/2.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Comparison of the M4-model (—) with experimental data (···O···): (a) uniaxial extension, (b) pure shear.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Evolution of the material parameters for the M4-model: (a) chain density (proportional), (b) chain length (proportional) and (c) number of active monomers in the network (proportional).</figcaption></figure>

## 4. Numerical results

The previous models are implemented in the finite element software Abaqus. Here, equations necessary for implementation are briefly recalled, then two examples are described.

### 4.1. Implementation

The previous constitutive equations assume the incompressibility of the material. This assumption greatly simplifies analytical solutions for homogeneous deformation, but it leads to major numerical difficulties such that a quasi-incompressible approach should be adopted (Bonet and Wood, 1997). In this context, the classical separation between isochoric and purely spherical deformation is considered (Ogden, 1984) and the strain energy density has to be decomposed into isochoric and spherical parts:

$$W = \bar{W}(\bar{I}_1, \bar{I}_2) + U(I_3). \qquad (16)$$

In this equation, the isochoric strain energy $\bar{W}$ depends on the two first isochoric strain invariants ($\bar{I}_1 = I_3^{1/3}I_1$ and $\bar{I}_2 = I_3^{2/3}I_2$) and the spherical strain energy $U$ is only expressed in terms of the third strain invariant $I_3$.

Finite element implementation necessitates the determination of the elasticity tensor. Thus, the Lagrangian elasticity tensor $\mathcal{E}^L$ is derived thanks to the differentiation of the second Piola-Kirchhoff stress tensor $\mathbf{S}$ with respect to the right Cauchy-Green strain tensor $\mathbf{C}$:
$$ \mathcal{C}^{\text{L}} = \mathcal{C}^{\text{L,iso}} + \mathcal{C}^{\text{L,sph}} = 2 \frac{\partial \mathbf{S}^{\text{iso}}}{\partial \mathbf{C}} + 2 \frac{\partial \mathbf{S}^{\text{sph}}}{\partial \mathbf{C}}, \qquad (17) $$

where $\mathbf{S}^{\text{iso}}$ and $\mathbf{S}^{\text{sph}}$ stand for the isochoric and spherical stress tensors given by

$$ \mathbf{S} = \mathbf{S}^{\text{iso}} + \mathbf{S}^{\text{sph}} = \frac{2}{\det(\mathbf{C})^{\frac{1}{2}}} \frac{\partial \bar{W}}{\partial \bar{I}_1} \left( \mathbf{I} - \frac{1}{3} I_1 \mathbf{C}^{-1} \right) + 2 I_3 \frac{\partial U}{\partial I_3} \mathbf{C}^{-1}. \qquad (18) $$

The two elasticity tensors $\mathcal{C}^{\text{L,iso}}$ and $\mathcal{C}^{\text{L,sph}}$ of Eq. (17) are nowadays well established (Bonet and Wood, 1997; Holzapfel, 2000). The isochoric elasticity tensor depends on the choice of the isochoric strain energy. As the M3- and M4-models only depend on the first strain invariant, the isochoric elasticity tensor reduces to

$$ \mathcal{C}^{\text{L,iso}} = \frac{4}{3 \det(\mathbf{C})^{\frac{1}{2}}} \frac{\partial \bar{W}}{\partial \bar{I}_1} \left( I_1 \frac{\partial \mathbf{C}^{-1}}{\partial \mathbf{C}} - \mathbf{C}^{-1} \otimes \mathbf{I} - \mathbf{I} \otimes \mathbf{C}^{-1} + \frac{1}{3} I_1 \mathbf{C}^{-1} \otimes \mathbf{C}^{-1} \right) \\ + \frac{4}{3 \det(\mathbf{C})^{\frac{1}{2}}} \frac{\partial^2 \bar{W}}{\partial \bar{I}_1^2} \left( 3 \mathbf{I} \otimes \mathbf{I} - I_1 \mathbf{I} \otimes \mathbf{C}^{-1} - I_1 \mathbf{C}^{-1} \otimes \mathbf{I} + \frac{1}{3} I_1^2 \mathbf{C}^{-1} \otimes \mathbf{C}^{-1} \right). \qquad (19) $$

Thus, only the two first derivatives of the strain energy density with respect to $\bar{I}_1$ have to be defined to compute this elasticity tensor. For our new models, they are given by

• M3-model:

$$ \frac{\partial \bar{W}}{\partial \bar{I}_1} = \frac{E J_m}{6} \frac{1}{J_m - \bar{I}_1 + 3}, \qquad (20) $$

$$ \frac{\partial^2 \bar{W}}{\partial \bar{I}_1^2} = \frac{J_m}{6} \frac{\partial E}{\partial \bar{I}_1} \frac{1}{J_m - \bar{I}_1 - 3} + \left[ \frac{E J_m}{6} - (\bar{I}_1 - 3) \frac{\partial J_m}{\partial \bar{I}_1} \right] \frac{1}{(J_m - \bar{I}_1 - 3)^2} \qquad (21) $$

with

$$ \frac{\partial E}{\partial \bar{I}_1} = \begin{cases} -E_0 E_1 \exp(-E_1 \bar{\alpha}) & \text{if } \alpha = \bar{\alpha}, \\ 0 & \text{otherwise,} \end{cases} \qquad (22) $$

$$ \frac{\partial J_m}{\partial \bar{I}_1} = \begin{cases} J_{m0} J_{m1} \exp(J_{m1} \bar{\alpha}) & \text{if } \alpha = \bar{\alpha}, \\ 0 & \text{otherwise.} \end{cases} \qquad (23) $$

• M4-model:

$$ \frac{\partial \bar{W}}{\partial \bar{I}_1} = C_1 \exp[C_3(\bar{I}_1 - 3)^2], \qquad (24) $$

$$ \frac{\partial^2 \bar{W}}{\partial \bar{I}_1^2} = \left[ \frac{\partial C_1}{\partial \bar{I}_1} + C_1 \frac{\partial C_3}{\partial \bar{I}_1} (\bar{I}_1 - 3)^2 + 2C_1 C_3 (\bar{I}_1 - 3) \right] \exp [C_3 (\bar{I}_1 - 3)^2] \qquad (25) $$

with

$$ \frac{\partial C_1}{\partial \bar{I}_1} = \begin{cases} -C_{10} C_{11} \exp(-C_{11} \bar{\alpha}) & \text{if } \alpha = \bar{\alpha}, \\ 0 & \text{otherwise,} \end{cases} \qquad (26) $$

$$ \frac{\partial C_3}{\partial \bar{I}_1} = \begin{cases} -C_{30} C_{31} \exp(-C_{31} \bar{\alpha}) & \text{if } \alpha = \bar{\alpha}, \\ 0 & \text{otherwise.} \end{cases} \qquad (27) $$

Similarly, the spherical part of the elasticity tensor is calculated after having chosen a spherical strain energy density. Here, we adopt $U(I_3) = 1/2K(I_3 - 1)^2$, $K$ being the compressibility modulus. For the details of the corresponding elasticity tensor, the reader can refer to Bonet and Wood (1997).
4.2. Simple examples

First, simple experiments, i.e. uniaxial extension and pure shear tests are considered. Simulations are con-
ducted using the two new stress-softening constitutive equations, i.e. M3- and M4-models.

For uniaxial extension, a classical bone specimen is considered. Its length is equal to 20 mm, its thickness is
2 mm and the central section is 4 mm large. It is meshed with linear elements in displacement. One of its
extremity is fixed and the other one is subjected to a prescribed displacement cycle. In regards with the thick-
ness, the plane stress assumption is adopted. Pure shear tests are also simulated, only a quarter of the specimen
is meshed because of the symmetries. Linear displacement elements are used. The parameter values of the pre-
vious section are adopted for the two models.

First, the strain-stress response predicted by the simulation in the central section is identical to the analyt-
ical response for both models. As demonstrated with uniaxial extension computations, greater is the stress-
softening effect, smaller is the number of active monomers in the network. Fig. 11 shows the evolution of
*n*N which is proportional to the number of active monomers. Obviously, stress-softening is revealed homoge-
neous in the centre of the sample and smaller in the neighbourhood of sample extremities. Similarly, numerical
predictions for pure shear exactly reproduce analytical results for both models. In Fig. 12, the number of
active monomers after a pure shear cycle is presented in a quarter of the specimen. In comparison with uni-
axial predictions, the whole sample is softened. The stiffness decrease is quasi-uniform, except near the free
edge where the deformation is different than the one in the rest of the specimen.

These two results validate the implementation of both models and illustrate the inhomogeneity of stress-
softening even in simple structures subjected to simple loading conditions. Such computations are necessary
to exhibit areas of structures in which the level of stress-softening is noticeable and then in which the stiffness
will be highly reduced for further deformation.

<figure><img src="image_12.png" /><figcaption>Fig. 11. Evolution of the number of active monomers in the rubber network for the M3- and M4-models.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12. Number of active monomers in the rubber network after a pure shear cycle for the (a) M3- (b) M4-models.</figcaption></figure>

### 4.3. An example of 3D simulation

The efficiency of the M3- and M4-models is evaluated with a more complicated problem. A diabolo sample of 60 mm height and 30 mm diameter is considered. It is meshed with hybrid elements linear in displacement and constant in pressure. The parameter values of the M3- and M4-models are those of the previous section. First the sample is softened in uniaxial extension and compression for different levels of stretch. Afterwards, torsion of 5 rad is prescribed. The aim of this example is to investigate the influence of stress-softening on subsequent material response.

Fig. 13 presents the torque-angle responses of the sample for the M3-model (the results obtained with the M4-model are similar and are not presented here). In the case of compressive pre-loading, there is no major differences with the case without pre-loading: even if the structure is deformed, local strain are not very large and thus the stress-softening phenomenon is not activated. For uniaxial extension pre-loading, it is shown that the level of extension highly influences the subsequent response of the sample. Finally, for the lateral extension which renders the pre-loaded structure non-symmetric, the structure is less softened than for uniaxial tensile

<figure><img src="image_14.png" /><figcaption>Fig. 13. Torque-angle responses of a diabolo sample after different softening conditions: (—) no previous loading, (· · ·) 10 mm uniaxial compression, (---) 30 mm uniaxial extension, (–) 45 mm uniaxial extension, (––) 60 mm uniaxial extension, (—) 50 mm lateral extension.</figcaption></figure>
loading conditions: a part of the central area of the sample is not softened and then the subsequent response is
revealed stiffer.

**5. Conclusion**

In this paper, the theory of network alteration proposed recently by Marckmann et al. (2002) has been
modified in order to take into account the occurrence of dangling chains in the network. As a consequence,
the number of monomers involved in the elastic response of the material is a decreasing function of the max-
imum deformation. This model is demonstrated to be efficient to predict stress-softening but its mathematical
formulation makes it difficult to use for numerical simulation. Then, considering that invariant formulations
are often preferred, the equivalence between some classical hyperelastic strain energy densities was used to
develop two new stress-softening constitutive equations. Their numerical abilities were demonstrated with
some examples and the influence of stress-softening on the subsequent response of rubber parts was
highlighted.

**References**

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

Beatty, M.F., Krishnaswamy, S., 2000. A theory of stress-softening in incompressible isotropic materials. J. Mech. Phys. Solids 48, 1931-1965.

Bergstrom, J.S., 1999. Large strain time-dependant behavior of elastomeric materials. Ph.D. thesis, Massachusetts Institute of Technology.

Bikard, J., Desoyer, T., 2001. Finite viscoelasticity, plasticity and damage of a class of filled elastomers: constitutive model. Mech. Res. Commun. 28, 693-702.

Bonet, J., Wood, R.D., 1997. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press.

Boyce, E.M., 1996. Direct comparison of the gent and the Arruda-Boyce constitutive models of rubber elasticity. Rubber Chem. Technol. 69, 781-785.

Bueche, F., 1960. Molecular basis for the Mullins effect. J. Appl. Polym. Sci. 3, 107-114.

Bueche, F., 1961. Mullins effect and rubber filler interaction. J. Appl. Polym. Sci. 5, 271-281.

Chagnon, G., 2003. Modélisation de l'effet Mullins dans les élastomères. Ph.D. thesis, Ecole Centrale de Nantes.

Chagnon, G., Marckmann, G., Verron, E., 2004a. A comparison of the physical model of Arruda-Boyce with the empirical Hart-Smith model and the Gent model. Rubber Chem. Technol. 77, 724-735.

Chagnon, G., Verron, E., Gornet, L., Marckmann, G., Charrier, P., 2004b. On the relevance of continuum damage mechanics as applied to the Mullins effect: theory, experiments and numerical implementation. J. Mech. Phys. Solids 52, 1627-1650.

DeSouzaNeto, E., Peric, D., Owen, D., 1994. A phenomenological three dimensional rate independent continuum damage model for highly filled polymers: formulation and computational aspects. J. Mech. Phys. Solids 42, 1533-1550.

Dorfmann, A., Ogden, R., 2003. A pseudo-elastic model for loading, partial unloading and reloading of particle-reinforced rubbers. Int. J. Solids Struct. 40, 2699-2714.

Drozdov, A.D., Dorfmann, A.I., 2001. Stress-strain relations in finite viscoelastoplasticity of rigid-rod networks: applications to the Mullins effect. Continuum Mech. Thermodyn. 13, 183-205.

Gent, A., 1996. A new constitutive relation for rubber. Rubber Chem. Technol. 69, 59-61.

Govindjee, S., Simo, J., 1991. A micro-mechanically continuum damage model for carbon black filled rubbers incorporating Mullins's effect. J. Mech. Phys. Solids 39, 87-112.

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

Holzapfel, G.A., 2000. Nonlinear Solid Mechanics—A Continuum Approach for Engineering. John Wiley.

Horgan, C.O., Ogden, R.W., Saccomandi, G., 2004. A theory of stress softening of elastomers based on finite chain extensibility. Proc. R. Soc. London A 460, 1737-1754.

Huntley, H.E., Wineman, A.S., Rajagopal, K.R., 1996. Chemorheological relaxation, residual stress and permanent set arising in radial deformation of an elastomeric hollow sphere. Math. Mech. Solids 1, 267-299.

Huntley, H.E., Wineman, A.S., Rajagopal, K.R., 1997. Stress softening, strain localization and permanent set in the circumferential shear of an incompressible elastomeric cylinder. IMA J. Appl. Math. 59, 309-338.

Kazakevičiute-Makovska, R., Kačiauskas, R., 2004. Modelling of stress softening in elastomeric materials: foundations of simple theories. Mech. Res. Comm. 31, 395-403.

Kilian, H.G., 1981. Equation of state of real networks. Polymer 22, 209-217.

Lion, A., 1996. A constitutive model for carbon black filled rubber: experimental investigations and mathematical representation. Continuum Mech. Thermodyn. 8, 153-169.

Marckmann, G., Verron, E., Gornet, L., Chagnon, G., Charrier, P., Fort, P., 2002. A theory of network alteration for the Mullins effect. J. Mech. Phys. Solids 50, 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.

Mooney, M., 1940. A theory of large elastic deformation. J. Appl. Phys. 11, 582-592.

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

Mullins, L., Tobin, N., 1957. Theoretical model for the elastic behavior of filler-reinforced vulcanized rubbers. Rubber Chem. Technol. 30, 551-571.

Ogden, R., 1972. Large deformation isotropic elasticity—on the correlation of theory and experiment for incompressible rubber like solids. Proc. R. Soc. Lond. A 326, 565-584.

Ogden, R., 1984. Recent advances in the phenomenological theory of rubber elasticity. Rubber Chem. Technol. 59, 361-383.

Perrin, G., 2000. Analytic stress-strain relationship for isotropic network model of rubber elasticity. C.R. Acad. Sci. Paris t328, 5-10.

Qi, H., Boyce, M., 2004. Constitutive model for stretch-induced softening of the stress-stretch behavior of elastomeric materials. J. Mech. Phys. Solids 52, 2187-2205.

Rivlin, R., Saunders, D., 1951. Large elastic deformations of isotropic materials—VII. Experiments on the deformation of rubber. Philos. Trans. R. Soc. A243, 251-288.

Simo, J., 1987. On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects. Comput. Methods Appl. Mech. Engng. 60, 153-173.

Wineman, A.S., Huntley, H.E., 1994. Numerical simulation of the effect of damaged induced softening on the inflation of a circular rubber membrane. Int. J. Solids Struct. 31, 3295-3313.

Zuniga, A.E., Beatty, M.F., 2002. A new phenomenological model for stress-softening in elastomers. Z. Angew. Math. Phys. 53, 794-814.