*2.4. Implementation*

The model was implemented in Matlab for the hypothetical condition of fully incompressible material, considering one material point loaded cyclically in the uniaxial direction and in a compressible full three dimensional form for the finite element analysis. For incompressible materials *J* = *I*3 = 1 and the Cauchy stress is related to the free energy density:

$$\boldsymbol{\sigma}\_{\rm el} = -p\mathbf{I} + 2\left(\frac{\partial W\_X}{\partial I\_1} + \frac{\partial W\_X}{\partial I^\*}\right)\mathbf{B} \tag{10}$$

The scalar *p* is the hydrostatic pressure (an indeterminate Lagrange multiplier) that can be determined only from the boundary conditions. **B** is the left Cauchy-Green deformation tensor. For compressible materials it is useful to split the deformation into distortional and volumetric (dilational) parts via a multiplicative decomposition.

$$\mathbf{F} = I^{1/3}\mathbf{\hat{F}}\tag{11}$$

$$\mathbf{B} = f^{2/3} \mathbf{B} \tag{12}$$

$$I\_1 = f^{2/3} I\_1 \tag{13}$$

The *J*1/3 and *J*2/3 terms are related to the volume changes, while the •¯ terms are a modified gradient and strain related to the distortional deformations. These decompositions express the stored energy as *<sup>W</sup>*(¯*I*1(*<sup>I</sup>*1, *J*), ¯*<sup>I</sup>*<sup>∗</sup>(*I*<sup>∗</sup>, *J*), *J*), that can be decomposed into volumetric and isochoric (distortional) responses (Equation (14)) and give the Cauchy stress [36] as in Equation (15)

$$\mathcal{W}(\mathbf{B}) = \mathcal{W}\_{\text{vol}}(I) + \mathcal{W}\_{\text{iso}}(\mathbf{B}). \tag{14}$$

$$\sigma\_{\rm el}(I\_1, I^\*, f) = \frac{2}{J} \left[ \frac{\partial \mathcal{W}\_X}{\partial \bar{I}\_1} + \frac{\partial \mathcal{W}\_X}{\partial \bar{I}^\*} \right] \mathbf{B} + k(J - 1). \tag{15}$$

The single Prony series element was integrated explicitly and the total stress is defined as:

$$
\sigma\_{\rm tot} = (1 - \phi)\,\sigma\_{\rm el} + \phi\,\sigma\_{\rm vis}.\tag{16}
$$

where *φ* is a fitting constant which is assumed to increase with filler volume fraction. Figure 1 shows the algorithm developed to implement the model as a user-defined material subroutine (VUMAT) coded in Fortran for Abaqus/Explicit.

The model requires 10 parameters, plus the bulk modulus in case of compressible material.

**Figure 1.** VUMAT Algorithm.

### **3. Materials and Experiments**

Seven compounds supplied by TARRC (Tun Abdul Razak Research Centre, Brickendonbury, Hertford, UK), were examined. Each compound contained a different volume fraction of carbon black (FEF N550) that had been mixed into natural rubber (NR, SMR CV60). The compound formulations in phr (part per hundred of rubber by mass) are given in Table 1.


**Table 1.** Compound formulation; phr (part per hundred of rubber by mass).

Characterisation of all this materials under cyclic tensile tests is available in [14] and it is briefly recalled in this section for the sake of completeness. Uniaxial tension tests were conducted using an Electropulse Instron Test Machine. Cyclic loading, unloading and reloading tests were performed using four different strain rates of 0.5/s, 1.5/s, 3/s and 6/s. All the specimens had a gauge length of 12 mm and a cross section of 2 mm × 0.5 mm. The nominal strain, *<sup>ε</sup>NOM*, was determined as the ratio of the axial displacement to the original length of a specimen (12 mm). This test was difficult to conduct because the large strain amplitude and rate approached the limit of the test machine. To reach the desired strain rates, the specimens had to be short. In addition, the strain rate was too fast to be reliably measured using the optical strain measuring device. The tensile force was measured using a 1 kN load cell.

A subset of the results are shown in Figures 2 and 3 for two of the materials. The results display strong nonlinearity, large hysteresis, complex Mullins behaviour (Figures 2 and 3b), cyclic stress relaxation (Figure 3a) and permanent set. Full detail of all the experimental results can be found in [14].

In the range of interest the strain rate effect is not a dominant feature, so in the following analysis it was ignored.

**Figure 2.** Stress–strain response of natural rubber filled with 20 phr of carbon black (NR20) when subjected to cyclic uniaxial tension. Different behaviour of the material after a different pre-strain is evident.

**Figure 3.** Stress–strain response of natural rubber filled with 30 phr of carbon black (NR30) submitted to cyclic uniaxial tension at 1.5/s. The solid curves are the trends for the first cycle at a given amplitude. The dotted curves are the trends after the initial cycle at given amplitude. (**a**) 5 cycles at the same amplitude, (**b**) 5 cycles at 4 different strain amplitude.
