A modified formulation of quasi-linear viscoelasticity for transversely isotropic materials under finite deformation.

The theory of quasi-linear viscoelasticity (QLV) is modified and developed for transversely isotropic (TI) materials under finite deformation. For the first time, distinct relaxation responses are incorporated into an integral formulation of nonlinear viscoelasticity, according to the physical mode of deformation. The theory is consistent with linear viscoelasticity in the small strain limit and makes use of relaxation functions that can be determined from small-strain experiments, given the time/deformation separability assumption. After considering the general constitutive form applicable to compressible materials, attention is restricted to incompressible media. This enables a compact form for the constitutive relation to be derived, which is used to illustrate the behaviour of the model under three key deformations: uniaxial extension, transverse shear and longitudinal shear. Finally, it is demonstrated that the Poynting effect is present in TI, neo-Hookean, modified QLV materials under transverse shear, in contrast to neo-Hookean elastic materials subjected to the same deformation. Its presence is explained by the anisotropic relaxation response of the medium.

Introduction

The ability to predict time-dependent deformation in soft, compliant solids is important when modelling a diverse range of materials, such as reinforced polymers, elastomers and rubbers, and is of increasing importance in the context of soft tissue mechanics. In many of these applications, the microstructure of the medium in question dictates that the material response is both strongly anisotropic and viscoelastic. Additionally, given the compliant nature of these materials, it is necessary to accommodate finite strains. The field of finite strain nonlinear viscoelasticity has a rich history and a huge range of alternative constitutive forms has been proposed. The review by Wineman [1] provides a comprehensive overview of the current state of the art of the field. In this Introduction, we provide a summary of some details of existing models in order to provide context and to motivate the present study, specifically with respect to the study of viscoelastic anisotropy.

In the most general viscoelastic setting, the constitutive relation relating stress to deformation has the Cauchy stress T(t), where t is time, written in the general form [1] where is known as a response functional, F is the deformation gradient and C = FTF is the right Cauchy–Green deformation tensor. The stress will also, in general, be a function of space but for the sake of succinctness this argument is omitted.

In order to make progress, the form of the functional in (1.1) clearly has to be specified. A fading memory hypothesis is generally assumed. This intuitive imposition simply states that more recent deformations or stresses are more important than those from the past. The most straightforward finite strain viscoelastic constitutive models are those of differential type, assuming that the response functional is dependent on time derivatives of the right stretch tensor evaluated at the current time [2-4]; however, more generality regarding the relaxation behaviour of the medium can be incorporated via integral forms. Single integral forms can be employed, which are essentially an extension of Boltzmann's superposition principle to finite deformations. Although such a principle clearly does not hold exactly in a nonlinear setting, it can often provide a reasonable approximation. Furthermore, such an approach is often much more amenable to implementation than multiple integral forms [5-8]. Coleman & Noll [9] introduced finite linear viscoelasticity where the response functional is linear in the Green strain. A model that has gained traction since its introduction, especially in recent times due to its flexibility, is the single integral Pipkin–Rogers model [10], which, if deformation is considered to begin at t = 0, takes the form [11] with the first term being associated with the instantaneous elastic response and where for some potential . This model has the advantage of allowing for strong nonlinearity and finite deformation. It also incorporates coupling between relaxation and strain, where necessary, although important decisions regarding the dependence of the potential on the explicit time term t − s must be made, motivated by experiments.

The theory of quasi-linear viscoelasticity (QLV), whose original form was proposed by Fung [12,13], is a special case of (1.2). It incorporates finite strains and assumes that Q has a time/deformation separation so that where e = 2∂W/∂C is interpreted as the elastic second Piola–Kirchhoff stress and W is then the usual elastic strain energy function employed for finite elasticity problems. The notation ‘:’ indicates the double contraction between a fourth-order and a second-order tensor, such that in Cartesian coordinates. The term involving e forms an auxiliary measure of the strain in the medium. The tensor is a time-dependent, reduced (non-dimensional) relaxation function tensor. In [14], the isotropic theory of QLV was revisited and the theory reformulated. The authors pointed out that a number of concerns raised recently as regards the efficacy of QLV were, in fact, unfounded.

