# Experimental investigation and theoretical modelling of induced anisotropy during stress-softening of rubber

Gilles Marckmann, Grégory Chagnon, Matthieu Le Saux, Pierre Charrier

► To cite this version:

Gilles Marckmann, Grégory Chagnon, Matthieu Le Saux, Pierre Charrier. Experimental investigation and theoretical modelling of induced anisotropy during stress-softening of rubber. International Journal of Solids and Structures, 2016, 97-98, pp.554-565. 10.1016/j.ijsolstr.2016.06.028. hal-01377314

HAL Id: hal-01377314

https://hal.science/hal-01377314

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>

Distributed under a Creative Commons Attribution 4.0 International License
Experimental investigation and theoretical modelling of
induced anisotropy during stress-softening of rubber

G. Marckmann<sup>a,*</sup>, G. Chagnon<sup>b</sup>, M. Le Saux<sup>a</sup>, P. Charrier<sup>c</sup>

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

<sup>b</sup>Univ. Grenoble Alpes, TIMC-IMAG, CNRS, F-38000 Grenoble, France

<sup>c</sup>Trelleborg Vibracoustic, CAE Durability Prediction department, Global Advanced Engineering, F-44474 Carquefou, France

Abstract

The Mullins effect refers to a stress-softening phenomenon of rubber-like
materials during cyclic loading. Anisotropy of the material behaviour is
generally observed after stretching. In this paper, a large set of original
suitable experiments are reported to characterise this effect under several
deformation conditions. Then, a phenomenological model is derived to cap-
ture the anisotropic distribution. For that, the affine micro-sphere model
(Miehe et al., 2004) is amended with a directional network alteration in or-
der to describe anisotropy. The alteration process, involving the breakage
and the slippage of the links embedded in the macromolecular network, is
modelled by the evolution of the average number of monomer segments per
chain during stretching. The average chain length and the chain density
are incrementally described by functions to allow both softening and stiffen-
ing, depending to the maximum and the minimum stretch rates and levels
endured in each direction. The good capacity of the model to reproduce

<footnote>*Corresponding author: gilles.marckmann@ec-nantes.fr</footnote>
experimental observations validates the above assumptions.

*Keywords:* Mullins effect, stress softening, constitutive equation, network alteration theory, anisotropy.

# 1. Introduction

Industrial competition in automotive industry has led to a growing use of numerical simulation in the design process of elastomeric parts such as tires, seals, engine mounts, shock absorbers... A good understanding and a well-established modelling of the material behaviour are essential prereq-uisites for efficient numerical modelling as Finite Element analysis. Elas-tomers are known to exhibit a strongly non-linear elastic behaviour under static loading conditions, characterised by the well-known S-shaped uniaxial response. Inelastic features such as hysteresis, stress-softening and residual strain are observed under cyclic loading. The strain-induced stress-softening phenomenon, widely known as the Mullins effect (Mullins, 1947, 1969), is characterised by a reduction of stiffness during the first loading cycles. More-over, the material response is perceptibly more compliant at strain smaller than the maximum strain endured during this first cycle, as pointed out by many authors (Mullins, 1947; Bonart, 1968; James and Green, 1975; Muhr et al., 1999; Robisson, 2000; Pawelski, 2001; Chagnon, 2003; Dorfmann and Pancheri, 2012; Machado et al., 2012; Merckel et al., 2012; Rebouah and Chagnon, 2014), the level of softening is not identical in all directions and present induced anisotropy. Moreover, rubber-like materials subjected to cy-cling loading conditions do not return to their initial state when unloaded, and residual strain is generally observed after unloading.
The Mullins effect does not find unanimous physical explanation at present (Diani et al., 2009) and occurs both in unfilled and filled elastomers, though its influence is more pronounced in filled rubber-like materials. To describe the strain-induced anisotropy of the Mullins effect, different approaches were developed. Besdo and Ihlemann (2005) introduced a tensorial history function based on second-order tensors evaluated using ellipsoid representations of all right Cauchy-Green tensors previously encountered during the deformation process. Shariff (2006) proposed a phenomenological damage anisotropic model based on the assumption that stress-softening is related to principal stretches, and to maximum and minimum strain values of the principal-direction line elements encountered during the deformation history. Göktepe and Miehe (2005) investigated the molecular damage concept developed by Govindjee and Simo (1991), and decomposed the network structure into a particle-to-particle and a crosslink-to-crosslink part. They exclusively associated the anisotropic stress-softening phenomenon to the particle-to-particle contribution, which involves breakage of bonds between chains and filler particles, and associated a scalar damage parameter with each direction of the hyperelastic micro-sphere model (Miehe et al., 2004). Thereby, the model inherently reproduces the strain-induced anisotropic distribution of the Mullins effect and naturally includes residual strain. Diani et al. (2006b) combined the Pawelski (2001) directional molecular-based hyperelastic model with the Marckmann et al. (2002) alteration theory and derived an anisotropic model for the Mullins effect and the inherent induced residual strain. Dargazany and Itskov (2009) proposed a similar approach by taking into account the existence of different chains with different lengths in each direction. In a
similar approach, Itskov et al. (2010) proposed three damage evolution func-
tions for the three principal stretches. Dorfmann and Pancheri (2012) used a
pseudo elastic approach to take into account the stress softening and changes
in material symmetry. Merckel et al. (2013) used the strain amplification fac-
tor concept introduced in the microsphere approach. Rebouah et al. (2013)
and Machado et al. (2014) proposed to develop a constitutive equation that
depends both on the maximum and actual deformation in a full network
approach.

A residual strain is generally observed after removal of the stress. The
magnitude of this residual strain increases with the amount of fillers in rubber
and with the maximum stretch endured before unloading (Dorfmann and Og-
den, 2004). To model both stress-softening and residual stretch, Ogden and
Roxburgh (1999) used the previous pseudo-elastic approach to describe the
Mullins effect, and incorporated a single additional damage variable function
of the maximum strain previously achieved and for which no restriction is
henceforth imposed. Holzapfel et al. (1999) extended this model to the inclu-
sion of an anisotropic description of the damage mechanism governed by three
damage variables related to the evolution of principal stretches. Dorfmann
and Ogden (2004) also investigated the theory of pseudo-elasticity and devel-
oped a model combining stress-softening and residual strain by including two
additional parameters respectively related to these inelastic phenomena. The
authors pointed out that the anisotropic material response under uniaxial de-
formation (transversal isotropy) may be described by a strain energy function
which is not symmetrical in terms of principal stretches (i.e. which does not
satisfy the well-known Valanis and Landel (1967) assumption). They as-
sociated the induced anisotropy only with the residual strain variable but
specified that a more general model, where the induced anisotropy is related
to the softening variable, can be derived. Recently, Rebouah and Chagnon
(2014) proposed a directional model were the stress becomes zero in some
directions before the zero deformation is obtained to generate residual elon-
gation.

