## Abstract

In this paper, we review critically some basic models derived from the statistics of long-chain molecules and then discuss the status of such models within the three-dimensional nonlinear theory of elasticity. We draw attention to some deficiencies of certain worm-like chain (WLC) models when viewed within the three-dimensional continuum mechanics framework. Modifications of the WLC models motivated by such considerations and avoiding these deficiencies are then discussed and shown to correspond well with data generated by the exact WLC model.

## 1. Introduction

Rubber elasticity has one of the longest and distinguished histories in all of polymer science (see, for example, the discussion in Saccomandi & Ogden 2004). However, the complex molecular nature of rubbers (and other elastomeric materials) and the need for the use of nonlinear theories to describe their mechanical behaviour have made difficult the development of realistic mathematical models. The theory of incompressible, isotropic hyperelasticity is a starting point for the modelling of elastomeric materials, and within this framework two main directions can be followed in developing specific constitutive models. These are the *phenomenological* approach based on continuum mechanics arguments (Saccomandi 2004), and the *statistical* approach based on molecular concepts (Erman 2004).

In following the molecular approach, the first step is to derive an elasticity model for a single polymeric chain within the rubber network and then to introduce a suitable averaging procedure to obtain a macroscopic constitutive equation for the material response. For example, in the elementary molecular theory of networks it is postulated that the elastic free energy of a network is equal to the sum of the elastic free energies of the individual chains (Erman 2004). On the other hand, the main tools of the phenomenological approach are the axioms of continuum mechanics together with a systematic and rigorous use of the restrictions imposed by thermodynamics, the use of representation formulae based on symmetry arguments and *a priori* restrictions induced by various mathematical requirements (Saccomandi 2004).

Rubber-like elasticity is also fundamental in the study of biological proteins and materials that are characterized by high resilience, large strains and low stiffness (Gosline *et al*. 2002). This is, for example, the case in the pioneering paper of Fung (1967) on mesentery tissue. In his paper, which may be considered as a starting point for the modern theory of cardiovascular solid mechanics (Humphrey 2002), an empirical relation for the highly nonlinear uniaxial stress-stretch behaviour of the tissue is proposed. This relation is the origin of the celebrated exponential Fung model for soft tissue mechanics (e.g. Humphrey 2002). However, this one-dimensional empirical relation is not consistent with the equations of nonlinear elasticity since the associated uniaxial stress does not contain the factor *λ*−1 that arises for simple tension in the general theory (*λ* being the stretch in the direction of extension). Thus, without modification, it cannot be extended to the general three-dimensional setting that is required to ensure its consistency within the more general framework. What is today referred to as the Fung model in this general setting was introduced, in the case of isotropy, by Demiray (1972). In recent years, there has been a steady stream of theoretical and experimental research aimed at constructing models compatible with the laws of continuum mechanics but also related, where possible, to the molecular approach (e.g. Holzapfel & Ogden 2003).

The situation that we have described above for soft tissue mechanics is similar to that arising in consideration of single-molecule studies within the context of DNA mechanics and the mechanics of other biological macromolecules and proteins. Stretching of individual biomolecules is now achieved by a variety of techniques that allow measurement of forces from 10 femtonewtons to hundreds of piconewtons (Bao & Suresh 2003). The entropic elasticity of these biomaterials is usually described by using one-dimensional models based on statistical mechanics expressions of the entropy that idealize each macromolecule as a freely jointed chain (FJC) or a worm-like chain (WLC).1 This is the case for DNA (Strick *et al*. 2003) and also for some structural proteins such as elastin (Urry *et al*. 2002) and collagen (Sun *et al*. 2004). A detailed introduction to the WLC may be found in appendix G of the book by Flory (1988).

The idea of using a general three-dimensional model based on a chain network in which the basic building block is the WLC has recently been proposed by Bischoff *et al*. (2002*a*,*b*) and Kuhl *et al*. (2004). The aim of the present paper is to build upon this work by considering a three-dimensional version of the WLC model within the context of the continuum theory of nonlinear elasticity.

