## Abstract

This paper models stress softening during cyclic loading and unloading of an elastomer. The paper begins by remodelling the primary loading curve to include a softening function and goes on to derive nonlinear transversely isotropic constitutive equations for the elastic response, stress relaxation, residual strain and creep of residual strain. These ideas are combined with a transversely isotropic version of the Arruda–Boyce eight-chain model to develop a constitutive relation that is capable of accurately representing the Mullins effect during cyclic stress softening for a transversely isotropic, hyperelastic material, in particular, a carbon-filled rubber vulcanizate.

## 1. Introduction

When a rubber specimen is loaded, unloaded and then reloaded, the subsequent load required to produce the same deformation is smaller than that required during primary loading. This stress-softening phenomenon is known as the Mullins effect, named after Mullins (1947), who conducted an extensive study of carbon-filled rubber vulcanizates. Diani *et al.* (2009) have written a recent review of this effect, detailing specific features associated with stress softening and providing a précis of models developed to represent this effect.

Many authors have modelled the Mullins effect since Mullins, for example, Ogden & Roxburgh (1999), Dorfmann & Ogden (2004), Diani *et al.* (2009) and Tommasi *et al.* (2006), who present an interesting micromechanical model. However, most authors model a simplified version of the Mullins effect, neglecting the following inelastic features: hysteresis, stress relaxation, residual strain and creep of residual strain.

Mullins (1947) observed experimentally that when a rubber vulcanizate sheet undergoes an equibiaxial tension, softening occurs in all three directions. The degree of softening is not the same in all three directions, and therefore anisotropic stress–strain properties are developed. We expect that any model capable of representing accurately the experimental data on stress softening would need to take this feature into consideration.

Not all inelastic features may be relevant for a particular application. Therefore, in order to develop a functional model, we require that specific parameters could be set to zero to exclude any particular inelastic feature yet still maintain the integrity of the model.

The time dependence of a rubber specimen that is cyclically stretched up to a particular value of the strain is as represented in figure 1. Initially, loading starts at point *P*_{0} at time *t*_{0} and the specimen is loaded to the particular strain at point *P*_{1} at time *t*_{1}, the material then being unloaded to zero stress at the point at time with corresponding strain , where . Further recovery, known as residual creep, then occurs at zero stress before reloading commences at time at point with strain , where . This reloading terminates at the same strain as before, but now at point *P*_{2} and time *t*_{2}. This pattern continues throughout the unloading/reloading process.

In this paper, we derive a transversely isotropic constitutive model to represent the Mullins effect for cyclic stress softening under uniaxial tension. In §2, we present the isotropic elastic model as developed by Rickaby & Scott (2012). Section 3 focuses on developing a stress-softening model for the primary loading path. Section 4 follows the work of Spencer (1984) and lays the foundations for a transversely isotropic model, which is then developed through §§5–7, where transversely isotropic models are presented for Arruda–Boyce eight-chain elasticity, stress relaxation and residual creep, respectively. In §8, we present a new transversely isotropic constitutive model and compare it with experimental data. Finally, in §9, we conclude that the present model of transverse isotropy and the introduction of stress softening on the primary loading path provide a much better fit to experimental data in comparison with the isotropic model of Rickaby & Scott (2012).

Preliminary results of the model were presented in Rickaby & Scott (2011).

## 2. Isotropic elastic response

In the reference configuration, at time *t*_{0}, a material particle is located at the position **X** with Cartesian components *X*_{1},*X*_{2},*X*_{3} relative to the orthonormal basis {**e**_{1},**e**_{2},**e**_{3}}. After deformation, at time *t*, the same particle is located at the position ** x**(

**X**,

*t*) with components

*x*

_{1},

*x*

_{2},

*x*

_{3}relative to the same orthonormal basis {

**e**

_{1},

**e**

_{2},

**e**

_{3}}. The deformation gradient is defined by An isochoric uniaxial strain is taken in the form 2.1The right Cauchy–Green strain tensor

**C**=

**F**

^{T}

**F**is given by and has principal invariants 2.2the last being a consequence of isochoricity.

An incompressible isotropic hyperelastic material possesses a strain energy function *W*(*I*_{1},*I*_{2}) in terms of which the Cauchy stress is given by
2.3where *p* is an arbitrary pressure, *I* is the unit tensor and **B**=**F****F**^{T} is the left Cauchy–Green strain tensor. We are concerned here only with uniaxial tension in the 1-direction and so may fix the value of *p* by the requirement
Using this value of *p* in equation (2.3) then gives the only non-zero component of stress to be the uniaxial tension
2.4in which equation (2.2)_{1} has been used. The uniaxial tension (2.4) vanishes in the reference configuration, where *λ*=1.

