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.

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 G is known as a response functional, F is the deformation gradient and C = F T F 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][3][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][6][7][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] T(t) = F(t) Q[C(t), 0] + with the first term being associated with the instantaneous elastic response and where Q = 2∂W/∂C for some potential W = W(C(s), t − s). 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 W 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 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 tensorG, 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 C with components C ijk in Cartesian coordinates. This tensor possesses the symmetries C ijk = C jik = C ij k = C k ij and the tensor relates elastic stress to strain in the form σ e ij = C ijk k . 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]): A tensor basis can be employed in order to write C in a convenient form. This choice is non-unique, but for the form (2.1), the modulus tensor can be written as where j 1 = λ, j 2 = j 3 = α, j 4 = β, j 5 = μ T , j 6 = μ L and the basis tensors J n , n = 1, 2, . . . , 6 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, tr → 0 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 R has the symmetries R ijk (t) = R jik (t) = R ij k (t), but importantly, it does not, in general, possess the major symmetry R ijk (t) = R k ij (t) [41]. In the context of transverse isotropy, the tensor R(t) can be written in terms of TI tensor bases in the form R n (t)J n , (2.5) where, with reference to (2.1), R n (t) are time-dependent relaxation functions (with dimensions of stress), chosen such that and lim where λ ∞ , α ∞ , β ∞ , μ T∞ and μ L∞ are the long-time moduli corresponding to λ, α, β, μ T and μ L , respectively. We note that these restrictions require R 2 (t) and R 3 (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 n (t) take the form of Prony series. A one-term Prony series for R 1 (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 (2.9) 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 To see this, we assume that deformation begins at t = 0, and therefore integrate (2.11) by parts to obtain : σ e (τ ) dτ . (2.12) In order for us to recover the correct elastic limit at t = 0 we must impose the condition G(0) = I, recalling that I is the fourth-order identity tensor, with components I ijk = (δ ik δ j + δ i δ jk )/2 in Cartesian coordinates. Assuming a tensor basis decomposition for G of the form for some new tensor basis K n , n = 1, 2, . . . , 6, with G n (0) = 1 for all n, this means that we must have 6 n=1 K n = I. (2.14) The basis J n 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 n and G n , which are given explicitly in (A 21)-(A 23). Finally, the basis K n decomposes a second-order tensor as follows: K n : σ e = σ e 1 + σ e 2 + σ e 3 + σ e 4 + σ e 5 + σ e 6 , (2.15) where the explicit forms of the terms σ e n 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.

(a) General compressible form
We begin by defining X = 3 i=1 X i E i and x(t) = 3 j=1 x j (t)e j as the position vectors that identify a point of the body in the initial configuration (at t = 0) and current configurations of the body, B 0 and B(t), respectively. The deformation gradient F(t) is defined as The left (right) Cauchy strain tensor is defined as B = FF T (C = F T F) and its three isotropic invariants are given by where the subscript i denotes differentiation with respect to the ith strain invariant. The notation T e 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 −1 T e F −T . Upon integrating by parts, the constitutive equation for QLV (see equation (1.3)) can be written as whereG is a fourth-order reduced relaxation function tensor. If we were to splitG 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 n (t) are the components of the relaxation function G introduced in §2, such that G n (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 n are given by where T e n are given by equations (A 16) and (A 17), but with the components σ e ij replaced by T e ij . We note that equation (3.6) cannot be written in the form of equation (3.5) for any choice ofG 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 n (t) are related to the functions R n (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 n (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 n with R n in (3.7) and the Cauchy stress can then be written as: where the P e n , 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.

(b) 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 J → 1 and λ → ∞. The elastic constitutive equation in (3.4) in this limit is (3.10) Moreover, from equations (A 24) and (A 25), where E L 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 timedependence of this relaxation function becomes instantaneous so that R 1 (t) = λ for all t, and therefore R 1 (t) = 0. Moreover, in (3.12), the relaxation function R(t) and the following terms have been introduced: ,  where, as we shall show in the next section, R(t) is associated with relaxation in the direction of the axis of anisotropy m, where E L 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 R 4 (t), R 5 (t) and R 6 (t) all appear in (2.10), but R 2 (t) does not. However, the elastic constant associated with R 2 (t), α in the incompressible limit reduces to lim λ→∞ α = μ T . We therefore take R 2 (t) = 1 2 R 5 (t) ∀t and we use the following constitutive equation: Under the assumptions mentioned above, the stress can be written in terms of three relaxation functions: R(t), R 5 (t) and R 6 (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 R 1 (t), R 2 (t) and R 3 (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.
(c) Linkage of relaxation functions to small-strain deformation modes In this section, we briefly summarize how, for an incompressible material, the three relaxation functions R(t), R 5 (t), R 6 (t) 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: inplane 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 E 3 -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 = 33 e 3 ⊗ e 3 , where 33 is the strain component along the e 3 axis. The associated stress response from (2.10) becomes (3. 16) We note that equation (3.16) depends on three relaxation functions (since R(t) = R 4 (t) − 1 2 R 5 (t) + 2R 6 (t)); therefore, from this test we are only able to obtain information about the composite relaxation function R(t − τ ), 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 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ), 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 R 5 (t). The third and last test required is a shear deformation in one of the two planes containing the axis of anisotropy (either e 1 -e 3 or e 2 -e 3 ). Upon choosing the plane e 1 -e 3 , the strain tensor can be written as = 13 (e 1 ⊗ e 3 + e 3 ⊗ e 1 ), where 13 is the amount of shear in the e 1 -e 3 plane. The shear stress response, then, is given by which allows us to determine the parameters appearing in R 6 (t). Once R(t), R 5 (t) and R 6 (t) are known, they can be used to calculate R 4 (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 R 1 (t), R 2 (t) and R 3 (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 = E 3 , so that for a fibre-reinforced composite, for example, the fibres are all aligned along the E 3 -direction in the undeformed configuration. All deformations begin at t = 0 and we use the notation X 1 , X 2 , X 3 and x 1 (t), x 2 (t), x 3 (t) for Cartesian coordinates in the undeformed and deformed configurations, B 0 and B(t), 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 U × [0, T ] as follows: f (y, s) = sup{sup{|f (y, s)| : y ∈ U} : s ∈ [0, T ]}.
We then set: where T * = Tn is the traction specified on the boundary ∂B t (t) of the deformed body with normal n, and u * is the displacement specified on the boundary ∂B u (t). 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 t ex imposed by the boundary condition on u * and the internal time t in , associated with the intrinsic relaxation time of the material. The two time scales associated with the characteristic times t ex and t in are independent, we can therefore assumex(t) =x(t ex , t in ). 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 n (t) (n = {1, . . . , 6}). An incompressible TI material with constitutive behaviour described by equation (3.12) will then have three relaxation times, namely τ R , τ 5 and τ 6 . We therefore define whereγ = V/L is the strain rate. For the sake of simplicity and clarity, we restrict the attention to the case t in = t R = t 5 = t 6 . As will be shown at the end of this section, the results will apply to the general case as well. Equation (4.1) then becomeŝ For the quasi-static approximation to be valid, we shall require Q 1 , Q 2 and Q 3 to be small. However, we note that Q 2 = Q 1 Q 3 , 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 ∂ 2x /∂t 2 ex and ∂ 2x /∂t 2 in 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-variablet defined as follows: so that the contributions of the terms ∂ 2x /∂t 2 ex and ∂ 2x /∂t 2 in are now O(1). By combining equations (4.4) and (4.8), the initial time interval where inertial effects cannot be neglected is therefore given byt ∈ [0, 1] or equivalently: We note that the upper bound of this time interval is independent on bothγ and τ n . The result in equation (4.9) is thus valid for the general case where t R = t 5 = t 6 . In conclusion, the quasi-static assumption is justified and therefore the proposed model can be used for t > L √ ρ/S. 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.

(a) 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 E 3 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 T 33 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 , Λ 2 (t) . (4.13)

From equation (3.15), the stress T(t) is given by
(4.14) The Lagrange multiplier p(t) can be calculated from the second equation of (4.11) which yields p(t) = 0. Furthermore, Π e L33 (τ ) =T e (τ ), (4.15) where the termT e is defined analogously toσ e 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 R 4 (t), R 5 (t) and R 6 (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).

(b) In-plane (transverse) shear
Let us now consider a homogeneous simple shear deformation in the plane of isotropy e 1 -e 2 , 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 where the Lagrange multiplier p(t) can be calculated from the last equation of (4.19) as follows: The corresponding elastic stresses T e 11 , T e 22 , T e 12 andT e , can be calculated from equations (3.10) and (A 18). We note that the shear stress T 12 (t) only depends on the relaxation function R 5 (t) as in the linear regime (see equation (3.17)). In the small-strain limit, T 11 and T 22 in (4.23) tend to zero and the equation for T 12 becomes identically equal to (3.17).

(c) 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: x 1 (t) = X 1 , x 2 (t) = X 2 and x 3 (t) = κ 3 (t)X 1 + X 3 , (4.24) where κ 3 (t) is the amount of shear in the e 1 -e 3 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 T 22 = 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 R 6 (t); however, in the small-strain limit, T 11 and T 33 in (4.30) tend to zero and the equation for T 13 becomes identically equal to (3.18). In the next section, we illustrate some key properties of the proposed 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.

Key features of the modified transversely isotropic quasi-linear viscoelasticity 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 timedependence 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.

(a) Comparison between linear viscoelasticity and modified transversely isotropic quasi-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 nondimensional parameter α MR . We further assume that the reduced relaxation functions take the form of classical one-term Prony series: with E L∞ , μ T∞ and μ L∞ indicating the long-time analogues of E L , μ T and μ L , respectively, and τ R , τ 5 and τ 6 being the associated relaxation times. We note that, by assuming that R(t) takes the form of a one-term Prony series, the relaxation function R 4 (t) = R(t) + 1 2 R 5 (t) − 2R 6 (t) will, in general, be a three-term Prony series. Alternatively, we could have chosen R 4 (t) to be a one-term Prony series, which would have led R(t) 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-andhold 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 T e → σ e . Moreover, the stress σ 33 (t) in (3.16) can be obtained by considering Λ → 1 + 33 . 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. T 12 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 L∞ /E L = 0.3, τ R = 1, μ T∞ /μ T = 0.9 , τ 5 = 2 and E L /μ T = 75. The yellow curves have been obtained from equation (3.12) by setting E L∞ /E L = μ T∞ /μ T = 0.9, τ R = τ 5 = 2, α = 0, R 2 (t) = 0, ∀t and E L /μ T = 3. The latter equation is a result of the assumption that the isotropic material is incompressible. The restriction R 2 (t) = 0, ∀t follows from the fact that for isotropic materials α = 0.

(b) 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 (T 22 ) and shear (T 12 ) 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 R 2 (t) = R 3 (t) = R 4 (t) = 0, R 6 (t) = R 5 (t) ∀ t, α = β = 0 and μ L = μ T . For the modified TI QLV model, we have set E L ∞/E L μ T ∞/μ T so that the function R(t) relaxes considerably more than R 5 (t).
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 T 22 is mostly determined by the integrals associated with R(t), whereas, for the modified isotropic QLV model, the relaxation behaviour is entirely dictated by the function R 5 (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 R(t), R 5 (t) and R 6 (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, T 22

(c) 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 T 22 from (4.23). Upon normalizing the stress on μ T , the parameters appearing in equation (4.23) are the mechanical parameters E L∞ /E L , μ T∞ /μ T and α MR , and the relaxation times, τ R and τ 5 . In figure 5, we plot curves for the stress T 22 /μ T 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.
As discussed at the beginning of the section, a hyperelastic model with a strain energy function given by (5.1) predicts a normal stress T e 22 = 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 T 22 in (4.23). We further note that when R(t)/E L = R 5 (t)/2μ T and R 2 (t) = R 3 (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 L ∞/E L > μ T∞ /μ T (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.

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.