In §2, we begin by reviewing critically some basic models derived from the statistics of long-chain molecules with particular attention to problems that may arise in fitting the models to the data. This is followed by discussion of the status of such models within the three-dimensional nonlinear theory of elasticity. In §3, we draw attention to deficiencies, from the continuum mechanics perspective, of a specific WLC model, with particular reference to simple tension and simple shear. Sections 4 and 5 are concerned with modifications of the Gent model of rubber elasticity (Gent 1996) that are motivated by the above considerations and avoid the deficiencies encountered in the WLC model. Some closing remarks are contained in §6.

## 2. Basic equations

### (a) Some long-chain molecular models

An introduction to the statistical properties of long, flexible objects such as polymer chains may be found in, for example, the book by de Gennes (1979). In the freely jointed chain model, a molecule is considered to be made up of rigid orientationally independent chains. Alignment of segments by a tension, denoted *f*, is described by using non-Gaussian statistics with the relative stretch *r*/*L* as a variable, where *r* is the end-to-end distance of the chain and *L* its *contour length*. Specifically,(2.1)where *k* is the Boltzmann constant, *T* the absolute temperature, *N* the numbers of freely jointed segments (each of length *l*, so that *L*=*Nl*), and 0<*r*≤*L*. The symbol ^{−1} stands for the inverse of the Langevin function, which is defined by(2.2)This is a classical model from which various three-dimensional micromechanical network chain models have been derived (Boyce & Arruda 2000).

In the WLC model, the chain is treated as a flexible beam that curves as a result of thermal fluctuations. An excellent treatment of the WLC may be found in the paper by Bouchiat *et al*. (1999), where this model is derived by using modern methods of statistical physics. Because, in the case of the WLC, an analytically exact expression for the force in a single chain is not available, it is usual to consider an interpolation formula of the form(2.3)as introduced by Marko & Siggia (1995), where *A*, with *l*≤*A*≤*L*, is a parameter referred to as the *persistence length* of the chain. A more complex form than (2.3) is also given in Bouchiat *et al*. (1999), where an amended version is proposed that takes account of some enthalpic corrections for large forces. This is written(2.4)where *α*_{i}, *i*=2, …, 7, are phenomenological parameters that account for the enthalpic corrections.

Here, we use the term ‘WLC model’ for these interpolation formulae instead of for the exact model, and it is therefore important to point out that our comments are relative to the interpolation formulae rather than the *real* WLC model.

The main characteristic of the functions (2.1), (2.3) and (2.4) is that when *r* approaches the fully extended length of the chain, i.e. *r*→*L*, the force *f* blows up, i.e. *f*→∞.

At this point some comments on the formulae (2.1) and (2.3) are called for, for which purpose we introduce the extension ratio *z*=*r*/*L*, such that 0<*z*≤1. Equation (2.3) may be rewritten as(2.5)This is the ratio of a third- and a second-order polynomial and has a singularity at *z*=1. Also *f*(0)=0. In (2.4), a polynomial term is added to the basic function (2.3) in order to improve the overall fitting of the model to the data, and we now discuss this ‘improvement’ in detail. Our aim is to show the importance of the singularity in (2.3) and for this reason we introduce the reduced WLC function(2.6)i.e. equation (2.5) with the term 4*z* omitted.

In figure 1, log plots of both (2.5) and (2.6) are shown for comparison with the ‘data’ taken from Bouchiat *et al*. (1999). These are data generated numerically from the exact WLC model. The least squares fitting problem can be expressed in general as min_{p} *F*(* p*), where the squared two-norm of residues (the ‘objective function’) is defined by(

*x*

_{i},

*y*

_{i}),

*i*=1, …,

*m*, are the data points and

*=[*

**p***p*

_{1}, …,

*p*

_{n}] represents the vector of parameters to be identified,