The aim of the present study is to experimentally exhibit and model the
evolution of stiffness of rubber in different directions depending on the his-
tory and the directions of first stretching. The residual strain and the loss of
isotropy inside the material will be analysed through a large panel of loading
conditions. In this way, we propose a phenomenological constitutive model
established under micro-mechanical considerations and that principally ac-
counts for the anisotropy induced by the Mullins effect. The first key con-
tribution is concerned with the investigation of a set of original appropriate
homogeneous experiments conducted under various deformation conditions.
Results obtained lead to a thorough understanding of the anisotropic na-
ture of stress-softening and the residual strain and highlight the occurrence
of a stiffening effect along contracted directions (stretch lower than 1) be-
yond the pre-softening stretch level. Pursuing the network alteration theory
proposed by Marckmann et al. (2002), these inelastic phenomena are asso-
ciated with micro-structural alteration mechanisms involving slippage and
breakage of links between polymer chains. This rearrangement is reflected
by the increase of chain contour length in tensile stretched regions and its
slow decrease along contracted directions.

The paper is organised as follows. Section 2 reports experimental inves-
## 2.1. Simple tension tests of previously softened specimens

The observation of the influence of anisotropic effects requires mechanical tests involving different stretching directions. One of them consists first in performing five pre-softening uniaxial tensile cycles up to a stretch equal to 3.25 on a virgin large flat dumbbell sample extracted from a calendered plate. The later sample is 110 mm long and 40 mm large in the useful area. The thickness is equal to 2 mm. Smaller dumbbell samples are then cut from this large specimen with different orientations with respect to the direction of pre-softening stretch (reference direction)<fnref n="1" />. The later are 10 mm long and 2.5 mm large in the useful area. Fig. 1(a) shows how the small dumbbell samples are extracted from the large one. Then, they are submitted to tensile stretch loading until breakage.

Small samples are also cut in a virgin large specimen and similar tensile tests were performed in order to evaluate the initial anisotropy of the plate. Fig. 1(b) shows the evolution of the nominal stress (first Piola-Kirchhoff stress) $\pi$ as a function of the principal stretch $\lambda$ along the loading direction. These experimental results exhibit an appreciable stiffer mechanical response when approaching the calendering direction.

As we only focus on the anisotropy induced by Mullins effect and not by the calendering process, the initial anisotropy observed above is taken into account by normalizing the mechanical response of specimens for each direction. The calendering direction is retained as the reference response (0°-direction) and a corrective factor is calculated for each angle. This factor is

<footnote n="1">¹All the samples cut in a specimen which was already stretched will be called “secondary specimens”.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Anisotropic distribution of stress-softening developed under tensile stretching. (a): geometry of the large and the secondary small dumbbells specimens. (b): initial anisotropy of the large *virgin* dumbbell samples extracted from a calendered plate. Stress-stretch responses of small test pieces cut at various angles from the reference direction are presented.</figcaption></figure>

evaluated so that all virgin curves are superimposed with the reference one (cf. Table 1).

We observe that this factor is linear with the angle. Then, mechanical responses of the pre-softened small samples are corrected by using the aforementioned corrective factor depending on the extraction angle.

Note furthermore that, due to residual strain, the dimensions of the small specimen are necessarily measured in a pre-softened configuration. Since the pre-softening tensile test generates an appreciable residual stretch $\lambda_{res}$, stress-free configurations before and after this pre-softening are different (see
<table><caption>Table 1: Value of corrective factors depending on the direction of extraction of the small specimen relatively to the calendered direction of the large specimen</caption><thead><tr><td>angle of extraction (°)</td><td>corrective factor</td></tr></thead><tbody><tr><td>0</td><td>1.0</td></tr><tr><td>22.5</td><td>1.036</td></tr><tr><td>45</td><td>1.072</td></tr><tr><td>67.5</td><td>1.108</td></tr><tr><td>90</td><td>1.144</td></tr></tbody></table>

Fig. 2).

<figure><img src="image_3.png" /><figcaption>Figure 2: Cinematics of the deformation</figcaption></figure>

Distribution of the residual stretches naturally depends of the angle of extraction $\theta$. Comparison of stress-stretch responses of small samples rigorously requires the use of a unique reference configuration. The total Lagrangian approach is used in the numerical experiments detailed in Section 3. Thus, the unstrained virgin material configuration is adopted as the reference one. Reference dimensions of the cut samples are determined with respect to this configuration. Thus, they are determined from the ratio of the small samples initial dimensions in the pre-softened configuration to the residual stretches
in the corresponding principal directions of the secondary tensile deformation.

<figure><img src="image_4.png" /><figcaption>Figure 3: Reference basis of the small specimen extracted from the large pre-softened specimen.</figcaption></figure>

Considering residual stretches, defined by the residual deformation gradient tensor $F_{res}$, and the principal basis of deformation $\{e_i\}_{i=1...3}$, residual stretches are defined in a new basis, $e_i^* = Re_i$, which is the strain principal directions of the cut samples during the succeeding stretches (see eg. 3). and where $R$ is the rotation tensor between $e$ and $e^*$. Then, the new reference deformation gradient tensor is given by:

$$F_{res}^* = R^T F_{res} R \quad (1)$$

Fig. 4 shows stress-stretch responses of secondary specimens extracted at different angles (0°, 22.5°, 45°, 67.5° and 90°) from the reference direction.

The nominal stress and stretch are both expressed in the virgin unstrained configuration and the stretch-stress curves are corrected according to the initial anisotropy of the *virgin* material as described above. Mechanical responses of the pre-deformed specimens appear to be lower than the virgin material one. This reflects the strain-induced stress-softening phenomenon. The anisotropic distribution of stress-stretch properties is also clearly ob-
<figure><img src="image_5.png" /><figcaption>Figure 4: Anisotropy induced by the Mullins effect under cyclic tensile conditions up to a stretch of 3.25. Comparison of tensile stress-stretch uniaxial responses of a virgin sample and pre-softened samples cut at different angles from the reference direction. Error bars shows experimental dispersal on four samples.</figcaption></figure>

servable. Although the prior tensile stretch produced stress-softening in all selected directions, the softening level significantly decreases as the direction deviates from the one of the initial stretch. It also appears to be nearly null in the orthogonal direction (90° in the figure). Furthermore, all secondary specimens are conditioned to the same stretch ratio. The material gets stiffer when the stretch level approaches the maximum stretch endured during the pre-softening stretch level (i.e. λ = 3.25). The material does not completely recover its virgin behaviour beyond the pre-softening stretch level (cf. table 2).