QLV offers an attractive approach to modelling nonlinear viscoelastic materials that can be implemented rather straightforwardly in computations. A particularly attractive aspect of QLV is that relaxation functions can be determined from experiments in the linear viscoelastic regime. Although this restricts the constitutive form to modelling materials that do not exhibit strain-dependent relaxation, it does immediately reduce the complexity of the model. Significant interest has focused on the case of transverse isotropy in the QLV context due to its importance in the application of modelling soft tissues. Small-strain QLV analyses have been conducted (see [15], for example), but the main focus for soft tissues has to be finite strains. Until now, the general trend in QLV theory has been to employ a scalar relaxation function implementation, i.e. letting , where G is a scalar function and is the fourth-order identity tensor. This has also generally been the case in the isotropic case, as was pointed out in [14]. As key examples in the transversely isotropic (TI) scalar relaxation function QLV context, Huyghe et al. [16] considered such an implementation for heart muscle tissue, Puso & Weiss [17] studied ligaments, Sahoo et al. [18] and Chatelin et al. [19] studied the brain, Motallebzadeh et al. [20], the eardrum and Jennesar et al. [21] focused on the spinal cord under tension. Vena et al. implemented incompressible QLV with a scalar relaxation function but with separate relaxation contributions from fibres and matrix [22]. It is important to stress that there is no reason to expect that the relaxation response of complex viscoelastic materials should be the same in all modes of deformation in general. Indeed, even in an isotropic scenario, the hydrostatic and deviatoric relaxation response are almost always very different in their nature, not only in their relaxation spectra but also in their functional form [23].

As pointed out by Weiss et al. [24], the in vivo distribution of stress and strain in ligaments and tendons is highly inhomogeneous. This is also true of many other deformed soft tissues in vivo. Simulations are important for many reasons, but in particular for simulating surgery [25]. An appropriate, fully three-dimensional constitutive model is therefore extremely important for accurate stress and deformation predictions. This is a fundamental motivation of the research carried out in the articles referred to above as well as of that presented here.

Work incorporating more than one relaxation function was carried out by Miller et al. who employed Ogden-type polynomial expansions for W, in a QLV framework where relaxation terms accommodate distinct relaxation times depending upon the order of the term in the expansion [26,27]. Since, however, these different orders are not associated with any specific physical deformations, they are incorporated purely to curve-fit to experiments. A large body of work that implicitly incorporates more than one relaxation function in anisotropic models is that associated with internal variable viscoelasticity theory, which was motivated by some of the earliest work on finite strain in viscoelastic isotropic solids [28-30]. This framework employs the uncoupled volumetric/deviatoric elasticity split dating back to Flory [31] and associates the time-dependent viscoelastic response to the deviatoric part only. The free energy function thus comprises a volumetric and isochoric elastic response, as well as a contribution due to configurational free energy associated with viscoelasticity. The decoupled stress then consists of equilibrium and non-equilibrium parts. The latter, which are described by evolving internal variables, dictate the viscoelastic response, and are governed by rate equations motivated by the linear theory [28]. This theory was driven forward by Holzapfel and colleagues, who extended the work to more realistic strain energy functions, beyond the Gaussian network theory [32] and to anisotropic solids [33,34], always providing highly detailed analysis of finite element implementations. Peña et al. developed this theory for TI materials in particular and applied it to ligaments and tendons [35]. More recently, anisotropic viscoelasticity based on internal variable theory has been employed to model the eye [36], and, more generally, in modelling soft tissues [37]. The internal variable approach has much in common with the QLV methodology, as was discussed in the isotropic case in [14], and in particular, it exhibits time/deformation separation. In the anisotropic setting, however, since distinct relaxation functions do not appear explicitly, it is non-trivial to link them to specific physical modes of deformation. Very recent work studied the time/deformation separability assumption with reference to experiments on filled rubbers and a wide range of models [38]. Interestingly, the internal variable model of Simo [30], which is equivalent to QLV for isotropic, incompressible materials appears to fit data fairly well across a variety of experiments. Apparently, no such experimental data are yet available for anisotropic materials.

The detailed discussion above regarding the state of the art on viscoelastic anisotropic theories motivates the work developed in this paper. A TI theory is proposed here, based on QLV, with a focus on the utility of the model, particularly in the incompressible regime, since this is a scenario of great practical importance. A tensor basis is employed for the relaxation tensor , which, in the incompressible limit, accommodates four relaxation functions. The proposed model is then used to predict the stress response in three common deformation modes: uni-axial extension along the fibres, and transverse and longitudinal shear. The results are presented for a specified strain energy function and specific relaxation functions in §5. Moreover, we compare the proposed model with a standard isotropic QLV model to highlight the importance of including more than one relaxation function when modelling TI soft tissues. We then use the TI QLV model to predict the Poynting effect in TI materials. Our results give new insights on the role played by viscoelasticity in determining the Poynting effect for such materials. Conclusions are drawn in §6.

Linear viscoelasticity

Anisotropic linear elastic materials are characterized by their tensor of elastic moduli with components in Cartesian coordinates. This tensor possesses the symmetries and the tensor relates elastic stress to strain in the form . In the case of TI materials, where the axis of anisotropy is in the direction of the unit vector M, a classical form is the following (used extensively in the context of fibre-reinforced materials [39]): where I∥ = M⊗M, ϵ∥ = M · M,  = (M)⊗M + M⊗(M) and λ, α, β μ and μ are the elastic constants. A tensor basis can be employed in order to write in a convenient form. This choice is non-unique, but for the form (2.1), the modulus tensor can be written as where j1 = λ, j2 = j3 = α, j4 = β, j5 = μ, j6 = μ and the basis tensors are given in appendix A. The choice of basis is motivated principally by the modes of deformation of interest and by which set of elastic moduli one wishes to work with.

