## Abstract

This paper offers a reappraisal of Fung's model for quasi-linear viscoelasticity. It is shown that a number of negative features exhibited in other works, commonly attributed to the Fung approach, are merely a consequence of the way it has been applied. The approach outlined herein is shown to yield improved behaviour and offers a straightforward scheme for solving a wide range of models. Results from the new model are contrasted with those in the literature for the case of uniaxial elongation of a bar: for an imposed stretch of an incompressible bar and for an imposed load. In the latter case, a numerical solution to a Volterra integral equation is required to obtain the results. This is achieved by a high-order discretization scheme. Finally, the stretch of a compressible viscoelastic bar is determined for two distinct materials: Horgan–Murphy and Gent.

## 1. Introduction

The study of nonlinear viscoelastic deformations of solid materials has a very long history, with a consequent proliferation of a diverse and extensive array of constitutive models. For a comprehensive overview of the topic, the reader is directed to the recent paper by Wineman [1] and related references therein. The subject is still active, with models continuing to be developed across the field, from highly mathematical approaches where implementation is not a concern to very applied studies where ease of application is essential. Each approach has its advantages: the models that more accurately describe the microphysics tend to prove difficult to employ in practical engineering or biomedical situations, whereas simpler approaches often miss crucial details. If the viscous model assumes a fading memory effect of the strain history, then this usually gives rise mathematically to Volterra-type integral equations, which in the linear case can generally be solved either analytically or numerically. However, in the nonlinear case (of interest to describe finite deformations), these problems are difficult to solve by any methods, even finite-element approaches. It is therefore crucial to develop constitutive models that are simple enough to be amenable to straightforward (and rapid) solution methods, yet include enough detail to capture the important physics underlying these relevant materials. One such ‘compromise’ approach was offered by Fung [2], in order to study the uniaxial elongation of biological soft tissue. His constitutive assumption, often called *quasi-linear viscoelasticity* (QLV) or *Fung's model* of viscoelasticity, assumes that the viscous relaxation rate is independent of the instantaneous local strain. This model, which is a special case of a more general Pipkin–Rogers constitutive model [3], predicts that at any time a stress that is equal to the instantaneous elastic stress response decreased by an amount depending on the past history, assuming that a Boltzmann superposition principle holds.

Fung's QLV model is perhaps the most widely used today, but has met with some criticism mainly because it has been suggested that it does not always yield physically ‘reasonable’ behaviour. Of course, as with any constitutive model for a complex nonlinear material, QLV has limitations. However, it can be expected to be appropriate for materials whose relaxation-rate coefficients are weakly dependent on deformation, or where the deformation comprises small perturbations about a large initial deformation. Various experimental results have appeared in the literature that appear to confirm that QLV is, in practice, a reasonable model for a range of materials; see for example [4–6].

The apparent deficiency of QLV will be discussed below, but the purpose of this paper is to reappraise Fung's approach, suggesting how it can be reformulated against existing interpretations in the literature. It will be seen that the form offered herein has similarities to other viscoelastic formulations, including that by Simo [7], and most importantly exhibits behaviour that makes it both (relatively) easy to use and physically useful.

Analysis will commence with the general tensorial form of (QLV) proposed by Fung [2]
**G**(*t*) is, in comparison with linear theory, the stress relaxation second-order tensor, **Π** denotes the second Piola–Kirchhoff stress, while **Π**^{e} is an instantaneous strain measure. The latter may be thought of as an equivalent (instantaneous) *elastic* stress, hence the superscript ‘e’. This tensorial integral identity is the natural generalization of the simple one-dimensional relationship proposed by Fung [2], which preserves *objectivity*.

As mentioned, thanks to the relative simplicity of Fung's approach over more general nonlinear viscoelastic models, it has proved extremely popular. This is especially true in the biomechanics field where it has been employed to predict the large deformations of soft tissues. However, despite the widespread use of QLV, published studies (especially in the incompressible limit) appear to interpret (1.1) in different ways, which then lead to quite different results. In the paragraphs below, a number of common approaches are discussed, where the interpretations of the model may be questioned. The main aim of this paper is then to re-derive QLV, starting from basic principles and also to ensure consistency in the limit of infinitesimal deformations, thus recovering the Boltzmann theory of linear viscoelasticity.

