On the thermodynamic consistency of Quasi-Linear Viscoelastic models for soft solids

Originating in the field of biomechanics, Fung's model of quasi-linear viscoelasticity (QLV) is one of the most popular constitutive theories employed to compute the time-dependent relationship between stress and deformation in soft solids. It is one of the simplest models of nonlinear viscoelasticity, based on a time-domain integral formulation. In the present study, we consider the QLV model incorporating a single scalar relaxation function. We provide natural internal variables of state, as well as a consistent expression of the free energy to illustrate the thermodynamic consistency of this version of the QLV model. The thermodynamic formulation highlights striking similarities between QLV and the internal-variable models introduced by Holzapfel and Simo. Finally, the dissipative features of compressible QLV materials are illustrated in simple tension.


Introduction
Nonlinear viscoelastic behaviour is observed in many soft solids, such as elastomers, gels, and biological materials.More specifically, the nonlinear mechanical response of viscous soft solids exhibits relaxation and creep phenomena in large deformation quasi-static tests.In dynamic tests, the response of such materials is sensitive to strain rates.Moreover, marked hysteresis loops in loading-unloading experiments are evident in such media [1,2].
Historically, the mechanical modelling of nonlinear viscoelastic solids has been approached in many different ways.Fundamentally, however, modelling viscoelastic effects amounts to a specification of the constitutive law to provide an accurate functional relationship between the instantaneous stress and the entire strain history [3,4].In the monograph by Truesdell and Noll [5], this definition corresponds to the concept of simple materials (Sec.29 therein), see also Sec. 6.7 of the book by Malvern [6].
If the stress depends only on a very short interval of the recent history of the deformation, then it can be expressed as a function of the time derivatives of the deformation gradient up to a finite or-der (cf.Truesdell and Noll [5] Sec.35).Often, time derivatives of the strain up to first order are considered (i.e., the stress is expressed in terms of strain and strain rates), which leads to Newtonian-type viscosity models [7,8].Similarly to the linear Kelvin-Voigt model, nonlinear strain-rate differential models fail to describe stress relaxation phenomena [9,10].
Differential models involving time derivatives of stress as well as strain provide a natural bridge to models involving longer durations of time history.Indeed, the stress can then be expressed via an integral representation or via models with memory variables.These models can be seen as more or less equivalent, with memory variables manifesting themselves as terms within a stress relaxation function that arises in integral models.There is a vast literature for both integral and memory variable approaches.Reviews of viscoelastic constitutive modelling contrast these different approaches [3,4,1], or introduce them in a disconnected way [2,11].
A popular and simple model of the integral type is Fung's quasi-linear viscoelasticity (QLV) [12,13] sor and the elastic stress response.Originating in biomechanics, Fung's QLV has been successfully employed in related applications, with the development of associated experimental techniques [14,15,16] and dedicated computational methods [17,18].It has also been employed to model polymers and rubbers [19,20].
The computation of the instantaneous QLV stress by means of the convolution product requires the storage of the whole strain history, since the constitutive law is non-local in time.In computational applications, this can be avoided by expanding the relaxation tensor as a Prony series, leading to the natural definition of memory variables [17,18].Suitable memory variables account for the deformation history in a time-local fashion [21], allowing efficient evaluation of the mechanical response.
As described extensively by Maugin [22,23], constitutive models with memory variables rely on thermodynamics with internal variables of state to ensure that the addition of new variables entails a dissipative contribution.Various viscoelastic constitutive models of the literature satisfy these principles by design [1].In particular, viscous strain variables may be consistently introduced within a multiplicative decomposition of the deformation gradient tensor, in combination with Maxwell-type rheologies [24,25].
To the present authors' knowledge, a thermodynamically-consistent expression of the energy in terms of the memory variables of QLV is yet to be presented.In this paper, we recall the equations governing the motion of QLV solids with a single scalar relaxation function (Sec.2).The main result is derived in Sec. 3, namely a thermodynamically consistent expression of the free energy with appropriate internal variables of state.This way, we establish links between Fung's QLV and the internal-variable models by Simo [26], and Holzapfel and Simo [27].
Using a neo-Hookean model, the dissipative features of compressible QLV solids are illustrated in Sec. 4. The incompressible case is addressed in A.