Indeed, due to residual strain, the secondary loading curves should rejoin
the virgin curve for stretch levels slightly higher than the maximal previously
<table><caption>Table 2: Residual strain of the small specimens in their principal axes</caption><thead><tr><th>angle of extraction (°)</th><th>λ<sub>1res</sub></th><th>λ<sub>2res</sub></th><th>λ<sub>3res</sub></th></tr></thead><tbody><tr><td>0</td><td>1.05</td><td>0.975</td><td>0.9768</td></tr><tr><td>22.5</td><td>1.0390</td><td>0.9859</td><td>0.9768</td></tr><tr><td>45</td><td>1.0125</td><td>1.0125</td><td>0.9768</td></tr><tr><td>67.5</td><td>0.9859</td><td>1.0390</td><td>0.9768</td></tr><tr><td>90</td><td>0.975</td><td>1.05</td><td>0.9768</td></tr></tbody></table>

endured one (Johnson and Beatty, 1993). As Muhr et al. (1999) pointed out, the softening generally remains above the prior maximum endured stretch level (Diani et al., 2009). The authors specified that application of the residual strain correction, usually adopted for the idealisation of the Mullins effect, would reduce this drop below the virgin curve.

Mechanical responses of the 0°, 22.5° and 45°-small samples exhibit a slope close to the virgin curve one beyond the maximum stretch previously reached. On the other hand, relatively to the virgin material response, significant stiffening appears for angles of 67.5° and 90° when the pre-softening stretch level is exceeded. The latest samples have submitted stretches lower than 1 during the pre-softening cyclic test. Consequently, it seems that stress-softening phenomenon is coupled with a stiffening effect in these directions. Even if the observed effect is rather small in comparison with the experimental dispersal, the same phenomenon is observed for each lot of samples.

We may presume that stiffening is linked to the compression loading conditions. This assumption is consistent with the common perception of the Mullins effect which is generally ascribed to the quasi-irreversible rearrange-
ment of the molecular network under stretching. In that case, the material is assumed to have a selective memory of the maximum stretch ratio encountered during its loading history and to adapt a new equilibrium configuration softer than the virgin one as long as this maximum stretch is not exceeded. Finally note that, in conformity with the incompressibility of the material, samples cut at 0° and 22.5° reveal a residual stretch greater than 1, those oriented at 67.5° and 90° manifest a contraction residual stretch (< 1) whilst the specimen extracted at 45° does not show any.

## 2.2. *Pure shear tests from previously softened specimens*

Similar experiments are conducted with a pure shear pre-softening test. Here, the specimen used for the pure shear test is moulded and is assumed to be initially homogeneous and isotropic unlike the calendered plate. The specimen is 300 mm long, 40 mm width and 2 mm thick. The width inside fixation rolls is equal to 35 mm. It is submitted to five cycles of pre-softening pure shear loading up to a stretch of 4 (Fig. 5: the width of the specimen is stretched of 400%).

<figure><img src="image_6.png" /><figcaption>Figure 5: The pure shear sample and the secondary small dumbbells samples used for the characterisation of the anisotropic distribution of stress-softening developed under pure shear stretching</figcaption></figure>

As shown in Fig. 5, small dumbbell samples are cut from this pre-softened specimen at different orientations (0°, 30°, 45°, 60° and 90°) in comparison with the reference direction (belonging to the pre-softening stretch one).
and undergo uniaxial stretch until breakage. The corresponding experimen-
tal results, expressed in the unstrained virgin configuration, are shown in
Fig. 6. The level of softening that occurs after the pre-softening pure shear
test is clearly as lower as the direction deviates from the pure shear load-
ing direction. Observations similar to those relative to the previous uniaxial
tensile tests can be made about the recovering of the virgin material be-
haviour beyond the pre-softening stretch level (λ = 4). However, a stiffening
effect which occurs in directions close to the direction orthogonal to the pre-
softening uniaxial tensile stretch loading is not clearly observable here. Since
the material is not stretched nor contracted in the direction orthogonal to the
pure shear loading one, this observation coincides with the previous analysis
where the stiffening phenomenon appears when stretches are lower than 1.
The smaller the angle is (in comparison with the reference one), the greater
the amount of residual stretch is. The sample oriented to 90° does not exhibit
any residual strain.

2.3. *Modified biaxial tension tests from previously softened specimens*

Tests inspired by the so-called modified biaxial tension test proposed by Pawelski (2001) are performed. These experiments consist beforehand in pre-softening successively in two directions and then in cutting small secondary samples at different angles in order to evaluate the anisotropic effects induced by stress-softening. As illustrated in Fig. 7, a cross-shaped sample cut in a calendered plate is considered and its *arms* are slitted to ensure the homogeneity of deformation in the central area. The central area is a square 40 mm long and 2 mm thick. The lengths of the *arms* are equal to 60 mm. This specimen is first submitted to five cycles perpendicularly to
<figure><img src="image_7.png" /><figcaption>Figure 6: (a) Anisotropy induced by the Mullins effect under cyclic pure shear conditions up to $\lambda = 4$. Comparison of tensile stress-stretch uniaxial responses of a virgin sample and samples cut at different angles from the reference direction, (b) close-up to residual stretches</figcaption></figure>

the calendering direction with a standard tensile test machine, while the two arms orthogonal to the loading direction remain free. A similar pre-softening loading is next applied in the calendering direction called reference direction. Afterwards, small secondary dumbbell samples are extracted in the central area of the pre-softened cross-shaped specimen at different angles (0°, 22.5°, 45°, 67.5° and 90°). Then, simple tensile tests are performed until breakage.

First, pre-softening tensile stretches of 3.25 are applied in the central area of the cross-shaped sample in both directions. Mechanical responses of the cut samples are compared to the unstrained virgin material response in Fig. 8. Note that the calendered plates used for this experiment are similar to the
<figure><img src="image_8.png" /><figcaption>Figure 7: The cross-shaped sample and the secondary small dumbbell samples used for the characterisation of the anisotropic distribution of stress-softening induced by modified *equibiaxial* stretching. (a) Undeformed, (b) deformed and (c) pre-softened cut configurations of the cross-shaped specimen.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Anisotropy induced by the Mullins effect under cyclic modified *equibiaxial* conditions up to a stretch of 3.25 in both directions. Comparison of tensile stress-stretch uniaxial responses of a virgin sample and samples cut at different angles from the reference direction.</figcaption></figure>
ones used for the uniaxial tensile test of section 2.1. A correction is performed
to take into account the initial anisotropy of the plates. The responses of
the material are reported in Fig. 8. As expected, the 45°-secondary sample
exhibits the smaller amount of softening since the maximum pre-stretch en-
dured in this direction is the lowest one. However, a remaining significant
softening occurs beyond the pre-softening stretch level. The 22.5° and 0°-
samples are slightly stiffer than the virgin sample beyond the prior maximum
stretch ratio, while 67.5° and 90°-samples exhibit a significant drop below the
virgin curve. Consequently, it seems that, under the maximum stretch ra-
tio endured during the pre-softening test, the succession of compressive and
tensile loading affects the material response. Indeed, due to the testing pro-
cedure, this experiment differs from a pure equi-biaxial test. Load histories
in the two pre-stretching directions are different. The 22.5° and 0°-directions
successively experience contraction followed by tensile stretch, whilst the
67.5° and 90°-directions undergo tensile stretch followed by compressive one.
These observations confirm the previously assumed existence of a directional
chain network rearrangement which takes place during stretching. Mecha-
nisms involved in this process may generate a new stress-free equilibrium
configuration stiffer than the virgin one after compression.