The simplest interpretation of Fung's relation is to assume a purely one-dimensional homogeneous deformation. This is attractive for its ease of use and has been employed extensively in biomedical applications [8–10], for example to model the shearing deformation of brain tissue [6]. However, in practical applications, especially for incompressible or near incompressible materials purely one-dimensional deformations, such as pure uniaxial extension, will not be realizable. Thus, a tensorial form of QLV must be employed. As mentioned, this relation must, at the minimum, satisfy objectivity, and hence the form chosen in (1.1). Nevertheless, a number of authors have erroneously expressed Fung's relation not for the second Piola–Kirchhoff stress, which guarantees objectivity, but in terms of other stress measures. For example, see Fung's original discussion on the subject in [2, page 253], in which the QLV relation is expressed in terms of the Kirchhoff stress! In the latest version [6.12] of ABAQUS, although the viscoelastic constitutive model is objective, the model is rather heuristic and as we shall see later, it assumes the same relaxation behaviour in both the shear and bulk parts of an isotropic material. Further technical aspects of objectivity are discussed in the paper by Liu [11].

Another approach to QLV sometimes employed in the literature (e.g. [12–14]) is, in an analogous fashion to nonlinear elasticity, to write the stress in terms of the instantaneous derivative (with respect to the principal stretch) of a strain energy function. Clearly, in this case the strain energy function must be dissipative and depend on the history of the strain, and so these authors express it as a fading-memory integral with integrand given as a hyperelastic strain energy function. This approach may be somewhat hard to justify and does not allow the user to consider the problem in terms of an auxiliary instantaneous measure of strain, i.e. the effective elastic stress term **Π**^{e} given above.

Perhaps, a more subtle point to those mentioned above is the choice of **Π**^{e} in (1.1). In this paper, the reasonable assertion is made that **Π**^{e} must be zero whenever the deformation is zero. However, this point is often not recognized, especially for incompressible materials, as can be seen, for example, by inspection of the integrands (setting the stretch equal to unity) in the articles by [15] (see eqns (8) and (12) therein), [16] (eqns (9) and (13)) and [17] (eqn (9)). Note that this requirement is satisfied in the one-dimensional models in [8,9] but not in [10]. If **Π**^{e} does not vanish for zero deformation then, in general, it can be expected that the actual stress field **Π** will be non-zero prior to the application of the imposed load or stretch, i.e. the solution will be non-causal! However, for incompressible materials, an arbitrary pressure term (Lagrange multiplier) can be adjusted so that the stress field is indeed causal, but then it must take a specific form at later times. That is, whenever the material is deformed, and then returns back to zero, the stress **Π** will instantaneously return to zero too. This is clearly an unrealistic consequence of the particular choice of model, for, as Fung states [2, page 228]: ‘the tensile stress at any time *t* is equal to the instantaneous [elastic] stress response … decreased by an amount depending on the past history’. This issue is examined further in §4*a*.

There are other variants of QLV employed in the literature that offer slightly modified forms of the governing equation to that suggested by Fung. The reader is referred, for example, to articles [18–20]. In this paper, the method employed follows closely that proposed by Fung but does not exhibit any of the limitations just discussed.

The paper is organized as follows. In the following section, the usual (Boltzmann) linear viscoelastic model is reviewed by the way of introduction to a new interpretation of Fung's QLV theory, presented in §3. In §4, results from the new model are contrasted with those in the literature for the case of uniaxial elongation of a bar: in §4*a* for an imposed stretch of an incompressible bar, in §4*b* for an imposed load, and in §4*c* for stretch of a compressible bar. A numerical solution to a Volterra integral equation is required for the results in §4*b*, and the discretization (time-stepping) procedure used is described in appendix A. Finally, concluding remarks are offered in §5.

## 2. Boltzmann's linear viscoelastic law

It is helpful to commence analysis by recapping the theory of linear viscoelasticity. Under the assumption of isotropy, infinitesimal elastic deformations can be described by the constitutive law
** σ** and

**are the second-order stress and strain tensors, respectively, and**

*ϵ**μ*and

*κ*represent the infinitesimal shear modulus and modulus of compression (or bulk modulus), respectively. To incorporate viscoelastic behaviour, the most natural extension of (2.1) is to assume that the stress

*remembers*the past history of the rate of strain with some

*fading memory*, and then apply the

*superposition principle*(or Boltzmann's principle); hence,

*μ*(

*t*),

*κ*(

*t*) are now time-dependent

*relaxation*functions, and the lower limit of the integral must be taken from

*t*=0, yields

Now, in the elastic case (2.1), incompressibility may be considered as the dual limit of *p*(*t*) may be considered as a Lagrange multiplier of incompressibility. For the viscoelastic counterpart, the limit of incompressibility can be considered in a similar fashion, noting first that *μ*,*κ* in the elastic case), respectively. Thus, taking *t* owing to the fading memory assumption) and *p*(*t*). Thus, in the limit of incompressibility, equations (2.2) and (2.3) become