*f*(

*,*

**p***x*

_{i}) is the force predicted by the constitutive model for the value

*x*

_{i}and

*y*

_{i}is the corresponding experimental value of the force. If

**p**^{*}is the optimum set of parameters obtained by the above fitting procedure the corresponding residual is

*S*=

*F*(

**p**^{*}) and the percentage relative error is defined byIn order to evaluate the ‘goodness of fit’ it is also usual to calculate the root mean square error (RMSE), here defined aswhere

*m*−

*n*is the number of degrees of freedom of the optimization problem and

*f*

_{ref}is the mean of the

*y*

_{i}data values. Note that, in contrast to

*S*,

*ϵ*is dimensionless. In the following, we also report this value for each example.

When we deal with equations (2.5) and (2.6), the only parameter to be identified is *kT*/4*A* and in this case the problem can easily be solved as a linear problem. To deal with the nonlinear least squares problems in the following sections we use the *lsqcurvefit* routine in Matlab 6.5 by selecting the Trust Region Method.

In the fitting reported in figure 1, we record that for (2.5) we have *kT*/4*A*=0.2486 and that the residual *S* is 2.3295 and the RMS error is *ϵ*=0.0106, while for (2.6) the corresponding values are *kT*/4*A*=0.2506, *S*=1.2018 and *ϵ*=0.0076. Thus, the residual *S* and the RMSE *ϵ* are less for the reduced form (2.6), although this is not apparent from inspection of the curves since the fit is clearly better for (2.5) in the lower range of extensions. The differences in *S* and *ϵ* are accounted for by the much better fit of (2.6) at larger extensions. This is clear from the numerical values but is not apparent from the figure because of the large gradients in this region. Our analysis shows that the function in (2.6) needs to be improved but that ‘improvement’ by addition of the term 4*z*, introduced to fit the data for small extensions, is not good. Indeed, the quality of the fitting deteriorates at the higher values of the extension although a better agreement for small extensions is achieved. For this reason, the formula (2.4) was proposed by Bouchiat *et al*. (1999). We note, however, that the large number of constitutive parameters introduced in (2.4) may be problematic. Indeed, in Ogden *et al*. (2004) it was shown that increasing the number of polynomial terms in phenomenological models may introduce numerical instabilities and non-uniqueness of the optimal sets of parameters. The problem of multiple solutions has not been apparent previously in the present context because DNA elongation data have been based on the formula (2.4) with the coefficients *α*_{i} set to the values identified by Bouchiat *et al*. (1999). In fact, (2.4) is normally used in a nonlinear fitting procedure to identify only the contour length *L* and the persistence length *A* of actual data. We deal with this point in a forthcoming paper (Ogden *et al*. submitted). It should be emphasized that it is the term with the singularity that is responsible for the very good representation of the non-Gaussian behaviour of the data at larger extensions.

We end this section by recalling that the elastic energy function *W* associated with (2.5) is(2.7)where *W*_{0} denotes a constant energy with respect to a reference level. This is a one-dimensional strain energy. In §2*b*, we consider the possibility of extending this model to three dimensions.

### (b) Nonlinear elasticity and molecular models

It is important to point out that the expressions for the tensile force in equations (2.1), (2.3) and (2.4) are *not* compatible with the three-dimensional theory of elasticity. This is not just an academic consideration since there are many situations where single molecule experiments are performed involving torsional deformations (e.g. Strick *et al*. 2003), and torsion requires three-dimensional analysis. Within the framework of the phenomenological theory of nonlinear elasticity the constitutive equations for the strain-energy function of an isotropic material are based on the choice of a scalar function of a suitable set of deformation invariants or, equivalently, on the principal stretches of the deformation. The usual invariants adopted are the principal invariants of the right Cauchy–Green deformation tensor * C*=

**F**^{T}

*, where*

**F***is the deformation gradient. These are defined by(2.8)In terms of the principal stretches*