In the incompressible limit, and such that (2.1) becomes where we note that the term αϵ∥ in (2.1) can be incorporated into the Lagrange multiplier term −p in (2.3). This means that there are only three independent elastic moduli for an incompressible, TI, linear elastic medium. Physically, this is explained by the fact that the restriction of zero volume change is not a purely hydrostatic condition (as is the case in an isotropic material)—the requirement means that the extension along the axis of anisotropy, M here, must also be constrained.

The constitutive form for anisotropic linear viscoelasticity theory for small strains, assuming fading memory and Boltzmann's superposition principle, takes the form [40] where the fourth-order relaxation tensor has the symmetries , but importantly, it does not, in general, possess the major symmetry [41]. In the context of transverse isotropy, the tensor can be written in terms of TI tensor bases in the form where, with reference to (2.1), R(t) are time-dependent relaxation functions (with dimensions of stress), chosen such that and where λ∞, α∞, β∞, μ and μ are the long-time moduli corresponding to λ, α, β, μ and μ, respectively. We note that these restrictions require R2(t) and R3(t) to be equal in the instantaneous and long-time limits; however, in general, they may relax at different rates and so cannot be assumed to be equal for all values of t. This further distinguishes the viscoelastic theory, which therefore has six independent relaxation functions, from the elastic theory, which only has five independent elastic constants.

A common choice is to let the R(t) take the form of Prony series. A one-term Prony series for R1(t), for example, would be given by where τ1 is the relaxation time associated with λ. The explicit form of (2.4) with (2.5) is and in the incompressible limit this becomes Next, with a view to the development of a modified quasi-linear theory of viscoelasticity for TI materials in the large deformation regime, i.e. in the form of (1.3), let us write the linear TI viscoelastic constitutive equation in the form noting that we have now written e, the elastic (instantaneous) stress introduced in (2.1), under the integral and therefore introduced the reduced (non-dimensional) relaxation tensor . The choice of basis employed for in (2.11) is important since the time derivative of the instantaneous elastic stress is present under the integral and so in the elastic limit one should recover the pure elastic stress. To see this, we assume that deformation begins at t = 0, and therefore integrate (2.11) by parts to obtain In order for us to recover the correct elastic limit at t = 0 we must impose the condition , recalling that is the fourth-order identity tensor, with components in Cartesian coordinates. Assuming a tensor basis decomposition for of the form for some new tensor basis , with G(0) = 1 for all n, this means that we must have The basis introduced above does not have this property, but a basis for transverse isotropy that does is given in (A 12) and (A 13) of appendix A. Using (2.13) and (2.1) in (2.11), and equating the resulting expression to (2.4), yields the connections between R and G, which are given explicitly in (A 21)–(A 23). Finally, the basis decomposes a second-order tensor as follows: where the explicit forms of the terms e are stated in (A 16)–(A 19). Expression (2.11) is now employed as the basis for a modified QLV theory for nonlinear materials subject to finite deformations. In this regime, appropriate stress measures that account for finite strains have to be used, as shall now be discussed.

Modified quasi-linear viscoelasticity

The viscoelasticity theory described above is now extended in order to deal with materials that are subject to finite deformation and whose constitutive response is nonlinear. In particular, a modified QLV theory is developed where relaxation is independent of deformation. Initially, the theory will be developed for general, compressible materials before the incompressible limit is taken. This yields a relatively compact constitutive model for incompressible TI viscoelastic materials that is suitable for use in computational models and for further development to model more complex materials.

General compressible form

We begin by defining and as the position vectors that identify a point of the body in the initial configuration (at t = 0) and current configurations of the body, and , respectively. The deformation gradient F(t) is defined as The left (right) Cauchy strain tensor is defined as B = FFT (C = FTF) and its three isotropic invariants are given by Let M and m = FM be the vectors along the principal axis of anisotropy of the TI material in question in the undeformed and deformed configurations, respectively. The anisotropy could be associated with, for example, the direction of the axis of aligned fibres in a medium, the direction orthogonal to parallel layers, or something more complex. Transverse isotropy requires an additional two anisotropic invariants By assuming the existence of a strain energy function W(I) with , the elastic Cauchy stress can be written in the general form [42] where the subscript i denotes differentiation with respect to the ith strain invariant. The notation Te distinguishes this nonlinear form from its linear counterpart (2.1), which we call e.