## 3. Quasi-linear viscoelasticity

When the strain is not infinitesimal, linear theory becomes inappropriate to describe deformations; hence, a nonlinear constitutive law has to be considered. As discussed in the Introduction, Fung's hypothesis (1.1) is examined here as a means of describing the motion of viscous nonlinearly elastic materials. In his renowned work [2], Fung introduced this quasi-linear constitutive model in order to capture the nonlinear stress–strain relationship of living tissues; however, it also has applicability to elastomeric materials.

Before deriving the quasi-linear theory, it is useful to introduce some standard definitions and notations. The deformation gradient tensor **F** is defined by
**x**(*s*) denoting the position of a generic particle *P* at time *s* ∈ [0,*t*], and **X** its position at the initial reference time. Note that the start time of the deformation, and any imposed tractions, will be taken as *t*=0. The quantity *J*=det **F**, expressing the local volume change, is a constant *J*=1 when the deformation is isochoric. Furthermore, from the deformation gradient tensor **F** the left Cauchy–Green tensor **B**=**F****F**^{T} is obtained, together with its principal invariants

Now, Fung [2] makes the assumption that the QLV stress depends linearly on the (superposed) time history of a related nonlinear elastic response (a nonlinear instantaneous measure of strain). This allows, for example, for incorporation of a finite hyperelastic theory in the analysis. In index notation, Fung's theory (1.1) can be written as^{1}
*G*_{ijkl} as a *reduced relaxation function tensor*. The ‘crucial’ simplification of Fung's theory is that this term is *independent* of the strain. Moreover, if the material is isotropic then **G**, a tensor of rank four, has just two independent components, being therefore consistent with linear theory.^{2}

Following the analysis of the previous section, for isotropic materials it is convenient to split the equivalent (instantaneous) *Cauchy stress* into two parts, one which accounts for *microscopic isochoric deformations* of the material and the other that measures purely *compressive deformations*. These two components can be expected to have different instantaneous elastic behaviours as well as distinct relaxation rates, associated for example with the unwinding/unravelling of polymeric fibres as opposed to their stretching. It is assumed that this decomposition can be achieved by taking the deviatoric and hydrostatic components of the equivalent elastic Cauchy stress:
*objectivity*. Instead, the second Piola–Kirchhoff stress tensor associated with (3.5) has to be introduced, which is defined by
**do not** refer to the deviatoric and hydrostatic parts of the second Piola–Kirchhoff stress, but correspond to the second Piola–Kirchhoff stress of the deviatoric and hydrostatic Cauchy stress components, respectively. Assuming a superposition principle as for the linear case, it is now possible to introduce an objective viscoelastic law, relating the second Piola–Kirchhoff stress to the past history of the nonlinear rate of strain measure. This is taken as
*J*^{−1}**F** and post-multiplying by **F**^{T} yield the Cauchy viscoelastic stress

As mentioned earlier, to allow for large (nonlinear) deformations, a hyperelastic theory can be employed assuming that the measure of the *effective elastic stress* **Π**^{e} is derived from an elastic potential. Let us then specialize equations (3.5)–(3.12), considering the existence of a strain energy function (SEF) *W*, say, which in the isotropic case is dependent on the principal invariants of the deformation *I*_{1},*I*_{2},*I*_{3}, or of the principal stretches λ_{1},λ_{2},λ_{3}:
*β*_{j}=*β*_{j}(*I*_{1},*I*_{2},*I*_{3}) are the so-called material (or elastic) response functions. In terms of the strain energy function, they are given by
**T**^{e},
**C**=**F**^{T}**F** is the right Cauchy–Green deformation tensor. The viscoelastic stress is obtained from (3.19) to (3.20) via (3.10) or for the Cauchy stress from (3.12).

Note that (as required) objectivity is now preserved for the viscoelastic model (3.11), but the first part is *not* deviatoric, and so the second term is not purely hydrostatic. In fact, in general both the compressive and shear components of the stress history contribute to the deviatoric and hydrostatic parts of **T**. However, the main point to note is that when there is no deformation, i.e. the principal stretches are unity, λ_{j}=1,*j*=1,2,3 and **B**≡**I**, then *I*_{1}=*I*_{2}=3,*I*_{3}=1, which reveals that the effective deviatoric elastic stress **T**^{e}_{D} vanishes. Similarly, **T**^{e}_{H} vanishes as the strain energy function has always to satisfy the additional relation (e.g. [22])
**Π**^{e}_{D} and **Π**^{e}_{H} vanish as *t*=0, justifies taking the integration range for the pair of integrals in (3.12) as 0 to *t*.

