A phenomenological constitutive model for the viscoelastic deformation of elastomers

This study proposes a one-dimensional constitutive model for elastomeric materials based on recent observations regarding the separation of elastic and viscous contributions in uniaxial cyclic tensile experiments on EPDM rubber. The focus is on capturing the changes in constitutive behavior and energy dissipation associated with the Mullins effect. In the model, this is achieved through the evolution of both permanent set and hyperelastic parameters of an Edwards-Vilgis function to account for the Mullins effect, and with a viscosity associated with the effective stretch rate of the network to describe the non-linear flow stress. The simulations are able to reproduce the observed constitutive response and its change with increasing levels of pre-deformation. The model is less able to accurately reproduce the virgin loading response, which is achieved via extrapolation to zero pre-strain. However, for practical purposes, where scragging of elastomeric products is the norm, the model is able to predict the cyclic response and the dissipated energy, and their change with different scragging levels in good agreement with experimental data.


Introduction
Elastomers are essential to a wide range of industrial applications, including tyres, dampers and seals to name a few. Constitutive models are a critical element in the design process to enable simulation of the in-service behavior of these products during their life cycle. Owing to stress softening, permanent set and viscoelasticity, amongst others, these materials exhibit a complex mechanical response. Furthermore, there is still much debate in the academic literature regarding the physical origins of several aspects of the underlying material behavior. A combination of the two factors described above has resulted in numerous modeling approaches being suggested by the research community.
A significant stress softening is noted for rubber-like materials after an initial loading. This phenomenon is referred to as the Mullins effect (Mullins 1948;Tobin 1965, 1957), and remains a major hurdle in the formulation of constitutive models, with the majority of rubber models still limited to behavior following a single prescribed preconditioning history. Several authors have attributed this phenomenon to bond rupture (Clough et al. 2016), molecular slippage (Houwink 1956), filler rupture (Kraus et al. 1966), molecular disentanglement (Hanson et al. 2005) and multiple networks (Fukahori 2005), although the cause of this phenomenon is still the subject of on-going research. Nevertheless, there are a limited number of physically motivated and phenomenological constitutive models that address the Mullins effect. Marckmann et al. (2002) implemented a network alteration model, which considers the rupture of filler-matrix and polymer chain bonds (i.e. weak interactions and cross-links). As the bonds break, the average number of monomer segments N in a polymer chain increases. The parameter N is described by a linear function of the maximum pre-deformation seen by the network and is determined empirically. As N increases, the average number of active chains per unit volume n decreases, causing a softening of the network. Mullins and Tobin (1957) suggested a two-phase model; the soft phase is the rubber matrix and the hard phase is the filler. The hard phase is initially surrounded by bound rubber from the soft phase, immobilized by means of entanglements and intermolecular forces. The softening effect arises from the release of bound rubber into the soft phase as a result of deformation and causes an increase in the effective volume of the soft phase. Models of this type were later implemented by Qi and Boyce (2004), and by Fernandes (2016). The former fits the initial and final volume fractions via optimization with respect to experimental data, whereas in the latter the effective volume of the soft phase was obtained by thermogravimetric analysis and swelling experiments.
Phenomenological descriptions of the Mullins phenomenon are broadly based on a concept introduced by Simo (1987), where strain energy density functions used to describe hyperelastic materials are multiplied by a reducing parameter. These models differ from each other in their definition of the damage variable, and they have been comprehensively reviewed by Diani et al. (2009).
Alongside stress softening, deformed elastomers exhibit a residual strain (or permanent set) on load removal. Simo's theory (Simo 1987) can be extended to account for permanent set by incorporating additional reducing parameters. An approach of this type has been followed by Dorfmann and Ogden (2004), and by Peña (2014). Maher et al. (2012) suggested an alternative approach, where an additive split of the stress tensor is used to capture the permanent set such that a negative stress is observed at a strain of zero; this is analogous to having a residual strain at a zero stress. This paper presents a newly developed constitutive model whose parameter evolution leads to a Mullins effect and permanent set. This model is inspired by a recent study by De Focatiis et al. (2009) involving cyclic tests following various levels of pre-deformation where the elastic and viscous contributions were extracted from the stress-stretch response. By isolating the elastic response, they demonstrated that the parameters of an Edwards-Vilgis (E-V) strain energy function (Edwards and Vilgis 1986) fitted to the elastic contribution evolve with pre-deformation. Once permanent set was accounted for, the viscosity obtained from the isolated viscous contribution was shown to be independent of predeformation and a unique function of network stretch. Although non-linear viscosities have been incorporated in previous constitutive models, to the best of the authors' knowledge, none have been linked to strain in elastomers.
This paper formulates a one-dimensional (1D) constitutive model based on the ideas of an evolving hyperelastic component and a non-linear viscoelastic Maxwell element. It then explores the extent to which such a model can capture the observed cyclic response and the dissipated energy of an EPDM rubber deformed in a series of complex strain histories. Although the implementation is currently only one-dimensional, the process of extending the model to 3D is discussed.