QLV requires a constitutive model of a form similar to (2.11); however, the theory cannot be directly formulated in terms of the Cauchy stress since objectivity must be ensured [14]. Instead, (2.11) can be written with respect to the second Piola–Kirchhoff stress tensor (t) and its elastic counterpart e = JF−1TeF−T. Upon integrating by parts, the constitutive equation for QLV (see equation (1.3)) can be written as where is a fourth-order reduced relaxation function tensor. If we were to split in terms of fundamental bases, then, in the isotropic case, this would correspond to splitting the elastic second Piola–Kirchhoff stress tensor into its hydrostatic and deviatoric parts, which do not have a clear physical interpretation. Therefore, instead of this approach, we choose to apply the split to the Cauchy stress in the following manner: which, by using the TI bases introduced in (A 12) and (A 13) amounts to writing where G(t) are the components of the relaxation function introduced in §2, such that G(0) = 1. This modified version of the QLV theory still preserves the property of the relaxation functions being independent of the deformation; however, the bases introduced in (A 12) and (A 13) do depend on the deformation through the vector m = FM. The terms e are given by where Te are given by equations (A 16) and (A 17), but with the components σe replaced by Te. We note that equation (3.6) cannot be written in the form of equation (3.5) for any choice of and, therefore, this constitutive expression is not the classical QLV approach; however, due to the similarities between the forms, we use the term modified QLV to describe this expression. It should be noted that this comment also applies to the compressible isotropic theory developed in [14].

The components G(t) are related to the functions R(t) (which are associated with the natural split of linear TI viscoelasticity in (2.9)) via equations (A 21)–(A 23). Note that the relaxation functions G(t) are independent of deformation, which is a fundamental and important assumption of both QLV and modified QLV and it means that the constitutive equation is restricted to materials for which this is a good approximation. It does mean, however, that such functions can be measured directly by small-strain tests on the medium in question, which is an attractive aspect of the theory, as we shall discuss later.

The connections (A 21)–(A 23) can now be employed to replace G with R in (3.7) and the Cauchy stress can then be written as: where the e, n = 1, 2, …, 6 are given in (A 26). Equation (3.9) can be used to model compressible materials. In the next section, we shall consider the incompressible limit, as this is of great utility in a number of important applications, including polymer composites and soft tissues.

The incompressible limit

We now have a model that can accommodate fully compressible TI behaviour in the large deformation regime. In order to illustrate the applicability of the model, let us consider the important case of incompressibility. The incompressible limit is recovered by setting and . The elastic constitutive equation in (3.4) in this limit is Moreover, from equations (A 24) and (A 25), where EL is the longitudinal Young modulus. For details on how to derive the last equality in (3.11), we refer to [43] and references therein. The incompressible limit of equation (3.9) then becomes where the Lagrange multiplier p(t) is given by The first integral in (3.12) vanishes if we assume that in the incompressible limit the time-dependence of this relaxation function becomes instantaneous so that R1(t) = λ for all t, and therefore R1′(t) = 0. Moreover, in (3.12), the relaxation function and the following terms have been introduced: where, as we shall show in the next section, is associated with relaxation in the direction of the axis of anisotropy m, where E is the axial Young's modulus. The subscripts T and A on e are associated with in-plane (transverse) shear in the plane of isotropy and anti-plane (longitudinal) shear, respectively. We note that R4(t), R5(t) and R6(t) all appear in (2.10), but R2(t) does not. However, the elastic constant associated with R2(t), α in the incompressible limit reduces to . We therefore take ∀t and we use the following constitutive equation: Under the assumptions mentioned above, the stress can be written in terms of three relaxation functions: and R6(t). These can be determined independently via three linear viscoelastic tests associated with uniaxial loading, in-plane shear (figure 1b) and anti-plane shear (figure 1c), respectively. This is an appeal of the model in the sense that the viscoelastic behaviour can be fully characterized by experiments in the linear viscoelastic regime (in the fully compressible case, three more experiments would be required to determine R1(t), R2(t) and R3(t)). Of course, in reality, not all materials respond viscoelastically in this manner and the relaxation functions can depend on strain amplitude, in principle. For now, we accept that the model cannot incorporate this effect but emphasize that not only is it capable of incorporating large deformations, but that it can also accommodate distinct relaxation behaviours associated with anisotropy, as opposed to the vast majority of existing models. We note further that most of these existing models also do not incorporate strain-dependent relaxation.

Figure 1.

(a) A block of a TI material with fibres pointing in the direction of the E3-axis, (b) under the simple transverse shear deformation in (4.17) and (c) under the simple longitudinal shear deformation in (4.24).

Linkage of relaxation functions to small-strain deformation modes

In this section, we briefly summarize how, for an incompressible material, the three relaxation functions can be associated with three independent small-strain regime tests: a simple extension in the direction of the axis of anisotropy and two shear deformations: in-plane and anti-plane shear. A range of small strain time-dependent experiments then permit the experimental determination of these relaxation functions and their associated relaxation spectra.