If it proves convenient to express the strain energy function in terms of the principal stretches **F** the quasi-linear viscoelastic stress relation (3.12) may be rewritten as

The final consideration of this section is the constraint of incompressibility for all possible deformations, i.e. *J*≡1. In this limit, the speed of propagation of compressive disturbances tends to infinity, and hence the relaxation time for viscous dilatational motions can be assumed to tend to zero. Therefore, the second term in (3.11) (or (3.12)), in an analogous fashion to that for linear elasticity (2.7), reduces to
*p*(*t*) can be considered as a Lagrange multiplier. Note that this expression in itself is not hydrostatic. It contains both a hydrostatic term, and a component that combines with the first term in (3.12) to make it deviatoric; hence,
*anisotropic* nonlinear viscoelastic materials rapidly becomes rather demanding, although strain measures for anisotropic QLV were discussed in [24]. It is expected that, despite the complexity, the present approach is extendable to certain classes of anisotropy.

## 4. Uniaxial elongation

One of the simplest, useful and hence most popular experiments to measure the properties of a material is to subject a specimen to a simple elongation test. There is a huge literature of results on such a deformation, and it is therefore useful to examine this homogeneous solution, comparing the present result with extant published work. To appreciate the general QLV model developed herein, and the approach needed to obtain a solution, three specific cases will be examined. First, in §4*a* a uniaxial stretch is imposed on the specimen, with the stress determined as a relaxing function of the history of the stretch. For comparison, both Yeoh and Mooney–Rivlin incompressible strain energy functions are considered. In §4*b*, a tensile load is imposed, and the stretch history (creep) determined from this. In this case, the integral equation must be solved numerically, which is the focus of discussion in appendix A. Finally, in §4*c*, the procedure in §4*a* is repeated for a compressible material, with the effect of compressibility on the viscoelastic behaviour highlighted.

### (a) Simple extension

In the case of simple extension, assumed homogeneous (in space) and incompressible, the principal stretches may be specified as
*X*_{1},*X*_{2},*X*_{3}) and (*x*_{1},*x*_{2},*x*_{3}) are the Cartesian coordinates in the undeformed and deformed state, respectively, and the stresses are
*F*_{ij}=∂*x*_{i}/∂*X*_{j} is of diagonal form
*J*=1, gives a relationship between the stretches λ_{1} and λ_{2}, in particular _{1}=λ, the deformation gradient tensor **F** and the left Cauchy–Green tensor **B**=**F****F**^{T} become
*p* can be eliminated, and so it is possible to rewrite the stress–strain relationship (4.6) and (4.7) as

Note that equations (4.9) and (4.10) depend on the specific choice of strain energy function *W*, and so it is useful to examine two specific examples. Let us start by assuming the instantaneous response is modelled by a two-term Yeoh strain energy function [25]
*α* is a positive constant and *μ* is the shear modulus from infinitesimal theory. Thus,
*γ* is a constant in the range −1/2≤*γ*≤1/2. Note that (4.17) reduces to the neo-Hookean model when *γ*=1/2. The partial derivatives of *W* are

