## Abstract

A great variety of models can describe the nonlinear response of rubber to uniaxial tension. Yet an in-depth understanding of the successive stages of large extension is still lacking. We show that the response can be broken down in three steps, which we delineate by relying on a simple formatting of the data, the so-called Mooney plot transform. First, the small-to-moderate regime, where the polymeric chains unfold easily and the Mooney plot is almost linear. Second, the strain-hardening regime, where blobs of bundled chains unfold to stiffen the response in correspondence to the ‘upturn’ of the Mooney plot. Third, the limiting-chain regime, with a sharp stiffening occurring as the chains extend towards their limit. We provide strain-energy functions with terms accounting for each stage that (i) give an accurate local and then global fitting of the data; (ii) are consistent with weak nonlinear elasticity theory and (iii) can be interpreted in the framework of statistical mechanics. We apply our method to Treloar's classical experimental data and also to some more recent data. Our method not only provides models that describe the experimental data with a very low quantitative relative error, but also shows that the theory of nonlinear elasticity is much more robust that seemed at first sight.

## 1. Introduction

The mechanics of rubber-like solids has a long and prolific history. Following World War II, a huge research effort was launched to find an explicit strain-energy function able to describe accurately the experimental data obtained from the testing of natural and synthetic rubbers. However, in spite of decades of intensive work in that area, to this day there is still no effective model able to perform this task in a satisfying and universal way.

This state of affairs is a plain fact, which cannot to be hidden by the countless and seemingly successful models and simulations to be found in the literature. These simulations may be concretely descriptive but in the end, they apply only to some special phenomena. From the point of view of physical sciences, constitutive models must be *universal*, not in the sense that a single model should describe the mechanical behaviour of all elastomers, but in the sense that, for a given soft material (e.g. a given sample of natural rubber), a given model should describe its mechanical response in a satisfactory manner for all deformation fields and all stretch ranges physically attainable. Here, a *satisfactory* model is defined as a model able to describe the experimental data first of all from a qualitative point of view and then from a quantitative point with acceptable relative errors of prediction with respect to the data. In the end, it must be kept in mind that a simulation in computational solid mechanics, irrespective of its degree of sophistication, is as good, *or as bad*, as the model it relies upon.

Over the years various mathematical models have been proposed and are now listed in the literature reviews, e.g. Boyce & Arruda [1], Marckmann & Verron [2], Beda [3]. By scanning all these constitutive models, we can identify *three fundamental breakthroughs*. First, the Mooney–Rivlin strain-energy density [4]: a purely phenomenological theory stemming from the early tremendous effort devoted to rewrite the theory of Continuum Mechanics using the language of Tensor Algebra. The Mooney–Rivlin model led to the exploration of the nonlinear theory of elasticity in deep and unexpected ways, yielded significative classes of non-homogeneous exact solutions and provided a new perspective to the interpretation of experimental data. The second breakthrough has been the Ogden strain-energy density function [5]: a rational re-elaboration of the Valanis–Landel hypothesis. For the first time, it became possible to fit accurately theoretical stress–strain curves to experimental data for a variety of deformations and a large range of strains. The third breakthrough is more complex to describe: it consists in the recent re-elaboration of the ideas underpinning the classical derivation of the neo-Hookean strain energy based on the basic tools of *statistical mechanics*. Here there are two possible approaches. One is based on micro-mechanical considerations, see e.g. [6] for a recent exploration in this direction and the review by Puglisi & Saccomandi [7]. The other is based on molecular considerations, see the paper by Rubinstein & Panyukov [8] on the elasticity of polymer networks for a survey. Below we summarize the micro-mechanical multi-scale approach.

The basic assumption used to derive, from microscopic considerations, the usual models of the mechanical behaviour of biological and polymeric networks is the *entropic nature* of their elasticity [9]. Because the mechanical response of these materials is due to the deformation of the individual chains or filaments composing the network, there is a strict relationship between the conformations of these macromolecules and the mechanical response of the full macroscopic network. This situation allows for a simple and direct method to determine the macroscopic strain energy starting from simple mesoscopic considerations. Hence the nonlinear force–deformation relationship of an ideal chain is easily obtained by considering that the chain's free energy is purely entropic; then the passage from the single chain model to the full network model is achieved by means of some *phenomenological* average procedures [9].