Take the E3-axis to be the axis of anisotropy. We first consider a scenario where a sample of TI material is stretched along this axis. The infinitesimal strain tensor is then  = ϵ33e3⊗e3, where ϵ33 is the strain component along the e3 axis. The associated stress response from (2.10) becomes We note that equation (3.16) depends on three relaxation functions (since ); therefore, from this test we are only able to obtain information about the composite relaxation function , which is associated with uniaxial deformation. To explicitly determine all three relaxation functions, two further tests are required. An in-plane shear test can be carried out in the plane of isotropy. This deformation is associated with a strain tensor of the form  = ϵ12(e1⊗e2 + e2⊗e1), where ϵ12 is the amount of shear in the isotropic plane. The resulting stress response is therefore, from an in-plane shear test we can extract the parameters appearing in the relaxation function R5(t). The third and last test required is a shear deformation in one of the two planes containing the axis of anisotropy (either e1-e3 or e2-e3). Upon choosing the plane e1–e3, the strain tensor can be written as  = ϵ13(e1⊗e3 + e3⊗e1), where ϵ13 is the amount of shear in the e1–e3 plane. The shear stress response, then, is given by which allows us to determine the parameters appearing in R6(t). Once , R5(t) and R6(t) are known, they can be used to calculate R4(t) from the first equation of (3.14). As mentioned above, to use the compressible theory, three additional, similar experiments would need to be carried out in order to determine the relaxation functions R1(t), R2(t) and R3(t).

Deformation

Let us now consider specific deformations of the medium in question. As above, let us take the axis of anisotropy to be M = E3, so that for a fibre-reinforced composite, for example, the fibres are all aligned along the E3-direction in the undeformed configuration. All deformations begin at t = 0 and we use the notation X1, X2, X3 and x1(t), x2(t), x3(t) for Cartesian coordinates in the undeformed and deformed configurations, and , respectively.

In all cases, it is assumed that the deformations are slow enough that inertial terms can be neglected (quasi-static assumption) and, therefore, since the deformations considered are homogeneous, they automatically satisfy the equations of motion: where ρ is the mass density in the deformed configuration. As pointed out in [44], in the context of finite visco-elasticity the existence of the quasi-static limit needs to be supported by a global existence result. To justify the quasi-static assumption, we derive a dimensionless form of equation (4.1) and we define a corresponding set of non-dimensional parameters which enables us to identify the quasi-static regime. We then consider three homogeneous deformations and calculate the quasi-static solution to equation (4.1) for each deformation mode.

By following the analysis carried in [45], let us first define the norm ∥f(y, s)∥ of a bounded function f defined on a set as follows: We then set: where T* = Tn is the traction specified on the boundary of the deformed body with normal n, and u* is the displacement specified on the boundary . Let L be the length-scale of the body, we can then define the non-dimensional quantities: Finally, we need a non-dimensional measure of the time variable t. For a viscoelastic material, there are two types of characteristic times, the external time tex imposed by the boundary condition on u* and the internal time tin, associated with the intrinsic relaxation time of the material. The two time scales associated with the characteristic times tex and tin are independent, we can therefore assume . Upon assuming that the relaxation functions R(t) take the form of a one-term Prony series, for a compressible TI material with constitutive equation in (3.9) we can identify six different relaxation times, each associated with a function R(t) (n = {1, …, 6}). An incompressible TI material with constitutive behaviour described by equation (3.12) will then have three relaxation times, namely and τ6. We therefore define where is the strain rate. For the sake of simplicity and clarity, we restrict the attention to the case . As will be shown at the end of this section, the results will apply to the general case as well. Equation (4.1) then becomes where For the quasi-static approximation to be valid, we shall require Q1, Q2 and Q3 to be small. However, we note that , therefore a necessary condition for the quasi-static assumption to be valid is The above equation is satisfied when for instance, the test is very slow, the tested sample is small and the internal relaxation time is long. However, as pointed out in [45] the condition (4.7) is not sufficient because the term and can in principle still be large. The problem is indeed singular for small times. This implies that there will always exists a small interval in the vicinity of t = 0 where the quasi-static approximation is not valid. This interval can be determined by re-scaling equation (4.5) with the new time-variable defined as follows: so that the contributions of the terms and are now . By combining equations (4.4) and (4.8), the initial time interval where inertial effects cannot be neglected is therefore given by or equivalently: We note that the upper bound of this time interval is independent on both and τ. The result in equation (4.9) is thus valid for the general case where . In conclusion, the quasi-static assumption is justified and therefore the proposed model can be used for . Furthermore, when designing an experimental test, one should aim at reducing as much as possible the upper bound in equation (4.9), to reduce the initial time where inertial effects are not negligible.

Let us now consider three specific deformation modes: simple extension in the fibres direction, in-plane and anti-plane shear.

Uniaxial deformation

