On nonlinear viscoelastic deformations: a reappraisal of Fung's quasi-linear viscoelastic model.

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.

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] where 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 §4a.

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 §4a for an imposed stretch of an incompressible bar, in §4b for an imposed load, and in §4c for stretch of a compressible bar. A numerical solution to a Volterra integral equation is required for the results in §4b, and the discretization (time-stepping) procedure used is described in appendix A. Finally, concluding remarks are offered in §5.

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 where 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, where μ(t), κ(t) are now time-dependent relaxation functions, and the lower limit of the integral must be taken from in order to correctly consider the initial deformation. Integrating (2.2) by parts, and assuming that the deformation commences at t=0, yields where the ′ denotes differentiation with respect to the argument of the function. Note that this expression incorporates any jump discontinuity when the motion starts. Moreover, the first term in (2.2) (or equivalently the first two terms in (2.3)) is trace free, or deviatoric, and accounts for shear deformations and loss, while the second term accounts for the hydrostatic part representing compressive deformations and loss.

Now, in the elastic case (2.1), incompressibility may be considered as the dual limit of and , which results in the stress–strain relationship i.e. the second term in (2.1) yields a finite non-zero limit, where 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 and are the long time shear and bulk moduli (the equivalent of μ,κ in the elastic case), respectively. Thus, taking (which implies for all t owing to the fading memory assumption) and in (2.2) (or in (2.3)), it is found that where the limit is again assumed to have a non-zero finite value, −p(t). Thus, in the limit of incompressibility, equations (2.2) and (2.3) become and respectively.

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 with 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=FFT is obtained, together with its principal invariants which, alternatively, can be expressed in terms of the principal stretches through

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] Fung refers to G 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: which can be expressed as The split of equation (3.5) into shear and dilatational components is the nonlinear (hyperelastic) analogue of linear stress–strain law (2.1). However, it cannot be generalized to viscoelasticity, as in (2.2), because it would not preserve objectivity. Instead, the second Piola–Kirchhoff stress tensor associated with (3.5) has to be introduced, which is defined by with in which It must be emphasized that the subscripts D and H 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 where now and are two scalar (independent)-reduced relaxation functions (with without loss of generality). The latter relaxation functions relate to the inherent viscous processes involved with instantaneous shear and compressional (volumetric) deformations, respectively. Clearly, if the material was anisotropic, then a more complex tensorial relaxation function would be required. Furthermore, pre-multiplying by J−1F and post-multiplying by FT yield the Cauchy viscoelastic stress and integrating by parts, following (2.7), gives

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 I1,I2,I3, or of the principal stretches λ1,λ2,λ3: The general form of elastic Cauchy stress may be written (e.g. [22,23]) as where β=β(I1,I2,I3) are the so-called material (or elastic) response functions. In terms of the strain energy function, they are given by where It is straightforward to calculate the trace of Te, and so the deviatoric and hydrostatic elastic Cauchy stress components in (3.6) become From this, the second Piola–Kirchhoff counterparts, (3.9), are given by and where it is recalled that C=FTF 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, λ=1,j=1,2,3 and B≡I, then I1=I2=3,I3=1, which reveals that the effective deviatoric elastic stress TeD vanishes. Similarly, TeH vanishes as the strain energy function has always to satisfy the additional relation (e.g. [22]) Thus, both terms ΠeD and ΠeH vanish as , and hence, as there is no applied stress until the initial time 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 , then for diagonal F the quasi-linear viscoelastic stress relation (3.12) may be rewritten as in which where now refers to the derivative .

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 where 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, where now from the first equation in (3.19), As described above, the present theory has been developed for isotropic materials. The study of 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.

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 §4a 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 §4b, 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 §4c, the procedure in §4a is repeated for a compressible material, with the effect of compressibility on the viscoelastic behaviour highlighted.

Simple extension