These two viscoelastic models can be compared by considering an ‘experiment’ where the stretch is imposed and the stress is measured. To specify matters, the stress relaxation function is chosen to be the classical one-term Prony series
*μ*, *τ* is the relaxation time. These are set as
*a*) is applied. The time variation is assumed slow so that inertial terms in the balance equations can be neglected. Figure 1*b* shows the resultant stress predictions *T*/*μ* for the Yeoh hyperelastic model (4.14) (dotted curves with *α*=1,2) and for the Mooney–Rivlin hyperelastic model (4.19) (dashed curves with *γ*=1/6,−1/3). It is interesting to observe that the two hyperelastic models depart from the solid curve, obtained when the strain energy function is of the neo-Hookean type, i.e. *α*=0,*γ*=1/2, in different ways. The Yeoh material is found to harden as the parameter *α* increases, whereas the effect of decreasing *γ* from 1/2 leads to a softening of the Mooney–Rivlin material's behaviour.

As mentioned in the Introduction, it is useful to contrast the results presented here with those discussed, for example, by Ciambella *et al.* [15]. In that article, the authors employed a form of QLV, which, with the incompressible Yeoh SEF, yields the stretch–stress law (see eqn (16) in [15]):
*et al.*'s article, their equation (4.23) was compared with that offered by the formulation used in the ABAQUS finite-element analysis (FEA) package (Version 6.7, see [26]), namely
*T*/*μ*. Employing the same stretch history as above (figure 1*a*), the relaxation function in (4.21) with parameters as in (4.22) and *α*=1.0, yields the curves in figure 1*c* predicted from (4.24) (ABAQUS: dotted), (4.23) (Ciambella *et al.*: dashed) and (4.14) (De Pascalis *et al.*: solid), respectively. An alternative way of representing this deformation is via a stretch–stress diagram (figure 1*d*). It is clear that the present QLV approach, illustrated by the solid curve in both figures, gives a result that lies remarkably close to that found from the ABAQUS model, whereas Ciambella *et al.*'s solution is quite distinct. The latter model was developed by the authors in order to address the deficiency they highlighted regarding v. 6.7 of ABAQUS. However, their result also exhibits a substantial limitation: it predicts instantaneous zero stress whenever λ recovers back to unity after some (arbitrary) deformation, thus showing no ‘memory’ of the history of the stress in this situation. This is at variance with physically observed behaviour, and in particular that found for linear viscoelasticity.

### (b) Simple tensile load

Rather than an imposed stretch, perhaps a more important experiment for measuring viscoelastic properties of materials is application of a simple tensile load: a slowly varying uniaxial stress is imposed on the body for which the resultant stretch is measured over time. For linear viscoelasticity, equation (2.2) can be inverted to yield a straightforward creep relation to determine the strain. However, for nonlinear viscoelasticity, the solution procedure is somewhat more difficult, as inversion is not possible. A typical QLV relation is given in (4.19), which may be considered in the general form
*f*_{j}(λ(*t*)) terms means that an inversion operator cannot be introduced. Instead, a numerical scheme must be employed that can evaluate the stretch as a function of time, subjected to a prescribed stress. A suitable numerical discretization procedure is offered in appendix A, which has error *O*(*δt*^{4}), where *δt* is the step size and so gives rapid convergence.

An example is illustrated in figure 2, with the imposed stress history shown on the left graph. The one-term Prony series relaxation function (4.21), with constant values as given in equation (4.22), is again employed. The resultant stretch is given in figure 2*b* for the Yeoh strain energy function (4.14) (dotted curves with *α*=1,2), and for the Mooney–Rivlin strain energy function (dashed curves with *α*=1,2). As before, the solid curve is the neo-Hookean prediction (*α*=0,*γ*=1/2) and it can be seen that increasing *α* in the Yeoh model leads to material hardening, while decreasing *γ* leads to softening of the Mooney–Rivlin material.

### (c) Simple extension of a compressible bar

The purpose of this section is to repeat that analysis in §4*a* for a compressible material. As before, a bar, with zero tractions on its lateral faces, is assumed to undergo simple homogeneous extension, caused by an imposed stress or stretch in the *X*_{1} direction, say. Then, the principal stretches may be given by
*X*_{2} and *X*_{3} directions is clear, and so the deformation gradient tensor is
**B** becomes

The above can be substituted into the governing viscoelastic equation (3.12), with effective elastic stresses (3.19) and (3.20), and simplified. The *X*_{1} component of the stress yields, after simplification,
*X*_{2} and *X*_{3} equations are both
*W*_{j} is the derivative of *W* with respect to the indicated invariant (3.16). Clearly, the particular choice of the strain energy function affects the results and so, as before, a couple of examples are presented. Assuming the Horgan–Murphy strain energy function [27], used for modelling a material with a small compressibility, the SEF is
*γ* is an arbitrary constant, *μ* is the infinitesimal shear modulus and *κ* is the infinitesimal bulk modulus. Then
*γ*=1/2, for ease of presentation, the stresses become
*J*_{m} is a constant limiting value for *I*_{1}−3, taking into account the finite extensibility of polymeric chains within the material. Hence,

## 5. Concluding remarks