First, consider the following deformation, which is associated with uniaxial deformation under tension with no lateral traction. Assume that the deformation begins at t = 0 and for t≥0, we have This corresponds to a simple extension, with stretch Λ(t), in the direction of the axis of anisotropy. If anisotropy is induced by fibres aligned along the E3 direction, equation (4.10) represents a simple extension in the direction of the fibres. Assuming that this deformation has been generated by a non-zero axial stress T33 with the lateral surfaces being free of traction, the stress state is assumed to take the form The deformation gradient associated with equation (4.10) is diagonal for t > 0, which gives rise to the following left Cauchy–Green strain tensor From equation (3.15), the stress T(t) is given by The Lagrange multiplier p(t) can be calculated from the second equation of (4.11) which yields p(t) = 0. Furthermore, where the term is defined analogously to in equation (A 18) and can be calculated by combining equations (3.10) and (4.13). Upon substituting equation (4.15) into equation (4.14), we obtain We note that the above equation depends only on R4(t), R5(t) and R6(t), as was the case for its linear counterpart in (3.16) and, in the small-strain limit, (4.16) becomes identical to (3.16).

In-plane (transverse) shear

Let us now consider a homogeneous simple shear deformation in the plane of isotropy e1–e2, as sketched in figure 1b. This type of shear is often called transverse shear [46,47]. The deformation is written as where κ(t) is the amount of shear. In this case, the deformation gradient F(t) and the left Cauchy–Green strain tensor B(t) are given by A traction-free boundary condition is imposed on the surface with normal N = (0, 0, 1), which implies that The non-zero components of the Cauchy stress tensor from equation (3.9) are then where the Lagrange multiplier p(t) can be calculated from the last equation of (4.19) as follows: Moreover, we have

Finally, substituting equations (4.21) and (4.22) into equation (4.20), we obtain the stress components T11, T22 and T12

The corresponding elastic stresses Te11, Te22, Te12 and can be calculated from equations (3.10) and (A 18). We note that the shear stress T12(t) only depends on the relaxation function R5(t) as in the linear regime (see equation (3.17)). In the small-strain limit, T11 and T22 in (4.23) tend to zero and the equation for T12 becomes identically equal to (3.17).

Anti-plane (longitudinal) shear

Finally, let us consider longitudinal shear—a simple shear along the fibre direction as depicted in figure 1c. This deformation can be written in the following form: where κ3(t) is the amount of shear in the e1–e3 plane. The deformation gradient and the left Cauchy–Green strain tensors associated with this anti-plane shear deformation are A traction-free boundary condition is imposed on the lateral surface with normal N = {0, 1, 0} in the undeformed configuration, which leads to Therefore, the non-zero components of the Cauchy stress in (3.15) are where the Lagrange multiplier p(t) can be calculated from the traction free condition T22 = 0: The terms along with (4.28), can be substituted into (4.27) to obtain the three stress components:

The shear stress in (4.30) depends on two relaxation functions, unlike its linear counterpart in (3.18), which depends only on R6(t); however, in the small-strain limit, T11 and T33 in (4.30) tend to zero and the equation for T13 becomes identically equal to (3.18). In the next section, we illustrate some key properties of the proposed model.

Key features of the modified transversely isotropic quasi-linearviscoelasticity model

A common procedure for investigating the time-dependent behaviour of a material is to perform a step-and-hold test. Usually, either the deformation or the load can be imposed on a sample of the material being tested. In the former case, the test is called a stress relaxation test, in the latter, a creep test. In this section, we shall focus on stress relaxation tests. In a stress relaxation test, the strain in the sample is increased very rapidly up to a maximum value, and is then held fixed during the holding phase. The response of the material is generally recorded in terms of forces or moments. Another type of test is the ramp-and-hold test, in which the sample is deformed over a finite time interval and is then, as with the step-and-hold test, held fixed for a prescribed time period. Some experiments also account for a recovery phase—a phase where the sample is allowed to return to its original state. Examples of these two types of tests are illustrated in figure 2. The step-and-hold test is impossible to achieve experimentally because the time-dependence of the strain is required to have the form of a Heaviside function; however, this test is very useful for theoretical purposes, especially for comparing the stress relaxation responses of different models. The ramp-and-hold test provides additional information when one is interested in studying the recovery behaviour of the sample as well as its relaxation behaviour. In the next section, these two types of test will be used to illustrate the main features of the proposed modified TI QLV model.

Figure 2.

Input time-dependence of the strain in (a) a step-and-hold stress relaxation test and (b) a ramp-and-hold relaxation test. Note that the strain in (a) is often modelled by a Heaviside function; however, in practice, a step-change in strain cannot be achieved experimentally as this would require an infinite strain rate; therefore, instead, a very rapid strain rate is used, as illustrated.

Comparison between linear viscoelasticity and modified transversely isotropicquasi-linear viscoelasticity