Indeed, when we model real macromolecules we rely on several idealized assumptions. For example, when we assume that there are no interactions between the monomers composing the molecule, we are considering in fact the mathematical model of an *ideal chain*. It is possible in principle to compute in a careful and detailed way the conformations of such chains in space [10], but these computations are cumbersome and usually very complex from a mathematical point of view. For this reason, further ad hoc approximations for describing the *end-to-end-distance* of a chain are introduced [11]. The most common of these approximations is given by the *Gaussian Chain Model* [9]. The output of this approximation is a linear relationship between the applied force magnitude, *f*, and the average distance between the chain ends, 〈*R*〉, along the direction of the applied force. This approximation is clearly valid in the small force limit, but the nonlinear regime calls for more sophisticated approximations. Among the various possibilities, two such nonlinear models are very popular: the *Freely Jointed Chain* model (FJC) for polymers and the *Worm-Like Chain* model (WLC) for stiff biological molecules [11]. Both the FJC and the WLC models introduce the concept of *contour length*. The contour length of a polymer chain is its length *R*_{max} at the maximum physically possible extension. Because only one configuration can be associated to the maximum extension of the chain in the entropic theory of elasticity, we must have *R*〉→*R*_{max}. This requirement cannot be captured by the Gaussian chain model. The main difference between the FJC and the WLC models lies in the divergence behaviour of the force as 〈*R*〉→*R*_{max}. For the FJC model, we have a first-order singularity: *f*∼(〈*R*〉−*R*_{max})^{−1}, whereas for the WLC model, we have a stronger divergence: *f*∼(〈*R*〉−*R*_{max})^{−2}.

In between the Gaussian Chain model (early stretch regime) and the FJC or WLC models (late stretch regime), we find a stage of *strain hardening* in tension. For instance Toan & Thirumalai [12] recently summarized this situation for polymers in a good solvent, where the polymeric chains decompose into a succession of independent *blobs*, so named by Pincus [13] who proposed a scaling law *f*∼〈*R*〉^{3/2}. We emphasize that in general the behaviour of the polymeric chain is different from the specific one considered by Pincus but that still, a strong strain-hardening effect manifests itself when we go from the small deformation regime to the moderate deformation regime.

Now consider the most basic, but still open, problem of rubber-like mechanics: *a meaningful fitting of experimental simple tension data.* The demand for a model able to describe the entire range of attainable uniaxial tension data is manifest in technological applications of rubber-like materials. There, finite-element computations play a fundamental role and almost all commercial codes rely on the tangent modulus method. To be effective, this method needs a model able to capture the *entire* data gleaned from a uniaxial tension test. We propose to rely on the ideas exposed above for a multiscale approach, by following a rational (and reasonable) procedure in the framework of the nonlinear theory of elasticity.

In the next section, we recall the tenets of nonlinear elasticity theory, with a focus on modelling the mechanical response of an isotropic incompressible hyperelastic material such as rubber. We see that the strain-energy density *W* is a function of two variables only, the first two principal invariants of strain and that a reasonable mathematical candidate for *W* must at least be consistent with fourth-order weak nonlinear elasticity. For uniaxial data, an historical formatting—the Mooney transform—reveals three distinct regimes of stretch: small-to-moderate, strain-hardening, limiting-chain, which we delineate rigorously. We also show that the fitting procedure should rely on minimizing the relative errors as opposed to the classical (absolute) residuals, because it provides the most consistent fitting across all regimes and all measures of stress and strain. In §3, we present two sets of data on the uniaxial extension of rubber: the canonical 1944 data of Treloar (see [9]) and a more recent set by Dobrynin & Carrillo [14]. Then we proceed to model each regime of stretch in turn. For the small-to-moderate regime, we show that *W* must be a function of the two invariants, which rules out the whole class of generalized neo-Hookean materials, for which *W*=*W*(*I*_{1}) only, including the Yeoh model. We study three representative models (Mooney–Rivlin, Gent–Thomas, Carroll) giving excellent fits in that regime. The conclusion here is that the specific choice of a model is not that important for the small-to-moderate regime, as long as it is not a generalized neo-Hookean material. For the next regime, which corresponds to an upturn in the Mooney plot, we add a strain-hardening term to those models and again the fit is excellent, now across the resulting wider range of data. For the range of extreme stretches, we show how to determine the order of singularity, i.e. how to find out whether the rubber stiffens according to the WLC or the FJC model. Finally, we provide models that fit the data over the entire range of experimental data, with errors below 4%. The conclusion is that the term used to capture the asymptotic limiting-chain behaviour also captures the earlier strain-hardening regime and that the number of fitting parameters can be kept to three only (figure 1). In §4, we confirm the robustness of our methodical modelling based on simple tension by looking at its predictions for two other types of deformation: simple shear and inhomogeneous rectilinear shear. In this paper, we focus on fitting models using uniaxial data only, because it is the most basic problem faced in the applications of the nonlinear theory of elasticity to rubber elasticity, and the most available experimental data. Only once this problem has been tackled and solved can we attack the question of fitting models for rubber using other types of data, for equi-biaxial or multiaxial tension, simple shear, pure shear, inflation, bending, torsion, etc. Section 5 is a recapitulation of the results and a reflection on their consequences.

## 2. Material modelling

Here we present some basic results of the axiomatic theory of continuum mechanics for hyperelastic materials and determine a consistent measure of the goodness of fit for simple tension data.

### (a) Choice of material parameters

We consider incompressible isotropic materials which are hyperelastic. For these it is possible to define a strain-energy density *W*=*W*(*I*_{1},*I*_{2}), where ** B**=