**F***λ*

_{1},

*λ*

_{2},

*λ*

_{3}, we have(2.9)When the internal constraint of incompressibility is in force all admissible deformations must satisfy the restriction . A basic result from the general theory of nonlinear elasticity (e.g. Ogden 1984) is that for an incompressible isotropic elastic material characterized by a strain energy

*W*(

*I*

_{1},

*I*

_{2}), the Cauchy stress tensor

*has the form(2.10)where*

**T***is the identity tensor,*

**I***=*

**B**

**FF**^{T}is the left Cauchy–Green deformation tensor and

*W*

_{i}=∂

*W*/∂

*I*

_{i},

*i*=1, 2.

It is clear that any possible incompressible isotropic strain energy derived by molecular arguments must be expressed as a function of the first two invariants in (2.9), and the corresponding Cauchy stress tensor must be compatible with formulae such as (2.10). An important example derived by using molecular arguments starting from the non-Gaussian FJC theory is that due to Arruda & Boyce (1993), while Beatty (2003) has provided continuum formulations of the various strain energies that are based on the FJC model.

Consider next a single perfectly flexible and freely jointed chain whose undeformed mean length is assumed to be . The relative chain stretch *λ*_{r} is the ratio of the current chain vector length *r* to its fully extended limiting length *Nl*. Thus,(2.11)where *λ*_{chain}=*r*/*r*_{0}. Note that *λ*_{r}∈[*N*^{−1/2},1]. Now suppose that the chain vector **r**_{0}=(*X*, *Y*, *Z*) in the reference configuration is directed along the line *X*=*Y*=*Z* through the origin of a rectangular Cartesian system. This has length and its direction is . If the network is subjected to an affine deformation in which the rectangular Cartesian axes coincide with the principal axes then the squared stretch of a chain initially oriented in an arbitrary referential direction * m* is given by(2.12)which specializes to

*I*

_{1}/3 for the specific

*above. The same result holds for a chain whose end point is oriented initially according to −*

**m***X*=

*Y*=

*Z*or −

*X*=−

*Y*=

*Z*or

*X*=−

*Y*=

*Z*. For the Arruda–Boyce eight-chain model, we then have(2.13)

Beatty (2003) showed that some of the constitutive equations proposed on the basis of the FJC model lead to a Cauchy stress tensor in the simple invariant form(2.14)where(2.15)For a single FJC, we have(2.16)

To avoid the use of the inverse Langevin function, it is usual to replace ^{−1}(*λ*_{r}) by a polynomial expansion (e.g. Boyce 1996), for which the crucial characteristic of the singularity in (2.16) is then lost. Another possibility is to replace the inverse Langevin function by its Padé approximation, as, for example, in Miehe *et al*. (2004). However, Horgan & Saccomandi (2002*a*) have shown that a very good global approximation is obtained by using the phenomenological model first proposed by Gent (1996). The Gent model has strain energy(2.17)where *μ* is the shear modulus for infinitesimal deformations and *J*_{m}(>0) is the limiting value of *I*_{1}−3 associated with limiting chain extensibility. (The upper bound on *I*_{1} is consistent with that associated with (2.13) if we take *J*_{m}=3*N*–3.) In the limit, as the chain extensibility parameter tends to infinity (*J*_{m}→∞), (2.17) reduces to the classical neo-Hookean form(2.18)

The model (2.17) is a fundamental one because it is associated with the *simplest* rational approximation of the response function in the Cauchy stress representation formula for an incompressible isotropic elastic material for which there is a fixed singularity (Horgan & Saccomandi 2003). The Cauchy stress associated with (2.17) is given by(2.19)so that the stress has a singularity as *I*_{1}→*J*_{m}+3. Note that this is a special case of (2.14).

## 3. An eight-chain WLC-based model

Kuhl *et al*. (2004), following Bischoff (2002*a*,*b*), have developed a three-dimensional version of the WLC model by considering a three-dimensional strain energy associated with (2.3). This may be obtained by establishing first a connection between *r* and the invariants in (2.9) and then relations similar to (2.11) and (2.13), as for an eight-chain WLC model. In remark 1 of Kuhl *et al*. (2004), the relationis considered, where *a* is a constant parameter that may be chosen so that *r*/*L* coincides with (2.13). It should be noted that the range of values of *r*/*L* in the one-dimensional WLC models is (0,1] for the extension ratio, whereas the stretch in the direction of extension is greater than one.

From equation (2.7), we obtain(3.1)Now, for a fixed temperature, the persistence length *A* is a constant proportional to the bending stiffness of a chain segment. Then, when *N*→∞ in (3.1), we find that this strain energy blows-up, which is a pathological and un-physical result. It might be expected that for *N*→∞ we recover the neo-Hookean energy function or at least a ‘classical’ strain energy that is well defined in this asymptotic limit (Horgan & Saccomandi 2003). It does not seem possible to circumvent this problem for this model, even if *A* is treated as a phenomenological parameter that depends on *N*.

In what follows it is convenient to introduce the notation *J*_{0}=3*N* and, for curve fitting purposes, to treat *kTNl*/4*A* as a single parameter. Then, we re-write (3.1) as(3.2)from which we deduce that(3.3)which will be required in §3*a*.