An analogous test is conducted by applying a pre-stretching of 3.25 in
the first direction of the cross-shaped specimen, perpendicular to the calen-
dering one, whereas the reference direction is stretched up to a factor 2.25.
Responses are depicted in Fig. 9. The respective contributions of the two
pre-softening contributions, as well as the induced anisotropy, are clearly
distinguishable. Remarks similar to those enunciated above may be stated
<figure><img src="image_10.png" /><figcaption>Figure 9: Anisotropy induced by the Mullins effect under cyclic modified biaxial conditions up to a stretch of 2.25 in the first direction and 3.25 in the second one. Comparison of tensile stress-stretch uniaxial responses of a virgin sample and samples cut at different angles from the reference direction.</figcaption></figure>

concerning the influence of the loading history undergone in each material direction.

## 2.4. Conclusions

All the experimental observations made in the present section clearly highlight the anisotropic nature of the Mullins effect. The angular distribution of stress-softening occurring under various homogeneous deformation modes apparently depends on the maximum stretches undergone by the material. However, an induced stiffening phenomenon along contracted directions was exhibited and the influence of the loading history is pointed out. This large variety of experimental results will be used to validate the model
derived in the next section.

### 3. Constitutive model

The aim of this section is to build a constitutive model for rubber which is able to take into account anisotropy induced by pre-stretches. Mathematical symbols are defined at the end of the paper.

The previous experimental observations are consistent with a network alteration as proposed by Marckmann et al. (2002) where alteration depends on the stretch endured by chains and thus on the direction of stretching.

#### 3.1. Description of the model

Rubber elasticity can be described through micro-mechanically based models. As they are explicitly expressed in terms of these two parameters, the chain density $n$ and the average number of monomers per chains $N$, micromechanical network models based on non-Gaussian chain statistics are well-adapted to the use of stress-softening.

In the paper, the affine “Micro-sphere” model of Miehe et al. (2004) is chosen as the reference elastic model with a set of 42 directions. For convenience, denote $\lambda^c = \lambda(\mathbf{d}^c) = \|\mathbf{F} \cdot \mathbf{d}^c\| = \sqrt{\mathbf{d}^c T \cdot \mathbf{C} \cdot \mathbf{d}^c}$ the stretch of a single chain oriented in the direction of the unit vector $\mathbf{d}^c$. The representation of the macroscopic Cauchy stress tensor under the incompressibility assumption ($\det \mathbf{F} = 1$) is given by :

$$ \sigma = C_R \sqrt{N} \sum_{c=1}^{r} \left[ \frac{\omega^c}{\lambda^c} \mathcal{L}^{-1} \left( \frac{\lambda^c}{\sqrt{N}} \right) (\mathbf{F} \cdot \mathbf{d}^c) \otimes (\mathbf{F} \cdot \mathbf{d}^c) \right] - p\mathbf{I} \quad (2) $$
where **d**^c represents the direction vector of the chain *c*, $\lambda(\mathbf{d}^c)$ is the stretch of the chain *c*, $\omega^c$ is the weight factor associated to the chain *c*, $C_R = nkT$ represents the material stiffness expressed in terms of the chain's density *n*, **I** is the unit tensor and *p* is a pressure classically introduced to impose incompressibility. This hyperelastic constitutive equation contains two physically-motivated material parameters *N* and $C_R$.

Besides the alteration processes previously stated by Marckmann et al. (2002), we can presume that reversible mechanisms should be involved. In particular, we can mention the slippage of chains over filler particles, drawing of chains from coils or from their layers coating the surface of filler particles, slippage of macromolecular chains at topological entanglement junctions (Svistkov and Komar, 2005), etc. This latest mechanism may explain the stiffening phenomenon observed in Section 2 along contracted directions. From a statistical point of view, chains embedded in the macromolecular network are improbably oriented along a straight line. Considering a chain passing through an entanglement loop formed by another molecule, as shown in Fig. 10(a), ends of this chain are certainly oriented in different directions. Application of an uniaxial tensile external force induces the slippage of the chain through this entanglement junction. As a chain is a sequence of monomer segments between two subsequent junctions, the length and consequently the number of monomers of the chain oriented along the tensile direction increases (i.e. the part of the chain situated above the junction point in Fig. 10(b)), whilst the one of the other chain (i.e. the part of the chain situated below the junction point) decreases. We can presume that interactions with the neighbouring molecular network or eventually the
quasi-irreversible tightening of the entanglement loop make this change of conformation not fully recoverable during unloading (Fig. 10(c)). Another consequence is that the number of monomers in each direction changes and the mass conservation cannot be considered in each direction, but globally (if any).

<figure><img src="image_11.png" /><figcaption>Figure 10: Schematic representation of a chain, depicted in bold line, slipping through an entanglement loop formed by another macromolecule (in grey line): (a) initial undeformed configuration, (b) deformed configuration under uniaxial tensile loading and (c) stress-free configuration after application of the stretching.</figcaption></figure>

As a result, the breakage and slippage mechanisms occurring upon stretching lead to a new equilibrium stress-free configuration that exhibits memory effects of the loading pre-history. This behaviour could be considered as a transformation of the network where internal variables are evolving under different paths depending on loading history. This brings to the association of stress-softening and stiffening phenomena respectively with stretched and contracted directions. The average chain length increases under tension and
decreases under contraction.

With a view to simplification, material directions of the idealized network are assumed to be independent one on another. Stress-softening and stiffening phenomena in a direction are supposed to only depend on the alteration process in this direction, involving both breakage and slippage mechanisms. The alteration process is modelled by considering an alteration function of both a stress-softening and a stiffening variable, respectively denoted $\alpha^c$ and $\beta^c$, which affect the chain oriented in the direction $d^c$ of the discrete network. The alteration function is identical for all directions (isotropy of the mechanism) and describes the evolution of the average number of monomers per chain upon loading. The material is assumed to have a directional selective memory of the alteration process such that for $c = 1, ..., r$:

$$N^c(t) = N^0 f(\dot{\alpha}^c, \dot{\beta}^c, \alpha^c, \beta^c), \quad (3)$$

where $N^0 = N^c(t = 0)$ denotes the average number of monomers per chain within the *virgin* material and where alteration variables $\alpha^c(t)$ and $\beta^c(t)$ are defined along the considered chain direction $c$ by:

$$\alpha^c(t) = \max_{0 \le \tau \le t} [\lambda^c(\tau)] \quad (4)$$

and

$$\beta^c(t) = \min_{0 \le \tau \le t} [\lambda^c(\tau)] \quad (5)$$

$\dot{x}$ denotes the rate of $x$.

The present model differs from those proposed by Diani et al. (2006a) or Göktepe and Miehe (2005) because the alteration function does not only depend on the maximum stretch. To simplify, we assume that the evolution
of the number of monomers by chain only depends on the history of the each direction. We also assume that the material presents a continuous and discontinuous evolution law as proposed by Miehe (1995). A term of continuous damage $D^c$ is introduced. $D^c$ is supposed to increase whatever the chain increases or decreases, and only depends on the stretch rate $|\dot{\lambda}^c|$. A term of discontinuous evolution is added to take into account the dependence upon the maximum and the minimum stretches (respectively $\alpha^c$ and $\beta^c$) and the associated rates (respectively $\dot{\alpha}^c$ and $\dot{\beta}^c$).

$$ (N^c \dot{C}_R^c) = N^0 C_R^0 (1 - D^c) g(\dot{\alpha}^c, \dot{\beta}^c, \alpha^c, \beta^c) \quad (6) $$

with $D^c$ defined by $\dot{D}^c = D(|\dot{\lambda}^c|)$.

### 3.2. Thermodynamical aspects

Within this framework, the amended macroscopic strain energy density is defined by

$$ W(\mathbf{F}, \alpha, \beta, \dot{\alpha}, \dot{\beta}, \mathbf{D}) \approx \sum_{c=1}^{r} \omega^c (1 - D^c) n^c w(\lambda^c, N^c, C_r^c) - p\mathbf{I} \quad (7) $$

The alteration variables $\alpha = (\alpha^c)_{c=1,...,r}$, $\beta = (\beta^c)_{c=1,...,r}$ and $\mathbf{D} = (D^c)_{c=1,...,r}$ may be considered as internal variables. Then the evolution of $N^c$ must be consistent with the Clausius-Planck inequality that guarantees a positive macroscopic dissipation. The dissipation $D_{int}$ is defined for isothermal deformation by:

$$ D_{int} = \mathbf{\Pi} : \dot{\mathbf{F}} - \dot{W} \ge 0 \quad \forall \mathbf{F} \text{ with } \det \mathbf{F} = 1, \quad (8) $$
According to the model of Miehe et al. (2004) and Eqs. (6), the strain energy density $W$ can be written as:

$$W = C_R^0 N^0 \sum_{c=1}^{r} \left[ \omega^c (1 - D^c) g \left( \dot{\alpha}^c, \dot{\beta}^c, \alpha^c, \beta^c \right) \left[ \frac{\lambda^c}{\sqrt{N^c}} \xi^c + \ln \left( \frac{\xi^c}{\sinh \xi^c} \right) \right] \right] + W_0 \quad (9)$$

with $\xi^c = \mathcal{L}^{-1}(\lambda^c/\sqrt{N^c})$. We note $\chi(\xi) = \frac{\lambda^c}{\sqrt{N^c}} \xi^c + \ln(\frac{\xi^c}{\sinh \xi^c})$.

$N^c$ depends on $\alpha^c$ and $\beta^c$ while $(N^c C_R^c)$ depends on $\alpha, \beta, \dot{\alpha}$ and $\dot{\beta}$, and finally $\dot{D}$ depends on $|\dot{\lambda}^c|$. Therefore the expression of the dissipation $\mathcal{D}_{int}$ is obtained by developing partial derivatives of each term.

$$ \begin{aligned} \mathcal{D}_{int}^c = &+ \frac{\partial g}{\partial \dot{\alpha}^c} \frac{\partial \dot{\alpha}^c}{\partial t} - \omega^c C_R^0 N^0 (1 - D^c) \chi(\xi) \left[ \frac{\partial g}{\partial \alpha^c} \dot{\alpha}^c + \frac{\partial g}{\partial \beta^c} \dot{\beta}^c + \frac{\partial g}{\partial \dot{\beta}^c} \frac{\partial \dot{\beta}^c}{\partial t} \right] \\ &- \omega^c C_R^0 N^0 (1 - D^c) g N^0 \frac{\partial \chi(\xi)}{\partial N^c} \left[ \frac{\partial f}{\partial \alpha^c} \dot{\alpha}^c + \frac{\partial f}{\partial \dot{\alpha}^c} \frac{\partial \dot{\alpha}^c}{\partial t} + \frac{\partial f}{\partial \beta^c} \dot{\beta}^c + \frac{\partial f}{\partial \dot{\beta}^c} \frac{\partial \dot{\beta}^c}{\partial t} \right] \\ &+ \omega^c C_R^0 N^0 g \chi(\xi) \dot{D}^c = \mathcal{D}_{int}^{c(g)} + \mathcal{D}_{int}^{c(f)} + \mathcal{D}_{int}^{c(D)} \end{aligned} \quad (10) $$

$$ \mathcal{D}_{int} = \sum_{c=1}^{r} (\mathcal{D}_{int}^{c(g)} + \mathcal{D}_{int}^{c(f)} + \mathcal{D}_{int}^{c(D)}) = \mathcal{D}_{int}^{(g)} + \mathcal{D}_{int}^{(f)} + \mathcal{D}_{int}^{(D)} \quad (11) $$

It is to note that the damage always increases, as the bounds of the stretches, thus $\dot{D}^c \ge 0$, $D^c \ge 1$, $\dot{\alpha}^c \ge 0$ and $\dot{\beta}^c \le 0$. Nevertheless the stiffening phenomenon imposes $\partial N^c / \partial \beta^c \ge 0$. So, the only way to satisfy the positivity of $\mathcal{D}_{int}$ is that positive terms in Eq. (10) balance negative ones.