*FF*^{T}, where

**is the gradient of deformation and the**

*F**λ*

_{i}are the principal stretches of the deformation (the square roots of the eigenvalues of

**.) Another representation of the principal invariants is**

*B*This form of the strain-energy density leads to the usual representation formula for the Cauchy stress tensor ** T**, due to Rivlin [17]

*p*is the Lagrange multiplier associated with the incompressibility constraint (

*W*

_{i}=∂

*W*/∂

*I*

_{i}, (

*i*=1,2).

Materials with strain-energy density functions such that *W*=*W*(*I*_{1}) only are denoted *generalized neo-Hookean materials*. A prime member of this class is the so-called *neo-Hookean material*
*μ*_{0}>0 is the infinitesimal shear modulus.

Now, establishing which restrictions should be imposed on the strain-energy function to guarantee reasonable physical behaviour is known as Truesdell's *Hauptproblem* of the elasticity theory [18]. In this connection, we recall the point of view of Ball & James [19]:

*At the end of the day, perhaps it would have been realized that Hadamard notions of well-posedness are far too restrictive in the nonlinear setting, that non-uniqueness and even non-existence comprise acceptable behaviour, and that there are probably no fundamental restrictions on the strain-energy function at all besides those arising from material symmetry and frame-indifference.*

So that we may not expect much from mathematical requirements. We must thus resort to *mechanical arguments*.

First of all, we point out that we dealt with frame-indifference and material symmetry when we chose to write *W* as a function of the principal invariants (2.1) [17]. Next, we require as a normalization condition that the strain energy be zero in the reference configuration: *W*(3,3)=0. Also, for incompressible materials, the requirement that the stress is zero in the reference configuration is always satisfied by a suitable choice of the value of *p* in that configuration.

The next fundamental requirement is linked to the notion of the *generalized shear modulus function* *μ*, defined as
*infinitesimal shear modulus* *μ* as
*μ*_{0}>0.

A further basic requirement that we may impose on *W* comes from the *weakly nonlinear theory of elasticity*, where the strain energy is expanded in powers of the strain. For instance, by choosing the Green–Lagrange tensor ** E**=(

*F*^{T}

**−**

*F***)/2 as a measure of strain, the most general fourth-order expansion of incompressible elasticity [20–22] can be written in the form**

*I**μ*

_{0}is the infinitesimal shear modulus (the only parameter of the linear theory),

*A*is the Landau third-order constant and

*D*is the fourth-order elastic constant. It is now known that in order to capture fully nonlinear effects in solids, the fourth-order theory is the

*minimal*model to consider. For example, in physical acoustics [23], it is necessary to carry out the expansion of

*W*to fourth order to ensure that nonlinear effects in shear waves emerge [21,24]; in the elastic bending of blocks made of incompressible soft solids, the onset of nonlinearity involves third- as well as fourth-order constants [25]; etc. In conclusion, we will impose a consistent compatibility with the fourth-order weakly nonlinear theory, instead of compatibility with linear theory only.

### (b) Fitting to uniaxial data

From now on we concentrate on the data given by the (idealized) homogeneous deformation of *uniaxial extension* resulting from the application of a unidirectional tension. We call *λ* the stretch in that direction. It is a simple matter to compute the first and second invariants of the Cauchy–Green deformation tensors, as
*tensile Cauchy stress component* [5] *t*=*T*_{11}, as *t*=*λ*∂*W*/∂*λ*. In experiments, the current force applied per unit length *f*_{1} is recorded, and by dividing it by the reference cross-sectional area of the sample, we arrive at the *engineering tensile stress* *σ*=*λ*^{−1}*t*. It is given by
*λ*−*σ* curve for *λ*≥1, we can perform a curve-fitting exercise for a candidate strain-energy density *W* and access the values of its derivatives.

Alternatively, we could plot the *t*−*λ* curve to access the derivatives of *W*.

Historically, Rivlin [17] divided the relationship (2.8) across by 2(*λ*−*λ*^{−2}) and plotted the curve obtained by having the resulting left-hand side as the vertical axis coordinate and *λ*^{−1} as the horizontal coordinate. This modified version of (2.8) is the so-called *Mooney plot*, given by the transform
*z*−*g* curve for the range 0<*z*≤1 to access the values of the derivatives of *W*. The practical effect of the Mooney transform is to map the range 1≤*λ*≤2 (where the measurements are the most reliable) over 50% of the *z*-domain.

Now, let us consider the experimental data of stretches and engineering stresses (*λ*_{i},*σ*_{i}) for *i*=1,…,*m* in the *engineering space* *λ*_{i},*t*_{i}) for *i*=1,…,*m* in the *Cauchy space* *Mooney space* *i*=1,…,*m*, where *m* is the number of measurements in the uniaxial test.

At this juncture, we emphasize that different fitting procedures are possible and that the choice of a given procedure may have an impact on our modelling considerations.

### (c) Goodness of fit