We first show that the modified QLV model proposed in (3.15) is equivalent to the linear model (2.10) in the small-strain regime for the three modes of deformation described in §4. To proceed, it is necessary to choose a specific form for the strain energy function W. We shall choose W to be of the form: where This expression was proposed for modelling the nonlinear behaviour of TI incompressible soft tissues [46], and it is consistent with the linear theory in the small strain limit. We note that the isotropic contribution to the strain energy takes the form of a Mooney–Rivlin function, with non-dimensional parameter αMR. We further assume that the reduced relaxation functions take the form of classical one-term Prony series: with E, μ and μ indicating the long-time analogues of E, μ and μ, respectively, and , τ5 and τ6 being the associated relaxation times. We note that, by assuming that takes the form of a one-term Prony series, the relaxation function will, in general, be a three-term Prony series. Alternatively, we could have chosen R4(t) to be a one-term Prony series, which would have led to be a three-term series; however, none of the results that follow would have changed qualitatively had we instead chosen to make that assumption.

We now consider the three modes of deformation illustrated in the previous section (uni-axial extension, transverse shear and longitudinal shear) and compare the stress responses predicted by the proposed model to those of its linear counterpart derived in §3c. We consider a step-and-hold test where the strain consists of a rapid ramp (of 0.1 s) followed by a holding phase (of 9.9 s), as depicted in figure 2a. Figure 3 shows the results in both the small (0.5%)-strain and the large (50%)-strain regimes. We note that in the small-strain regime the predictions of the two models are in agreement in all the three modes of deformation, thus, the modified TI QLV model is able to recover the linear limit correctly. Indeed, the shear stresses σ12(t) and σ13(t) in (3.17) and (3.18), respectively, are recovered by taking the limit . Moreover, the stress σ33(t) in (3.16) can be obtained by considering . In the large-strain regime, the results for the two models vary considerably. In particular, by comparing figure 3b,c with 3e,f , we observe that, although the shear stress responses are the same for both models, there is a large discrepancy for the normal stresses. The linear TI model always predict zero normal stresses for all t both in the small- and the large-strain regimes, whereas the modified TI QLV predicts non-zero normal stresses. This feature of the modified TI QLV model will be further analysed in the following sections.

Figure 3.

Comparison between the linear TI model (Lin VE) (2.10) and the proposed modified TI QLV model (Mod QLV) (3.15) of the resulting stresses induced for the three modes of deformation: (a,d) uni-axial extension, (b,e) transverse shear and (c,f) longitudinal shear. Stress responses to step-and-hold tests in the small-strain regime are plotted in (a–c) and in the large-strain regime in (d–f). The curves are obtained by setting: E/E = 0.3, , μ/μ = 0.9, τ5 = 2, μ/μ = 0.8 , τ6 = 1.5, E/μ = 75 and μ/μ = 5. All stresses shown are normalized by μ.

Comparison between modified isotropic quasi-linear viscoelasticity and modified transversely isotropic quasi-linear viscoelasticity

We now consider a step-and-hold test in transverse shear and calculate the stress response predicted by the proposed model. We then implement the model proposed by De Pascalis et al. [14], where a single relaxation function is used to account for the behaviour of an incompressible, isotropic material, and compare the two results in order to highlight the main differences between the models.

We consider the strain to consist of a rapid ramp (of 0.1 s) where the amount of shear κ increases up to 0.5 followed by a holding phase (of 9.9 s), as depicted in figure 2a. We then calculate the normal (T22) and shear (T12) stress components from equation (4.23). In figure 4, the results for the two models are presented. The stress components are normalized on their respective maximum values at 0.1 s. We recall that the modified isotropic QLV model [14] can be recovered from the proposed model by setting R2(t) = R3(t) = R4(t) = 0, R6(t) = R5(t) ∀ t, α = β = 0 and μ = μ. For the modified TI QLV model, we have set E∞/E≪μ∞/μ so that the function relaxes considerably more than R5(t).

Figure 4.

The predictions of the proposed modified TI QLV model (blue) and the modified isotropic QLV model from [14] (yellow) for a material under transverse shear. The blue curves have been obtained by setting αMR = 0.25, E/E = 0.3, , μ/μ = 0.9 , τ5 = 2 and E/μ = 75. The yellow curves have been obtained from equation (3.12) by setting E/E = μ/μ = 0.9, , α = 0, R2(t) = 0, ∀t and E/μ = 3. The latter equation is a result of the assumption that the isotropic material is incompressible. The restriction R2(t) = 0, ∀t follows from the fact that for isotropic materials α = 0.

Although the predictions of the two models are in agreement for the shear stress response, their normal stress responses differ. For the modified TI QLV model, the relaxation curve of T22 is mostly determined by the integrals associated with , whereas, for the modified isotropic QLV model, the relaxation behaviour is entirely dictated by the function R5(t), as shown in figure 4, and the normalized normal and shear stress relaxation behaviours are identical. This is an important and unique property displayed by the proposed model, and it arises from the presence of a tensorial relaxation function with distinct components. In general, the stress response predicted by the proposed modified TI QLV model will depend on the competition between the integrals associated with each of the three relaxation functions , R5(t) and R6(t). The modified isotropic QLV model lacks this property and so do all QLV models that incorporate a single scalar relaxation function. It is therefore of utmost importance to include more than a single relaxation function when modelling TI materials that exhibit direction-dependent relaxation behaviours. Otherwise, the risk of running into errors when measuring the mechanical parameters can dramatically increase.