The Arruda & Boyce (1993) eight-chain model was developed to model nonlinear isotropic rubber elasticity by considering the properties of the polymer chains of which rubber is composed. It is characterized by the strain energy function
2.5where
2.6Here, *μ* is the shear modulus and *N* is the number of links forming a single polymer chain. The Langevin function is defined by
with inverse denoted by . Upon substituting for *W* from equation (2.5) into equation (2.3), we obtain the elastic stress in the Arruda–Boyce model,
2.7

## 3. Stress softening

### (a) Softening on the unloading and reloading paths

In order to model stress softening on the unloading and reloading paths, i.e. paths *B* and *C* of figure 1, Rickaby & Scott (2012) introduced the following softening function:
3.1where is the maximum strain energy that is achieved on the primary loading path at the maximum stretch before unloading commences. *W* is the strain energy value at the intermediate stretch *λ*, so that when . The quantities *b*_{ω}, *r*_{ω} are positive dimensionless material constants where

The softening function (3.1) has the property that
3.2thus providing a relationship between the Cauchy stress **T** in the unloading and reloading of the material, after the primary loading has ceased, and the Cauchy stress in the primary loading phase of an isotropic elastic parent material. In the model of Rickaby & Scott (2012), the stress is a purely elastic response, not subject to any softening. Softening functions are discussed in more detail by Rickaby & Scott (2012) and Dorfmann & Ogden (2003, 2004). None of these models consider the possibility of softening on the primary loading path, i.e. path *A* of figure 1.

### (b) Softening on the primary loading path

Mullins (1969) observed that in filled rubber vulcanizates, pronounced softening occurs during primary loading, but only at very small deformations. He conjectured that this was due to the breakdown of clusters of filler particles. For intermediate and large deformations, however, there was no noticeable softening during primary loading.

This property of softening on the primary loading path, that is, path *A* of figure 1, may be included within the model developed here by introducing a new softening function, *ζ*_{0}(*λ*), which is similar in form to the softening functions (3.1) previously defined, except that in (3.1) is replaced by . We therefore define *ζ*_{0}(*λ*) by
3.3where *r*_{0}, *b*_{0} and *ϑ*_{0} are positive constants. We model the empirical fact that softening on the primary loading path occurs only for small strains by requiring *ζ*_{0}(*λ*) to be close to 1 when *λ* is close to .

Upon combining *ζ*_{0}(*λ*) with equation (2.7) for an incompressible, isotropic material, the initial primary loading path can be modelled by
3.4Eliminating the pressure *p* by the requirement that , as before, gives the uniaxial tension
3.5

As *ζ*_{0}(*λ*) has been modelled to be very close to 1 on the primary loading path when *λ* is close to , it will not affect the modelling of any subsequent primary loading, or unloading and reloading, paths.

The inclusion of a softening function on the primary loading path has not previously been discussed in the literature in relation to the Mullins effect.

## 4. Transversely isotropic elastic response

Dorfmann & Ogden (2004) found that after a uniaxial deformation such as (2.1), there is an induced change in the material symmetry because some of the damage caused by the stretch is irreversible. The material symmetry therefore changes from being fully isotropic to being transversely isotropic, with preferred direction in the direction of uniaxial stretch. This change of material symmetry influences all of the subsequent response of the material. Horgan *et al.* (2004) conjectured that if loading is terminated at the stretch on the primary loading path, then the damage caused is dependent on the value of and that this must be reflected in the subsequent response of the material upon unloading and reloading. They too concluded that the material response must become transversely isotropic.

Diani *et al.* (2006) have also experimentally observed the transition from an isotropic to an anisotropic material for carbon-filled elastomers under uniaxial testing. Strain-induced anisotropy has been studied by other authors, including Park & Hamed (2000) and Dorfmann & Pancheri (2012).

Spencer (1984) characterized a transversely isotropic elastic solid by the existence of a single preferred direction, denoted by the unit vector field **A**(**X**). After deformation, the preferred direction **A**(**X**) becomes parallel to
which is not, in general, a unit vector.

The strain energy function *W* in a transversely isotropic material is a function of five invariants, namely, the three defined by (2.2) and the further two defined by
4.1