In general, let ** p**=[

*p*

_{1},…,

*p*

_{n}] be the set of material parameters involved in the definition of the strain energy

*W*, that have to be identified by

*matching*, as best as we can, the data with the desired mathematical model for

*W*. For example, for the neo-Hookean solid (2.3) we have

**=[**

*p**μ*

_{0}] and for the fourth-order elasticity model (2.6) we have

**=[**

*p**μ*

_{0},

*A*,

*D*]. To evaluate the best-fit parameters

**we solve a least squares (LS) problem, which can be linear or nonlinear depending on the functional dependence of the strain-energy densities**

*p**W*from

**.**

*p*To quantify the goodness of fit, we define on the one hand the following (*absolute*) *residuals* in the spaces *i*=1,…,*m*), respectively, and on the other hand, the following *relative residuals*,
*relative residuals are the same in the Engineering, Cauchy and Mooney spaces.* Moreover, (2.11) shows that the *relative errors* are non-dimensional quantities and can be expressed in percentage, as opposed to the absolute residuals. Then in the classical LS setting, the strategy is to minimize either of the following Euclidean norms ** p**,

In this paper, we make the choice of *minimizing the relative errors* in the two norm, that is

We shall leave aside minimization of classical absolute residuals from now on because we found that it gives large relative errors in the small-to-moderate range of extension, which is not desirable from experimental and modelling points of view. Moreover, the best-fit parameters obtained by minimizing absolute residuals are not the same when we use the engineering stress data, the Cauchy stress data, or the Mooney plot data, which is a problem as the parameters in a strain-energy density should be independent of the choice of stress measure.

To quantify the goodness of our fits, we will record the maximal relative error

## 3. Numerical results

### (a) Experimental data

To test our models, we will consider two comprehensive sets of data recorded for the uniaxial extension of rubber samples. The first set is due to Treloar [9], dating back to his canonical 1944 experiments. For that set we use the original tabulated data, with the engineering stress measured in Newton per square millimetre units, see raw data consisting of 24 points below. The second set is due to Becker & Kruger [26], as collected by Dobrynin & Carrillo [14]: we will use the data collected in their set labelled ‘DC9’; it consists of 48 data points.

Plotting the data in the Mooney space *three regimes of extension* for the elastomers (figure 2). Starting from *z*=*λ*^{−1} at value 1 (no extension) and decreasing (extension), we see at first a linear decrease of *g* with *z* for *small-to-moderate deformations*. Then, around a value for *z* of 0.3 (stretch of 230% or so), an *upturn* occurs as *g* goes through a minimum. Finally, we enter the *large deformation regime*, where the reduced tensile stress rises sharply and considerably.

Before we move on to model each regime in turn, we note that we had to exclude the first 14 points from the original DC9 data (greyed in the Mooney plot of figure 2), because they prevented the DC9 Mooney plot from having a linear decrease in the small-to-moderate regime. For reasons explained later, every soft solid *must* have such a linear decrease. The reason for the experimental discrepancy presented by the first 14 points is not known, but is probably due to non-optimal accuracy regime for the load-cell in the low load regime and/or the presence of a slight slack of the sample in its starting position.

### (b) From small to moderate strains

In this section, we focus on the term in the strain energy that will model the small-to-moderate regime of deformation for the samples. We identify this region by the *linear part* of the Mooney plot occurring before the upturn (figure 2). Hence we perform the curve-fitting exercise on the first *N*_{0} data points only, where we determine *N*_{0} by conducting a model-dependent sensitivity analysis (to be detailed later).

First, we see that the neo-Hookean model (2.3) is incapable of reproducing the simple extension data in any meaningful way. That is because, for *W*_{nH}, the Mooney plot formula (2.9) yields

For a polynomial dependence of *W* on *I*_{1}, we can consider the Yeoh strain-energy density [27]
*c*_{i} are constitutive parameters (by (2.5), the infinitesimal shear modulus of this material is *μ*=2*c*_{1}). This model gives of course a better fit in the small-to-moderate region than the neo-Hookean potential, figure 1, lower panels, and table 1.

However, as is evident from the formula (2.9), we must include a dependence on the second principal invariant *I*_{2} in order to make meaningful progress.

Hence we will consider, in turn, the following two classical strain energies: the *Mooney–Rivlin model* [4],
*Gent–Thomas model* [28],
*Carroll model* [16]
*C*_{1} and *C*_{2} are the material constants to be found from the best-fit procedure. As noted by Carroll [16], all three belong to the class of strain energies proposed by Klingbeil & Shield [29]. When these strain energies are expanded up to third order in the powers of the Green–Lagrange strain tensor ** E**, we find that they all give

To determine rigorously how many data points we should take for the small-to-moderate region, we performed a sensitivity analysis for all three models. We found that the maximal relative error is at its lowest value when *N*_{0}=7 for the Treloar dataset and *N*_{0}=26 for the DC9 dataset (figure 4). We see from the resulting best-fit curves of figure 5 that, as expected, the *I*_{2}-dependence is precisely the missing ingredient to obtain excellent agreement in the small-to-moderate regime. The figure is eloquent on two features.