Materials and methods
The experimental data employed in this study have been described previously (De Focatiis et al. 2009), and only a brief overview of the experimental procedure is given here.

Materials and manufacturing
A carbon-black filled (50 phr) EPDM rubber was compression molded for 13 minutes at 160 • C into approximately 0.5 mm thick sheets. Uniaxial tensile test specimens were cut from the sheet using a hand operated Wallace specimen cutting press fitted with a dumbbell shaped cutter 1BA according to BS ISO 527-2.

Mechanical testing
Uniaxial tensile cyclic tests were performed on an Instron 4204 tensile testing machine fitted with a counterbalanced traveling extensometer to record the strain at a constant nominal strain rate of 0.03 s −1 at room temperature, (24 ± 1) • C. The specimens were subjected to four load-unload cycles through to a specified pre-deformation λ max , as illustrated in the inset of Fig. 1, and ten separate λ max increments between 1 and 6 inclusive were explored. To prevent buckling, the specimens were always unloaded to a stretch λ (0.1N) corresponding to a small tensile force of 0.1 N.
Representative stress-stretch responses to the cyclic loading are illustrated in Fig. 1. Many of the characteristic features of the deformation of filled elastomers can be observed, including stress softening, permanent set and hysteresis. The stress softening and permanent set increase with increasing stretch λ.

Extraction of the elastic and viscous contributions
Decompositions into elastic and viscous contributions are commonly used in constitutive models describing elastomers and polymers, as first suggested by Haward and Thackray (1968), and as implemented by Bergström and Boyce (1998), and several others. In this study, it is assumed that the changes to the constitutive behavior of the rubber due to the Mullins effect arise as a result of the maximum level of deformation experienced, and hence simply depend on λ max . As subsequent unloading-reloading loops are similar (see Fig. 1), any unloading-reloading loops following λ max can be used in the decoupling of the mechanical response. The third unload-reload loop has been chosen in this instance since the constitutive response shows the most dramatic change in the first few cycles.
The assumption is that the steady-state response at a given strain level can be decomposed as the sum of an elastic (equilibrium) σ e and a dissipative, or viscous σ v stress. σ v changes sign depending on the sign of the strain rate, positive on loading and negative on unloading. Fig. 1 The stress-stretch response of selected EPDM specimens subjected to cyclic loading as shown in the inset to maximum stretches λ max of 2, 3, 4, 5 and 6. The specimens are unloaded to a stretch corresponding to a load of 0.1 N λ (0.1N) to prevent buckling. A virgin specimen was used for each test Fig. 2 The experimental stress-stretch response (symbols) to cyclic tension to a maximum stretch λ max of 3, showing only the pre-deformation and the third unload-reload loop. The elastic response σ e (dashed line) and viscous response σ v (dot-dash line) are obtained from the mean and from half of the difference between the third loop unloading and reloading stress, respectively. The transient portion of the data (here indicated as a strain of 0.33 at each end of the deformation) and the data used in the procedure are highlighted  Prisacariu et al. (2005), σ e and σ v are computed from the experimental measurements as the mean and as half of the difference between the loading and the unloading stress, respectively. Here, to ensure that the viscoelastic transients have saturated, the initial and final 0.33 strains of each cycle are discarded. The elastic and viscous contributions obtained in this way are highlighted in Fig. 2.