For an incompressible material, *I*_{3}=1, and so the strain energy takes the form
The elastic stress in an incompressible transversely isotropic elastic material is then given by
4.2where *p* is an arbitrary pressure and ⊗ denotes a dyadic product. Equation (4.2) is equivalent to that presented by Spencer (1984, eqn (67)).

The preferred direction **A** lies in the direction of the uniaxial tension (2.1), so that, in components,
4.3From equations (4.1) and (4.3), the invariants *I*_{4} and *I*_{5} are given by
The stress (4.2) reduces to

For the stress to vanish in the reference configuration, we require 4.4

The invariant *I*_{4}=** a**⋅

**is clearly a measure of stretch in the direction of transverse isotropy. Merodio & Ogden (2005) observed that the invariant**

*a**I*

_{5}is more connected with shear stresses acting normally to the preferred direction. For the uniaxial deformation (2.1), there are no such shear stresses and so in our elastic model, we take the strain energy to be independent of

*I*

_{5}. A strain energy function that is independent of

*I*

_{5}, vanishes in the reference configuration, and satisfies the derivative condition (4.4)

_{1}, is given by 4.5where

*s*

_{1}and

*s*

_{2}are constants. The transversely isotropic strain energy function

*W*

_{ti}above is found in §8 to fit the experimental data extremely well.

## 5. The eight-chain model in transverse isotropy

We follow Kuhl *et al.* (2005) in developing a model of transversely isotropic elasticity based on the original (Arruda & Boyce 1993) eight-chain model of isotropic elasticity. Rubber is regarded as being composed of cross-linked polymer chains, each chain consisting of *N* links, with each link being of length *l*. We introduce the two lengths
5.1The locking length *r*_{L} is the length of the polymer chain when fully extended. The chain vector length *r*_{0} is the distance between the two ends of the chain in the undeformed configuration. Because of significant coiling of the polymer chains, this length is considerably less than the locking length. The value is derived by statistical considerations.

In this extension of the Arruda–Boyce model, we consider a cuboid aligned with its edges parallel to the coordinate axes, as in figure 2*a*. The edges parallel to the *x*_{1}-axis, which is the preferred direction of transverse isotropy, have length *a*, and the remaining edges all have length *b*. Each of the eight vertices of the cuboid is attached to the centre point of the cuboid by a polymer chain, as depicted in figure 2*a*. Each of these eight chains is of the same length, which we take to be the vector chain length *r*_{0}.

Using figure 2*a* and equation (5.1)_{2}, we see that the chain vector length may be written as
5.2where *α*=*b*/*a* is the aspect ratio of the cuboid, with *α*=1 corresponding to material isotropy. The special case *α*=0 is illustrated in figure 2*b* and corresponds to *b*=0. Kuhl *et al.* (2005) discuss this special case and regard it as representing unidirectional fibre reinforcement.

Now suppose that the material undergoes triaxial extension in directions parallel to the cuboid edges, so that the new dimensions of the cuboid are (*aλ*_{1},*bλ*_{2},*bλ*_{3}). Each of the eight chains has the same length, and this new length is given by
which can be rewritten in terms of the invariants *I*_{1} and *I*_{4} as
5.3

The argument of the inverse Langevin function is, as in Arruda & Boyce (1993),
where *r*_{L} is given by equation (5.1)_{1}. We have, using also equations (5.1)_{2} and (5.2),
The quantities *γ* and *β* are defined by
5.4where, as before, *α*=*b*/*a* is the aspect ratio of the cuboid in this extension of the Arruda–Boyce model. We recall that *α*=1 corresponds to material isotropy so that *I*_{4} then cancels out of equation (5.4), reducing it to
5.5which is consistent with equation (2.6) of the isotropic Arruda–Boyce model.

By substituting equations (5.4) into the strain energy (2.5) of the isotropic Arruda–Boyce model, we obtain the following expression for the strain energy in the transversely isotropic Arruda–Boyce model:
5.6where *h*_{4} is a constant chosen so that the stress vanishes in the undeformed state, i.e. chosen so that equation (4.4)_{1} is satisfied.

The stress resulting from equation (5.6) is
For this stress to vanish in the reference configuration, where *I*_{1}=3,*I*_{4}=1 and , we must take
For an isotropic material, *α*=1 and we find that *h*_{4}=0, as expected.

The total strain energy is obtained by adding contributions from (5.6) and (4.5), From (4.2), the total elastic stress in the transversely isotropic Arruda–Boyce model is 5.7

## 6. Stress relaxation in transverse isotropy

Bernstein *et al.* (1963) developed a model for nonlinear stress relaxation that has been found to represent accurately experimental data for stress relaxation (see Tanner (1988) and references therein).