First, the relative errors are dramatically reduced over that regime, to less than 1.8% for Treloar's data and less than 3.6% for the DC9 data (table 2). This improvement is achieved with a set of only two material parameters ** p**=[

*C*

_{1},

*C*

_{2}] at our disposal, in contrast with the Yeoh model (3.3) which gave much higher errors with a set of three fitting parameters

**=[**

*p**c*

_{1},

*c*

_{2},

*c*

_{3}].

Second (as noted earlier by Carroll [16]), the actual dependence of *W* on *I*_{2} does not matter much, with all three models (3.4), (3.5), (3.6) performing equally well. This is not altogether surprising in view of the third-order expansion (3.7) and its universality and relevance to the small-to-moderate range of extension. In effect, the equivalence of the Mooney–Rivlin strain energy and the third-order expansion shows that every incompressible isotropic hyperelastic solid must present a linear decrease at first in the Mooney space.

In the next section, we denote the functional dependence of *W* on *I*_{2} generically by *f*(*I*_{2}) to save space.

### (c) Strain-hardening regime

We now move on to the *strain-hardening regime*, which is how we called the range of data corresponding to the *upturn* in the Mooney plots (figure 2). Thus, we will perform our curve fitting exercises on the *N*_{0}+*N*_{p} first data points, where *N*_{p} is the number of points required to capture that the Mooney plot has gone through a minimum. We take *N*_{p}=5.

Keeping in line with our methodical approach to modelling we now add a power-law term to the strain-energy density of the previous section and consider *W* in the form
*C*_{3}>0 (N mm^{−2}) and *n*>0 (non-dimensional) are constants [30]. As noted by Carroll [16], as the tensile stretch *λ* increases, the principal force associated with an energy term *λ*^{2n−1} (Note that we could have chosen to add a *λ*^{4m−1}.) Hence, for an ideal (no solvent) polymer, *n*=1 and we recover the Gaussian Chain model; for a polymer in a good solvent, *n*=5/4 and we recover the Pincus correction of the force behaving as *λ*^{3/2}.

Here we wish to model the strain stiffening as illustrated by the upturn in the Mooney plot. As *λ* increases, *λ*^{2n} and *σ* as *λ*^{2(n−1)}. Hence in the Mooney plot, as *z* goes towards zero, *g*(*z*) behaves as *z*^{−2(n−1)}, indicating that *n*>1 is required for an upturn. We point out that the optimization of the material parameter set ** p**=[

*C*

_{1},

*C*

_{2},

*C*

_{3},

*n*] is a

*nonlinear*procedure. In general, nonlinear curve fitting exercises lead to some serious computational problems, the most common being the emergence of multiple minima for the same level of error, see the discussion by Ogden

*et al.*[31] for rubber models. Also, the choice of an adequate initial guess is a crucial issue in order to avoid false minima [32]. To circumvent these problems we propose the following ad hoc procedure. For definiteness we consider that the optimal value of

*n*should be found in the range 1.0<

*n*<2.5: the lower bound indicates a departure from the linear fit of the small-to-moderate Mooney plot region and, at twice the Pincus value, the upper bound is a reasonable ceiling to capture the upturn of the Mooney plot prior to the large strain-stiffening regime. We limit the range of data to the

*N*

_{0}+

*N*

_{p}first points, spanning the small-to-moderate regime (linear variation in the Mooney plot) and the strain-hardening regime (upturn in the Mooney plot). Then we fix

*n*at the beginning of the 1<

*n*<2.5 range, at

*n*=1.01 say, and perform the fit for the reduced set

**=[**

*p**C*

_{1},

*C*

_{2},

*C*

_{3}]. This fit is

*linear*and thus removes the risk of multiple minima and the requirement of a good initial guess. Then we increase

*n*to span the range and record the corresponding maximum relative errors err(

*n*). Finally we keep the optimal value of

*n*, corresponding to

For Treloar's data we fix *N*_{0}+*N*_{p}=7+5, and for the DC9 data we take *N*_{0}+*N*_{p}=24+5. For the fitting of the former set, the Gent–Thomas model with a strain-hardening correction gives the lowest error, while for the later set it is the Mooney–Rivlin model with a strain-hardening correction which performs best. However, similarly to the previous section, the differences in performance between the three models are negligible, see table 3 for a summary. Overall the fit is excellent over the small-to-moderate and strain-hardening regimes, as attested by the low values of the relative errors, around or less than 2%. In table 3 we also report the corresponding constants of weakly nonlinear elasticity *μ*_{0}, *A*, and *D*. They are obtained by expanding *W* in (3.9) up to fourth order in the Green–Lagrange strain tensor. We find that now
*D* in (2.6) we have
*C*_{1} and *C*_{2} parameters from tables 2–3, reflecting that these parameters are simply the result of a curve-fitting exercise. The weakly nonlinear elasticity constants *μ*_{0} and *A* however carry consistently from one model to the other.