For a given load (for instance a general biaxial extension), stretches $\lambda_{i \in 1,2,3}$ are linked by the kinematic constraint of incompressibility $(\lambda_1 \lambda_2 \lambda_3 = 1)$ unlike $(\alpha^c)_{c=1,...,r}$ and $(\beta^c)_{c=1,...,r}$. Reported on the Micro-sphere representation, $\alpha$ and $\beta$ form yield surfaces whose evolution are directional and history dependent (see Fig. 11). Unfortunately, these surfaces are not regular and not convex in the undeformed configuration. We can only notice that the lower
<figure><img src="image_12.png" /><figcaption>Figure 11: Representation of the internal state variables in the Micro-sphere. Vectors $d^c$ are scaled respectively with 1 (white), with $\alpha^c$ (red) and with $\beta^c$ (green) for a biaxial extension (a) and an uniaxial extension (b). The three surfaces are concentric; offsets on axis 3 have been applied for illustration.</figcaption></figure>

value of $\beta$ is generally reached in a direction orthogonal to a direction where $\alpha$ is the greater. The model induces a complex expression of $D_{int}$ which does not allow a simple relation between terms to enforce its positiveness. We propose to verify it a posteriori in section 4.2.

# 4. Comparison of the model with experimental data

## 4.1. Proposed model and identification

The aim of this section is to determine the alteration function involved in the model and to estimate the material parameters in order to reproduce
the trend of experimental data. The ability of the formulation to predict stress-softening, stiffening and residual strain phenomena is assessed.

Since evolution of the average chain contour length cannot be experimentally observed, the evolution of alteration during stretching is described by the following empirical alteration functions:

$$f(\dot{\alpha}^c, \dot{\beta}^c, \alpha^c, \beta^c) = \exp(a_1(\alpha^c - 1)) + b_1(\beta^c - 1), \quad (12)$$

$$g(\dot{\alpha}^c, \dot{\beta}^c, \alpha^c, \beta^c) = a_3 \exp(a_3(\alpha^c - 1))\dot{\alpha}^c + b_3 \exp(b_3(\alpha^c - 1))\dot{\beta}^c \quad (13)$$

where $a_1, b_1, a_3, b_3$, and $b_4$ are material parameters. These equations must satisfy conditions Eq. (10). Note that, insofar as it clearly depends on the nature of the material and the forming process, this form is specific to the material used in the present study. The current average number of monomers per chain $N^c(t)$ then follows Eq. (3). The larger is the value of $N^c$, the greater is the chain locking stretch value.

The damage function is chosen as:

$$\dot{D} = \nu|\dot{\lambda}| \quad (14)$$

where $\nu$ is a material parameter. That means that less and less macro-molecular chains are concerned by the mechanical response of the material. Nevertheless, by assuming $\nu = 0$, the density of active macromolecules stay constant. The material parameters were fitted on all the experimental tests, the values are summarized in Table 3.

An illustration of the abilities of the model is presented in Fig. 12. The model demonstrates its possibilities to describe cyclic loadings, i.e., the difference between the first and second loadings and the difference with the next
<table><caption>Table 3: Values of the material parameters fitted on the experimental results of Section 2</caption><thead><tr><th>Parameters</th><th>C<sup>c</sup><sub>R0</sub>N<sup>c</sup><sub>0</sub></th><th>N<sup>c</sup><sub>0</sub></th><th>a<sub>1</sub></th><th>a<sub>3</sub></th><th>b<sub>1</sub></th><th>b<sub>3</sub></th><th>b<sub>4</sub></th><th>ν</th></tr></thead><tbody><tr><td>Values</td><td>0.239</td><td>5.65</td><td>0.410</td><td>0.370</td><td>0.578</td><td>0.669</td><td>0.796</td><td>0.00590</td></tr></tbody></table>

loading. Moreover by changing the material parameters along the direction, the model can be initially anisotropic, and the different evolutions of the parameters according to the direction induced anisotropy.

<figure><img src="image_13.png" /><figcaption>Figure 12: Example of a cyclic loading with the "extension" component and the "path" component of the stress softening model.</figcaption></figure>

First, the uniaxial tension test with secondary specimens is numerically reproduced by firstly simulating kinematics of the preconditioning tensile
test in order to evaluate the new stress-free configuration (average chain length and density, prior maximum and minimum chain stretches) in each discretized direction. While considering these characteristics as initial ones, analytical simulation of the subsequent tensile tests on cut samples is then achieved for the retained set of loading directions (oriented at 0°, 22.5°, 45°, 67.5° and 90° from the pre-softening loading direction). As shown in Fig. 13, the proposed model provides a good agreement with experimental data depicted in Fig. 4, more precisely, it successfully predicts the anisotropic distribution of stress-softening, stiffening and residual strain phenomena under uniaxial tensile stretching.

<figure><img src="image_14.png" /><figcaption>Figure 13: Model prediction (right) of the anisotropy induced by the Mullins effect under cyclic tensile conditions up to a stretch of 3.25, compared with experiments (left).</figcaption></figure>

Fig. 14 presents the evolution of variables involved in the proposed formulation and highlights the influence of history of each chain in the global response of the specimen. Loading is decomposed in time ranges. From $t = 0$
<figure><img src="image_15.png" /><figcaption>Figure 14: Evolution of internal variables during uniaxial stretch for chains oriented at 0° and 90° respectively.</figcaption></figure>

to $t_1$, the micro-sphere is stretched in direction 90° to $\lambda = 3.25$. Then from $t_1$ to $t_2$ the microsphere is unstained, then stretched in direction 0° to $\lambda = 4.5$ (reaching value $\lambda = 3.25$ at time $t_3$). Finally, the microsphere is unstrained. When the both chains reach the same value of $\alpha_c$ and $\beta_c$ (i.e. $\alpha_{90} = \alpha_0$ and $\beta_{90} = \beta_0$) the parameters $N_c$ are identical ($N_{90} = N_0 = 12.21$) but the damage parameters are different ($D_0 = 0.029$ and $D_{90} = 0.018$) as they depends on the loading paths. That explains the difference of the responses observed in Fig. 13.

Pure shear tests with secondary specimens are then simulated in a similar way. As illustrated in Fig. 15, model predictions correctly describe the
dependence of the material response on the direction, this dependence being
induced by the prior pure shear loading. Nevertheless, the remaining soft-
ening effect, significant on experimental results above the prior maximum
stretch level endured, is not reproduced by the constitutive equation.

<figure><img src="image_16.png" /><figcaption>Figure 15: Model predictions (right) of the anisotropy induced by the Mullins effect under cyclic pure shear conditions up to a stretch of 4, compared with experiments (left).</figcaption></figure>

Results of calculations for the modified *equibiaxial* test with secondary
samples are shown in Fig. 16. The response of the model is corrected by the
initial anisotropy in order to be compared with experimental data of Fig. 8.
Despite a small discrepancy about the amount of softening, the general trend
of data is fairly well described below and beyond the pre-softening stretch
level. While comparing experimental and theoretical results for the modified
biaxial respectively shown in Fig. 9 and 17, a similar conclusion can be drawn.
<figure><img src="image_17.png" /></figure>