In the paper by Kuhl *et al*. (2004), the basic model (3.1) is considered in a more general framework in which anisotropy, compressibility and growth are included, and their paper focuses more on the description of remodelling of biological tissue than on a detailed study of the mechanical properties of the strain energy associated with the WLC model. In the subsequent paper by Kuhl *et al*. (in press), the model derived within their framework is analysed mathematically with a view to establishing various convexity properties.

Here, we study some basic mechanical properties of (3.1) and, to this end, we start by considering the mechanical response of this model for some simple homogeneous deformations.

### (a) Simple tension

For an incompressible material, the principal stretches in simple tension are (e.g. Ogden 1984)(3.4)where *λ* is the stretch in the direction of extension. It follows that(3.5)and the principal components of the Cauchy stress tensor (2.10) are given by(3.6)The requirement that the lateral surfaces of the specimen undergoing the simple extension are traction free, *t*_{2}=*t*_{3}=0, yields(3.7)and the tensile force *f* per unit of undeformed cross-sectional area necessary to achieve the stretch is given by(3.8)We emphasize that here and subsequently *f* is, therefore, the nominal (or engineering) *stress*, whereas in the one-dimensional models it is the tensile *force*.

In respect of (3.3), we obtain(3.9)The tensile stress *f* in (3.9) blows up when *λ*→*λ*^{*}, where(3.10)and we also note that *f*→∞ as *N*→∞.

It is interesting to compare (3.9) with the tensile stress associated with the Gent material (2.17), namely(3.11)In the case of (3.11), the limit *J*_{m}→∞ enables us to recover the tensile force(3.12)associated with the Gaussian (neo-Hookean) model (2.18).

Again it is possible to propose a reduced form of the WLC model by omitting the first and third terms in (3.9), which then reduces to(3.13)

Figure 2 shows a comparison of (3.9) and (3.13), which parallels that shown in figure 1 in respect of (2.5) and (2.6). In the case of (3.9), we have *kTNl*/4*A*=0.0875 and *J*_{0}=5 with a residual *S*=0.041 and RMSE *ϵ*=0.0014, while for (3.13) *kTNl*/4*A*=0.0878, *J*_{0}=5 and the residual is *S*=0.1316 and *ϵ*=0.0026 (figure 2*c* shows the relative errors). There is, therefore, no essential difference between the reduced and the full versions of the three-dimensional WLC, which can be appreciated from the log plot in figure 2*b*. In contrast to figure 1, however, the full version has the smaller residual and RMSE. A better agreement between the theoretical model and the data is obtained than for the one-dimensional model. This is related, in particular, to the fact that since the one-dimensional tension functions proposed in order to give an explicit representation of the WLC were not derived in a rigorous way from the three-dimensional theory they do not contain the important factor (*λ*−*λ*^{−2}).