In this section we showed that a strain energy with three terms and four material parameters can cover the range of 500% extension for Treloar's data (1≤*λ*≤6) and 300% extension for the DC9 data (1≤*λ*≤4). It is worth noting that for each model the solution for the optimal set of parameters is unique and the problems are not ill-conditioned.

If we wish to go further and also capture the behaviour of soft matter in a regime of extreme extension (6≤*λ*≤8 and 4≤*λ*≤6, respectively), we have to recognize that the corresponding stiffening of the curve is associated with a singularity of the fitting function. In other words, *limiting-chain effects* come to dominate in the latter regime and impose an asymptotic barrier. In the next section we present a method to determine the order of the singularity by examination of the experimental data.

### (d) Estimate of the order of singularity: Freely Jointed Chain versus Worm-Like Chain model

In this section and the next, we focus on the *extreme stretch range*, where the material stiffens rapidly with strain. First we consider only the last *N*_{f} points of the force–extension data, and estimate how fast the data curve ‘blows up’.

Assume that in the large stretch range, the tensile force behaves as
*A*_{0} is the cross-section area of the sample in the undeformed state, *c* is a constant, *λ*=*λ*_{m} indicates the location of the unknown vertical asymptote, and *k* is the *order of singularity*: *k*=1 for the FJC model and *k*=2 for the WLC model. In logarithmic coordinates we have
*k*, as follows. For the last *N*_{f} data points, we can calculate *λ*_{m} is greater than *λ*_{f}, the largest stretch reported in the experiments (hence *λ*_{f}=7.6 for the Treloar set, see (3.1), and *λ*_{f}=5.54 for the DC9 set). Then we start by fixing *λ*_{m}=*λ*_{f}+0.001, say, and plot *λ*_{m}, and so on. Here for both datasets, we find that the maximal relative error goes through a minimum as *λ*_{m} increases. When it is at its minimal value, we can say that we have identified the actual order of singularity *k** and the actual limiting stretch

In figure 7, we report these results for the Treloar's and DC9 data, in the upper and lower panels, respectively, where the optimal solutions yielding *k** and *N*_{f}=9, and find *k**=0.962, with a maximal relative error of 2.1%, indicating a *first-order singularity* location for large deformations. For the DC9 rubber dataset, we take *N*_{f}=5 and find *k**=2.03 with a maximal relative error of 0.6%, clearly pointing out to a *second-order singularity*.

If our goal is to provide a one-dimensional interpretation of the data at extreme stretches, then our method shows that the rubber studied by Treloar behaves according to the FJC model, while the DC9 rubber behaves according to the WLC model.

The FJC model is easily extended to a three-dimensional strain-energy function, as is often captured by the Arruda & Boyce [33] 8-chain model, itself represented well by the phenomenological Gent model [34],
*C*_{1}>0 is the initial shear modulus and *J*_{m}>0 is a stiffening (limiting-chain) parameter.

The WLC model is not easily translated into a three-dimensional strain-energy function, see Ogden *et al.* [35]. Dobrynin & Carrillo [14] proposed to mimic the WLC behaviour with the following phenomenological model
*G*>0 and is the ‘network infinitesimal shear modulus’ and *β* is a stiffening parameter. In passing, note that by (2.5) the infinitesimal shear modulus is

Guided by these models, the next section will try to capture the nonlinear behaviour of the Treloar data and its first-order singularity by invoking the first model above, and the behaviour of the DC9 data and its second-order singularity with the second model above.

### (e) Limiting-chain effect

In this section, we bring together all our accumulated knowledge to reproduce the nonlinear behaviour of the rubber tested by Treloar over the *entire range of experimental stretches*.

The natural summation of our analysis so far is that the whole model to be studied should be given, for the Treloar data, by following form
*W*(3,3)=0.)

Here the main difficulty is to identify the optimal value of the limiting-chain parameters *J*_{m} or *β*, because their presence implies that the fitting procedure is a *nonlinear* least-square optimization problem. A possible consequence is that a set ** p**=[

*C*

_{1},

*C*

_{2},

*J*

_{m}] or [

*C*

_{1},

*C*

_{2},

*β*] could be wrongly identified as the optimal one, while there exist other minima giving a better or similar error [31,32].

To investigate fully the nature of our optimal set, first we solved the nonlinear least-squares (NLS) problem by applying the `Lsqcurvefit` routine in Matlab and the `NonlinearFit` routine in Maple, both based on the Levenberg–Marquardt method. Second, since the models contain only one nonlinear term, we were also able to apply the VARPRO method devised by Golub & Pereyra [37]. This ‘Variable Projection method’ algorithm is designed for nonlinear LS problems where some of the parameters to be identified appear in linear terms, as here. We used a recent Matlab routine [38], updated and improved with respect to the original Fortran routine. Third, we used a ‘semi-linear’ approach. In that case we identified the values of the constants *C*_{1}, *C*_{2} by a linear fit on the first *N*_{0}+*N*_{1} data points. Then we varied the value of the stiffening parameter *J*_{m} or *β* and recorded the corresponding maximum relative error. Then we kept the value of *J*_{m} or *β* giving the lowest relative error over the range. In the end, these three strategies yielded very similar results, with small differences in the values of the constants and with low error, comparable with the experimental error that can be expected from uniaxial tension experiments. There is little value in presenting all the results of all three procedures and instead we just present those of the NLS routines, which are the simplest to implement.