Equilibrium contributions
The evolution of σ e as a function λ, with increasing λ max , is illustrated in Fig. 3. The elastic response is strongly affected by the increase in historical maximum strain, and this implies that at least some of the Mullins effect arises from an evolution of the underlying elastic network, here simply dependent on λ max . Various hyperelastic models can be used to capture the elastomer's equilibrium response arising from the entropic nature of the network. In Fig. 3 The experimental σ e (symbols) and model σ E−V (lines) equilibrium contributions as a function of stretch for five pre-deformation λ max levels between 2 and 6. The model is an Edwards-Vilgis function accounting for the permanent set λ set this instance, an Edwards-Vilgis (E-V) hyperelastic function (Edwards and Vilgis 1986) is chosen due to the physically inspired model parameters. The E-V strain energy density W function is given by where k B is the Boltzmann constant, T is the absolute temperature, N C is the cross-link density, N S is the slip-link density, η is a measure of slip-link mobility, α is a measure of finite chain inextensibility and λ i are the stretches in the principal directions. The stress in the E-V spring can be determined by differentiation of Eq. (1) with respect to the stretch direction, accounting for isochoric deformation. E-V parameters were obtained by optimization on each third loop elastic contribution, and the process was repeated for all values of λ max . The optimization algorithm was coded in Matlab using the 'lsqcurvefit' tool to obtain values of N C , N S , α and η by minimization of the error, defined here as the root-mean-square (rms) of the difference between experimental and numerical stress values within one unload-reload loop. It was found that an E-V function on its own was unable to produce satisfactory fits to the data, and in order to achieve good agreement, it was necessary to incorporate a measure of permanent set λ set , achieved via a multiplicative decomposition of λ into where λ eff is the effective stretch, i.e. the stretch actually experienced by the network, and λ set is a permanent plastic stretch arising from irreversible effects. In addition, N C was always several orders of magnitude smaller than N S , and hence insignificant as far as the total elastic stress is concerned. This is most likely due to the difficulty in separating the effects Fig. 4 (a) Viscous stress σ V as a function of stretch λ, and (b) the viscosity η V as a function of effective stretch λ eff . When expressed as a function of λ eff , the viscosity data overlay to form a master curve. A simple function is fitted to describe the viscosity master curve of cross-links and slip-links. Hence, for the purpose of modeling the elastic response, it was sufficient to express the E-V function as This combination produced excellent agreement between the elastic data and the E-V functions, as shown in Fig. 3, with rms errors typically of the order of 0.01 MPa, and always less than 0.02 MPa.

Viscous contributions
The viscous contribution to the stress σ V is illustrated in Fig. 4a for all values of λ max as a function of λ. The varying levels of λ max lead to different σ V vs. λ curves. This viscous contribution is associated with flow of the network, and, as such, a viscosity η V is defined as where σ V,true is the true flow stress andε true is the true strain rate. It was found that σ V differed with varying levels of λ max , however, once λ set obtained from the elastic optimization is accounted for, and η V is plotted as a function of λ eff , the curves overlap to form a single master curve as shown in Fig. 4b. The implication of this is that η V may be modeled as a Mullins-independent, but λ eff -dependent quantity. A plausible explanation for this effect is that molecular alignment leads to an increasingly anisotropic flow process. Similar effects have been noted in polymer glasses (De Focatiis and Buckley 2006;Senden et al. 2010).
Starting from the definition of true strain, the time derivative of Eq. (5) is calculated to determineε true and is expressed aṡ whereλ eff is the effective stretch rate andλ set is the rate of change of permanent set. The data used to determine the viscosity is obtained from the third loop. The evolution of λ set predominantly takes place during the first loading, and it is assumed thatλ set = 0 for subsequent loadings. Equation (6) is thus reduced tȯ The relation between true σ V,true and nominal σ V flow stresses can be expressed as σ V,true = λσ V , and hence, by substitution of Eq. (7) into (4), η V can be expressed as