In figure 3, a comparison between the one-dimensional formula (2.5), the three-dimensional WLC (3.9) and the Gent formula in (3.11) is shown, with the previous data shifted from the range (0,1) to the range (1,2), corresponding to the stretch *λ*. Note that for *J*_{0}=5, the asymptotic value of the stretch calculated from (3.10) is *λ*^{*}=2. These plots enable the differences between the models to be appreciated. We note, in particular, that the Gent model is unable to fit the data and that the stress derived from the three-dimensional model gives the better results for the considered data. For (2.5) and (3.9), the same material constants were used as in figures 1 and 2, respectively, while for (3.11) we found the values *μ*=1.2872, *J*_{m}=1.9188, with residual *S*=580.83 and RMSE *ϵ*=0.1709.

### (b) Simple shear

Although simple shear deformations are not, *per se*, important in the context of biomolecular experiments, torsional deformations are important, and torsion is *locally* a simple shear. It is, therefore, of interest to investigate the response of the various models for a deformation that involves simple shear, which we do in this section. A simple shear deformation in the (*X*_{2}, *X*_{3}) plane is defined by(3.14)where *κ*>0 is the (constant) amount of shear (e.g. Ogden 1984). It follows that(3.15)and the squared principal stretches are computed as(3.16)

The shear stress is obtained from the specialization of (2.10) as(3.17)where the secant modulus of shear is defined as(3.18)The infinitesimal shear modulus is defined as(3.19)and we note that in the *linear* theory of incompressible isotropic elasticity this is the only constitutive parameter that arises.

In the case of (3.2), we have(3.20)and(3.21)

From (3.21), which, by (3.19) and the definition *J*_{0}=3*N*, defines the infinitesimal shear modulus *μ*(*J*_{0}) as a function of *J*_{0}∈[3,∞), we obtain and . These limits are surely unrealistic.

For the Gent model, we have(3.22)and therefore , a result that is independent of *J*_{m}∈[3,∞).

## 4. A WLC modification of the Gent model

Horgan & Saccomandi (2003) have shown that the model (2.17) proposed by Gent is indeed obtainable by a systematic procedure if *rational* approximations of the response functions for isotropic incompressible materials are considered. This result was obtained by generalizing the usual approximation method developed by Rivlin and others (e.g. Rivlin 1960) based on the expansion of the free energy as a truncated power series in the principal extensions *e*_{i}, *i*=1, 2, 3, which are related to the principal stretches by *λ _{i}*=1+

*e*

_{i}.

In terms of the principal extensions, the principal invariants (2.8) read(4.1)and the constraint is written as(4.2)In order to compute the power series it is convenient to introduce the equivalent set of invariants(4.3)from which it is easily shown that, when (4.2) holds,(4.4)Here, we have used the usual notation *F*=*O*(*a*) as *a*→0, if there exist constants *δ*>0 and *ϵ*>0 such that |*F*|<*δ*|*a*| whenever |*a*|<*ϵ*.

The method of Rivlin involves expanding the strain-energy function in the form(4.5)where *B*_{αβ} are suitable constitutive parameters and *α* and *β* may be integers or rational exponents. The technique proposed by Horgan & Saccomandi (2003) is based on approximating the response coefficients *W*_{1} and *W*_{2} in the representation formula (2.10) with the approximations based on rational functions instead of polynomials. The need for approximation is dictated by the fact that in experiments it is not the energy that is measured directly but quantities related to the first derivatives of the energy (forces or stress components). It is also, in general, necessary to consider supplementary conditions to ensure the existence of a strain-energy function, but such considerations are not needed here. We remark that it is well known that rational approximations may give more accurate results than polynomial approximations, in particular when singularities such as those associated with limiting chain extensibility models arise. Approximations using rational functions such as Padé approximants in this context have been described in Horgan & Saccomandi (2002*b*), where the relationships between the Gent model (2.17) and molecular models based on the inverse Langevin function are established. It is worth mentioning that in the recent review by Shaqfeh (2005) on the dynamics of the single-molecule DNA in flow, it was observed that the interpolation formula (2.3) is indeed a Padé approximant of the numerically simulated force obtained from the real WLC model.