For completeness and comparison, we used both models above on both the Treloar and the DC9 data. We called them GMR, GG, GC and DCMR, DCG, DCC, corresponding to the Gent (3.16) or Dobrynin and Carrillo (3.17) models, respectively, with an *I*_{2} dependence of the Mooney–Rivlin, Gent, and Carroll type, in turn. Notice that the GG model first appeared in Pucci & Saccomandi [39] and the DCMR model in Dobrynin & Carrillo [14].

We also calculated the corresponding constants of second- and third-order elasticity, as

In table 4 we collected the results of the curve fitting exercises over the whole range of available data for the Treloar data, and similarly in table 5 for the DC9 data.

For all models, the fit is excellent. Hence with *only three parameters*, we are able to cover the full reported data with an error that compares favourably with the experimental error (recall that Treloar reported values with only three significant digits, see (3.1)). In particular, the ‘Gent–Gent model’ [31,39]
*N*_{0}+*N*_{p} data points), so that the model (3.9) with a power-law term can be advantageously replaced with the Gent–Gent model (figure 1).

## 4. Simple and rectilinear shear

Now it remains to figure out how robust and predictive our results are with respect to other deformations. For instance, Ogden *et al.* [31] showed that two different Ogden models could give the same excellent fit for uniaxial tension and widely different predictions for rectilinear shear, thereby invalidating them. Here we first consider the *simple shear deformation*,
*K* is the amount of shear. This deformation is used to idealize some experimental test protocols such as the *double sandwich test* and the *quadruple shear test* [40]. However, it is not easy to find simple shear test data in the literature on nonlinear elastic solids, for the reasons discussed by Destrade *et al.* [41]: in nonlinear elasticity theory, simple shear is not due to a pure shear stress field, and normal stress differences have to be taken into account.

In the literature, simple shear data are often inferred from a pure shear experiment, see Treloar [9]. Yeoh [27] does present some simple shear data obtained with a quadruple test. He observes that ‘These normal components of stress do not have much effect on the overall load/deformation behaviour and are ignored.’ However, we note there exists no protocol to measure normal stress components during a simple shear test, and therefore we cannot quantify and decide whether these stress quantities can indeed be neglected. In fact, normal stress effects (Poynting effect) can be observed directly in experiments, by drilling holes in the shearing plates and observing the bulging of the material [42,43].

Now to draw meaning comparisons using the Mooney–Rivlin, Gent–Thomas and Carroll models (3.4), (3.5) and (3.6), we must stay in the small-to-moderate range of stretches, up to a maximal stretch of extension *λ*_{1}≃2, indicating that the material cannot be contracted to a stretch lower than *K*=0.71. In figure 9*a*, we plot the predictions in simple shear for Treloar's rubber, using the material constants of table 2. Clearly all three models are in good accord, with the Carroll and Gent–Thomas models indicating a slight softening for shears *K*≥0.5.

Moving on the models (3.9), including a strain-hardening term to model the upturn, we now limit the simple shear predictions to the range 0≤*K*≤1.0, as the upper bound reflects a tilting angle of 45^{°} for the slanted faces, a clear *physical* limit for a realistic simple shear test. Figure 9*b* displays the comparison between each model, using the material parameters of table 3, and again shows good agreement in the predictions. These results confirm that the *robustness* already noted in the fitting of simple extension data is conserved when considering a deformation with a completely different character (recall that simple shear is a plane deformation, with changing orientation of the principal axes with the deformation, in contrast with the three-dimensional uniaxial tension with conserved principal orientations).

Next, we consider the following *rectilinear shear deformation*
*h*. This inhomogeneous deformation leads to a simple boundary-value problem (BVP) which has to be solved numerically for each different model. The BVP is as follows: a slab is located in the region −*H*≤*Y* ≤*H*, and its deformation is driven by a gradient of pressure in the *X*-direction, with *h*(−*H*)=*h*(*H*)=0, see Ogden *et al.* [31] for details. The aim here is again to check whether the various models with material parameters listed in tables 2 and 3 make comparable predictions of the deformation field. The answer is found in figure 10, showing *h*(*Y*) through the slab when *H*=1, for increasing levels of pressure gradient: the variation between the different models is hardly notable, showing again the robustness of the methodical modelling presented in the paper. This result is in sharp contrast with the simulations conducted by Ogden *et al.* [31] for the same BVP: there the same Ogden models, characterized by different values of the material parameters, gave the same close fitting to uniaxial tension and yet predicted different deformations.