Preliminaries
In what follows, we present the basic equations of Lagrangian solid dynamics [1,28].We consider a homogeneous and isotropic solid continuum on which no external volume force is applied.A particle initially at position X in the reference configuration moves to position x in the current configuration.The deformation gradient tensor is the second-order tensor where u = x − X is the displacement field, I is the identity tensor, and Grad denotes the gradient operator with respect to the material coordinates X (Lagrangian gradient).If the Euclidean space is described by an orthonormal basis {e 1 , e 2 , e 3 } and a Cartesian coordinate system, then I = [δ i j ], where δ i j is the Kronecker delta.The volume dilatation equals the ratio ρ 0 /ρ of the mass densities in the reference (undeformed) and deformed configurations.One can define various strain tensors as functions of F , such as the right Cauchy-Green deformation tensor C = F F , and the Green-Lagrange strain tensor E = 1 2 (C − I ).Frequently, principal stretches λ i are introduced, whose squares λ 2 i correspond to the eigenvalues of C .Thus, the principal invariants I i of C are given by ( The dynamics of the continuum in question are governed by conservation of momentum, which involves the divergence of a stress tensor.Typically, in the Lagrangian description, this equation of motion involves the first Piola-Kirchhoff stress tensor P , and the Eulerian version involves the Cauchy stress tensor σ = J −1 P F .Specified by the constitutive law, these stress measures are also related to the second Piola-Kirchhoff stress tensor S = F −1 P .

Fung's quasi-linear viscoelasticity
Fung's quasi-linear viscoelasticity (QLV) is presented below (see Sec. 7.13 of [12]).This model is based on the assumption that the stress is linearly dependent on the history of the elastic stress response, and a Boltzmann superposition principle between both quantities is assumed.The second Piola-Kirchhoff stress is given by where the elastic response [13] is derived from a strain energy density function W (I 1 , I 2 , I 3 ), and G is a fourth-order relaxation tensor.Here the colon denotes the double contraction where Einstein notation is used and the dot denotes the material time derivative.The notation W i is shorthand for the derivative ∂W /∂I i .
Also, in a similar fashion to Taylor et al. [21], we introduce the Flory decomposition of the deformation into volumetric and deviatoric parts (see Holzapfel [1] Sec. 6.4).Thus, we introduce the volumepreserving Cauchy-Green strain tensor C = J −2/3 C and its volume-changing counterpart J 2/3 I .We perform the change of variable W = W ( Ĩ1 , Ĩ2 , J ) in the expression of the strain energy, where The invariants Ĩi in Eq. ( 6) describe volumepreserving deformation, while the dilatation J describes volume-changing deformation.The third invariant Ĩ3 of C is equal to one, as deduced from the definitions in Eq. (3).
To compute the elastic response (5) in terms of the new variables ( Ĩ1 , Ĩ2 , J ), let us recall expressions for the , where the fourth-order unit tensor is . Thus, we introduce the decomposition where Dev(•) = (•) −1 3 (• : C )C −1 denotes the deviatoric operator in the Lagrangian description [1].The expression of Se is deduced from Eq. ( 5) with W3 = 0. Converting back to the variables (I 1 , I 2 , I 3 ), the chain rule yields which are the same expressions as in De Pascalis et al. [13] (Eqs.(3.19)-(3.20)therein).In a standard fashion [1], Taylor et al. [21] assumes the separability of isochoric and volumetric deformations W = W iso ( Ĩ1 , Ĩ2 ) + W vol (J ), which is a particular case of the present expressions.
In a similar fashion to related works [3,21], we assume that the relaxation is the same in all directions, i.e.G = G I s where we have defined the fourth-order symmetric identity tensor and G is a scalar function.Thus, the constitutive law (4) yields Fung [12] initially proposed an integral expression of G with a continuous spectrum of relaxation.This expression leads to high computational costs as it requires the storage of the whole deformation history [9].This drawback can be avoided by approximating G as an exponential series [17,18] 1 of the form with an arbitrary number n of relaxation mechanisms [21,15].The Heaviside step function H is included in Eq. ( 10) for convenience.Such a Prony series with magnitudes g k > 0 and characteristic relaxation times τ k > 0 can be linked to generalised Maxwell-type rheologies.

Generalities
The consequences of the first and second principles of thermodynamics are summarized below.We con-sider deformable solids whose associated constitutive law involves n second-order tensorial internal variables of state α 1 , . . ., α n [22,23].

Isentropic modelling
We consider the set of variables of state {η, E , α 1 , . . ., α n }, where η denotes the entropy per unit mass.The state is assumed local in time, i.e., only its instantaneous value is considered.The first principle of thermodynamics is reflected in the conservation of energy ρ ė = σ : D, where ė is the material time-derivative of the internal energy e per unit mass, and D = F − Ė F −1 denotes the strain-rate tensor (i.e., the symmetric part of the Eulerian velocity gradient Ḟ F −1 ).The second principle of thermodynamics imposes the increase of entropy ρ η ≥ 0. Assuming an adiabatic process, the dissipation per unit of reference volume reads D = ρ 0 T η (W/m 3 ), where T > 0 is the absolute temperature.Thus, combining the local equations of thermodynamics with the Gibbs identity, the Clausius-Duhem inequality is obtained: Here, U = ρ 0 e is the internal energy per unit reference volume.Since the inequality (11) must be satisfied for all states and all evolutions, the coefficient T − ∂e/∂η must equal zero.Using the identity J σ : D = S : Ė , see e.g.Ref. [29], the Clausius-Duhem inequality ( 11) is rewritten as We call this framework isentropic because the partial derivatives ∂/∂E , ∂/∂α k are evaluated at constant entropy [1].
Isothermal modelling This approach involves the variables of state {T, E , α 1 , . . ., α n }.It is linked to the above expressions by introducing the partial Legendre transform ψ = e −T η of e, which is Helmholtz' free energy per unit mass.We have where Ψ = ρ 0 ψ is the Helmholtz free energy per unit of reference volume.The Clausius-Duhem inequality ( 12) is re-expressed as We observe that thermodynamic restrictions have the same form in the isentropic and in the isothermal frameworks.In the isothermal framework described here, proving thermodynamic consistency of Fung's QLV (9) amounts to finding Ψ, α k such that the Clausius-Duhem inequality ( 14) is always satisfied.

Fung's quasi-linear viscoelasticity
Memory variables Using the expression of the relaxation function (10), the constitutive law ( 9) is rewritten as [21] where the viscous stress can thus be interpreted as a memory variable.Computing its material time-derivative, one shows that S v k satisfies the linear evolution equation [21] Thus, the convolution product ( 9) is replaced by a sum of n memory variables, which satisfy a linear differential equation.By construction, Fung's QLV model reduces to hyperelasticity for particular relaxation functions involving certain limits of relaxation times: • Relaxed elastic solid.The relaxed elastic limit corresponds to infinite durations, i.e. to short relaxation times τ k → 0. Hence, the evolution equation ( 17) produces S v k = g k S e .If the motion is causal, the convolution product (9) reduces to S = (1 − k g k ) S e where the coefficient 1 − k g k defines the relaxed elastic modulus.
• Unrelaxed elastic solid.The unrelaxed elastic limit corresponds to infinitesimal durations, i.e. to long relaxation times τ k → +∞.Hence, the evolution equation (17) gives S v k = 0 for causal signal.The convolution product (9) reduces to S = S e .With respect to the relaxed elastic solid, the effective elastic moduli differ by a scalar coefficient.
These elastic limits correspond to zero dissipation [30].
Dissipation Consider a (presumably convex) strain energy function W from which the elastic response S e is obtained by differentiation.We define the free energy in such a way that S = ∂Ψ/∂E is satisfied: The arbitrary functions Φ k are (presumably convex) potentials whose dependency on the variables S v k of Eq. ( 16) needs to be specified.If the viscous stresses S v k governed by Eq. ( 17) are internal variables of state α k , then the dissipation ( 14) reads Let us introduce the Legendre transform The free energy (18) and the dissipation (19) become To ensure the positivity of the dissipation, a compatible choice of potentials is Consequently, the scaling property of the Legendre transformation imposes where the potential Φ known as complementary energy density defines the Legendre transform S e : E − W (E ) of W , with E = ∂Φ/∂S e [31,28].Finally, we have By virtue of the convexity inequality [28], the dissipation D is non-negative for any convex strain energy density function W of E .Therefore, the present compressible QLV model is thermodynamically admissible.

Connections with other models
The derivation of the model introduced by Simo [26] is very similar.In fact, this model is based on the free energy where the separability of isochoric and volumetric deformations W = W iso ( Ĩ1 , Ĩ2 ) + W vol (J ) is assumed.Similarly to Eq. ( 17), a consistent linear evolution equation for the viscous stresses is proposed.Thus, the main difference between Eqs. ( 22) and ( 18) lies in the presence of the volume-preserving version Ẽ = 1 2 ( C − I ) of E .Note that the "Simo" model described in the MSC Nastran Implicit Nonlinear user guide is actually a QLV model (Ref.[32], Chap.10, Eq. (10-88)).Nevertheless, that implementation assumes "that the viscoelastic behavior . . .acts only on the deviatoric behavior".While it remains unclear how this is done, the implementation may therefore be consistent with the original Simo model [26].The free energy (18) of Fung's QLV can be rewritten as The symbol Ψ ∞ denotes the free energy of the relaxed elastic solid, such that S ∞ = ∂Ψ ∞ /∂E is the corresponding Piola-Kirchhoff stress.Introducing the variable yields the constitutive law S = S ∞ + k Q k and the rate equation Qk = g k Ṡe − Q k /τ k .These expressions are very similar to Sec. 6.10 of the monograph by Holzapfel [1].Indeed, the latter follow from the works of Simo [26] and Holzapfel and Simo [27], as well as Govindjee and Simo [33].Thus, the above nonlinear viscoelasticity theories are strongly related to QLV.
Moreover, they are equivalent in the incompressible limit (see A), as observed in previous work [20,34].
The expression of the free energy may also be rewritten as a hereditary integral.Miller and Chinzei [35] and related studies [21,11] introduce the convolution product Ψ * = G * Ẇ from which the stress S = G * Ṡe is derived.While the final expression of the stress can easily be related to the present study (Eq.( 9)), it is not as straightforward to explain the form of the free energy Ψ * .This may be a consequence of a daring differentiation of the free energy with respect to the instantaneous strain inside the convolution integral (see Ref. [21] Eq. ( 8)).Note that in the infinitesimal strain limit, the consistent potential energy reads as a double convolution of the relaxation function with the strain rates (see e.g.Carcione [10] Chap.2).This remark further questions the expression Ψ * of the free energy as a single convolution product.

Illustrations
The theoretical analysis of Sec. 3 is now illustrated by means of simple deformations.We consider compressible neo-Hookean QLV solids described by the strain energy [28] where the Lamé parameters µ , µ are positive.The corresponding bulk modulus reads µ + 2 3 µ.The elastic response S e = ∂W /∂E deduced from Eqs. ( 7)- (8) reads Due to consistency with linear elasticity in the infinitesimal strain limit, this constitutive law is at least locally invertible.Thus, we introduce the complementary energy density Φ(S e ) such that E = ∂Φ/∂S e .The components of F and C for uniform extension along the X -direction are of the form [28] where λ > 0 is the tensile stretch.In simple tension, the tractions transverse to the X -direction vanish, so Fig. 2a shows the evolution of s e with respect to the stretch λ, as well as the evolution of the Cauchy stress component λ 2 s e /J .Here, we have chosen ϑ = 1 3 × 10 −3 , which is a typical value for nearlyincompressible rubber-like soft solids.One observes that the stress-stretch relationship is one-to-one over the range displayed in the figure.Within this range, we can deduce the stretch λ > 0 from s e using Eqs.( 28)-( 29), e.g. by means of a root-finding algorithm.Thus, we can retrieve the deformation in Eq. ( 27) from the stress.
The viscoelastic stress S is given by the convolution product in Eq. ( 9).This constitutive law is rewritten as S = S e − k S v k in terms of the memory variables S v k , which depend on the whole history t → S e (t ) of the elastic stress.Here, the viscous stresses are of the form S v k = s v k e 1 ⊗ e 1 for all times, where 1 − e −(t −s)/τ k ṡe (s) ds (30) is obtained by componentwise integration of Eq. ( 16).The viscous strains E v k are deduced from the viscous stresses, as specified in Sec. 3. Thus, these memory variables are also functions of the entire stress history.In fact, one shows that k /g k ) by using the definition of Φ k with respect to Φ.In other words, the mappings E → S e and E v k → S v k /g k based on the strain energy W are the same, and they admit the same inverse based on the complementary energy Φ.In the present one-dimensional case, it suffices to replace s e , λ by s v k /g k , λ v k in Eqs. ( 28)-( 29) to retrieve the stretches λ v k from the viscous stresses s v k .This way, the corresponding deformations E v k are obtained.Note in passing that besides being quite involved, the computation of E v k is not needed in most practical applications.
The dissipation D of Eq. ( 21) involves the sum of terms of the form whose evaluation follows from Eqs. ( 28)-( 29) in the present one-dimensional case.The positivity of D k for all deformations and all evolutions is illustrated in Fig. 2b.The black line marks the locus E = E v k of the relaxed elastic solid limit where D k vanishes.Finally, the compressible neo-Hookean QLV solid model is thermodynamically admissible in simple tension over the range of stretches considered here.