In the case of simple extension, assumed homogeneous (in space) and incompressible, the principal stretches may be specified as where (X1,X2,X3) and (x1,x2,x3) are the Cartesian coordinates in the undeformed and deformed state, respectively, and the stresses are having assumed that the lateral surfaces are free of stress. The deformation gradient F=∂x/∂X is of diagonal form and the constraint of incompressibility, J=1, gives a relationship between the stretches λ1 and λ2, in particular . Setting λ1=λ, the deformation gradient tensor F and the left Cauchy–Green tensor B=FFT become The principal invariants are therefore and hence taking the diagonal terms of equation (3.25) yields the principal Cauchy stresses and respectively. Thus, by subtraction, the Lagrange multiplier p can be eliminated, and so it is possible to rewrite the stress–strain relationship (4.6) and (4.7) as where from (3.26) and

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] which yields the stress–strain relationship where α is a positive constant and μ is the shear modulus from infinitesimal theory. Thus, It is straightforward to obtain from (4.8) the relation where The second example is the instantaneous response modelled by a Mooney–Rivlin strain energy function which yields the stress–strain relationship where γ 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 and so (4.8) becomes, after simplification, in which

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 where μ, are the infinitesimal shear modulus and the long-time infinitesimal shear modulus, respectively, and τ is the relaxation time. These are set as A dynamic imposed stretch history (shown in figure 1a) is applied. The time variation is assumed slow so that inertial terms in the balance equations can be neglected. Figure 1b 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.

Figure 1.

The imposed stretch history is shown in graph (a). The resultant dimensionless stress T/μ is plotted in graph (b), found for the Yeoh model (4.14) (dotted) and for the Mooney–Rivlin material (4.19) (dashed) where the solid curve is the neo-Hookean limit α=0 (or γ=1/2). (c) T/μ is plotted from the predictions of (4.24) (dotted), (4.23) (dashed) and (4.14) (solid), respectively. (d) The dimensionless stress, T/μ, is plotted against stretch, λ, from the predictions of (4.24) (dotted), (4.23) (dashed) and (4.14) (solid).

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]): Similar equations can be derived from the constitutive models in [16] or [17]. In Ciambella 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 The three alternative viscoelastic expressions (4.24), (4.23) and the present result (4.14) may be conveniently compared by again imposing a stretch and determining the resultant dimensionless stress T/μ. Employing the same stretch history as above (figure 1a), the relaxation function in (4.21) with parameters as in (4.22) and α=1.0, yields the curves in figure 1c 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 1d). 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.

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 Although the integrand is of separable type, the presence of the various f(λ(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(δt4), 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 2b 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.

Figure 2.

Plot of the dimensionless stress history T/μ (a) and the resultant stretch λ (b), from the Yeoh model predictions (4.14) (dotted) and that from Mooney–Rivlin predictions (4.19) (dashed). The solid curve is the neo-Hookean limit α=0 (or γ=1/2).

Simple extension of a compressible bar

The purpose of this section is to repeat that analysis in §4a 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 X1 direction, say. Then, the principal stretches may be given by where symmetry between the X2 and X3 directions is clear, and so the deformation gradient tensor is and the left Cauchy–Green tensor B becomes The principal invariants are therefore

The above can be substituted into the governing viscoelastic equation (3.12), with effective elastic stresses (3.19) and (3.20), and simplified. The X1 component of the stress yields, after simplification, and the normal stresses in the X2 and X3 equations are both where as before W 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 where γ is an arbitrary constant, μ is the infinitesimal shear modulus and κ is the infinitesimal bulk modulus. Then and which are substituted in (4.30) and (4.31). Setting γ=1/2, for ease of presentation, the stresses become and respectively. By contrast, for the Gent model (e.g. [28]), where J is a constant limiting value for I1−3, taking into account the finite extensibility of polymeric chains within the material. Hence, which can then substituted into (4.30) and (4.31) to yield the equivalent results for the stresses. These are omitted for brevity.

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.