In the next section, we will show that the modified TI QLV model predicts the so-called Poynting effect. This is a commonly observed phenomenon in nonlinear elastic materials undergoing simple shear or torsion. Such materials tend to expand or contract in the direction perpendicular to the shear direction. To prevent such an expansion or contraction, a normal stress is required. For the transverse shear deformation illustrated in §4b, for instance, T22 can be either positive or negative when the material is sheared. The general convention is that a negative (positive) normal stress is associated with the positive (negative) Poynting effect.

The Poynting effect

It is well known that nonlinear elastic materials display the Poynting effect. Since 1909, when the phenomenon was first discovered by Poynting [48], many studies have focused on theoretical and experimental aspects of this effect. In particular, within the class of incompressible, hyperelastic materials under simple shear, the problem has been studied for isotropic [49,50] and TI media [47,51]. Specifically, for isotropic materials, it has been shown that no Poynting effect occurs in neo-Hookean materials. Similarly, fibre-reinforced neo-Hookean materials exhibit no Poynting effect in transverse shear. Interestingly, however, the Poynting effect in viscoelastic materials seems to have received little attention. This section highlights how the proposed modified TI QLV model is capable of giving new theoretical insights into the Poynting effect in viscoelastic TI materials.

We now consider a ramp-and-hold test as in figure 2b for the transverse shear deformation in equation (4.17). The amount of shear κ increases up to 0.5 in 10 s, is then held constant for 10 s before decreasing back to 0 in 10 s and being held there for a further 10s. We calculate the corresponding normal stress response T22 from (4.23). Upon normalizing the stress on μ, the parameters appearing in equation (4.23) are the mechanical parameters E/E, μ/μ and αMR, and the relaxation times, and τ5. In figure 5, we plot curves for the stress T22/μ for different values of αMR. The Poynting effect decreases with increasing αMR, but is still non-zero even for αMR = 1. We recall that αMR = 1 is associated with a neo-Hookean TI material, whereas αMR = 0 indicates a pure Mooney–Rivlin TI material.

Figure 5.

The normal stress response of the modified TI QLV material under transverse shear. The normal stress T22/μ is plotted against the time t for different values of α = {0, 0.25, 0.5, 0.75, 1}, for fixed E/E = 0.66, μ/μ = 0.9, and τ5 = 2.

As discussed at the beginning of the section, a hyperelastic model with a strain energy function given by (5.1) predicts a normal stress Te22 = 0 for fibre-reinforced neo-Hookean materials under simple shear, indicating no Poynting effect. In contrast to the hyperelastic model, the modified TI QLV model does predict a (positive) Poynting effect, as illustrated by the purple curve (αMR = 1) in figure 5. This interesting behaviour suggested by the numerical results will of course need to be verified by bespoke experiments; however, the curves depicted in figure 6 remarkably indicate that the Poynting effect is dictated by the competition between the two relaxation functions appearing in the equation for T22 in (4.23). We further note that when and R2(t) = R3(t) = 0 for all t, i.e. in the isotropic case, the Poynting effect vanishes (green, dotted curve) in agreement with the elastic behaviour of neo-Hookean TI materials. Finally, we note that when E∞/E > μ/μ (purple curve in figure 6), the proposed model predicts a negative Poynting effect. This behaviour might be displayed by materials with soft fibres, i.e. where the matrix is stiffer than the embedded fibres and therefore the relaxation effects due to the matrix are stronger than those arising from the fibres.

Figure 6.

The Poynting effect predicted by the modified TI QLV model. The normal stress component T22/μ is plotted against time t for a neo-Hookean viscoelastic TI material (αMR = 1). The following parameters have been set: , while E/E spans over {0.1, 0.5, 0.8, 0.9, 0.95}.

Conclusion

In this paper, a modified TI QLV theory for finite deformations has been developed. Transverse isotropy is accommodated both in terms of elastic anisotropy and relaxation functions, thus improving on existing scalar relaxation function TI QLV models. The numerical results presented in §5 have shown that incorporating distinct relaxation functions is crucial when modelling TI materials, and that simplified models with only one relaxation function would fail to capture the Poynting effect, for example. Another appeal of the proposed model is that the relaxation functions can be determined from small-strain mechanical tests. Moreover, the formulation in terms of tensor bases motivates similar analyses for other important viscoelastic anisotropies, such as orthotropy. Finally, the theory developed here can be used as a starting point for more complex, fully three-dimensional nonlinear viscoelastic theories that are able to incorporate strain-dependent relaxation.