Constitutive modeling
This section describes the implementation of a constitutive model incorporating the experimental observations stated above: (1) an additive decomposition of stress into elastic and viscous components; (2) the need for an effective stretch to account for permanent set; (3) network elasticity evolving with maximum stretch, and (4) a stretch-dependent viscosity.
In addition, the behavior within viscoelastic transients will be considered to complete the model. The proposed model, illustrated in Fig. 5, is a modified viscoelastic standard linear solid (SLS) model in series with a custom 'slider' element representing the permanent set. The modified SLS model consists of an E-V hyperelastic spring in parallel with a non-linear Maxwell element. The non-linearity of the Maxwell element arises from the non-linearity of the dashpot viscosity.

Slider
In the proposed model, λ is decomposed into λ eff and λ set via a multiplicative decomposition as described in (2). The significance of λ eff in this model is that it refers to the stretch that the E-V spring and the Maxwell element actually experience. The permanent set here represents the combined effects of irreversible deformation experienced by rubber and filler combined. It is modeled by a slider as shown on Fig. 5 whose stretch is, for the case of simple uniaxial deformation dealt with here, dependent on a monotonic function of λ max . To cover more general cases this may need to be a more complex tensorial function of time, temperature and deformation history.

Non-linear Maxwell element
The viscoelastic behavior exhibited by the elastomer arises from local bond stiffness and viscous flow of monomer segments past each other and past the filler. Here it is modeled by a non-linear Maxwell element consisting of a linear elastic spring in series with a non-linear dashpot, as shown in Fig. 5. The governing differential equation for this element iṡ where σ M is the stress in the Maxwell element,σ M is the stress rate, τ is a relaxation time, λ eff is the effective stretch rate, and E is the stiffness of the linear spring. The relaxation time τ is given by η V /λλ eff E, where η V (λ eff ) is the viscosity of the dashpot, which itself is a function of the effective stretch λ eff , i.e. the viscosity master curve.

Parameter evolution
The values and uncertainties (expressed as 95% confidence intervals) of N S , η, α and λ set obtained from the optimizations carried out on the equilibrium contributions are illustrated on Fig. 6, for the full range of experiments carried out to different λ max . At small values of λ max , there is considerable uncertainty on these parameters due to the limited strain range of equilibrium data away from transients. All parameters exhibit a dependence on λ max , suggesting that the maximum level of deformation influences all aspects of the elastic network. To describe this evolution, simple mathematical functions were fitted to the parameters as shown in Fig. 6, for values of λ max ≥ 2.5 only. Several forms of these functions were explored, and the functions selected are by no means unique, but provide simple and numerically stable representations of the parameter evolutions. The forms of the functions selected to describe the evolution of the E-V parameters and of λ set are log 10 (N S ) = C 1,N S λ max + C 2,N S (10) α = C 1,α λ 2 max + C 2,α λ max + C 3,α (11) η = C 1,η exp −C 2,η (λ max − 1) λ set = 1 + C 1,set exp C 2,set λ max − exp C 2,set (13) Fig. 6 The evolution of (a) slip-link density N S , (b) slip-link mobility η, (c) chain inextensibility α and (d) permanent set λ set with pre-deformation λ max . The parameters extracted from the elastic contributions are shown as circles. The corresponding uncertainties (95% confidence interval) are also shown. To describe the evolution of these parameters, simple functions (lines) are fitted to the extracted parameter values as functions of λ max for λ max ≥ 2.5 (larger circles) To ensure that the permanent deformation is zero for virgin specimens, the condition λ set = 1 at λ max = 1 is applied to an exponential equation of the form λ set = C 3,set + C 1,set exp C 2,set λ max , resulting in expression (13). The estimates for the coefficients C j,k fitted to the data using linear regression for the functions are reproduced in Table 1.

Viscosity master curve
To describe the viscosity master curve as illustrated in Fig. 4b, a function of the form is proposed, and parameters are obtained by minimizing the normalized error and provided in Table 1. Different approaches have been proposed to describe viscosity within a rubber constitutive model. These approaches often modify the dashpot element within a generalized Maxwell framework to capture inelastic effects, for example see the work of Rendek and Lion (2010), and Jalocha et al. (2015). Currently, to the best of the authors' knowledge, despite the experimental evidence suggesting a Mullins-independent but stretch-dependent viscosity, such an implementation has not been previously attempted for elastomers.