For a transversely isotropic incompressible viscoelastic solid, we can build on the work of Lockett (1972, pp. 114116) and Wineman (2009, §12) to write down the following version of the Bernstein *et al.* (1963) model for the relaxation stress in transverse isotropy:
6.1for *t*>*t*_{0}. The first line of (6.1) is that derived by Lockett (1972, pp. 114116) for full isotropy.

Eliminating the pressure *p* from equation (6.1) by the requirement that gives the uniaxial tension
6.2with vanishing for *t*≤*t*_{0}. In (6.2), *A*_{0} is a material constant and *A* _{l}(*t*), where *l*∈{1,2,4,5}, are material functions that vanish for *t*≤*t*_{0} and are continuous for all *t*.

Figure 1 represents a cyclically loaded and unloaded rubber specimen with stress relaxation occurring along path *P*_{0}*P*_{1}, from the point *P*_{0} at time *t*_{0} up to the point *P*_{1} where , which is reached at time *t*_{1}. Stress relaxation then follows the unloading path down to the position of zero stress, which is reached at time *t**_{1} and stretch . Stress relaxation continues at zero stress and decreasing strain along the path , whose end-point is reached at time and stretch . On the reloading path stress-relaxation proceeds until point *P*_{2} is reached, at time *t*_{2}, where once again . This pattern then continues throughout the unloading and reloading process.

For cyclic stress relaxation, the material functions *A* _{l}(*t*) are replaced by
6.3where *A*_{l}(*t*), *l*∈{1,2,4,5}, are continuous material functions that vanish for *t*≤*t*_{0}. In equation (6.3), *ϕ*_{0}, *ϕ*_{1} and *ϕ*_{2} are continuous functions of time. For simplicity, on the unloading paths and on the stress-free paths, we employ the same function *ϕ*_{1} as the argument for *A*_{l}(*t*).

In the present model, we assume that stress relaxation commences from the point of initial loading at time *t*_{0}. Stress relaxation may proceed at different rates in unloading and reloading. This is governed by the functions *ϕ*_{1} and *ϕ*_{2} in equation (6.3). Separate functions *ϕ*_{1} and *ϕ*_{2} are needed for the unloading and reloading phases, respectively, in order to model better the experimental data in §8.

Employing equation (6.3), from equation (6.1), we can derive the following isotropic and transversely isotropic relaxation stresses, respectively:
6.4and
6.5for *t*>*t*_{0}.

The total stress for a transversely isotropic relaxing material is then given by 6.6where , , and are the stresses (3.4), (5.7), (6.4) and (6.5), respectively.

The total stress (6.6) falls to zero in *t*>*t*_{0} and so we must have for *t*>*t*_{0}, implying that for *λ*>1. Each of the quantities *A*_{0}, *A*_{1}(*t*), *A*_{2}(*t*), *A*_{4}(*t*), *A*_{5}(*t*) occurring in equation (6.5) has positive coefficient for *λ*>1 and so at least one of them must be negative to maintain the requirement for *λ*>1.

## 7. Creep of residual strain in transverse isotropy

We postulate that the residual strain that is apparent after a loading and unloading cycle is caused by creep during that cycle and any previous cycles. The creep of residual strain may proceed at different rates in unloading and reloading. We consider that the creep of residual strain does not operate during primary loading.

Following on from the work of Bergström & Boyce (1998), Rickaby & Scott (2012) showed that during cyclic unloading and reloading, the creep causing residual strain can be modelled, in the case of isotropy, as a stress of the form
7.1for *t*>*t*_{1} and *λ*>1 with vanishing for *t*≤*t*_{1}. Here, in the case of isotropy. In equation (7.1), *a*_{1} and *d*_{ω} are material constants, *d*_{1} for unloading and *d*_{2} for reloading, with *d*_{2}≤*d*_{1}. The function *a*(*t*) is defined by
7.2where *Φ*_{1} and *Φ*_{2} are continuous functions of time for the unloading and reloading phases, respectively. For simplicity, on the stress-free paths, we also employ *Φ*_{1} as the argument for *a*(*t*).

The polymer chain length extension ratio is denoted by *λ*_{chain} and defined by
7.3in which equations (5.2)–(5.4) have been used. In the isotropic case, *α*=1, equation (7.3) reduces to , as already recorded after equation (7.1).

The transversely isotropic residual strain model is obtained by substituting *λ*_{chain} from equation (7.3) into equation (7.1) to give the creep stress
7.4for *t*>*t*_{1} and *λ*>1.