### (a) A 0/4 approximation

Our analysis of the interpolation form of the WLC suggests that a stronger form of singularity than that reported by Gent needs to be considered. We therefore take(4.6)where *a*_{1} and *b*_{i}, *i*=1, 2, 3, are constants. This is a 0/4 approximation in the principal extensions. The constants in (4.6) are not arbitrary since they must satisfy some standard requirements of the nonlinear theory of elasticity. To reflect the finite chain length assumption a parameter *J*^{*} is chosen so that as *J*_{1}→*J*^{*} both the stress and the strain energy associated with (4.6) become singular. It is also desirable that, as *J*^{*}→∞, (4.6) and the associated strain energy should reduce to a standard polynomial form. Without loss of generality, the constants can be normalized by, for example, setting *b*_{1}=1. We then have the approximation(4.7)where *a* and *b* are the two remaining (arbitrary) parameters.

The response coefficient (4.7) is associated with the strain-energy function(4.8)

For this model, we have(4.9)and, recalling (3.19), we deduce the connection(4.10)Equation (4.8) may now be written(4.11)and for *J*^{*}→∞, we recover the neo-Hookean energy (2.18), now written(4.12)

In simple extension, the tensile stress calculated from (4.11) is(4.13)wherein *J*_{1}=*λ*^{2}+2*λ*^{−1}−3. The factor *J ^{*}–J*

_{1}has now been isolated to emphasize that the asymptote is fixed by the value of

*J*

^{*}and is unaffected by the value of

*b*. While

*μ*has the interpretation as the infinitesimal shear modulus, the role of

*b*is to control the rate at which the curves approach the asymptote and the flatness of the curve for small and moderate values of the stretch.

The fitting of this model to the data is considered in figure 4*a* as a log plot. Figure 4*b* shows the relative errors. We note that the infinitesimal shear modulus *μ* is very small (approximately 0.0052), while *b*=−0.4942 and *J*^{*}=1.9803. The residual is *S*=0.1456 and the RMSE is *ϵ*=0.0028, so that it is clear that the fit for (4.13) is better overall than for the WLC model.

## 5. An improvement of the WLC modification

It is well known that generalized neo-Hookean strain-energy functions of the form *W*=*W*(*I*_{1}) cannot describe completely the elastic properties of elastomeric materials. The criticisms of such models date back to Rivlin's investigation of natural rubber specimens (e.g. Saccomandi 2004), and have recently been reiterated by Horgan & Saccomandi (1999) and Wineman (2005) in considering the implications of a global universal relation for torsional deformations. Indeed, for generalized neo-Hookean materials the simplified representation formula for the Cauchy stress tensor involves only one response function, and this allows a new universal relation to be determined (i.e. a relation between the stress and strain components that is independent of the special form of the strain energy). This universal relation does not apply to a general incompressible isotropic energy function *W*=*W*(*I*_{1}, *I*_{2}). On comparison of experimental data with this special universal relation it can be deduced that the invariant *I*_{2} cannot be excluded if the energy function is to be regarded as realistic. For this reason, we examine here an *enthalpic* correction to the proposed strain energy involving *I*_{2}.

For simplicity, we consider the two examples(5.1)and(5.2)where is defined in (4.11) and, in each case, *C* is an additional material constant. The first of these two models represents a generalization of the classical Mooney–Rivlin model, which is linear in *I*_{2}, while the second is a generalization of the Gent–Thomas model (Gent & Thomas 1958). In Pucci & Saccomandi (2002), it was shown that for *b*=0 the model (5.2) allows the whole range of the experimental data in simple extension of a particular rubber to be described to an accuracy within 5%.