Determination of the linear spring modulus
The non-linear Maxwell element includes a linear elastic spring that, together with the nonlinear dashpot, describes the transient response of the viscoelastic arm of the model. In order to determine E, the response within the transient regions must be considered. One can determine E via the numerical evaluation of Eq. (9) for the case of constant effective rate of deformation, where E is the only unknown.
The viscous contribution is obtained by subtracting the extrapolated elastic contribution predicted by the E-V model from the experimental third loop response within the transient, i.e. σ − σ E−V . This is necessary because it is not possible to use the mean equilibrium experimental data within a transient. Fig. 7a illustrates an example of the third loop experimental data and the model elastic response extrapolated into the transients. The stress resulting from the subtraction process is shown in Fig. 7b.
In order to obtain a suitable value of E, the loading portion of σ M is used in a least squares algorithm against a numerical solution for Eq. (9). This procedure was implemented on all experimental data sets to obtain E = 3.07± 1.09 MPa. Fig. 7b shows an example of the experimental and model transient responses for the case λ max = 2.5. Although the agreement Fig. 8 Comparison of experimental (symbols) and model (lines) third loop response for five pre-deformation λ max levels between 2 and 6. The inset shows the deformation history imposed on the specimens and the third unload-reload loop between model and experiment is far from perfect, it should be noted that experimental data within the transients is noisy, and that for an accurate reproduction of the transients it is likely that a spectrum of relaxation times would be required.
A summary of the steps required in order to extract model parameters and their evolution is provided in Appendix A.

Simulation of the third loop response
Once all model parameters have been obtained, they can be employed in generic simulations of the response. A comparison between the experimental and model third loop responses is shown on Fig. 8 for five λ max levels ranging from 2-6. Using a single set of model parameters, some of which evolve with λ max , simulation and experimental data are generally in good agreement, with the shape of the curves and the degree of permanent set matching reasonably well for all levels of stretch.

Simulation of more complex deformation histories
To probe the capability of the model, it was subjected to complex loading conditions probing cycles within a previously defined maximum stretch. Two protocols are employed, both of which condition the material by subjecting it to three cycles through to λ = 4. The first protocol (TP1) is followed by load/unload loops, where the specimen is loaded to increasing stretch levels and always unloaded (to 0.1N). The second protocol (TP2) starts at the maximum λ of 4, and unloads to progressively smaller λ, always reloading to λ = 4. The two protocols are shown in the insets of Figs. 9 and 10.
Comparisons between the experimental and the model responses for TP1 and TP2 are shown on Figs. 9 and 10, respectively. The primary loading is also shown here, and it can be seen that the model underestimates the stress in the latter stages of the loading. Beyond the Fig. 9 A comparison between the experimental and model responses for TP1. The protocol consists of three uniaxial load/unload cycles, followed by five reload/unload loops and a reload/unload cycle, as shown in the inset. The dashed line in the inset represents the data that has been omitted from the plot Fig. 10 A comparison between the experimental and model responses for TP2. The protocol consists of two uniaxial load/unload cycles, a reload step, four unload/reload loops, and lastly an unload/reload cycles and an unload step, as shown in the inset. The dashed line in the inset represents the data that has been omitted from the plot initial loading, during the cycles, the model reproduces the stress well during the remaining deformation histories, although the transient responses saturate somewhat faster in the simulations compared to the experiments.