This paper has focused on reappraising Fung's method for QLV. It has been shown that some of the negative features commonly associated with the approach are merely a consequence of the way it has been applied elsewhere. The approach outlined herein was shown in §4 to yield ‘sensible’ results and offers a straightforward approach in solving a wide range of models. The present method exhibits similarities with Simo's approach to nonlinear viscoelasticity [7], although the latter assumed a scalar relaxation function acting on the stress. Therefore, it is noted that Simo's method would not reduce in the linear limit to the usual relation (2.3) if the latter shear and bulk moduli have different temporal behaviour. (Note that ABAQUS v. 6.12 employs the same relaxation time constants in both the shear and bulk parts of the field.) Herein, the relaxation function is a tensor, which for isotropic materials reduces to two distinct scalar relaxation functions, one acting on the compressive part and the other on the shear component of the stress. This therefore is consistent with that found for linear theory.

It was shown that for an imposed deformation (e.g. stretch) it is simple matter to solve the QLV equation directly. For imposed stress, the stretch has to be deduced by solving (inverting) the integral equation, and this is achieved here using a discretized scheme accurate to an order of the cube of the time-step size. Muliana *et al.* [29] have recently offered a QLV model where the strain is expressed as a function of the stress, which may be viewed as a dual model to (1.1). However, it is not clear as yet whether their approach is as effective at modelling viscoelasticity as Fung's scheme.

The authors are currently applying the present approach to a range of problems in viscoelasticity. The numerical procedure described in appendix A can be employed in a straightforward manner for any deformation that is incompressible and homogeneous, for example simple shear. It is presently being extended so that it can used for compressible materials undergoing simple extension or shear. The QLV method is also adaptable to inhomogeneous problems, for example those studied in [30], large amplitude radial deformations in two and three dimensions, and simple torsion. Such models lead to partial differential nonlinear Volterra integral equations and extending the present numerical scheme to these should offer improvement over current methods. It is also anticipated that the present approach will be more suited to satisfying the equations of equilibrium than, say, the more heuristic QLV relation employed in ABAQUS (v. 6.12). Finally, the authors are using the model for several biomechanics problems, to derive a perturbation theory for the viscoelastic evolution of small deformations on the top of large ones and to derive effective properties of ligaments and other tissues.

## Acknowledgements

The authors are grateful to the Engineering and Physical Research Council for the award (grant no. EP/H050779/1) of a postdoctoral research assistantship for R.DeP.

## Appendix A. A numerical high-order solution procedure for Volterra integral equations

The purpose of this appendix is to offer a simple and rapid numerical scheme for solving Volterra integral equations of the type obtained by the QLV approach (see e.g. (4.14)). It is assumed that they have the separable form
,*h*_{j}(*X*) are all known functions of their respective arguments. Clearly, all the equations offered in §4*a* lie within this class but for more general problems, with inhomogeneous deformations, the analysis described below must be amended.

Let us now discretize the system and evaluate the equation at each time *t*_{n}=*nδt*, with *n*=0,1,2,3,…,*n*_{max}, where *δt* is a *small* time step. On writing the stretch at each time step as
*s*=*mδt*+*u* is used. The integrand functions can be further expanded as
*g* and *f*_{j} as
*A*_{n−1},*B*_{n−1},*C*_{n−1} are the following expressions:
*n* and at each level of approximation. Hence, keeping terms up to *O*(*δt*) in (A 10), and rearranging, gives
*O*(*δt*^{2}) and *O*(*δt*^{3})) to zero yields, respectively,

Now, the procedure for solving the discretized equations (A 10), for specified stress, is clear. Deformation commences at *t*=0 (or *n*=0) and so from (A 1):
*T*(0)=0, would mean that λ_{0}=1. The derivatives of λ_{0} are derived from (A 14) to (A 16) (setting *n*=1 in these equations) and hence λ_{1} is found from (A 4). Then *n* is incremented, the values of λ′,λ′′,λ′′′ determined, and these are again used to determine the stretch at the next time step. The process is repeated at each following step. Note that size of the sums in *A*_{n−1},*B*_{n−1},*C*_{n−1} grows with *n*, and so for large times, or using small time steps, evaluation of the system can become rather slow. However, this difficulty can be alleviated, with the present system, by obtaining a relationship between *A*_{n} and its previous increment *A*_{n−1}, and similarly for *B*_{n},*C*_{n}. The precise forms of these relationships are determined from the given *f*_{j}(λ) and the relaxation function *O*(*δt*^{4}), i.e. for a time step of *δt*=0.01 then the error is *O*(10^{−8}).

## Footnotes

↵1 Note, in order to clarify exposition, the space variable in (1.1) is omitted; the same will be done henceforth, as long as no confusion occurs.

↵2 This is not always taken in account, for example in [1,21] the authors introduce an isotropic model which contains only one scalar relaxation function

*G*.

- Received January 22, 2014.
- Accepted March 4, 2014.

© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.