For simple tension, *I*_{2}=2*λ*+*λ*^{−2} and the stresses associated with (5.1) and (5.2) are, respectively(5.3)and(5.4)where is given by (4.13).

In figure 5, the stresses associated with (5.3) and (5.4) are fitted to the data on the same basis as for figure 4. Each model gives better results than (4.13). The values of the parameters are, for (5.3), *μ*=0.0020, *b*=−0.4977, *J*^{*}=1.9912, *C*=0.0730, and, for (5.4), *μ*=0.0007, *b*=−0.4991, *J*^{*}=1.9966, *C*=0.2880. The sums of residues are *S*=0.0179 and 0.0074, respectively; the RMS errors are *ϵ*=0.9998×10^{−3} and 0.6428×10^{−3}, respectively.

## 6. Concluding remarks

In this paper we have considered the development of a phenomenological version of the WLC model. This is of interest for several reasons. First, because the WLC models used in the literature are in fact interpolation formulae of the real WLC model and therefore have a parallel phenomenological status to the continuum mechanics models. Second, the WLC models are essentially one-dimensional, whereas the models discussed in the present paper are three-dimensional and without the drawbacks that arise in the direct translation of the WLC classical interpolation formulae to three dimensions. Our methodology enables us to find three-dimensional phenomenological models that admit the accurate fitting of the ‘data’ generated numerically from the real WLC model. To appreciate how such models work for real experimental data, we consider finally an application to data from an actual single molecule experiment.

In figure 6, we show a set of data taken from Baumann *et al*. (2000) compared with the predictions of the model (4.13). These data have also been used recently by Rosa *et al*. (2005) to test a new interpolation formula for semi-flexible polymers. The data have been re-scaled by dividing the abscissa values by the first abscissa value so that the stretch *λ* can be used for the abscissa. Some general comments are appropriate at this point. Before approaching the vertical asymptote the data are very flat, which means that a very small value of the infinitesimal shear modulus is required, as was also the case with the ‘data’ used in §5. Moreover, from the data we can estimate the value of the limiting chain parameter. This information suggests possible starting values of the fitting parameters. In particular, we have linked *J*^{*} to the location of the asymptote, while *b* is related to the flatness of the curve for *f*. This is an important issue since it is known that in nonlinear least squares fitting it is possible to have several optimal sets of parameters (Ogden *et al*. 2004) if random starting guesses in a wide range of values are considered. The parameters obtained here are *μ*=1.9466×10^{−5}, *b*=−0.4995, *J*^{*}=38.7247, the residual is *S*=0.8374 and the RMSE is *ϵ*=0.0587.

The approach that we have adopted allows the WLC model to be set within a three-dimensional continuum mechanical framework, and has been illustrated for some simple example models. Although there is clearly scope for refining the model the fit obtained using either of (5.3) or (5.4) is very satisfactory. We hope that our method will allow the possibility of describing accurately the mechanical properties of biological and other polymers with a small number of phenomenological parameters. Moreover, our model can be implemented within standard finite element codes for the analysis of, for example, torsional stresses in DNA loops (e.g. White & Bauer 2004).

## Acknowledgements

We thank Angelo Rosa for providing us with the digitized data used in §6. We are very grateful to Gerhard Holzapfel for his careful reading of the manuscript and for his valuable suggestions. The work of G.S. and I.S. was partially supported by GNFM of Italian INDAM and by PRIN2004 Modelli Matematici per la Dinamica del DNA. RWO is also grateful to GNFM of Italian INDAM for support.

## Footnotes

↵The original WLC model is due to Kratky & Porod (1949) and is often referred to as the Kratky–Porod chain.

- Received May 13, 2005.
- Accepted October 12, 2005.

- © 2005 The Royal Society