Figure 16: Model predictions (right) of the anisotropy induced by the Mullins effect under cyclic modified "equibiaxial" conditions up to a stretch of 3.25 in both directions, compared with experiments (left).

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

Figure 17: Model predictions (right) of the anisotropy induced by the Mullins effect under cyclic modified "biaxial" conditions up to a stretch of 2.25 in the first direction and 3.25 in the second one, compared with experiments (left).
## 4.2. Thermodynamics

As seen in section 3.2, the only way to verify the positiveness of the energy dissipation is to verify it a posteriori. For the proposed equation, the derivatives to introduce in eq. 10 are:

$$ 
\begin{aligned}
\frac{\partial f}{\partial \alpha^c} &= a_1 \exp(a_1(\alpha^c - 1)), \quad \frac{\partial g}{\partial \alpha^c} = a_3 \exp(a_3(\alpha^c - 1))\dot{\alpha}^c + b_3 b_4 \exp(b_4(\alpha^c - 1))\dot{\beta}^c \\
\frac{\partial f}{\partial \dot{\alpha}^c} &= 0, \quad \frac{\partial g}{\partial \dot{\alpha}^c} = a_3 \exp(a_3(\alpha^c - 1)) \\
\frac{\partial f}{\partial \beta^c} &= b_1, \quad \frac{\partial g}{\partial \beta^c} = 0 \\
\frac{\partial f}{\partial \dot{\beta}^c} &= 0, \quad \frac{\partial g}{\partial \dot{\beta}^c} = b_3 \exp(b_4(\alpha^c - 1))
\end{aligned}
\qquad (15) $$

We can notice that, for the present model, the positiveness of the energy dissipation can not be ensured mathematically. This invalidates it for further predictive simulations. But, in this paper, the model is only used to validate the assumption of a stiffening phenomenon involved in the stretching of rubber. In the following, we only try to validate the positiveness of the dissipation of the present simulation, with the material parameters given in table 3.

We first verified the positiveness of the dissipation for the above experiments. Fig. 18 shows the dissipation of the sample as a function of time during the test for the tensile test (corresponding to Fig. s 13 and 14). We can notice that the dissipation $D_{int}$ is always positive. The same verifications are made for all the test-cases with the same conclusion.

To support this verification we propose to use the above fitted model and to impose a quiet complex load history. Load histories are built including cyclic loads of four different types among uniaxial extension, equi-biaxial
<figure><img src="image_19.png" /><figcaption>Figure 18: Dissipation inside the material during the uniaxial tensile test. A first extension to $\lambda = 3.25$ is applied in direction 0° then a second extension of $\lambda = 6$ in a different direction.</figcaption></figure>

extension, general biaxial extension, and pure shear deformations. Then the dissipation $D_{int}$ obtained by the development of Eq. 8 is computed.

An example of arbitrary biaxial extension is given in Fig. 19(a) according to the plan $\lambda_1$ and $\lambda_2$ (because of the incompressibility the state of deformation can be resumed by both parameters $\lambda_1$ and $\lambda_2$). The positiveness of the corresponding dissipation is observed as shown in Fig. 19(b).

This process is repeated for many different load paths and the second principle of thermodynamics is systematically verified.

Fig. 19(b) shows that of negative part $D_g$ of the dissipation due to function g of equation (9) is counterbalanced by the other terms of the dissipation. It is also to note the part $D_D$ of the damage function is negligible with respect to the other terms.
<figure><img src="image_20.png" /><figcaption>Figure 19: (a) Arbitrary biaxial load path given as function of λ₁ and λ₂, (b) corresponding dissipation given as function of λ₂</figcaption></figure>

**5. Concluding remarks**

A relevant experimental procedure was proposed to characterise the anisotropic distribution of stress-softening under various homogeneous deformation modes. Besides the generally observed stress-softening and residual strain phenomena, a stiffening effect has been stated in contracted regions beyond the maximum stretch level previously undergone by the material. Within this framework, we proposed a new network mechanism for the anisotropic description of these inelastic effects, which are considered to result from the breakdown and slippage of junctions bridging polymer chains and/or filler particles. This alteration scenario leads to an increase of the average chain contour length in tensile loaded regions and a decrease in contracted zones. This approach is implemented in a hyperelastic network-based formulation
by describing the evolution of the material parameters through an alteration
function involving both a softening and a stiffening variable. Since the struc-
tural rearrangements are assumed to generate a new equilibrium stress-free
configuration, this evolution is described in an incremental way. Discretized
integration formula, well-adapted to the representation of initial or induced
anisotropy, is retained to numerically perform the homogenisation procedure.
It was shown that, our model reproduces the residual stretch as a natural
consequence of the anisotropic character of the alteration process. Suitability
of the resulting formulation is demonstrated on the basis of available experi-
mental results highlighting the anisotropic material properties generated by
prior loading. Finally, the proposed formulation permits to fit globally a
large set of experimental data that are very different in term of strain states
with few material parameters.

Acknowledgment

We would like to thank Odile Duroc of the Advanced Engineering De-
partment of Trelleborg Automotive for her enthusiasm and help to carry out
the unconventional experiments reported in of this paper.

We also thanks Pr. Erwan Verron and Pr. Laurent Gornet for their help
during the Master-thesis intership of M. Le Saux.
## Notations

$x$ scalar $x$

$f(x)$ $f$ as function of $x$

$\frac{\partial f}{\partial x}$ partial derivative of $f$ with respect to $x$

$\dot{x}$ derivative of $x$ with respect to time

$\mathbf{b}$ vector $\mathbf{b}$

$\mathbf{A}$ tensor $\mathbf{A}$

$\mathbf{b}^T$ transpose of $\mathbf{b}$

$\mathbf{A} \cdot \mathbf{b}$ dot product of $\mathbf{A}$ by $\mathbf{b}$

$\mathbf{A} : \mathbf{B}$ double-dot product of $\mathbf{A}$ by $\mathbf{B}$

$\mathbf{A} \otimes \mathbf{B}$ dyadic product of $\mathbf{A}$ by $\mathbf{B}$

$\approx$ is approximated by...

## References

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

Bažant, Z. P., Oh, B. H., 1986. Efficient numerical integration on the surface of a sphere. Z. Angew. Math. Mech. 66, 37–49.

Besdo, D., Ihlemann, J., 2005. Directional sensitivity of Mullins effect constitutive models for rubber. In: Balkema (Ed.), Constitutive Models for Rubber IV. ECCMR05. pp. 229–235.
Bonart, R., 1968. X-ray investigations concerning the physical structure of cross-linking in segmented urethane elastomers. J. Macromol. Sci.: Physics Ed. B2, 115-138.

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