## 5. Conclusion and discussion

The fitting of constitutive parameters in the mathematical models of rubber-like materials is often considered to be a trivial task, although in fact this is not the case at all. Indeed more than 10 years ago, Ogden *et al.* [31] had already pointed out how the fitting procedure is a very delicate aspect of the modelling procedure. Here we went one step further and proposed a possible solution to the main issues. Our solution is not a new form of the strain energy; rather it is intrinsic to the theory of nonlinear elasticity and proves its robustness.

We used experimental data already published in the literature to show that our method is quite general and does not rely too heavily on the quality of the experimental data to capture the salient features of the behaviour of rubber. Clearly, it is possible to collect data in a more careful and uniform way for our treatment, for example, by fixing *a priori* the range of extension of interest and by fixing a uniform incremental step in the acquisition of the data for all the samples.

The starting point of our method was to delineate three ranges of deformations in the nonlinear framework: the small-to-moderate range, the strain-hardening range (the zone of the upturn point in the Mooney plot) and the region of limiting-chain singularity. The aim of the procedure was to find a mathematical model able to fit the data within a reasonable error in part or all of the range. The main focus was more on the computational aspect of the fitting procedure than on the corresponding potential experimental issues. Using a step-by-step method, we came to the conclusion that the theory of nonlinear elasticity is a more robust theory than was previously thought. We saw that minimizing the relative errors of (2.13), as opposed to minimizing the classical (absolute) residuals of (2.12), ensures consistency of the fitting procedure across the various measures of stress and strain. We demonstrated quantitatively that for incompressible materials, the strain energy must depend on both principal invariants *I*_{1} and *I*_{2}. Moreover, we found that any standard linear combination of function of these two invariants fits well the data in the small-to-moderate range.

Then in order to model the upturn zone, we proposed an additional term to incorporate strain-stiffening effects. We also required full compatibility with the weakly nonlinear theory of fourth-order elasticity. This compatibility is not a mathematical whim but is dictated by the generality of the response of rubber-like materials to large shear and large bending deformations. Again the resulting model provides a robust mathematical behaviour able to capture the data with low relative errors.

For the last step, we needed to understand what is going on at very large deformations. This range of deformation has been studied in details in the last two decades and two reliable models have emerged: the Worm Like Chain model [35] and the Freely Jointed Chain model, itself accurately modelled by the Gent model [34,44]. We proposed a practical way to investigate the mathematical behaviour of the vertical asymptote associated with the limiting-chain effects, and thus determine the order of the singularity for a given set of data. Hence we were able to confirm the insights coming from statistical mechanics computations on ideal macromolecular chains.

Finally, we showed that the three-parameter Gent–Gent model proposed by Pucci & Saccomandi [39] is capable of fitting the entire deformation range with low relative error for both sets of experimental data studied. This versatility singles out the Gent–Gent strain-energy density as a powerful model for all nonlinear aspects of rubber-like behaviour (see also [31,39,45]), although we note that the other five models investigated also performed well at the curve-fitting tasks, including the Carillo–Dobryin model proposed in 2011.

By deconstructing, decrypting and then recomposing all the puzzle of mathematical models for rubber-like materials, we provided evidence that the theory of nonlinear elasticity can describe accurately uniaxial tension data with a strictly bounded error and a low number of unique parameters. The theory of nonlinear elasticity is a basic, well-grounded theory of continuum mechanics. It is a mandatory passage required to explain the mechanical behaviour of many real-world materials. Therefore, it was fundamental to confirm that this theory is not only mathematically well founded but also provides a robust description of experimental data, which is not contingent to the choice of a specific strain-energy function, provided it satisfies some basic requirements and is adequate for (a) given regime(s). Using one of the other models proposed in the literature instead of the ones constructed here cannot change our main findings, because the heart of the modelling problem is independent on the functional form of the strain energy but relies instead on the methodology of fitting and on the way we read the data. This is essentially good news for the theory of nonlinear elasticity. Indeed if this theory was dependent on the choice of a given functional form, we would have to declare it a poor and restricted mathematical theory of mechanics.

Finally, because one of the most important elements of a computer simulation for deformable solid mechanics is the actual model of a material, our findings on the mathematical structure of the mathematical models for rubber-like materials is bound to have non-trivial consequences on their numerical implementation.

## Data accessibility

See (3.1).

## Authors' contributions

All the authors have contributed equally to all parts of the paper.

## Competing interests

There is no conflict of interest.

## Funding

G.S. is partially supported by GNFM of Istituto Nazionale di Alta Matematica, Italy.

## Acknowledgements

We are most grateful to Dr Carrillo for providing us with the tabulated version of the DC9 data used in his paper Dobrynin & Carrillo [14].

## Footnotes

This paper is dedicated to the memory of Michael A. Hayes (1936–2017).

- Received October 31, 2016.
- Accepted January 10, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.