The total stress for a transversely isotropic relaxing material is then given by 7.5in which, for notational convenience, we have defined stresses and where , , , and are given by equations (3.4), (5.7), (6.4), (6.5) and (7.4), respectively.

## 8. Constitutive model and comparison with experiment

### (a) Constitutive model

On substituting the individual stresses given by equations (5.7), (6.5) and (7.4) into equation (7.5), we obtain the following constitutive model for the transversely isotropic material: 8.1

### (b) Comparison with experimental data

In modelling the Mullins effect, we have used the Biot stress **T**_{B}, defined by
in order to compare our theoretical results with experiment.

Figures 3 and 4 provide a comparison of the constitutive model we have developed with experimental data, which came courtesy of Dorfmann & Ogden (2004) and was presented in their paper.

In comparing our model with experimental data, we have employed the Heaviside step function *H*(*t*) defined by
and we have approximated the inverse Langevin function by its Padé approximant derived by Cohen (1991), namely,
8.2This is a very good approximation, even close to the singularity at *x*=1, see Rickaby & Scott (2012, fig. 8). The close relation of the simple model of Gent (1996) to the approximate equation (8.2) has been made clear by Horgan & Saccomandi (2002).

Figure 3 has been obtained by applying the following constants and functions: We see in figure 3 that the transversely isotropic model provides a good fit with experimental data and is a significant improvement on the isotropic model of Rickaby & Scott (2012, fig. 15). Figure 4 has been obtained by applying the following constants and functions: Figure 4 shows that the transversely isotropic model provides a good fit with experimental data and is a significant improvement on the isotropic model of Rickaby & Scott (2012, fig. 16). The better fit here is partly due to the modelling of stress softening on the primary loading path.

## 9. Conclusions

This model appears to be the first appearance in the literature of a transversely isotropic stress softening and residual strain model which has been combined with a transversely isotropic version of the Arruda–Boyce eight-chain constitutive model of elasticity in order to develop a model that is capable of accurately representing the Mullins effect in uniaxial tension when compared with experimental data.

The model has been developed in such a way that any of the salient inelastic features could be excluded and the integrity of the model would still be maintained. The proposed model should prove extremely effective in the modelling of many practical applications.

Figures 3 and 4 provide a comparison between the experimental data of Dorfmann & Ogden (2004) and the transversely isotropic model presented here. By comparing figure 3 and fig. 15 of Rickaby & Scott (2012), it can be seen that the present transversely isotropic model provides a much better fit to the data than does the original isotropic model of Rickaby & Scott (2012). Similarly, for the higher concentration of carbon black of figure 4, comparing with fig. 16 of Rickaby & Scott (2012) shows that the transversely isotropic model provides an equally better fit. This is partly due to the new feature of the modelling of stress softening on the primary loading path.

After an applied uniaxial deformation, the induced transverse isotropy means that the directions perpendicular and parallel to the deformation have sustained different degrees of damage, which is borne out by the experimental data of Diani *et al.* (2006, fig. 3). Thus, while some of the material parameters of our model remain the same along the direction of uniaxial tension and perpendicular to it, the anisotropic terms, in particular, the stress-relaxation terms, may not necessarily be equal for two perpendicular directions.

Our present model can be modified to cope with multiple stress/strain cycles, with increasing values of maximum stretch, and we hope to present a comparison between theory and experiment at a later date.

The version of the model developed here is for uniaxial tension. We expect that the results presented here could be extended to include equibiaxial tension, pure or simple shear, and a general three-dimensional model. These ideas will be developed in later papers.

A further application of the model could be to the mechanics of soft biological tissue. The comparison between the stress softening associated with soft biological tissue, in particular muscle, and filled vulcanizated rubber has been discussed in detail by Dorfmann *et al.* (2007, §2). For cyclic stress softening, both materials exhibit stress relaxation, hysteresis, creep and creep of residual strain. Similar observations have been made for arterial material, see, for example, Holzapfel *et al.* (2000). After preliminary investigations, it is apparent that the model presented here could be extended to include biological soft tissue, though the inherent anisotropy of biological tissue would need to be taken into consideration.

## Acknowledgements

One of us (S.R.R.) is grateful to the University of East Anglia for the award of a PhD studentship. The authors thank Prof. Luis Dorfmann for most kindly supplying experimental data. Furthermore, we thank the reviewers for their constructive comments and suggestions.

- Received August 2, 2012.
- Accepted September 4, 2012.

- This journal is © 2012 The Royal Society