Conclusion
In the framework of Fung's quasi-linear viscoelasticity model, the reinterpretation of the natural memory variables as internal variables of state provides an expression of the free energy.The dissipation is shown to be positive provided that the strain energy function is convex.We note that the model equations soobtained can be linked to other models in the literature.
We should be aware of the limitations of the QLV model, which is known to hardly capture the discrepancy between creep and relaxation time scales [18] and which in general does not exhibit straindependent relaxation effects.Moreover, QLV may only be valid at moderate deformations [11].However, its strong similarities to other models mean that its regimes of validity are more-or-less equivalent.The above results could be extended and employed in the modelling of anisotropic materials [34] and thermoelastic materials [1].The results may also be extended to compressible materials with distinct relaxation functions in shear and in compression [13].More general results could be obtained by exploiting the notion of fading memory, and the similarity between QLV and the linear viscoelasticity formalism [36,37].

A Incompressible case
In the case of incompressible materials, the volume dilatation J ≡ 1 is prescribed at all times, and the mass density ρ = ρ 0 is constant.The deformation E is equal to its volume-preserving counterpart Ẽ = 1 2 ( C − I ).Moreover, the third invariant I 3 ≡ 1 is prescribed too, so that the strain energy can be reduced to a function W (I 1 , I 2 ).The constitutive law (9) becomes S = −pC −1 + G * Ṡe = −pC −1 + Ġ * S e , (32) where S e = ∂W /∂E is deduced from Eq. ( 8) with W 3 = 0, and p is an indeterminate Lagrange multiplier due to the incompressibility constraint.
In a similar fashion to the compressible case (15), the constitutive law (32) is expressed as where the viscous stresses S v k are governed by the linear evolution equation (17).We define the free energy Ψ using the same expression as in Eq. ( 18).This way, the Piola-Kirchhoff stress of Eq. ( 32) satisfies S = −pC −1 + ∂Ψ/∂E .The incompressibility constraint is included in the Clausius-Duhem inequality ( 14) by introducing the Lagrange multiplier p as follows: see Holzapfel [1] Sec. 6.3.The next steps of the derivation are analogous to the compressible case, and finally, we obtain the same expression of the dissipation as in Eq. (21).Therefore, the positivity of the dissipation is guaranteed for convex strain energy functions W .
Remark.In Eq. ( 32), the elastic stress S e may include hydrostatic stress contributions -in other words, we have S e : C ≡ 0. Assuming the elastic stress S e purely deviatoric is rather restrictive.Indeed, the corresponding Lagrange-Charpit equations yield W = f (I 1 / I 2 ) where f is an arbitrary function.Sometimes, incompressible QLV is formulated as follows [13] where S e D = Dev(S e ) is purely deviatoric.To proceed as above, one may replace the strain energy W by the strain energy W = W (I 1 , I 2 ) − 1 6 (S e : C )(I 3 − 1) (36) in Eq. ( 18), which is defined in such a way that S e D = ∂W /∂E under the incompressibility constraint.Alternatively, one may extend W ( Ĩ1 , Ĩ2 ) to non-unimodular deformation gradients, such that S e D = ∂ W /∂E is satisfied when incompressibility is assumed [38].Thus, the analysis of the dissipation's sign for Eq. ( 35) seems more intricate.

Figure 1 :
Figure 1: Number of QLV records by year in the Web of Science (WoS) database for the very specific query: AB=(visco* AND quasi* AND (Fung OR QLV*)).