Post-Mullins response
The model parameters were obtained by the use of cyclic experimental data sets, each of which remains within a pre-defined value of λ max . Thus, it is perhaps not surprising that the third unload-reload loops are reproduced to a good degree of accuracy. It should be noted, however, that the model does overestimate the stress for large values of λ max . This overestimation is attributed primarily to the combined inaccuracies arising from the functions describing the evolution of the E-V parameters. It is possible that selection of more complex functions may lead to a better overall fit, but at the expense of an increase in the number of constants. What is particularly encouraging is the agreement in the energy dissipated per cycle between simulation E mod C and experiment E exp C , illustrated in Fig. 11. This is of interest in damping applications, and it shows that the model can accurately account for changes to dissipation arising as a result of changes to the constitutive behavior associated with the Mullins effect.
A further analysis of the dissipated energy can be obtained by considering the more complex deformation histories explored in Fig. 9 and Fig. 10. E mod C and E exp C are compared in Fig. 12, and again show generally a good level of agreement. The model overpredicts the energy dissipated by a small fraction, and this is attributed to the discrepancy in the transient behavior noted in Figs. 9 and 10.

Pre-Mullins response
Although the virgin loading response is never employed in obtaining the model parameters, it can be obtained from the model by extrapolation of the E-V parameters to zero strain. This allows for a comparison between the model and experimental first loading as highlighted on Figs. 9 and 10. It is apparent that there is an increasing discrepancy in the stress for λ > 2, and the model underpredicts the stress during the first loading. It is plausible that the changes occurring when the material is exceeding its previous λ max (i.e. whenλ set > 0 and model parameters are evolving) might themselves be associated with a dissipative and timedependent process. This would not be surprising since there is evidence in the literature of Mullins healing (Corby and De Focatiis 2019).
In line with the above idea, one approach to achieving a larger stress during the first loading could be to include an additional viscosity dependent onλ set . A stress defined in this way would only be present during the initial loading stage, and hence not impact the rest of the simulations.

Extension to a fully three-dimensional model
In principle, the only major obstacles to a fully three-dimensional (3D) implementation of the proposed model concern the description of permanent set, of the elastic parameter evolution, and of a suitable effective stretch scalar, all of which are possible. It is well known, however, that both processing and deformation of elastomers can lead to anisotropic behavior (Diani et al. 2009;Fernandes and De Focatiis 2015;Mullins 1949). Constitutive models ought to accommodate this induced anisotropy, which will undoubtedly affect not only the elastic part but also the viscosity and the permanent set. This would require a much more extensive experimental data set involving sequential straining in different directions to identify where simplifications might be made in order to reduce the number of material constants. Typically, the notion of material directions is used to tackle anisotropy, for example see work by Göktepe and Miehe (2005), Diani et al. (2006), Itskov et al. (2010), Merckel et al. (2013) and Rebouah et al. (2013).
It is worth noting that, in an early study by Mullins (1949), when sheets of elastomers were subjected to sequential deformations in two perpendicular directions it was found that: (1) the ratio between the permanent sets in the alternating directions remained constant, largest in the direction of stretch; and, (2) the permanent deformation remained isochoric. It is therefore plausible that a significant part of the observed anisotropy could be attributed to permanent deformation alone, and that an isotropic (but evolving) hyperelastic formulation coupled with a simple 3D representation of permanent set might suffice in a generic 3D implementation. Work is on-going in our laboratory to assemble a 3D model based on this hypothesis.

Conclusions
This study has described how a constitutive model can be assembled on the basis of observations made from the separation of elastic and viscous contributions in uniaxial cyclic tensile experiments on EPDM rubber. The model includes a description of the evolution of both permanent set and hyperelastic parameters of an Edwards-Vilgis function to account for the Mullins effect. The novelty of this model arises from the use of a viscosity master curve to describe the non-linear Mullins-independent but strain-dependent flow stress, which was obtained directly from the viscous contribution from the experimental data. The approach described also presents a straightforward route to obtaining the relevant parameters from cyclic experimental data.
The model reproduces the cyclic response of the deformation histories explored in this study to a good level of accuracy. The model is also able to accurately predict the dissipated energy per cycle and its change with increasing levels of pre-deformation. Where cycles are smaller in amplitude, and the dissipated energy is more influenced by transients, the accuracy of the predictions is reduced since the model employs a single relaxation time.
When the prediction of the initial loading is compared to experimental data, the model underestimates the stress at larger strain levels. In practical applications, where elastomeric products are often scragged (i.e., subjected to pre-deformation), this limitation is less relevant, and the focus is always on the cyclic response, and how it changes with different scragging levels, something which the model predicts well.