Chagnon, G., Verron, E., Gornet, L., Marckmann, G., Charrier, P., 2004. 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.

Chagnon, G., Verron, E., Marckmann, G., Gornet, L., 2006. Development of new constitutive equations for Mullins effect in rubber using the network alteration theory. Int. J. Solids Struct. 43, 6817-6831.

Dargazany, R., Itskov, M., 2009. A network evolution model for the anisotropic Mullins effect in carbon black filled rubbers. Int. J. Solids Struct. 46, 2967-2977.

Diani, J., Brieu, M., Gilormini, P., 2006a. Observation and modeling of the anisotropic visco-hyperelastic behavior of a rubberlike material. Int. J. Solids Struct. 43, 3044-3056.

Diani, J., Brieu, M., Vacherand, J. M., 2006b. A damage directional constitutive model for the Mullins effect with permanent set and induced anisotropy. Eur. J. Mech. A/Solids 25, 483-496.

Diani, J., Fayolle, B., Gilormini, P., 2009. A review on the Mullins effect. Eur. Polym. J. 45, 601-612.
Dorfmann, L., Ogden, R. W., 2004. A constitutive model for the Mullins effect with permanent set in particule-reinforced rubber. Int. J. Solids Struct. 41, 1855-1878.

Dorfmann, L., Pancheri, F., 2012. A constitutive model for the Mullins effect with changes in material symmetry. Int. J. Non-Linear Mech. 47, 874-887.

Gillibert, J., Brieu, M., Diani, J., 2010. Anisotropy of direction-based constitutive models for rubber-like materials. Int. J. Solid Struct. 47, 640-646.

Göktepe, S., Miehe, C., 2005. A micro-macro approach to rubber-like materials. Part III: The micro-sphere model of anisotropic Mullins-type damage. J. Mech. Phys. Solids 53, 2259-2283.

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

Holzapfel, G. A., Stadler, M., Ogden, R. W., 1999. Aspects of stress softening in filled rubbers incorporating residual strains. Constitutive models for rubber, Dorfmann & Muhr (eds), Balkema, Rotterdam.

Itskov, M., Ehret, A., Kazakeviciute-Makovska, R., Weinhold, G., 2010. A thermodynamically consistent phenomenological model of the anisotropic Mullins effect. ZAMM - J. Appl. Math. Mech. 90, 370-386.

James, H. M., Green, A., 1975. Strain energy functions of rubber. II- The characterization of filled vulcanizates. J. Appl. Polym. Sci. 19, 2319-2330.
James, H. M., Guth, E., 1943. Theory of the elastic properties of rubber. J. Chem. Phys. 11 (10), 455-481.

Johnson, M. A., Beatty, M. F., 1993. The Mullins effect in uniaxial extension and its influence on transverse vibration of rubber string. Continuum Mech. Thermodyn. 5, 83-115.

Johnson, M. A., Beatty, M. F., 1995. The Mullins effect in equibiaxial extension and its influence on the inflation of a balloon. Int. J. Eng. Sci. 33, 223-245.

Kuhn, W., Grün, F., 1942. Beziehunger zwischen elastischen konstanten und dehnungsdoppelbrechung hochelastischerstoffe. Kolloideitschrift 101, 248-271.

Le Cam, J. B., Huneau, B., Verron, E., Gornet, L., 2004. Mechanism of fatigue crack growth in carbon black filled natural rubber. Macromol. 37, 5011-5017.

Machado, G., Chagnon, G., Favier, D., 2012. Induced anisotropy by the Mullins effect in filled silicone rubber. Mech. Mater. 50, 70-80.

Machado, G., Chagnon, G., Favier, D., 2014. Theory and identification of a constitutive model of induced anisotropy by the Mullins effect. J. Mech. Phys. Solids 63, 29-39.

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.
Merckel, Y., Brieu, M., Diani, J., Caillard, J., 2012. A Mullins softening criterion for general loading conditions. J Mech. Phys Solids 60 (7), 1257 – 1264.

Merckel, Y., Diani, J., Brieu, M., Caillard, J., 2013. Constitutive modeling of the anisotropic behavior of Mullins softened filled rubbers. Mech. Mater. 57, 30–41.

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

Miehe, C., Göktepe, S., Lulei, F., 2004. A micro-macro approach to rubber-like materials - Part I: the non-affine micro-sphere model of rubber elasticity. J. Mech. Phys. Solids 52, 2617–2660.

Muhr, A. H., Gough, J., Gregory, I. H., 1999. Experimental determination of model for liquid silicone rubber: Hyperelasticity and Mullins effect. In: Proceedings of the First European Conference on Constitutive Models for Rubber. pp. 181–187.

Mullins, L., 1947. Effect of stretching on the properties of rubber. J. Rubber Res. 16, 275–289.

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

Ogden, R. W., Roxburgh, D. G., 1999. A pseudo-elastic model for the Mullins effect in filled rubber. Proc. R. Soc. Lond. A 455, 2861–2877.
Pawelski, H., 2001. Softening behaviour of elastomeric media after loading in changing directions. Constitutive models for rubber, Besdo, Schuster & Ihleman (eds), 27–34.

Rebouah, M., Chagnon, G., 2014. Permanent set and stress softening constitutive equation applied to rubber like materials and soft tissues. Acta Mech. 225,, 1685–1698.

Rebouah, M., Machado, G., Chagnon, G., Favier, D., 2013. Anisotropic Mullins softening of a deformed silicone holey plate. Mech. Res. Comm. 49, 36–43.

Robisson, A., 2000. Comportement mécanique d'un élastomère chargé en silice. Etude de l'influence des charges et modélisation par une loi visco-hyperélastique endommageable. Ph.D. thesis, Ecole Nationale Supérieure des Mines de Paris.

Shariff, M. H. B. M., 2006. An anisotropic model of the Mullins effect. J. Eng. Math. 56, 415–435.

Svistkov, A. L., Komar, L. A., 2005. Modeling of relaxation processes in filled elastomeric materials. Polym. Sci. Ser. A 47 (4), 370–375.

Treloar, L. R. G., 1944. Stress-strain data for vulcanised rubber under various types of deformation. Trans. Faraday Soc. 46, 59–70.

Valanis, K. C., Landel, R. F., 1967. The strain-energy function of a hyper-elastic material in terms of the extension ratios. J. Appl. Phys. 38 (7), 2997–3002.
Wu, P. D., Giessen, E. V. D., 1992. On improved 3-D non-gaussian models for rubber elasticity. Mech. Res. Commun. 19 (5), 427-433.