A continuum and computational framework for viscoelastodynamics: finite deformation linear models
AA continuum and computational framework for viscoelastodynamics: finitedeformation linear models
Ju Liu a,b,c , Marcos Latorre d , Alison L. Marsden c,e a Department of Mechanics and Aerospace Engineering, Southern University of Science and Technology, Shenzhen,Guangdong 518055, P.R.China b Guangdong-Hong Kong-Macao Joint Laboratory for Data-Driven Fluid Mechanics and Engineering Applications, SouthernUniversity of Science and Technology, Shenzhen, Guangdong, 518055, P.R.China c Department of Pediatrics (Cardiology) and Institute for Computational and Mathematical Engineering, Stanford University,Clark Center E1.3, 318 Campus Drive, Stanford, CA 94305, USA d Department of Biomedical Engineering, Yale University, New Haven, CT 06520, USA e Department of Bioengineering, Stanford University, Clark Center E1.3, 318 Campus Drive, Stanford, CA 94305, USA
Abstract
This work concerns the continuum basis and numerical formulation for deformable materials with viscousdissipative mechanisms. We derive a viscohyperelastic modeling framework based on fundamental thermo-mechanical principles. Since most large deformation problems exhibit isochoric properties, our modelingwork is constructed based on the Gibbs free energy in order to develop a continuum theory using pressure-primitive variables, which is known to be well-behaved in the incompressible limit. A set of general evolutionequations for the internal state variables is derived. With that, we focus on a family of free energies thatleads to the so-called finite deformation linear model. Our derivation elucidates the origin of the evolutionequations of that model, which was originally proposed heuristically and thus lacked formal compatibilitywith the underlying thermodynamics. In our derivation, the thermodynamic inconsistency is clarified andrectified. A classical model based on the identical polymer chain assumption is revisited and is found tohave non-vanishing viscous stresses in the equilibrium limit, which is counter-intuitive in the physical sense.Because of that, we then discuss the relaxation property of the non-equilibrium stress in the thermodynamicequilibrium limit and its implication on the form of free energy. A modified version of the identical polymerchain model is then proposed, with a special case being the model proposed by G. Holzapfel and J. Simo.Based on the consistent modeling framework, a provably energy stable numerical scheme is constructed forincompressible viscohyperelasticity using inf-sup stable elements. In particular, we adopt a suite of smoothgeneralization of the Taylor-Hood element based on Non-Uniform Rational B-Splines (NURBS) for spatialdiscretization. The temporal discretization is performed via the generalized- α scheme. We present a suite ofnumerical results to corroborate the proposed numerical properties, including the nonlinear stability, robust-ness under large deformation, and the stress accuracy resolved by the higher-order elements. Additionally,the pathological behavior of the original identical polymer chain model is numerically identified with anunbounded energy decaying. This again underlines the importance of demanding vanishing non-equilibriumstress in the equilibrium limit. Keywords:
Continuum mechanics, Gibbs free energy, Viscoelasticity, Incompressible solids, Isogeometricanalysis, Nonlinear stability
1. Introduction
Many polymeric and biological materials may undergo large deformations, with mechanical behaviorcharacterized by hyperelasticity with intrinsic viscous dissipative mechanisms [1, 2, 3, 4, 5]. To incorporate
Email addresses: [email protected],[email protected] (Ju Liu), [email protected] (Marcos Latorre), [email protected] (Alison L. Marsden)
Preprint submitted to Computer Methods in Applied Mechanics and Engineering February 12, 2021 a r X i v : . [ m a t h . NA ] F e b he viscous dissipative mechanism into the solid model, there have been several different modeling approachesdeveloped to date. From a big picture perspective, prior modeling work can be categorized into at leasttwo different groups: one based on the concept of internal state variables [6, 7, 8] and one based on thehereditary integral [9, 10, 11, 12]. The former approach invokes a set of internal state variables, which isassociated with the physical process occurring at the microscopic level and is manifested at the macroscopicscale. The evolution of the internal state variables is constrained by the second law of thermodynamics andgenerally results in dissipative behavior. This concept has been successfully applied to modeling inelasticityand phase transitions in general [7, 8]. In the second approach, the model works directly on the stress-strainrelationship. The viscous effect is modeled via convolution of a relaxation function with the strain historyto describe the fading memory effect on the stress [13], and the convolution is also known as the hereditaryintegral. This approach includes the notable quasi-linear viscoelasticity theory that has been implementedby the finite element method [14] and applied to investigating soft tissues [15, Chapter 7]. It has also beengeneralized by utilizing fractional-order derivatives to capture the continuous relaxation spectrum, givingrise to a propitious alternative candidate for viscoelasticity modeling [16, 17, 18]. The constrained mixturetheory, which has been developed and applied for vascular growth and remodeling, can be categorized intothis general modeling framework, in which the fading memory effect due to the continual turnover is fittedinto the hereditary integral [19, 20]. We also mention that some formulations based on internal state variablescan be equivalently written into the hereditary integral formulation. The model discussed in this article,as will be shown, can be represented in terms of the hereditary integral with exponential relaxation kernels[21, 22].Focusing on the internal state variable approach, different underlying modeling assumptions have givenrise to different models and computational procedures for viscoelasticity at finite strains. In recent works,it has been common to employ a multiplicative decomposition of the deformation gradient to characterizeelastic and viscous deformations [23]. With that, viscoelastic materials at finite strains may be established[24, 25, 26, 27, 28, 29, 30], and when combined with Ogden’s hyperelastic function [31], are able to describeresponses that may adequately deviate from the thermodynamic equilibrium. Among those works, the modelconstructed by Reese and Govindjee [25, 26] is rather representative. In their works, the non-equilibriumstresses are consistently derived from a Helmholtz free energy, that can be additively split into equilibriumand non-equilibrium parts [32]. The derivation is similar to that of the finite strain elastoplasticity theory[33], and material isotropy was assumed in their derivation. Recently, with the isotropy assumption released,the modeling approach of [25] was further extended to account for material anisotropy [34, 35].Alternatively, a model proposed by Simo in [36] gained popularity over the years and inspires numeroussubsequent modeling works [37, 38, 39, 40, 41]. In his approach, the non-equilibrium stresses due to theviscous effect, which are also termed the “over-stress”, are governed by a set of linear evolution equations,with the elastic stress derived from strain energy. The evolution equations were proposed in a heuristic manner as a straightforward generalization of the standard Zener solid model. Due to the linear nature, themodeling approach is often termed as finite deformation linear viscoelasticity , or finite linear viscoelasticity for short. The linear evolution equations also enable one to express the non-equilibrium stresses in terms ofa simple convolution integral, which can then be computed via a one-step second-order accurate recurrenceformula. The algorithm incrementally integrates the constitutive laws by only using the information fromthe past one step [36, 42, 43]. Thanks to this recurrence formula, the viscoelasticity model becomes amenableto finite element implementations [42, 44, 45]. Another appealing feature is that the strain energy functionmay conveniently account for material anisotropy [39, 46, 47, 48], as the over-stresses are governed by aset of linear evolution equations. On the other side, although large deformation is allowed for such models,the inherent linearity assumed for the evolution equation is regarded as a drawback of this approach, asit is regarded to be suited only for small deviations from thermodynamic equilibrium [44, Section 6.10](see also [49, Chapter 10]). Yet, at least under physiological settings, we may reasonably expect the smalldeviation assumption to remain sound, and thus the model may still be well-suited to many biomechanicalapplications. A more critical issue is its lack of a rational thermodynamic foundation. Indeed, the linearevolution equations have not been explained by a rational thermomenchanical theory, although some relateddiscussions were made in [46]. In recent work, a finite-time blow-up solution has been identified for thistype of model [50], signifying its inconsistent nature. This alerting evidence hinders further adoption of this2odel for viscoelastic materials.In this work, we consider a general continuum formulation for fully coupled thermomechanical modelswith viscous-type dissipation. The viscous deformation is characterized by a set of internal state variables.Unlike the prior approach [36], the non-equilibrium stresses are distinguished from the variables conjugatingto the internal state variables. A set of fully nonlinear evolution equations for those conjugate variables arealso obtained following the standard Coleman-Noll argument [51]. The origin and relations of the conjugatevariables with the non-equilibrium stresses are elucidated through a careful derivation. In particular, itis shown that the non-equilibrium stresses and the conjugate variables can be identical only under a veryspecial circumstance, an issue that has long been ignored in the literature [36, 42, 44]. Next, we considera special form of the non-equilibrium part of the free energy, or the configurational free energy, whichgoverns the viscous responses. The particular form of the energy can be viewed as a generalization of theform proposed in [46]. If one further demands the energy to be quadratic in terms of the internal statevariables, a set of linear evolution equations for the conjugate variables is obtained, without invoking anylinearization technique [25]. To the best of our knowledge, this is the first time that the evolution equations’thermodynamic origin gets elucidated. We also mention that the right-hand side of the evolution equationsis governed by the fictitious second Piola-Kirchhoff stress, which is slightly different from the classical model[42, Chapter 10]. As the first instantiation of this framework, we considered the identical polymer chainmodel proposed in [46]. Interestingly, it can be easily seen that the non-equilibrium stresses of this model donot relax to zero in the thermodynamic equilibrium limit, which is physically counter-intuitive. It is indeed arequirement that the viscous stresses should vanish for static processes in general [49, Chapter 10]. Withouta prior multiplicative decomposition of the deformation, the non-equilibrium stresses are not guaranteed tofully relax in the limit. Therefore, to characterize the relaxation, we discuss the necessary and sufficientconditions from the perspective of the free energy design. Following that, we analyze first a simple modelthat has been considered in [46, Section 4.2]. This model can be viewed as an extension of the St. Venant-Kirchhoff model to the viscous part and is shown to be the only model in which the non-equilibrium stressesequal the corresponding conjugate variables. Furthermore, the configurational free energy of this modelsuggests that there is an additive split of the quadratic strain, making it parallel to the elastoplasticitytheory developed by Green and Naghdi [52, 53]. The last example considered is a modified version ofthe identical polymer chain model, inspired from the aforesaid one [40, 46] and is guaranteed to have thenon-equilibrium stresses fully relaxed in the thermodynamic equilibrium limit.For finite deformation problems, the deformation is typically highly isochoric. If one considers theHelmholtz free energy as the thermodynamic potential, the resulting system becomes singular in the incom-pressible limit (see Figure 1). The pressure-volume curve, as part of the stress-strain curve, has a very largeslope and becomes a vertical line in the incompressible limit. The value of the slope enters into the elasticitytensor and engenders an ill-conditioned stiffness matrix in the incompressible limit, a well-known issue inthe displacement-based elasticity formulation [54]. To handle this issue, one may consider performing a Leg-endre transformation for the free energy and switch the independent variable from the specific volume to itsconjugate counterpart, the pressure [55]. In doing so, the problem enjoys a saddle-point nature with the pres-sure acting as an independent variable. Correspondingly, in the incompressible limit, the volume-pressurecurve becomes horizontal with a slope approaching zero. With the Legendre transformation performed,the resulting thermodynamic potential becomes the Gibbs free energy, and it is, therefore, a well-suitedpotential for thermomechanical analysis. This argument is supported by an analysis made in the realm ofcomputational fluid dynamics (CFD). It has been shown that the pressure primitive variable is among thetwo sets of variables for interpolation if one wants to have a compressible CFD code that survives in the lowMach number limit [56, 57]. The aforementioned Gibbs free energy can be viewed as the thermodynamicexplanation of the effectiveness of the pressure primitive variables. A theory for finite elasticity based on theGibbs free energy was established [55]; the saddle-point nature has also been exploited to design precondi-tioners [58, 59]; it can also be shown to enjoy a provably nonlinearly stable semi-discrete formulation usinginf-sup stable elements [60]; it can be conveniently utilized to construct a strongly-coupled fluid-structureinteraction formulation [55, 61]. In this work, we follow our prior approach and use the Gibbs free energyto construct a mixed formulation for viscoelastodynamics. The thermodynamic consistent nature of thecontinuum model naturally allows us to devise a discretization that inherits the energy stability property.3 ! 𝜂, 𝐽, 𝑁= Θ𝜂 − 𝑃𝐽 + 𝜇𝑁𝐻 𝜂, 𝑃, 𝑁= Θ𝜂 + 𝜇𝑁 𝐺 Θ, 𝑃, 𝑁 = 𝜇𝑁𝐴 Θ, 𝐽, 𝑁= −𝑃𝐽 + 𝜇𝑁 Figure 1: Illustration of the Legendre transformation of thermodynamic potentials: the internal energy ι R , Helmholtz freeenergy A , Gibbs free energy G , and enthalpy H . In particular, the Helmholtz free energy A (Θ , J, N ) can be transformed to G (Θ , P, N ) assuming convexity of A with respect to J . The resulting constitutive laws based on A and G are given in theright two figures. In the theory based on the Helmholtz free energy, one has P = − ∂A/∂J . In the incompressible limit, thepressure-volume (i.e., stress-strain) curve becomes a vertical line, and the corresponding components in the elasticity tensorblow up to infinity, which in turn causes numerical issues. This partly explains that, for most displacement-based formulations,the tangent matrix is ill-conditioned and demands a direct solver oftentimes. In the constitutive laws based on G , the volume isgiven by J = ∂G/∂P , which becomes a horizontal line in the thermodynamic limit. The resulting continuum model constitutesa saddle-point problem and requires inf-sup stable or stabilized discretization techniques. The resulting mixed formulation necessitates inf-sup stable element pairs, and we adopt a suite of smoothgeneralization of the Taylor-Hood element by NURBS that has been numerically demonstrated to be stable[60, 62, 63, 64]. In particular, the same set of basis functions is utilized to represent the geometry and toconstruct the spaces for the displacement and velocity fields, which fits well into the paradigm of NURBS-based isogeometric analysis [65, 66]. One appealing feature of this choice is that it achieves higher-orderaccuracy [67] without sacrificing robustness [68], in contrast to conventional C -continuous finite elementsthat often becomes “fragile” when raising the polynomial order [69]. Lastly, we mention other promisingcandidates for the spatial discretization, such as the smooth generalization of the Raviart-Thomas elements[62, 70] as well as divergence-free discontinuous Galerkin methods [71].The body of this work starts in Section 2 where we derive a general continuum theory for viscoelasticitywith the Gibbs free energy as the thermodynamic potential. Following, we restrict the discussion to aspecial form of the free energy and provide a definition for the finite linear viscoelasticity. After revisitingthe identical polymer chain model, the relaxation condition for the viscous stress is discussed. Two materialmodels are then presented, which completes the continuum modeling section. In Section 3, we consider thespatial and temporal discretization of the constructed continuum model. In particular, we show the energystability of the proposed numerical formulation. In Section 4, a suite of three numerical tests is performedas an examination of the continuum model and verification of the numerical scheme. We give concludingremarks in Section 5. In Appendices, the implementation details of the three material models are presented.
2. Continuum Basis
Let Ω X and Ω t x be bounded open sets in R with Lipschitz boundaries. The motion of the body isdescribed by a family of smooth mappings parameterized by the time coordinate t , ϕ t ( · ) = ϕ ( · , t ) : Ω X → Ω t x = ϕ (Ω X , t ) = ϕ t (Ω X ) , ∀ t ≥ , X (cid:55)→ x = ϕ ( X , t ) = ϕ t ( X ) , ∀ X ∈ Ω X .
4n the above, x represents the current position of a material particle originally located at X , which implies ϕ ( X ,
0) = X . The displacement and velocity of the material particle are defined as U := ϕ ( X , t ) − ϕ ( X ,
0) = ϕ ( X , t ) − X , V := ∂ ϕ ∂t (cid:12)(cid:12)(cid:12)(cid:12) X = ∂ U ∂t (cid:12)(cid:12)(cid:12)(cid:12) X = d U dt . In this work, we use d ( · ) /dt to denote a total time derivative. The spatial velocity is defined as v := V ◦ ϕ − t .Analogously, we define u := U ◦ ϕ − t . The deformation gradient, the Jacobian determinant, and the rightCauchy-Green tensor are defined as F := ∂ ϕ ∂ X , J := det ( F ) , C := F T F . Since most materials of interest behaves differently in bulk and shear under large strains, the deformationis multiplicatively decomposed into a volumetric part J / I and an isochoric part ˜ F := J − / F [72, 73].Clearly, by construction one has the multiplicative decomposition of the deformation gradient as F = (cid:16) J I (cid:17) ˜ F . (2.1)The corresponding modified right Cauchy-Green tensor ˜ C is defined as˜ C := J − C . To facilitate the following discussion, we note the differentiation relation ∂ ˜ C ∂ C = J − P T with P := I − C − ⊗ C , wherein I is the fourth-order identity tensor. The projection tensor P furnishes deviatoric behaviors in theLagrangian description [44]. Furthermore, in this work, the magnitude of a second-order tensor A is denotedas | A | := tr[ AA T ] . (2.2) The motion of the continuum body has to satisfy the local balance equations as well as the second lawof thermodynamics [74, 75]. The advective form of the mass balance equation can be written as dρdt + ρ ∇ x · v = 0 , (2.3)wherein ρ is the mass density in the current configuration. This mass balance equation is equivalent to theequation for the volumetric strain J , dJdt = J ∇ x · v , or equivalently, dJdt = J ∇ X V : F − T . (2.4)The balance of linear momentum can be written as ρ d v dt = ∇ x · σ + ρ b , or equivalently, ρ d V dt = ∇ X · P + ρ B , (2.5)wherein σ denote the Cauchy stress, b represents the body force per unit mass, P := J σF − T is the firstPiola-Kirchhoff stress, and B := b ◦ ϕ t . The balance of angular momentum is satisfied by imposing symmetryon the Cauchy stress, i.e., σ = σ T . The balance of internal energy is stated as ρ dιdt = σ : ∇ x v − ∇ x · q + ρr, (2.6)5n which ι is the internal energy per unit mass, q denote the heat flux, and r is the heat source per unitmass. Further, we introduce the following quantities: the internal energy ι R defined with respect to thereference volume is related to ι by ι R := ρ ι ; the heat flux defined per unit referential surface area is denotedas Q , and it is related to q by the Piola transformation Q = J F − q ; the heat supply per unit referentialvolume R := ρ r . With these definitions, the balance of the internal energy can also be expressed as dι R dt = P : d F dt − ∇ X · Q + R. (2.7)Let us introduce S as the entropy per unit volume in the referential configuration, Θ as the absolutetemperature field in the reference configuration, and D as the dissipation. The second law of thermodynamicsstates the dissipation of the system is non-negative, i.e., D := d S dt + ∇ X · (cid:18) Q Θ (cid:19) − R Θ ≥ . (2.8)The Gibbs free energy per unit volume in the referential configuration is defined as [76, Chapter 5] G := ι R − Θ S + P J, (2.9)wherein P is the thermodynamic pressure defined on the referential configuration, which is the conjugatevariable to J . The thermodynamic pressure defined on the current configuration is denoted by p := P ◦ ϕ − t .Taking material time derivatives at both sides of (2.9) results in dGdt + S d Θ dt − J dPdt = dι R dt − Θ d S dt + P dJdt .
Substituting the balance equation of the internal energy (2.7), the second law of thermodynamics (2.8), andthe mass balance equation (2.4) into the above equation, one readily obtainsΘ D = J σ : ∇ x v − Q · ∇ X ΘΘ +
P J ∇ x · v − S d Θ dt + J dPdt − dGdt . (2.10)We may additively split the Cauchy stress into deviatoric and hydrostatic parts, σ = σ dev + 13 (tr [ σ ]) I , with which one may show that J σ : ∇ x v = 12 J ˜ F − σ dev ˜ F − T : ddt ˜ C + 13 J tr [ σ ] ∇ x · v . Consequently, the relation (2.10) can be rewritten asΘ D = 12 J ˜ F − σ dev ˜ F − T : ddt ˜ C + J (cid:18)
13 tr [ σ ] + p (cid:19) ∇ x · v − Q · ∇ X ΘΘ − S d Θ dt + J dPdt − dGdt . (2.11)From the above, we postulate that the Gibbs free energy G is a function of ˜ C , P , and Θ by invokingTruesdell’s principle of equipresence [77]. Additioanlly we assume that G also depends on a set of strain-likeinternal state variables { Γ α } mα =1 , G = G ( ˜ C , P, Θ , Γ , · · · , Γ m ) . (2.12)Here m is the number of relaxation processes characterizing the viscous property of the material. Oftentimes,the internal state variables Γ α is viewed as a strain tensor akin to the right Cauchy-Green strain tensor[37, 46]. We demand that for a given homogeneous reference temperature Θ > G ( I , , Θ , I , · · · , I ) = 0 , (2.13)6hich is commonly known as the normalization condition [44]. With the above function form of the Gibbsfree energy, material time derivative of the Gibbs free energy can be written explicitly as dGdt = 12 ˜ S : ddt ˜ C + ∂G∂P dPdt + ∂G∂ Θ d Θ dt + m (cid:88) α =1 ∂G∂ Γ α : ddt Γ α , (2.14)in which we introduced the fictitious second Piola-Kirchhoff stress˜ S := 2 ∂G ( ˜ C , P, Θ , Γ , · · · , Γ m ) ∂ ˜ C . Substituting (2.14) into (2.10) leads toΘ D = (cid:18) J ˜ F − σ dev ˜ F − T −
12 ˜ S (cid:19) : ddt ˜ C + J (cid:18)
13 tr [ σ ] + p (cid:19) ∇ x · v − Q · ∇ X ΘΘ − (cid:18) S + ∂G∂ Θ (cid:19) d Θ dt + (cid:18) J − ∂G∂P (cid:19) dPdt − m (cid:88) α =1 ∂G∂ Γ α : 12 ddt Γ α . (2.15)We introduce Q α as the conjugate variables to the internal state variables Γ α , Q α := − ∂G ( ˜ C , P, Θ , Γ , · · · , Γ m ) ∂ Γ α . (2.16)The relation (2.15) furnishes a natural means of choosing constitutive relations that ensure the satisfactionof the second law of thermodynamics. We may thereby make the following choices for the constitutiverelations, σ dev = J − ˜ F (cid:16) P : ˜ S (cid:17) ˜ F T = J − ˜ F (cid:18) P : (cid:18) ∂G∂ ˜ C (cid:19)(cid:19) ˜ F T , (2.17)13 tr [ σ ] = − p, (2.18) Q = − ¯ κ ∇ X Θ , (2.19) S = − ∂G∂ Θ , (2.20) ρ = ρ (cid:18) ∂G∂P (cid:19) − , (2.21) Q α = V α : (cid:18) ddt Γ α (cid:19) . (2.22)In the above, ¯ κ is the thermal conductivity and V α is a positive definite fourth-order viscosity tensor. Here,we also assume that there exists a fourth-order tensor ( V α ) − such that an inverse relation for (2.22) holds,12 ddt Γ α = ( V α ) − : Q α . (2.23)Based on (2.17) and (2.18), the Cauchy stress σ can be represented as σ := σ dev + 13 tr [ σ ] I = J − ˜ F (cid:16) P : ˜ S (cid:17) ˜ F T − p I = J − ˜ F (cid:18) P : (cid:18) ∂G∂ ˜ C (cid:19)(cid:19) ˜ F T − p I . S , as the pull-back operation performed on J σ , can be writtenas S := J F − σF − T = J F − (cid:18) σ dev + 13 tr [ σ ] I (cid:19) F − T = S iso + S vol , S iso := J F − σ dev F − T = J − P : ˜ S = J − P : 2 ∂G∂ ˜ C , S vol := 13 tr [ σ ] J F − F − T = − JP C − . With the above choices, it can be shown that the dissipation relation (2.15) reduces toΘ D = ¯ κ |∇ X Θ | − m (cid:88) α =1 Q α : 12 (cid:18) ddt Γ α (cid:19) = ¯ κ |∇ X Θ | + m (cid:88) α =1 (cid:18) ddt Γ α (cid:19) : V α : (cid:18) ddt Γ α (cid:19) = ¯ κ |∇ X Θ | + m (cid:88) α =1 Q α : ( V α ) − : Q α , (2.24)which remains non-negative for arbitrary kinematic processes, thereby automatically satisfying the secondlaw of thermodynamics. The choices of the constitutive relations (2.17)-(2.22) are certainly not a uniqueway for ensuring the second law of thermodynamics. Yet, as will be shown, these choices are general enoughto describe the mechanical behavior of a variety of materials. The definition of Q α (2.16) suggests that it isa function of Γ α , ˜ C , P , and Θ. Thereby, taking material time derivatives of Q α yields ddt Q α = 2 ∂ Q α ∂ Γ α : 12 ddt Γ α + 2 ∂ Q α ∂ ˜ C : 12 ddt ˜ C + ∂ Q α ∂ Θ ddt Θ + ∂ Q α ∂P ddt P. Invoking the constitutive relation (2.23) to replace the rate of Γ α in the above relation, we obtain a set ofevolution equations for Q α after rearranging terms, ddt Q α − ∂ Q α ∂ Γ α : ( V α ) − : Q α = 2 ∂ Q α ∂ ˜ C : 12 ddt ˜ C + ∂ Q α ∂ Θ ddt Θ + ∂ Q α ∂P ddt P. (2.25)The above constitutes a set of general evolution equations for Q α , and in the latter part of this work we willreveal how it may recover the familiar linear evolution equations [44, Chapter 6]. In the last, we characterizethe thermodynamic equilibrium state by the following two conditions, ∇ X Θ (cid:12)(cid:12)(cid:12) eq = , and ddt Γ α (cid:12)(cid:12)(cid:12) eq = O , for α = 1 , · · · , m. (2.26)In the above, the notation (cid:12)(cid:12) eq indicates that equality holds at the thermodynamic equilibrium limit. Theconditions (2.26) suggest that, in the equilibrium limit, there is no more dissipation in the system. Inparticular, (2.26) implies that there is no more change of the internal state variables in the thermodynamicequilibrium limit. According to the constitutive relation (2.22), (2.26) is equivalent to Q α (cid:12)(cid:12)(cid:12) eq = O , for α = 1 , · · · , m. (2.27) Here we consider a more refined structure of the thermodynamic potential G . First, we note that thedensity ρ has to be independent of ˜ C , since the modified right Cauchy-Green tensor describes volume-preserving deformations. Therefore, one may reasonably demand that the density is a function of thepressure and temperature. Then the constitutive relation (2.21) implies that ∂G ( ˜ C , P, Θ , Γ , · · · , Γ m ) ∂P = ρ ρ − ( p, θ, Γ , · · · , Γ m ) , (2.28)8ith the equality holds by recalling that P ( X ) = p ( ϕ t ( X ) , t ) and Θ( X ) = θ ( ϕ t ( X ) , t ). The relation (2.28)suggests that the free energy can be additively split as G = G ( ˜ C , P, Θ , Γ , · · · , Γ m ) = G ∞ vol ( P, Θ , Γ , · · · , Γ m ) + G iso ( ˜ C , Θ , Γ , · · · , Γ m ) , (2.29)by performing a partial integration of (2.28) with respect to P (see also [55, p. 559]). This split structureof the free energy has also been justified based on physical observations [78, 79]. It needs to be pointed outthat the split is a logical consequence of the multiplicative decomposition (2.1). Second, there is a widely-adopted assumption with experimental justifications that viscous effects only affect the isochoric motion,which suggests that the volumetric energy G ∞ vol is independent of the internal state variables. Third, asa commonly adopted approach [44], the isochoric part of the free energy G iso can be further split into anequilibrium part as well as a dissipative, or non-equilibrium, part. The last two assumptions further lead tothe following form of the Gibbs free energy, G = G ( ˜ C , P, Θ , Γ , · · · , Γ m ) = G ∞ vol ( P, Θ) + G ∞ iso ( ˜ C , Θ) + m (cid:88) α =1 Υ α (cid:16) ˜ C , Θ , Γ , · · · , Γ m (cid:17) , (2.30)The two terms G ∞ vol and G ∞ iso characterize the hyperelastic material behavior at the equilibrium state astime approaches infinity. Given the explicit structure of the Gibbs free energy (2.30), we may introduce theequilibrium and non-equilibrium parts of the fictitious second Piola-Kirchhoff stress as˜ S = ˜ S ∞ iso + m (cid:88) α =1 ˜ S α neq , ˜ S ∞ iso := 2 ∂G ∞ iso ( ˜ C , Θ) ∂ ˜ C , ˜ S α neq := 2 ∂ Υ α ( ˜ C , Θ , Γ , · · · , Γ m ) ∂ ˜ C . (2.31)Correspondingly, the isochoric part of the second Piola-Kirchhoff stress can be expressed as S iso = S ∞ iso + m (cid:88) α =1 S α neq , S ∞ iso := J − P : ˜ S ∞ iso , S α neq := J − P : ˜ S α neq . Due to (2.30), Q α can be defined in terms of the configurational free energy Υ α as Q α := − ∂G ( ˜ C , P, Θ , Γ , · · · , Γ m ) ∂ Γ α = − ∂ Υ α ( ˜ C , Θ , Γ , · · · , Γ m ) ∂ Γ α . (2.32)Correspondingly, the evolution equations (2.25) can be reduced to ddt Q α − ∂ Q α ∂ Γ α : ( V α ) − : Q α = 2 ∂ Q α ∂ ˜ C : 12 ddt ˜ C + ∂ Q α ∂ Θ ddt Θ , (2.33)since Q α is independent of P from (2.32). Here the evolution equations (2.33) together with (2.32) can beviewed as a set of nonlinear rate equations characterizing the evolution of the internal state variables Γ α aswell. Finally, to complete the thermomechanical theory, the entropy per unit reference volume adopts theform, S = S ∞ + m (cid:88) α =1 S α , S ∞ := − ∂G ∞ vol ∂ Θ − ∂G ∞ iso ∂ Θ , S α := − ∂ Υ α ∂ Θ . Remark 1.
In our derivation, we start from the classical definition of the Gibbs free energy (2.9) withoutassuming which variables it should depend on. Subsequently, the derived relation (2.11) together with theTruesdell’s principle of equipresence suggests that it should be a function of ˜ C , P , and Θ . With the additivesplit of the free energy into isochoric and volumetric parts, we may observe from (2.31) that the isochoricpart of the Gibbs free energy plays the same role as the isochoric part of the conventional strain energy ofthe Helmholtz type [44, Chapter 6]. Therefore, we may view that the Gibbs free energy is related to theHelmholtz free energy by performing the Legendre transformation on the volumetric part of the free energy nly [55, Section 2.4]. Indeed, the volumetric energy is commonly convex with respect to J , unless there aremultiple stable phases (e.g. the van der Waals model [80, 81]), and this convexity property guarantees thevalidity of the Legendre transformation. We note that some adopt a different notion of the Gibbs free energywhich depends on the whole stress [82, 83]. That type of theory requires the convexity of the free energy withrespect to the whole strain, which does not hold for many nonlinear materials [84]. Remark 2.
There is a requirement for the non-equilibrium stresses based on physical intuitions, that is S α neq = J − P : ˜ S α neq should fully relax, or vanish, in the thermodynamic equilibrium limit. This conditionapparently poses an additional constraint for the configurational free energy Υ α , and we will revisit this pointin Section 2.4.1. Remark 3.
It has recently been observed that the additive split of the energy (2.29) may lead to non-physical responses for anisotropic materials [79, 85, 86]. In this work, we restrict our discussion to isotropicmaterials, and the split (2.29) should be sound in both mathematics and physics [78, 79].
Remark 4.
If we assume the viscosity tensors V α are isotropic, they adopt explicit forms as V α = 2 η αD (Θ) (cid:18) I − I ⊗ I (cid:19) + 23 η αV (Θ) I ⊗ I , (2.34) where η αD (Θ) and η αV (Θ) are non-negative and represent the deviatoric and volumetric viscosities respectively[25, 87, 88]. The conjugate variables Q α are related to the internal state variables Γ α by a linear relation, Q α = η αD (Θ) (cid:18) ddt Γ α −
13 tr[ ddt Γ α ] I (cid:19) + 13 η αV (Θ)tr[ ddt Γ α ] I . (2.35) Consequently, the dissipation relation (2.24) reduces to Θ D = ¯ κ |∇ X Θ | + m (cid:88) α =1 (cid:32) η αD (Θ)2 (cid:12)(cid:12)(cid:12)(cid:12) ddt Γ α −
13 tr[ ddt Γ α ] I (cid:12)(cid:12)(cid:12)(cid:12) + η αV (Θ)6 tr[ ddt Γ α ] (cid:33) . In this section, we consider specifically the finite deformation linear viscoelasticity, the precise definitionof which will be given momentarily. We start by restricting the configurational free energy to the followingform, Υ α (cid:16) ˜ C , Θ , Γ , · · · , Γ m (cid:17) = H α ( Γ α , Θ) + (cid:32) ˆ S α − ∂G α ( ˜ C , Θ) ∂ ˜ C (cid:33) : Γ α − I F α ( ˜ C , Θ) . (2.36)In the above choice, the α -th configurational free energy Υ α depends on Γ α only; H α , G α , and F α are scalar-valued functions to be provided for material modeling; ˆ S α is a constant stress-like tensor whose physicalsignificance will be revealed in Section 2.4.1. A similar form of the configurational free energy was proposedin [46] with minor differences. To satisfy the normalization condition for Υ α , we demand that H α ( I , Θ ) = 0 , F α ( I , Θ ) = 0 , for a given homogeneous temperature Θ >
0. It is worth pointing out that the essence of the particularconfigurational free energy is that the elastic strain ˜ C and the internal state variable Γ α are coupled throughthe second term in (2.36) only, where the internal state variables Γ α are present in linear forms. Also,recalling that Γ α are akin to the right Cauchy-Green strain tensor, the terms ( Γ α − I ) / S α iso ( ˜ C , Θ) := 2 ∂G α ( ˜ C , Θ) ∂ ˜ C , α (cid:16) ˜ C , Θ , Γ , · · · , Γ m (cid:17) = H α ( Γ α , Θ) + (cid:16) ˆ S α − ˜ S α iso ( ˜ C , Θ) (cid:17) : Γ α − I F α ( ˜ C , Θ) . Due to the definition (2.32) and the form of the configurational free energy (2.36), we now have Q α ( ˜ C , Θ , Γ , · · · , Γ m ) = − ∂ Υ α ( ˜ C , Θ , Γ , · · · , Γ m ) ∂ Γ α = 2 (cid:32) ∂G α ( ˜ C , Θ) ∂ ˜ C − ∂H α ( Γ α , Θ) ∂ Γ α (cid:33) − ˆ S α = ˜ S α iso ( ˜ C , Θ) − ˆ S α − ∂H α ( Γ α , Θ) ∂ Γ α . (2.37)From the above relation, the material time derivative of ˜ S α iso can be expressed as ddt ˜ S α iso = ∂ Q α ∂ ˜ C : ddt ˜ C + ∂ Q α ∂ Θ ddt Θ + 2 ∂ H α ∂ Γ α ∂ Θ ddt Θ , with which the right-hand side of the evolution equations (2.33) can be written compactly as ddt Q α − ∂ Q α ∂ Γ α : ( V α ) − : Q α = ddt ˜ S α iso − ∂ H α ∂ Γ α ∂ Θ ddt Θ . (2.38)The non-equilibrium stresses (2.31) now can be represented as˜ S α neq = 2 ∂ Υ α ∂ ˜ C = 2 (cid:32) ∂F α ( ˜ C , Θ) ∂ ˜ C − ∂ ˜ S α iso ∂ ˜ C : Γ α − I (cid:33) , (2.39)in which we have made use of the major symmetry of the fourth-order tensor ∂ ˜ S α iso /∂ ˜ C . Based on the above discussion, we now provide the definition of finite linear viscoelasticity here.
Definition 1.
A finite linear viscoelastic material is described by the configurational free energy (2.36) with V α = 2 η α (Θ) I and a quadratic form for H α in terms of the internal strain variables, H α ( Γ α , Θ) = µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12) Γ α − I (cid:12)(cid:12)(cid:12)(cid:12) , (2.40) where µ α (Θ) is a temperature dependent shear modulus associated with the α -th relaxation process. The form V α = 2 η α (Θ) I is a special case for an isotropic tensor of order four by taking η α = η αD = η αV in (2.34). According to (2.37), one can show that the quadratic form (2.40) implies Q α = ˜ S α iso − ˆ S α − µ α (Θ) Γ α − I , and ∂ Q α ∂ Γ α = − µ α (Θ) I . (2.41)Apparently, the choice of the quadratic form for H α leads to the above linear relation between Q α and Γ α ,which further simplifies the evolution equations (2.38) as ddt Q α + Q α τ α = ddt ˜ S α iso − dµ α dt ( Γ α − I ) , (2.42)with τ α (Θ) := η α (Θ) /µ α (Θ) . τ α has the dimension of time and is commonly referred to as the relaxation time for the α -thprocess. For most polymers, the relaxation time τ α can be modeled by the Arrhenius equation, whichcharacterizes faster viscoelastic effects with the increase of temperature [89]. The second term on the right-hand side of (2.42) arises due to the temperature dependent material parameters [44, p. 365]. Given aninitial condition Q α := Q α | t =0 , the solution of (2.42) can be obtained in a hereditary integral form, Q α = exp ( − t/τ α ) Q α + (cid:90) t + exp ( − ( t − s ) /τ α ) (cid:18) dds ˜ S α iso − ( Γ α − I ) dµ α ds (cid:19) ds. (2.43)We recall that the general evolution equations (2.38) is nonlinear for Q α . For finite linear viscoelasticmaterials, the evolution equations reduce to a linear system, whose solution can be conveniently representedin terms of the hereditary integral form (2.43). Within this work, we assume that Γ α = I for α = 1 , · · · , m at time t = 0. Consequently, the relation (2.41) implies Q α = ˜ S α iso 0 − ˆ S α , with ˜ S α iso 0 := ˜ S α iso | t =0 . Therefore, the constant tensor ˆ S α is determined by the initial values of the stresses,ˆ S α := ˜ S α iso 0 − Q α , (2.44)which ensures the consistency of the initial condition. With (2.39) and (2.41), the non-equilibrium stressesfor finite linear viscoelasticity are given by˜ S α neq = 2 (cid:32) ∂F α ( ˜ C , Θ) ∂ ˜ C − ∂ ˜ S α iso ∂ ˜ C : Γ α − I (cid:33) = 2 ∂F α ( ˜ C , Θ) ∂ ˜ C − µ α (Θ) ∂ ˜ S α iso ∂ ˜ C : (cid:16) ˜ S α iso − ˆ S α − Q α (cid:17) . Remark 5.
The evolution equations (2.42) derived here are similar to those introduced and used in [36, 37,38, 42, 46, 90], which can be expressed as ddt Q α + Q α τ α = ddt (cid:16) J − P : ˜ S α iso (cid:17) − dµ α dt Γ α . We emphasize that the evolution equations in those works were proposed based on a purely heuristic argumentwith inspiration coming from the standard linear solid model [91]. The right-hand side terms of the aboveevolution equations are different from those in (2.42) , in which a deviatoric projected stress from ˜ S α iso wasutilized to drive the evolution of Q α .2.4.2. A model based on the identical polymer chain assumption Now for a finite linear viscoelasticity model, considering the form for H α given in Definition 1, themodeling work reduces to designing the form of G α and F α in (2.36). In [40, 46], a model was proposed bytaking G α ( ˜ C , Θ) = F α ( ˜ C , Θ) = β α ∞ G ∞ iso ( ˜ C , Θ) , with β ∞ α ∈ (0 , ∞ ) being non-dimensional constants. This model is motivated by the observation that thephenomenological viscoelastic behavior is induced by a medium composed of identical polymer chains. Withthis assumption, the viscoelastic material is completely characterized by the potential G ∞ iso together withnon-dimensional parameters β ∞ α . The configurational free energy for this model can be rewritten asΥ α (cid:16) ˜ C , Θ , Γ , · · · , Γ m (cid:17) = µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12) Γ α − I (cid:12)(cid:12)(cid:12)(cid:12) + (cid:16) ˆ S α − β α ∞ ˜ S ∞ iso (cid:17) : Γ α − I β α ∞ G ∞ iso ( ˜ C , Θ) . (2.45)The evolution equations can be expressed as ddt Q α + Q α τ α = β α ∞ ddt ˜ S ∞ iso − dµ α dt ( Γ α − I ) . (2.46)12ith Q α calculated, one may obtain the non-equilibrium stresses as˜ S α neq = β α ∞ ˜ S ∞ iso − β α ∞ µ α (Θ) ∂ ˜ S ∞ iso ∂ ˜ C : (cid:16) β α ∞ ˜ S ∞ iso − ˆ S α − Q α (cid:17) , (2.47)with ˆ S α = β α ∞ ˜ S ∞ iso 0 − Q α and ˜ S ∞ iso 0 := ˜ S ∞ iso | t =0 . For convenience, we introduce the fictitious elasticity tensor˜ C ∞ iso defined as ˜ C ∞ iso := 2 J − ∂ ˜ S ∞ iso ∂ ˜ C = 4 J − ∂ G ∞ iso ( ˜ C , Θ) ∂ ˜ C ∂ ˜ C , (2.48)with which the non-equilibrium stresses (2.47) can be expressed as˜ S α neq = β α ∞ ˜ S ∞ iso − β α ∞ µ α (Θ) J ˜ C ∞ iso : (cid:16) β α ∞ ˜ S ∞ iso − ˆ S α − Q α (cid:17) . As will be revealed in Section 2.4.3, in the thermodynamic equilibrium limit, the non-equilibrium stressesgiven by this model will not vanish in general. This behavior is counter-intuitive as the stress from thedissipative potential is expected to be fully relaxed in the limit from the material modeling perspective. Inour numerical experiences, a non-vanishing S α neq often leads the body to deform to an unexpected state andsometimes gives unstable material behavior (see Figures 6 and 7). It is therefore necessary to discuss thecondition characterizing the relaxation of the non-equilibrium stresses. S α neq in the thermodynamic equilibriumlimit Intuitively, the non-equilibrium stress S α neq (2.39) shall vanish in the thermodynamic equilibrium limit.This property can be conveniently achieved in the multiplicative viscoelasticity theory [25]. Without havingthe multiplicative decomposition of the deformation gradient, we need to characterize this property byanalyzing necessary and sufficient conditions for the vanishment of S α neq in the equilibrium limit. First, therelaxation of the non-equilibrium stresses in the limit is given by O = S α neq (cid:12)(cid:12)(cid:12) eq = P : ˜ S α neq (cid:12)(cid:12)(cid:12) eq . It can be proved that the above relation is equivalent to the vanishment of the fictitious stresses ˜ S α neq , O = ˜ S α neq (cid:12)(cid:12)(cid:12) eq , with the proof detailed in Appendix A. Recall that the equilibrium [42] is characterized by, Q α (cid:12)(cid:12)(cid:12) eq = O , or equivalently ddt Γ α (cid:12)(cid:12)(cid:12) eq = O , for α = 1 , · · · , m. For finite linear viscoelastic materials, the above conditions imply that O = Q α (cid:12)(cid:12)(cid:12) eq = ˜ S α iso (cid:12)(cid:12)(cid:12) eq − ˆ S α − µ α (Θ) Γ α − I (cid:12)(cid:12)(cid:12) eq , which is equivalent to Γ α − I (cid:12)(cid:12)(cid:12) eq = 12 µ α (Θ) (cid:18) ˜ S α iso (cid:12)(cid:12)(cid:12) eq − ˆ S α (cid:19) . S α neq (cid:12)(cid:12)(cid:12) eq = 2 (cid:32) ∂F α ( ˜ C , Θ) ∂ ˜ C (cid:12)(cid:12)(cid:12) eq − ∂ ˜ S α iso ∂ ˜ C (cid:12)(cid:12)(cid:12) eq : Γ α − I (cid:12)(cid:12)(cid:12) eq (cid:33) = 2 (cid:32) ∂F α ( ˜ C , Θ) ∂ ˜ C (cid:12)(cid:12)(cid:12) eq − µ α (Θ) ∂ ˜ S α iso ∂ ˜ C (cid:12)(cid:12)(cid:12) eq : (cid:16) ˜ S α iso − ˆ S α (cid:17) (cid:12)(cid:12)(cid:12) eq (cid:33) = 2 ∂∂ ˜ C (cid:18) F α ( ˜ C , Θ) − µ α (Θ) (cid:12)(cid:12)(cid:12) ˜ S α iso − ˆ S α (cid:12)(cid:12)(cid:12) (cid:19) (cid:12)(cid:12)(cid:12) eq . The last equality in the above derivation makes use of the symmetry property of ∂ ˜ S α iso /∂ ˜ C . Therefore, the necessary condition for ensuring ˜ S α neq (cid:12)(cid:12)(cid:12) eq = O is O = ∂∂ ˜ C (cid:18) F α ( ˜ C , Θ) − µ α (Θ) (cid:12)(cid:12)(cid:12) ˜ S α iso − ˆ S α (cid:12)(cid:12)(cid:12) (cid:19) (cid:12)(cid:12)(cid:12) eq , or equivalently, after a partial integration with respect to ˜ C , F α ( ˜ C , Θ) (cid:12)(cid:12)(cid:12) eq = (cid:18) µ α (Θ) (cid:12)(cid:12)(cid:12) ˜ S α iso − ˆ S α (cid:12)(cid:12)(cid:12) + T α (Θ) (cid:19) (cid:12)(cid:12)(cid:12) eq , wherein T α is an arbitrary function of the temperature. The configurational free energy Υ α (2.36) for finitelinear viscoelasticity in the equilibrium limit isΥ α (cid:16) ˜ C , Θ , Γ , · · · , Γ m (cid:17) (cid:12)(cid:12)(cid:12) eq = H α ( Γ α , Θ) (cid:12)(cid:12)(cid:12) eq + (cid:16) ˆ S α − ˜ S α iso (cid:17) : Γ α − I (cid:12)(cid:12)(cid:12) eq + F α ( ˜ C , Θ) (cid:12)(cid:12)(cid:12) eq = (cid:32) µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12) Γ α − I (cid:12)(cid:12)(cid:12)(cid:12) − µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12) Γ α − I (cid:12)(cid:12)(cid:12)(cid:12) (cid:33) (cid:12)(cid:12)(cid:12) eq + (cid:18) µ α (Θ) (cid:12)(cid:12)(cid:12) ˜ S α iso − ˆ S α (cid:12)(cid:12)(cid:12) + T α (Θ) (cid:19) (cid:12)(cid:12)(cid:12) eq = − µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12) Γ α − I (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) eq + µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12) Γ α − I (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) eq + T α (Θ) (cid:12)(cid:12)(cid:12) eq = T α (Θ) (cid:12)(cid:12)(cid:12) eq . (2.49)This suggests that, for the finite linear viscoelastic model considered, the mechanical part of the config-urational free energy will eventually vanish in the thermodynamic equilibrium limit, if we demand thenon-equilibrium stresses to vanish in the limit. Based on the above discussion, a convenient modeling choicecan be made by demanding that F α always satisfies the following relation, F α ( ˜ C , Θ) = 14 µ α (Θ) (cid:12)(cid:12)(cid:12) ˜ S α iso − ˆ S α (cid:12)(cid:12)(cid:12) + T α (Θ) , (2.50)from which we readily have˜ S α neq = 2 ∂ ˜ S α iso ∂ ˜ C : (cid:18) µ α (Θ) (cid:16) ˜ S α iso − ˆ S α (cid:17) − Γ α − I (cid:19) = 1 µ α (Θ) ∂ ˜ S α iso ∂ ˜ C : Q α . (2.51)Therefore, the non-equilibrium stresses S α neq automatically vanish when reaching the equilibrium limit, andthe relation (2.50) serves as a sufficient condition for the configurational free energy that guarantees thestress relaxation. In this work, we adopt (2.50) for the analytic form of F α , thereby leaving G α as the onlyundetermined part of the configurational free energy for modelers.14 emark 6. In the literature, the configurational free energy is often proposed with F α = G α in (2.36) [42, 46]. Based on the above analysis, the non-equilibrium stress arising from that configurational freeenergy will not vanish in the thermodynamic equilibrium limit, except for one special case to be discussed inthe next section.2.4.4. The Holzapfel-Simo-Saint Venant-Kirchhoff model For finite linear viscoelastic materials that satisfy the condition (2.50), the stresses are completely de-termined in terms of G α . Here, as an illustrative example, we consider a special form of G α , G α ( ˜ C , Θ) = µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ C − I (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (2.52)Based on (2.37), one readily has˜ S α iso = 2 ∂G α ( ˜ C , Θ) ∂ ˜ C = µ α (Θ) (cid:16) ˜ C − I (cid:17) , and Q α = ˜ S α iso − ˆ S α − µ α (Θ) ( Γ α − I ) = µ α (Θ) (cid:16) ˜ C − Γ α (cid:17) + Q α . The evolution eqation (2.42) can be written as ddt Q α + Q α τ α = ddt ˜ S α iso − dµ α (Θ) dt ( Γ α − I ) = µ α (Θ) ddt ˜ C + dµ α (Θ) dt (cid:16) ˜ C − I (cid:17) − dµ α (Θ) dt ( Γ α − I )= µ α (Θ) ddt ˜ C + dµ α (Θ) dt (cid:16) ˜ C − Γ α (cid:17) . The fictitious non-equilibrium stress ˜ S α neq can be represented as˜ S α neq = 1 µ α (Θ) ∂ ˜ S α iso ∂ ˜ C : Q α = 1 µ α (Θ) µ α (Θ) Q α = Q α . Consequently, we have the fictitious second Piola-Kirchhoff stress represented as˜ S = ˜ S ∞ iso + m (cid:88) α =1 ˜ S α neq = ˜ S ∞ iso + m (cid:88) α =1 Q α . This recovers one specific model proposed by Holzapfel and Simo in [46, Section 4.2], which can be viewedas a generalization of the Saint Venant-Kirchhoff model to the viscous regime. This is the reason for itsname used here. We should note that directly using Q α in place of ˜ S α neq is fairly common in the literature[37, 41, 42, 44], and the above analysis shows that this choice implicitly implies a special form of theconfigurational free energy. We note that, to ensure ˜ S α neq = Q α , one needs ∂ ˜ S α iso ∂ ˜ C = 2 ∂G α ( ˜ C , Θ) ∂ ˜ C ⊗ ˜ C = µ α (Θ) I , which implies a quadratic form of G α in terms of ˜ C . Therefore, enforcing ˜ S α neq = Q α implicitly requires G α to be quadratic in terms of ˜ C . This also suggests a linear relationship between ˜ S α iso and ˜ C . For generalnonlinear models, apparently one should distinguish Q α from ˜ S α neq . Furthermore, according to (2.50), wehave F α ( ˜ C , Θ) = 14 µ α (Θ) (cid:12)(cid:12)(cid:12) ˜ S α iso − ˆ S α (cid:12)(cid:12)(cid:12) + T α (Θ) = 14 µ α (Θ) (cid:12)(cid:12)(cid:12) µ α (Θ) (cid:16) ˜ C − I (cid:17) − ˆ S α (cid:12)(cid:12)(cid:12) + T α (Θ) . The configurational free energy (2.36) can be written asΥ α (cid:16) ˜ C , Θ , Γ , · · · , Γ m (cid:17) = H α ( Γ α , Θ) + (cid:32) ˆ S α − ∂G α ( ˜ C , Θ) ∂ ˜ C (cid:33) : Γ α − I F α ( ˜ C , Θ)15 µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12) Γ α − I (cid:12)(cid:12)(cid:12)(cid:12) + (cid:32) ˆ S α − µ α (Θ) ˜ C − I (cid:33) : Γ α − I
2+ 14 µ α (Θ) (cid:12)(cid:12)(cid:12) µ α (Θ) (cid:16) ˜ C − I (cid:17) − ˆ S α (cid:12)(cid:12)(cid:12) + T α (Θ)= µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ C − Γ α (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ˆ S α : ˜ C − Γ α µ α (Θ) (cid:12)(cid:12)(cid:12) ˆ S α (cid:12)(cid:12)(cid:12) + T α (Θ)= 14 µ α (Θ) (cid:12)(cid:12)(cid:12) µ α (Θ) (cid:16) ˜ C − Γ α (cid:17) − ˆ S α (cid:12)(cid:12)(cid:12) + T α (Θ) . (2.53)One can see that for the model considered above, the configurational free energy is a function of ( ˜ C − Γ α ) / Remark 7.
It is convenient to make a rheological interpretation of this model. It can be seen that in thismodel ˜ S α iso , = O as ˜ C | t =0 = I . If we choose Q α = O , then ˆ S α = O according to (2.44) . We may introducean ininitesimal strain ε as an approximation of (cid:16) ˜ C − I (cid:17) / , ε ≈ ˜ C − I . Similarily, there is a small-strain approximation of ( Γ α − I ) / , ε α v ≈ ˜ Γ α − I . If we denote ε αe := ε − ε α v , it can be shown that under the infinitesimal strain approximation, the configura-tional free energy (2.53) can be approximated by Υ α ≈ µ α (Θ) | ε αe | + T α (Θ) . This can be interpreted as the elastic energy of the springs within the Maxwell element of the standard Zenermodel.2.4.5. A modified model based on the identical-polymer-chain assumption
From the analysis in Section 2.4.3, we can see that the identical polymer chain model stated in Section2.4.2 will engender counter-intuitive behavior of the non-equilibrium stress. In this section, we revisit theidentical polymer chain model, by adopting the form (2.50) for F α . The remaining component for theconfigurational free energy is given by G α ( ˜ C , Θ) = β ∞ α G ∞ iso ( ˜ C , Θ) , (2.54)with β ∞ α ∈ (0 , ∞ ) being non-dimensional constants [36, 40]. With this modeling assumption, the materialbehavior is completely characterized by the form of G ∞ iso , and we have˜ S α iso = 2 ∂G α ( ˜ C , Θ) ∂ ˜ C = β ∞ α ∂G ∞ iso ( ˜ C , Θ) ∂ ˜ C = β ∞ α ˜ S ∞ iso . With this elasticity tensor, the non-equilibrium stress (2.51) can be represented as˜ S α neq = 1 µ α ∂ ˜ S α iso ∂ ˜ C : Q α = J β ∞ α µ α ˜ C ∞ iso : Q α = 2 β ∞ α µ α ∂ G ∞ iso ( ˜ C , Θ) ∂ ˜ C ∂ ˜ C : Q α . emark 8. In the literature, the terminology “ground-stress” is often used to refer the fictitious stress S ∞ iso [97, 98]; “over-stress” is used to refer Q α [25, 39] or S α neq [98]. The prevailing custom is to use Q α and ˜ S α neq interchangeably [36, 40, 44]. However, the above analysis indicates that, for a general nonlinear formof G α like (2.54) , ˜ S α neq (cid:54) = Q α . Remark 9.
Interestingly, if the functional form of G ∞ iso is given by the Neo-Hookean model, we have ∂ G ∞ iso ( ˜ C , Θ) ∂ ˜ C ∂ ˜ C = O , which indicates that S α neq = O . In other words, under the identical-polymer-chain assumption, the Neo-Hookean material cannot have viscous stress in finite linear viscoelasticity. We may summarize the thermodynamic potential and the resulting constitutive relations for finite linearthermal-visco-hyperelastic materials discussed in Section 2 as follows.
Constitutive laws for finite linear thermal-visco-elasticity
Gibbs free energy: G ( ˜ C , P, Θ , Γ , · · · , Γ m ) = G ∞ vol ( P, Θ) + G ∞ iso ( ˜ C , Θ) + m (cid:88) α =1 Υ α (cid:16) ˜ C , Θ , Γ , · · · , Γ m (cid:17) , Υ α (cid:16) ˜ C , Θ , Γ , · · · , Γ m (cid:17) = H α ( Γ α , Θ) − ∂G α ( ˜ C , Θ) ∂ ˜ C : Γ α − I S α : Γ α − I F α ( ˜ C , Θ) ,H α ( Γ α , Θ) = µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12) Γ α − I (cid:12)(cid:12)(cid:12)(cid:12) , F α ( ˜ C , Θ) = 14 µ α (Θ) (cid:12)(cid:12)(cid:12) ˜ S α iso − ˆ S α (cid:12)(cid:12)(cid:12) + T α (Θ) . Second Piola-Kirchhoff stress: S = S iso + S vol , S iso = J − P : ˜ S , S vol = − JP C − , ˜ S = ˜ S ∞ iso + m (cid:88) α =1 ˜ S α neq , ˜ S ∞ iso = 2 ∂G ∞ iso ∂ ˜ C , ˜ S α neq = 1 µ α (Θ) ∂ ˜ S α iso ∂ ˜ C : Q α , ˜ S α iso = 2 ∂G α ∂ ˜ C , ˜ C ∞ iso = 4 J − ∂ G ∞ iso ∂ ˜ C ∂ ˜ C . Cauchy stress: σ = σ dev − p I , σ dev = 1 J F S iso F T . Entropy: S = S ∞ + m (cid:88) α =1 S α , S ∞ := − ∂G ∞ vol ∂ Θ − ∂G ∞ iso ∂ Θ , S α := − ∂ Υ α ∂ Θ . Dissipation: D = ¯ κ Θ |∇ X Θ | + m (cid:88) α =1 η α (cid:12)(cid:12)(cid:12)(cid:12) ddt Γ α (cid:12)(cid:12)(cid:12)(cid:12) = ¯ κ Θ |∇ X Θ | + m (cid:88) α =1 η α | Q α | . Q α : ddt Q α + 1 τ α Q α = ddt ˜ S α iso − dµ α (Θ) dt ( Γ α − I ) , Q α | t =0 = Q α . Hereditary integral form for Q α : Q α = exp ( − t/τ α ) Q α + (cid:90) t + exp ( − ( t − s ) /τ α ) (cid:18) dds ˜ S α iso − ( Γ α − I ) dµ α (Θ) ds (cid:19) ds.
3. Numerical formulation
In this section, we develop a suite of discretization methods for the material model developed in Section2 by restricting our discussion to fully incompressible materials under the isothermal condition. We willprove that the proposed numerical formulation is stable in energy. Henceforth, we use the acronyms IPC,HS, and MIPC models to stand for the identical polymer chain model stated in Section 2.4.2, the Holzapfel-Simo-Saint Venant-Kirchhoff model stated in Section 2.4.4, and the modified identical polymer chain modelstated in Section 2.4.5, respectively.
Under isothermal conditions, the energy equation is decoupled from the mechanical system, and themotion of an incompressible continuum body is governed by the following system of equations, = d u dt − v , in Ω t x , (3.1)0 = ∇ x · v , in Ω t x , (3.2) = ρ ( J ) d v dt − ∇ x · σ dev + ∇ x p − ρ ( J ) b , in Ω t x . (3.3)The boundary Γ t x := ∂ Ω t x enjoys a non-overlapping subdivision: Γ t x = Γ g,t x ∪ Γ h,t x , wherein Γ g,t x and Γ h,t x arethe Dirichlet and Neumann parts of the boundary, respectively. The boundary conditions can be stated as u = g , on Γ g,t x , v = d g dt , on Γ g,t x , ( σ dev − p I ) n = h , on Γ h,t x . (3.4)Given the initial data u , p , and v , the initial-boundary value problem is to solve for u , p , and v thatsatisfy (3.1)-(3.4), and u ( x ,
0) = u ( x ) , p ( x ,
0) = p ( x ) , v ( x ,
0) = v ( x ) . (3.5)The equations (3.1)-(3.5) constitute an initial-boundary value problem. The constitutive relations for thematerial has been given in Section 2. For fully incompressible materials, it can be shown that the volumetricpart of the free energy adopts the form G ∞ vol = P/ρ [55], which leads to ρ = ρ . In the formulation, we adopta modified constitutive relation for the density, that is ρ ( J ) = ρ /J [60]. Apparently, at the continuum level,this relation is equivalent to ρ ( J ) = ρ , as the divergence-free condition for the velocity guarantees J = 1.Yet, we note that at the discrete level, the condition J = 1 in general rarely holds in the pointwise sense[54]. The adoption of the relation ρ ( J ) = ρ /J facilitates our discussion of the nonlinear numerical stability,as will be shown in the next section. 18 .2. Spline spaces Before introducing the semi-discrete formulation, we state the construction of B-splines and NURBSbasis functions, which are utilized to build the discrete function spaces. Given the polynomial degree p andthe dimensionality of the B-spline space n , the knot vector can be represented by Ξ := { ξ , ξ · · · , ξ n + p +1 } ,wherein 0 = ξ ≤ ξ ≤ · · · ≤ ξ n + p +1 = 1. The B-spline basis functions of degree p , denoted as N p i ( ξ ), for i = 1 , · · · , n , can then be defined recursively from the knot vector using the Cox-de Boor recursion formula[65, Chapter 2]. The definition starts with the case of p = 0, where the basis functions are defined aspiecewise constants, N i ( ξ ) = (cid:40) ξ i ≤ ξ < ξ i +1 , . For p ≥
1, the basis functions are defined recursively as N p i ( ξ ) = ξ − ξ i ξ i + p − ξ i N p − i ( ξ ) + ξ i + p +1 − ξξ i + p +1 − ξ i +1 N p − i +1 ( ξ ) . Given a set of weights { w , w , · · · , w n } , the NURBS basis functions of degree p can be defined as R p i ( ξ ) := w i N p i ( ξ ) W ( ξ ) , and W ( ξ ) := n (cid:88) j =1 w j N p j ( ξ ) . Importantly, the knots can also be represented with two vectors, one of the unique knots { ζ , ζ , · · · , ζ m } and another of the corresponding knot multiplicities { r , r , · · · , r m } . As is standard in the literature ofcomputer-aided design, we consider open knot vectors in this work, meaning that r = r m = p + 1. Wefurther assume that r i ≤ p for i = 2 , · · · , m −
1. Across any given knot ζ i , the B-spline basis functions have α i := p − r i continuous derivatives. The vector α := { α , α , · · · , α m − , α m } = {− , α , · · · , α m − , − } isreferred to as the regularity vector. A value of − α i indicates discontinuity of the basis functions at ζ i . We introduce the function space R p α := span { R p i } ni =1 , where the notation R p α is used to indicate that α i = α for i = 2 , · · · , m −
1, suggesting continuity C α for the spline function spaces. The construction ofmultivariate B-spline and NURBS basis functions follows a tensor-product manner. For l = 1 , ,
3, given p l , n l , and the knot vectors Ξ l = { ξ ,l , ξ ,l , · · · , ξ n l + p l +1 ,l } , the univariate B-spline basis functions N p l i l ,l arewell-defined. Consequently, the multivariate B-spline basis functions can be defined by exploiting the tensorproduct structure, N p , p p i ,i ,i ( ξ , ξ , ξ ) := N p i , ( ξ ) ⊗ N p i , ( ξ ) ⊗ N p i , ( ξ ) , for i l = 1 , , · · · , n l and l = 1 , , . Given a set of weights { w i ,i ,i } , the NURBS basis functions are defined by R p , p , p i ,i ,i ( ξ , ξ , ξ ) := w i ,i ,i N p , p p i ,i ,i ( ξ , ξ , ξ ) W ( ξ , ξ , ξ ) , W ( ξ , ξ , ξ ) := (cid:88) l =1 n l (cid:88) i l =1 w i ,i ,i N p , p p i ,i ,i ( ξ , ξ , ξ ) . Correspondingly, the NURBS function space is defined as R p , p , p α , α , α := span { R p , p , p i ,i ,i } n , n , n i =1 ,i =1 ,i =1 . We first define two discrete function spaces on ˆΩ := (0 , ,ˆ S h := R p + a , p + a , p + a α + b , α + b , α + b × R p + a , p + a , p + a α + b , α + b , α + b × R p + a , p + a , p + a α + b , α + b , α + b , ˆ P h := R p , p , p α , α , α , X ! GX ! HX + tx X = ' ! t ( x ) x = ' t ( X ) ^ + ' t ' ! t ! g;tx ! h;tx A ' t / A Figure 2: Illustration of the referential and current configurations with the boundary subdivisions. The spline functions aredefined on ˆΩ and are pushed forward to the two configurations via the mapping ψ and ϕ t ◦ ψ . with integer parameters 1 ≤ a and 0 ≤ b ≤ a . Assuming the referential configuration of the body can beparametrized by a geometrical mapping ψ : ˆΩ → Ω X . The boundary of Ω X can be partitioned into twonon-overlapping subdivisions as ∂ Ω x := Γ G X ∪ Γ H X , and the two subdivisions satisfy Γ G X = ϕ − t (Γ g,t x ) andΓ H X = ϕ − t (cid:0) Γ h,t x (cid:1) . The relation between the two configurations and the boundary subdivisions are illustratedin Figure 2. The discrete function spaces on Ω X can be defined through the pull-back operation, S h := { w : w ◦ ψ ∈ ˆ S h } , P h := { q : q ◦ ψ ∈ ˆ P h } . With the spaces S h and P h defined, we may specify the trial solution spaces on the referential configurationas S U h = (cid:110) U h : U h ( · , t ) ∈ S h , t ∈ [0 , T ] , U h ( · , t ) = G on Γ G X (cid:111) , S P h = (cid:110) P h : P h ( · , t ) ∈ P h , t ∈ [0 , T ] (cid:111) , S V h = (cid:26) V h : V h ( · , t ) ∈ S h , t ∈ [0 , T ] , V h ( · , t ) = d G dt on Γ G X (cid:27) , and the corresponding test function spaces are defined as V P h = (cid:110) Q h : Q h ( · , t ) ∈ P h , t ∈ [0 , T ] (cid:111) , V V h = (cid:8) W h : W h ( · , t ) ∈ S h , t ∈ [0 , T ] , W h ( · , t ) = on Γ G X (cid:9) . Given the displacement U h , the placement field is given by ϕ h = U h + X . Consequently, we may also statethe trial solution space defined on the current configuration as S u h = (cid:110) u h : u h ◦ ϕ h ∈ S h , t ∈ [0 , T ] , u h ( · , t ) = g on Γ g,t x (cid:111) , S p h = (cid:110) p h : p h ◦ ϕ h ∈ P h , t ∈ [0 , T ] (cid:111) , S v h = (cid:26) v h : v h ◦ ϕ h ∈ S h , t ∈ [0 , T ] , v h ( · , t ) = d g dt on Γ g,t x (cid:27) , V p h = (cid:110) q h : q h ◦ ϕ h ∈ P h , t ∈ [0 , T ] (cid:111) , V v h = (cid:110) w h : w h ◦ ϕ h ∈ S h , t ∈ [0 , T ] , w h ( · , t ) = on Γ g,t x (cid:111) . With the discrete function defined above, the semi-discrete formulation on the current configuration can bestated as follows. Find y h ( t ) := { u h ( t ) , p h ( t ) , v h ( t ) } T ∈ S u h × S p h × S v h such that for t ∈ [0 , T ], = B k ( ˙ y h , y h ) := d u h dt − v h , (3.6)0 = B p ( q h ; ˙ y h , y h ) := (cid:90) Ω t x q h ∇ x · v h d Ω x , (3.7)0 = B m ( w h ; ˙ y h , y h ) := (cid:90) Ω t x w h · ρ ( J h ) d v h dt + ∇ x w h : σ dev − ∇ x · w h p h − w h · ρ ( J h ) b d Ω x − (cid:90) Γ h,t x w h · h d Γ x , (3.8)for ∀ { q h , w h } ∈ V p h × V v h , with y h (0) := { u h , p h , v h } T . Here u h , p h , and v h are the L projec-tions of the initial data onto the finite dimensional trial solution spaces. Alternatively, the semi-discreteformulation can be pulled back to the referential configuration, which can be stated as follows. Find Y h ( t ) := { U h ( t ) , P h ( t ) , V h ( t ) } T ∈ S U h × S P h × S V h such that for t ∈ [0 , T ], = B k (cid:16) ˙ Y h , Y h (cid:17) := d U h dt − V h , (3.9)0 = B p (cid:16) Q h ; ˙ Y h , Y h (cid:17) := (cid:90) Ω X Q h J h ∇ X V h : F − Th d Ω X , (3.10)0 = B m (cid:16) W h ; ˙ Y h , Y h (cid:17) := (cid:90) Ω X (cid:16) W h · ρ d V h dt + ∇ X W h : (cid:0) J h σ dev F − Th (cid:1) − J h P h ∇ X W h : F − Th − W h · ρ B (cid:17) d Ω X − (cid:90) Γ H X W h · H d Γ X , (3.11)for ∀ { Q h , W h } ∈ V P h × V V h , with Y h (0) = { U h , P h , V h } T . Here the initial data are related by U h = u h ◦ ϕ t , P h = p h ◦ ϕ t , and V h = v h ◦ ϕ t . It can be shown that the above semi-discrete formulationinherits the dissipation property from the continuum model from the following theorem. Theorem 1 (A priori energy stability) . Assuming the Dirichlet boundary data G is time independent, thesolutions of the semi-discrete problem (3.9) - (3.11) satisfy ddt (cid:90) Ω X ρ (cid:107) V h (cid:107) + G ∞ iso ( ˜ C h ) + m (cid:88) α =1 Υ α ( ˜ C h , Γ αh ) d Ω X = (cid:90) Ω X ρ V h · B d Ω X + (cid:90) Γ X V h · H d Γ X − m (cid:88) α =1 (cid:90) Ω X (cid:18) ddt Γ αh (cid:19) : V α : (cid:18) ddt Γ αh (cid:19) d Ω X . (3.12) Proof.
Since the Dirichlet boundary data G is independent of time, one is allowed to choose Q h = P h and W h = V h , which leads to0 = B p (cid:16) P h ; ˙ Y h , Y h (cid:17) + B m (cid:16) V h ; ˙ Y h , Y h (cid:17) = (cid:90) Ω X P h J h ∇ X V h : F − Th d Ω X + (cid:90) Ω X V h · ρ d V h dt + ∇ X V h : (cid:0) J h σ dev F − Th (cid:1) − J h P h ∇ X V h : F − Th − V h · ρ B d Ω x − (cid:90) Γ H X V h · H d Γ x ddt (cid:90) Ω X ρ (cid:107) V h (cid:107) d Ω X + (cid:90) Ω X ddt F h : ∂ (cid:32) G ∞ iso ( ˜ C h ) + m (cid:88) α =1 Υ α ( ˜ C h , Γ αh ) (cid:33) ∂ F h − V h · ρ B d Ω X − (cid:90) Γ X V h · H d Γ X = ddt (cid:90) Ω X ρ (cid:107) V h (cid:107) d Ω X + ddt (cid:90) Ω X (cid:32) G ∞ iso ( ˜ C h ) + m (cid:88) α =1 Υ α ( ˜ C h , Γ αh ) (cid:33) d Ω X + (cid:90) Ω X m (cid:88) α =1 Q αh : 12 ddt Γ αh − V h · ρ B d Ω X − (cid:90) Γ X V h · H d Γ X = ddt (cid:90) Ω X ρ (cid:107) V h (cid:107) + G ∞ iso ( ˜ C h ) + m (cid:88) α =1 Υ α ( ˜ C h , Γ αh ) d Ω X + (cid:90) Ω X m (cid:88) α =1 (cid:18) ddt Γ αh (cid:19) : V α : (cid:18) ddt Γ αh (cid:19) − V h · ρ B d Ω X − (cid:90) Γ X V h · H d Γ X . Rearranging terms in the above equality results in ddt (cid:90) Ω X ρ (cid:107) V h (cid:107) + G ∞ iso ( ˜ C h ) + m (cid:88) α =1 Υ α ( ˜ C h , Γ αh ) d Ω X = (cid:90) Ω X ρ V h · B d Ω X + (cid:90) Γ X V h · H d Γ X (3.13) − m (cid:88) α =1 (cid:90) Ω X (cid:18) ddt Γ αh (cid:19) : V α : (cid:18) ddt Γ αh (cid:19) d Ω X , (3.14)which completes the proof.This numerical stability property guarantees that the proposed modeling and computational frameworkpreserve critical structures of the physical system. It is worth emphasizing that the energy stable scheme isconstructed based on the thermodynamically consistent continuum model derived in Section 2. To the bestof our knowledge, we are unaware of any other numerical scheme with a priori energy stability proved forincompressible viscoelasticity. A straightforward consequence of this stability property is that the energywill be monotonically decreasing for unforced mechanical systems (i.e. B = and H = ).For compressible materials, one may show that a pressure-squared term enters into the definition of theenergy, providing boundedness of the pressure field. This justifies the use of equal-order interpolation forcompressible materials. For fully incompressible materials, the pressure force has no contribution to theenergy, and its stability comes from the inf-sup condition [60, Section 3.2]. The inf-sup stability of thediscrete spaces given in Section 3.3 has been numerically examined in [60, Section 4.1], where it was foundthat b + 1 ≤ a guarantees the element pair to be inf-sup stable. In this work, we use a = 1 and b = 0 forthe discrete function spaces.In addition to the energy stability, one may also conveniently show the momentum conservation of thesemidiscrete formulatoin, which is stated in the following theorem. Theorem 2 (Semidiscrete momentum conservation) . If Γ t x = Γ h,t x , the following conservation propertieshold for the semidiscrete formulation (3.6) - (3.8) , ddt (cid:90) Ω X ρ V h d Ω X = (cid:90) Ω X ρ B d Ω X + (cid:90) Γ X H d Γ X ,ddt (cid:90) Ω X ρ ϕ h × V h d Ω X = (cid:90) Ω X ρ ϕ h × B d Ω X + (cid:90) Γ X ϕ h × H d Γ X . Proof.
The above conservation properties are direct consequences of choosing w h = e i and w h = e i × ϕ h respectively in (3.8), where e i is a unit vector in the i -th direction.22 emark 10. It can be shown that the dissipation term in the stability estimate (3.12) can be equivalentlywritten as m (cid:88) α =1 (cid:90) Ω X (cid:18) ddt Γ αh (cid:19) : V α : (cid:18) ddt Γ αh (cid:19) d Ω X = m (cid:88) α =1 (cid:90) Ω X Q α : ( V α ) − : Q α d Ω X . It should be pointed out that the stability is analyzed for the semi-discrete scheme. It remains an intriguingtopic to further extend this estimate to the fully discrete regime. The energy-momentum scheme [99, 100,101] is a promising candidate for this goal.3.4. Time integration algorithm
In this section, we first introduce a discrete algorithm for updating the stresses. Following that, we statethe fully discrete scheme using the generalized- α method. To obtain the stress, we need to perform time integration for the constitutive laws. Let the time interval(0 + , T ] be divided into N subintervals of size ∆ t n := t n +1 − t n delimited by a discrete time vector { t n } Nn =0 .The approximations to the velocity, pressure, and displacement and their first time derivatives at time t n are denoted as Y n := { V n , P n , U n } T and ˙ Y n := (cid:110) ˙ V n , ˙ P n , ˙ U n (cid:111) T , respectively. Correspondingly, the approximations to the deformation gradient and strain measures at time t n are represented as F n = I + ∇ X U n , J n = det( F n ) , C n = F Tn F n , ˜ C n = J − / n C n . The approximated projection tensor and elasticity tensor are given by P n +1 = I − C − n +1 ⊗ C n +1 , ˜ C ∞ iso n +1 = 4 J − n +1 ∂ G ∞ iso ( ˜ C n +1 ) ∂ ˜ C ∂ ˜ C The algorithmic stresses at time t n +1 read as S n +1 = S iso n +1 + S vol n +1 , S iso n +1 = J − n +1 P n +1 : ˜ S n +1 , S vol n +1 = − J n +1 P n +1 C − n +1 , ˜ S n +1 = ˜ S ∞ iso n +1 + m (cid:88) α =1 ˜ S α neq n +1 , ˜ S ∞ iso n +1 = 2 (cid:18) ∂G ∞ iso ∂ ˜ C (cid:19) n +1 , ˜ S α neq n +1 = J n +1 β ∞ α µ α ˜ C ∞ iso n +1 : Q αn +1 . To evaluate the stresses, we need to provide an algorithmic way to evaluate Q αn +1 based on the hereditaryintegral (2.43). Following the notation introduced in [42, 44], we first introduce a dimensionless parameter ξ α := − ∆ t n / τ α . The approximation to the variable Q α at time t n +1 is given by Q αn +1 = exp ( − t n +1 /τ α ) Q α + (cid:90) t n +1 + exp ( − ( t n +1 − s ) /τ α ) β ∞ α dds ˜ S ∞ iso ds = exp(2 ξ α ) exp( − t n /τ α ) Q α + exp(2 ξ α ) (cid:90) t n + exp ( − ( t n − s ) /τ α ) β ∞ α dds ˜ S ∞ iso ds + (cid:90) t n +1 t n exp ( − ( t n +1 − s ) /τ α ) β ∞ α dds ˜ S ∞ iso ds = exp(2 ξ α ) Q αn + (cid:90) t n +1 t n exp ( − ( t n +1 − s ) /τ α ) β ∞ α dds ˜ S ∞ iso ds exp(2 ξ α ) Q αn + exp( ξ α ) β ∞ α (cid:90) t n +1 t n dds ˜ S ∞ iso ds = β ∞ α exp( ξ α ) ˜ S ∞ iso n +1 + exp( ξ α ) (cid:16) exp( ξ α ) Q αn − β ∞ α ˜ S ∞ iso n (cid:17) . (3.15)In the second-to-last step of (3.15), a mid-point rule is applied for the exponential term in the time integralto obtain an approximation [42, 44]. Error analysis shows that the above approximation is second-orderaccurate, and the relation (3.15) is often referred to as the recurrence update formula for Q α . There existsan alternate second-order accurate recurrence update formula by applying the mid-point rule to the stressrate, rather than the exponential kernel term, in the time integral. Interested readers are referred to [42,p. 355] for details. Remark 11.
We mention that the recursive formula (3.15) is obtained by making use of the semigroupproperty of the kernel of the hereditary integral, and there exist other recursive formulas [42, Chapter 10].The recursive formula (3.15) is second-order accurate, and it is feasible to achieve higher-order accuracywith more involved update formulas for the internal state variables [94]. In fractional-order viscoelasticitymodels, however, a different recursive formula can be derived by invoking the fast convolution method [16].
Remark 12.
In the analysis calculations, the values of Q αn +1 can be conveniently initialized and stored atquadrature points. In the constitutive routine, their values at the quadrature points are updated by (3.15) ,which are utilized as the input for the calculation of stresses at the quadrature points. Yet, in the postpro-cessing (such as visualization), the stresses are typically not sampled or interpolated at the quadrature pointsused for analysis. The global smoothing procedure [102, 103] can be invoked to recover the values of Q αn +1 and subsequently the values of stresses.3.4.2. Fully discrete scheme With the time discrete stress given above, we may state the fully discrete algorithm by invoking thegeneralized- α method [104, 105]. At time t n , given Y n and ˙ Y n , the time step size ∆ t n , and the parameters α m , α f , and γ , find Y n +1 and ˙ Y n +1 , such that for ∀ { Q h , W h } ∈ V P h × V V h , = B k (cid:16) ˙ Y n + α m , Y n + α f (cid:17) , (3.16)0 = B p (cid:16) Q h ; ˙ Y n + α m , Y n + α f (cid:17) , (3.17)0 = B m (cid:16) W h ; ˙ Y n + α m , Y n + α f (cid:17) , (3.18) Y n +1 = Y n + ∆ t n ˙ Y n + γ ∆ t n (cid:16) ˙ Y n +1 − ˙ Y n (cid:17) , (3.19)˙ Y n + α m = ˙ Y n + α m (cid:16) ˙ Y n +1 − ˙ Y n (cid:17) , (3.20) Y n + α f = Y n + α f ( Y n +1 − Y n ) . (3.21)Let (cid:37) ∞ ∈ [0 ,
1] denote the spectral radius of the amplification matrix at the highest mode. The fol-lowing choice of the parameters ensures second-order accuracy, unconditional stability, and controllablehigh-frequency dissipation for linear first-order ordinary differential equations [105], α m = 12 (cid:18) − (cid:37) ∞ (cid:37) ∞ (cid:19) , α f = 11 + (cid:37) ∞ , γ = 11 + (cid:37) ∞ . Remark 13.
It is known that the many different temporal schemes may be recovered by the generalized- α scheme via distinct choices of the parameters. For example, choosing (cid:37) ∞ = 0 . renders a scheme that isspectrally equivalent to the second-order backward difference method [105]; choosing (cid:37) ∞ = 1 . recovers themid-point rule. It is also worth pointing out that the generalized- α scheme has been conventionally appliedto second-order structural dynamics. Recent work shows that the generalized- α method applied to a first-order structural dynamic system does not suffer from the ‘overshoot’ phenomenon [106], and thus possesses any desirable properties of implicit schemes, as noted by Hilber and Hughes [107]. Writing the structuraldynamics problem as a first-order system introduces three additional velocity degrees of freedom per node. Itcan be shown that these additional degrees of freedom can be solved in a segregated manner in a consistentNewton-Raphson algorithm [55, 58, 60, 108]. The additional cost induced by the velocity degrees of freedomis merely the memory for storing them and an explicit update formula, which are thus marginal. Remark 14.
The stability of the generalized- α scheme was analyzed based on linear problems [104, 105]and remains unclear for nonlinear problems. In fact, it is known that the mid-point rule, an instantiationof the generalized- α scheme with (cid:37) ∞ = 1 . , is often energetically unstable for nonlinear structural dynamics[109, 110]. It remains an interesting topic to construct fully discrete schemes that are provably stable inenergy, and the family of energy-momentum methods serves as a promising candidate in this role [99, 109,110, 111] Remark 15.
It is worth pointing out that the generalized- α scheme may achieve the claimed second-orderaccuracy if all unknowns are collocated at the intermediate time step, as is done in (3.16) - (3.21) . In a verypopular approach, the pressure is collocated at the time step t n +1 with the rest variables collocated at theintermediate time step following the rule of the generalized- α scheme. It was found recently that the claimedsecond-order accuracy is lost in that approach [64].3.5. Elasticity tensor Here we provide the isochoric elasticity tensors that are used in the consistent linearization of the model.We first define the elasticity tensors approximated at time t n +1 as C iso n +1 := (cid:18) ∂ S iso ∂ C (cid:19) n +1 , C ∞ iso n +1 := (cid:18) ∂ S ∞ iso ∂ C (cid:19) n +1 , and C α neq n +1 := (cid:18) ∂ S α neq ∂ C (cid:19) n +1 . Based on the additive split structure of the stress, the isochoric elasticity tensor C iso n +1 can be expressedas C iso n +1 = (cid:18) ∂ S ∞ iso ∂ C (cid:19) n +1 + m (cid:88) α =1 (cid:18) ∂ S α neq ∂ C (cid:19) n +1 = C ∞ iso n +1 + m (cid:88) α =1 C α neq n +1 . It can be shown that C ∞ iso n +1 = P n +1 : ˜ C ∞ iso n +1 : P Tn +1 + 23 Tr (cid:16) J − n +1 ˜ S ∞ iso n +1 (cid:17) ˜ P n +1 − (cid:0) C − n +1 ⊗ S ∞ iso n +1 + S ∞ iso n +1 ⊗ C − n +1 (cid:1) , ˜ C ∞ iso n +1 := 4 J − n +1 (cid:18) ∂ G ∞ iso ∂ ˜ C ∂ ˜ C (cid:19) n +1 , Tr ( · ) = ( · ) : C n +1 , ˜ P n +1 := C − n +1 (cid:12) C − n +1 − C − n +1 ⊗ C − n +1 . The derivation of the formula can be found in [44, p. 255]. Following similar steps, we have C α neq n +1 = P n +1 : ˜ C α neq n +1 : P Tn +1 + 23 Tr (cid:16) J − ˜ S α neq n +1 (cid:17) ˜ P − (cid:0) C − n +1 ⊗ S α neq n +1 + S α neq n +1 ⊗ C − n +1 (cid:1) , where ˜ C α neq n +1 := 2 J − n +1 (cid:32) ∂ ˜ S α neq ∂ ˜ C (cid:33) n +1 . Based on (3.15), we have 2 (cid:18) ∂ Q α ∂ C (cid:19) n +1 = δ α J n +1 ˜ C ∞ iso n +1 , δ α := β ∞ α exp( ξ α ). We may now express ˜ C α neq n +1 as˜ C α neq n +1 = 4 β ∞ α J − n +1 µ α (cid:18) ∂ G ∞ iso ∂ ˜ C ∂ ˜ C ∂ ˜ C (cid:19) n +1 : Q αn +1 + 2 δ α β ∞ α µ α (cid:18) ∂ G ∞ iso ∂ ˜ C ∂ ˜ C (cid:19) n +1 : ˜ C ∞ iso n +1 = 4 β ∞ α J − n +1 µ α (cid:18) ∂ G ∞ iso ∂ ˜ C ∂ ˜ C ∂ ˜ C (cid:19) n +1 : Q αn +1 + δ α β ∞ α J n +1 µ α ˜ C ∞ iso n +1 : ˜ C ∞ iso n +1 . The algorithms for calculating the stresses and elasticity tensors for the IPC, HS, and MIPC models aredocumented in Appendix B, Appendix C, and Appendix D, respectively.
4. Numerical results
In this section, we investigate the proposed viscoelastic model by a suite of numerical examples usingthe numerical scheme proposed in the previous section. Unless otherwise specified, we use p + a + 1 Gaussquadrature points in each direction. Also recall that we have fixed a = 1 and b = 0 in the construction ofdiscrete function spaces. We take Q α = O in all numerical studies, which implies ˆ S α = ˜ S α iso 0 . Y XZ V Material properties: G ∞ iso = c (cid:16) ˜ I − (cid:17) + c (cid:16) ˜ I − (cid:17) , ρ = 1 . × kg/m , E = 1 . × Pa, c = c = E/ β ∞ = 1 . µ = 10 c , τ = 1 s.Reference scales: L = 1 m, M = 1 kg, T = 1 s. Table 1: Three-dimensional beam bending: problem setting, boundary conditions, initial conditions, and material properties.Notice that the parameter β ∞ is only used in the IPC and MIPC models. In this example, we consider a three-dimensional beam vibration problem with bending dominated de-formation. The problem setting and the material properties are defined in Table 1. In particular, thebackground elastic material is characterized by the Mooney-Rivlin model. The bottom surface of the beamis fully clamped while the rest boundary surfaces are specified with zero traction boundary conditions. Thebody is initially in a stress-free condition with zero displacements. The initial velocity is given by V ( X ,
0) = (cid:18) V ZL , , (cid:19) T , V = 53 m / s , (a) Total energies (b) Kinetic energies (c) Potential energies Figure 3: The total, kinetic, and potential energies (i.e., G ∞ iso ( ˜ C h )+Υ ( ˜ C h , Γ h )) of the IPC (red), HS (blue), and MIPC (black)models over time. The simulations are performed with a fixed time step size ∆ t/T = 1 × − . The solid lines illustrate resultsobtained from a spatial mesh with p = 2, a = 1, b = 0, and 5 × ×
30 elements; the dashed lines illustrate results obtainedfrom a spatial mesh with p = 1, a = 1, b = 0, and 1 × × E is chosen tobe the total energy at time t = 0, which is 1 . × kg m /s .
27S MIPC IPC
Figure 4: The snapshots of the pressure fields at time t/T = 2 . t/T = 1 × − plotted on the deformed configuration. Two sets of spatial meshes are used: a spatial mesh with p = 2, a = 1, b = 0, and5 × ×
30 elements and a spatial mesh with p = 1, a = 1, b = 0, and 1 × × which initiates the vibration. For this specific problem, the energy stability given by Theorem 1 suggeststhat the total energy monotonically decreases with respect to time, i.e., ddt (cid:90) Ω X ρ (cid:107) V h (cid:107) + G ∞ iso ( ˜ C h ) + Υ ( ˜ C h , Γ h ) d Ω X = − (cid:90) Ω X (cid:18) ddt Γ h (cid:19) : V : (cid:18) ddt Γ h (cid:19) d Ω X = − (cid:90) Ω X Q : (cid:0) V (cid:1) − : Q d Ω X ≤ . In the first set of simulations, we investigate the problem with two spatial meshes and integrate in timeuntill 5 T . The generalized- α method is utilized with a fixed time step size and (cid:37) ∞ = 1 .
0, which recoversthe mid-point rule. The evolutions of the total, kinetic, and potential energies are illustrated in Figure 3.We can observe the energy decay of all three material models with different dissipation rates. The HS modelleads to the fastest dissipation of the energy, while the energies of the IPC and MIPC models dissipaterelatively slower. In Figure 4, the snapshots of the pressure field on the deformed configuration at time t/T = 2 . α scheme is illustrated. The generalized- α schemes with ∆ t/T = 1 . × − and (cid:37) ∞ = 0 .
0, 0 .
5, and 1 . t/T = 1 . × − .The numerical dissipation introduced by the time-stepping algorithm is negligible in comparison with thephysical dissipation generated by the model from the energy evolution. In the detailed view shown in thethree bottom figures of Figure 5, we see that, as expected, choosing a smaller value of (cid:37) ∞ induces morenumerical dissipation since this parameter dictates the dissipation on the high-frequency modes according tothe analysis of linear problems. The differences between the two simulations using (cid:37) ∞ = 1 . µ = 0 . c ,with all other settings equal to those in Table 1. The snapshots of the deformation states and the pressure28 IPC HS MIPC
Figure 5: The total energies of the IPC (red), HS (blue), and MIPC (black) models over time. The simulations are performedwith a spatial mesh with p = 2, a = 1, b = 0, and 5 × ×
30 elements. The reference value E is chosen to be the total energyat time t/T = 0, which is 1 . × kg m /s . Detailed view of the energies in the vicinity of t/T = 5 is depicted in thebottom. T = 0 .
25 0 .
50 0 .
75 1 .
00 1 .
50 2 .
00 2 .
50 3 .
00 3 . Figure 6: The snapshots of the pressure fields at nice time instances of the IPC (top) and MIPC (bottom) models with afixed time step size ∆ t/T = 1 × − plotted on the deformed configuration. The material parameters are the same as thosereported in Table 1, except µ = 0 . c here. The simulations are performed using a spatial mesh with p = 1, a = 1, b = 0, and1 × × Figure 7: The total, kinetic, and potential energies of the IPC (red) and MIPC (black) models over time. The simulations areperformed with µ = 0 . c , a fixed time step size ∆ t/T = 1 × − , and a spatial mesh with p = 1, a = 1, b = 0, and 1 × × E is chosen to be the total energy at time t = 0, which is 1 . × kg m /s . t/T = 1 . t/T = 2 .
2, after which the potential and total energies quickly become negative. Indeed,the IPC model cannot guarantee the boundedness of the configurational free energy (see, e.g. (2.49)). Thus,the evolution to unbounded negative energy is indeed possible for this very model, as is discovered in thiscase. This unstable behavior of the IPC model again suggests the original IPC model, which results innon-vanishing non-equilibrium stresses, is materially unstable and may produce non-physical results.
Y XZ V Material properties: G ∞ iso = c (cid:16) ˜ I − (cid:17) + c (cid:16) ˜ I − (cid:17) , ρ = 1 . × kg/m , E = 1 . × Pa, c = c = E/ β ∞ = 0 . µ = c , τ = 0 . β ∞ = 0 . µ = c , τ = 0 . β ∞ = 0 . µ = c , τ = 1 . L = 1 m, M = 1 kg, T = 1 s. Table 2: Three-dimensional beam torsion: problem setting, boundary conditions, initial conditions, and material properties.Notice that the parameters β ∞ , β ∞ , and β ∞ are only used in the IPC and MIPC models. In this second example, we consider the torsion of a three-dimensional beam. The problem setting andmaterial properties are illustrated in Table 2. In particular, the boundary conditions are identical to those ofthe beam bending example, and the body is initially stress-free with zero displacements. The initial velocityis given by V ( X ,
0) = V (cid:18) − YL , XL , (cid:19) T , V = 100 sin (cid:18) πZ L (cid:19) m / s . Again, the problem constitutes an unforced mechanical system with the total energy monotonically decreas-ing with time. We numerically investigate this problem using two sets of spatiotemporal discretizations. Ina coarse discretization, we use a fixed time step size ∆ t/T = 5 . × − and a spatial mesh with p = 1 and2 × ×
12 elements; in a fine discretization, we use a fixed time step size ∆ t/T = 1 . × − and a spatialmesh with p = 2 and 5 × ×
30 elements. The time integration is performed with the mid-point rule (i.e., (cid:37) ∞ = 1 . xy -component of the Cauchy stress are plotted on the current configuration at time t/T = 0 .
1. Forcomparison purposes, the results of the three models calculated by the two discretizations are illustrated.31 (a) Total energies (b) Kinetic energies (c) Potential energies
Figure 8: The total, kinetic, and potential energies (i.e., G ∞ iso ( ˜ C h ) + Υ ( ˜ C h , Γ h )) of the IPC (red), HS (blue), and MIPC(black) models over time for the beam torsion problem. The solid and dashed lines illustrate results obtained from the fineand coarse spatiotemporal discretizations. The reference value of the total energy E is chosen to be the total energy at time t = 0, which is 2 . × kg m /s .
32S MIPC IPC
Figure 9: The snapshots of the pressure fields plotted on the deformed configuration at time t/T = 0 . HS MIPC IPC
Figure 10: The snapshots of the stress component σ xy plotted on the deformed configuration at time t/T = 0 .
33e observe that both the pressure and the Cauchy stress are well resolved by the coarse mesh and are indis-tinguishable. This again demonstrates the superior capability of the chosen element technology in resolvingthe stresses.
In the third example, we consider a clamped cylindrical support characterized by viscoelastic materialbehaviors with its inner surface connected to a vibrating device. This is a benchmark problem designedfor examing the rate-dependent behavior with hysteresis loops [22, 112]. The geometrical setting of thecylindrical support together with its material properties are illustrated in Table 3. The vibration of thedevice is described by U ( X , t ) = ( U sin ( ωt ) , , T , U = 5 . × − m , which represents a translation along the x -direction, and, in this study, the strain rate ω takes the values of 5s − , 10 s − , and 20 s − , respectively. The outer surface of the support is fixed while traction-free boundaryconditions are prescribed on the two end annular surfaces. The geometry of this cylindrical support canbe exactly represented via NURBS. We use 24 elements in the circumferential direction, 6 elements in theradial direction, and 5 elements in the axial direction. We choose p = 2 and 3 for the discrete pressurefunction space, respectively, to generate two sets of meshes. For the given strain rate, the period of onecycle is 2 π/ω , and we simulate the problems for three cycles. For the mesh with p = 2, the time step sizeis ∆ t = 2 . × − s; for the mesh with p = 3, the time step size is ∆ t = 1 . × − . In this problem, wechoose (cid:37) ∞ = 0 . Z YX . . . Material properties: G ∞ iso = c (cid:16) ˜ I − (cid:17) + c (cid:16) ˜ I − (cid:17) , ρ = 1 . × kg/m , E = 1 . × Pa, c = c = E/ β ∞ = 1 . µ = 4 c , τ = 0 . L = 1 m, M = 1 kg, T = 1 s. Table 3: Three-dimensional clamped cylindrical support: problem setting, boundary conditions, initial conditions, and materialproperties. Notice that the parameter β ∞ is only used in the MIPC models. To investigate the viscous dissipative effect, we calculate the force on the inner surface Γ inner as F := (cid:90) Γ inner P N d Γ . The lateral component in the x -direction F x is plotted against the lateral displacement U x for the HS andMIPC models in Figure 11 for both meshes. The hysteretical curves are almost indistinguishable, suggestingthe results are mesh independent. Comparing the results between the HS and MIPC models, we may observethat the difference in the configurational free energy engenders differences in the hysteretical loop, with theMIPC model giving higher values of the maximum force. With the increase of the strain rate ω , bothmaterial models need more cycles to reach equilibrium; the increase in strain rate leads to an increase in34 HS MIPC -6 -4 -2 0 2 4 6-2000-1000010002000
Figure 11: The lateral force F x on Γ inner plotted against the lateral displacement U x for ω = 5s − (top row), ω = 10s − (middle row), ω = 20s − (bottom row). In the left column, the configurational free energy is given by the HS model while onthe right column it is given by the MIPC model.
35S MIPC ω = 5 s − , p = 2 ω = 5 s − , p = 3 ω = 5 s − , p = 2 ω = 5 s − , p = 3 ω = 10 s − , p = 2 ω = 10 s − , p = 3 ω = 10 s − , p = 2 ω = 10 s − , p = 3 ω = 20 s − , p = 2 ω = 20 s − , p = 3 ω = 20 s − , p = 2 ω = 20 s − , p = 3 Figure 12: The snapshots of the pressure plotted on the deformed configuration at time t = π/ ω of the HS and MIPC modelsfor the clamped cylindrical support problem. The results from the meshes with p = 2 and 3 are shown for comparison. p = 2 and 3 show good agreement for both models under the three strain rates. Thisagain demonstrates the superior stress resolving property of the proposed discretization method [58].
5. Conclusion
In this work, we start by considering a general continuum theory for viscoelasticity with the viscousdeformation characterized by a set of internal state variables. With the incompressible constraint as acommon material property in mind, we choose to develop a theory based on the Gibbs free energy [55, 60],which leads to a pressure primitive variable formulation. A set of nonlinear evolution equations is derivedfor the viscous deformation. With that, we consider a special form of the configurational free energy, whichleads to the definition of the finite linear viscoelasticity. This viscoelasticity theory is endowed with a setof linear evolution equations, which is amenable to finite element implementation and anisotropic materialmodeling. It is revealed through the derivation that the original viscoelasticity model developed in [36, 42]can be rectified to maintain thermodynamic consistency. In particular, the following rectifications need tobe pointed out.1. The right-hand side of the evolution equations (2.38) is driven by the fictitious second Piloa-Kirchhoffstress ˜ S α iso . In contrast, in the original model, the right-hand side is driven by J − P : ˜ S α iso [36, 42].2. The non-equilibrium stresses S α neq and the conjugate variables Q α are two different quantities. Theybecome identical only when the configurational free energy takes a special form (2.52).3. The configurational free energy needs to satisfy the normalization condition (2.13); the relaxation ofthe non-equilibrium stresses poses an additional constraint on the form of the configurational freeenergy (see Section 2.4.3). In particular, the identical polymer chain model considered in [46] doesnot guarantee the vanishment of the viscous stress in the equilibrium limit, and we have identified anunstable solution of that model (see Figures 6 and 7).Indeed, the original finite linear viscoelasticity model has been criticized for the lack of a thermomechanicalfoundation, and a finite time blow-up solution has been previously identified. The above three points couldbe the potential source of the blow-up phenomenon observed by S. Govindjee, et al. in [50]. Additionally,in one instantiation of the model, the free energy indicates that there is an additive split of the elasticand viscous strain (see (2.53)), which makes this theory analogous to the elastoplasticity theory proposedby A. Green and P. Naghdi [52, 53]. Based on the consistent continuum theory, we construct a numericalformulation that inherits the stability to the semidiscrete formulation. We invoke the smooth generalizationof the Taylor-Hood element based on NURBS for the spatial discretization. The well-known recurrenceformula is utilized to integrate the constitutive relation, and the dynamic integration is performed by thegeneralized- α scheme. A variety of benchmark examples are presented and corroborate the properties of thecontinuum and numerical formulations in different deformation states. The superior stress accuracy of theadopted NURBS basis function is demonstrated as well.Based on this work, we will extend the proposed viscoelasticity theory to fiber-reinforced materials witha particular focus on arterial wall modeling and vascular fluid-structure interaction. It is also of interest toinvestigate the approximation of the viscoelastic constitutive relation directly from the hereditary integral,which may open the door for the numerical modeling of biological tissue growth and remodeling [19]. Acknowledgements
We want to thank Prof. Jay D. Humphrey at Yale University for many helpful discussions. This work issupported by the National Institutes of Health under the award numbers 1R01HL121754, 1R01HL123689,R01EB01830204, the startup grant provided by the Southern University of Science and Technology under37he award number Y01326127, the Guangdong-Hong Kong-Macao Joint Laboratory for Data-Driven FluidMechanics and Engineering Applications under the award number 2020B1212030001, the computationalresources from the Center for Computational Science and Engineering at Southern University of Science andTechnology, the Stanford Research Computing Center, and the Extreme Science and Engineering DiscoveryEnvironment supported by the National Science Foundation grant ACI-1053575.
Appendix A. An analysis of the null space of P In this section, we show that a fictitious stress living in the null space of P has to be a zero stress. Bydefinition, this stress satisfies O = P : ˜ S = ˜ S − (cid:16) ˜ S : C (cid:17) C − . From the above definition, one has˜ S = 13 (cid:16) ˜ S : C (cid:17) C − = 13 (cid:16) ˜ S : ˜ C (cid:17) ˜ C − = a ˜ C − , with a := 13 (cid:16) ˜ S : ˜ C (cid:17) . (A.1)Recall that the fictitious stress is given by an isochoric energy G iso ( ˜ C , Θ , Γ , · · · , Γ m ),˜ S = 2 ∂G iso ∂ ˜ C . Now, let C ( t ) be an arbitrary Green-Lagrange tensor that smoothly varies over time t . We consider a freeenergy G defined for C ( t ) with a conjugate stress satisfying the relation (A.1). Therefore, one has G ( C ( t )) = G ( C ( t )) − G ( I ) = 12 (cid:90) t S ( C ( s )) : ˙ C ( s ) ds = 12 (cid:90) t a C − ( s ) : ˙ C ( s ) ds = (cid:90) t a C − ( s ) : F T dF ds = (cid:90) t a tr[ d ] ds, with d being the rate of deformation tensor. Now we consider that the Green-Lagrange tensor characterizesan isochoric motion, implying tr[ d ] = 0. Therefore, we have G ( ˜ C ) = 0 for ˜ C characterizing volume-preserving deformations. This suggests that for a fictitious stress satisfying the relation (A.1), it has to bea zero stress, ˜ S = 2 ∂G iso ∂ ˜ C = O . In other words, an isochoric free energy cannot induce a hydrostatic fictitious stress.
Appendix B. Algorithm for the stress and elasticity tensor in the identical polymer chainmodel with G α = F α In the original model based on the identical polymer chain assumption, it is taken that G α ( ˜ C , Θ) = F α ( ˜ C , Θ) = β ∞ α G ∞ iso ( ˜ C , Θ) [46]. Then the configurational free energy takes the following form.Υ α (cid:16) ˜ C , Θ , Γ , · · · , Γ m (cid:17) = H α ( Γ α , Θ) + (cid:32) ˆ S α − ∂G α ( ˜ C , Θ) ∂ ˜ C (cid:33) : Γ α − I G α ( ˜ C , Θ) , = µ α (Θ) (cid:12)(cid:12)(cid:12)(cid:12) Γ α − I (cid:12)(cid:12)(cid:12)(cid:12) + (cid:16) ˆ S α − β ∞ α ˜ S ∞ iso (cid:17) : Γ α − I β ∞ α G ∞ iso ( ˜ C , Θ) . To the best of our knowledge, numerical analysis for that model has not been performed, probably due tothe appearance of a six-order tensor in the definition of the elasticity tensor. Here, we provide the algorithmfor computing the stress as well as the elasticity tensor of this model. Given the time step ∆ t n , and the valueof ˜ S ∞ iso n and Q αn at each quadrature points, proceed through the following steps to compute the isochoricpart of the second Piola-Kirchhoff stress and the isochoric part of the elasticity tensor.38. Calculate the deformation gradient and the strain measures based on the displacement U n +1 , F n +1 = I + ∇ X U n +1 , J n +1 = det ( F n +1 ) , C n +1 = F Tn +1 F n +1 , ˜ C n +1 = J − / n +1 C n +1 .
2. Calculate the fictitious second Piola-Kirchhoff stress˜ S ∞ iso n +1 = 2 (cid:18) ∂G ∞ iso ∂ ˜ C (cid:19) n +1 .
3. Calculate the fictitious elasticity tensor˜ C ∞ iso n +1 = 4 J − n +1 (cid:18) ∂ G ∞ iso ∂ ˜ C ∂ ˜ C (cid:19) n +1 .
4. For α = 1 , · · · , m , calculate Q αn +1 = β ∞ α exp ( ξ α ) ˜ S ∞ iso n +1 + exp ( ξ α ) (cid:16) exp ( ξ α ) Q αn − β ∞ α ˜ S ∞ iso n (cid:17) .
5. For α = 1 , · · · , m , calculate˜ S α neq n +1 = β ∞ α ˜ S ∞ iso n +1 − β ∞ α µ α J n +1 ˜ C ∞ iso n +1 : (cid:16) β ∞ α ˜ S ∞ iso n +1 − ˆ S α − Q αn +1 (cid:17) .
6. Calculate the fictitious stress ˜ S n +1 = ˜ S ∞ iso n +1 + m (cid:88) α =1 ˜ S α neq n +1 .
7. Calculate the projection tensor P n +1 = I − C − n +1 ⊗ C n +1 .
8. Calculate the isochoric part of the second Piola-Kirchhoff stress S iso n +1 = J − n +1 P n +1 : ˜ S n +1 .
9. Calculate the fictitious elasticity tensors for the non-equilibrium part,˜ C α neq n +1 = β ∞ α ˜ C ∞ iso n +1 − β ∞ α µ α J − n +1 (cid:18) ∂ G ∞ iso ∂ ˜ C ∂ ˜ C ∂ ˜ C (cid:19) n +1 : (cid:16) β ∞ α ˜ S ∞ iso n +1 − ˆ S α − Q αn +1 (cid:17) − β ∞ α µ α J n +1 (1 − δ α ) ˜ C ∞ iso n +1 : ˜ C ∞ iso n +1 .
10. Calculate the fictitious elasticity tensor˜ C n +1 = ˜ C ∞ iso n +1 + m (cid:88) α =1 ˜ C α neq n +1 .
11. Calculate the fourth-order tensor˜ P n +1 = C − n +1 (cid:12) C − n +1 − C − n +1 ⊗ C − n +1 .
12. Calculate the isochoric part of the elasticity tensor C iso n +1 = P n +1 : ˜ C n +1 : P Tn +1 + 23 Tr (cid:16) J − n +1 ˜ S n +1 (cid:17) ˜ P n +1 − (cid:0) C − n +1 ⊗ S iso n +1 + S iso n +1 ⊗ C − n +1 (cid:1) . ppendix C. Algorithm for the stress and elasticity tensor in the Holzapfel-Simo model The Holzapfel-Simo model presented in Section 2.4.4 is different from the classical one documented inclassical works [44, Chapter 6.10]. Here, we provide the algorithm for computing the stresses as well asthe elasticity tensor. Given the time step ∆ t n , and the value of ˜ S ∞ iso n and Q αn at each quadrature points,proceed through the following steps to compute the isochoric part of the second Piola-Kirchhoff stress andthe isochoric part of the elasticity tensor.1. Calculate the deformation gradient and the strain measures based on the displacement U n +1 , F n +1 = I + ∇ X U n +1 , J n +1 = det ( F n +1 ) , C n +1 = F Tn +1 F n +1 , ˜ C n +1 = J − / n +1 C n +1 .
2. Calculate ˜ S ∞ iso n +1 = 2 (cid:18) ∂G ∞ iso ∂ ˜ C (cid:19) n +1 .
3. For α = 1 , · · · , m , calculate Q αn +1 = µ α exp ( ξ α ) ˜ C n +1 + exp ( ξ α ) (cid:16) exp ( ξ α ) Q αn − µ α ˜ C n (cid:17) .
4. Calculate the fictitious stress ˜ S n +1 = ˜ S ∞ iso n +1 + m (cid:88) α =1 Q αn +1 .
5. Calculate the projection tensor P n +1 = I − C − n +1 ⊗ C n +1 .
6. Calculate the isochoric part of the second Piola-Kirchhoff stress S iso n +1 = J − n +1 P n +1 : ˜ S n +1 .
7. Calculate the fictitious elasticity tensors for the equilibrium and non-equilibrium parts,˜ C ∞ iso n +1 = 4 J − n +1 (cid:18) ∂ G ∞ iso ∂ ˜ C ∂ ˜ C (cid:19) n +1 , and ˜ C α neq n +1 = 2 µ α exp( ξ α ) J − n +1 I , for α = 1 , · · · , m.
8. Calculate the fictitious elasticity tensor˜ C n +1 = ˜ C ∞ iso n +1 + m (cid:88) α =1 ˜ C α neq n +1 .
9. Calculate the fourth-order tensor˜ P n +1 = C − n +1 (cid:12) C − n +1 − C − n +1 ⊗ C − n +1 .
10. Calculate the isochoric part of the elasticity tensor C iso n +1 = P n +1 : ˜ C n +1 : P Tn +1 + 23 Tr (cid:16) J − n +1 ˜ S n +1 (cid:17) ˜ P n +1 − (cid:0) C − n +1 ⊗ S iso n +1 + S iso n +1 ⊗ C − n +1 (cid:1) . ppendix D. Algorithm for the stress and elasticity tensor in the modified identical polymerchain model The modified identical polymer chain model presented in Section 2.4.5 rectifies the non-physical behaviorof the original identical polymer chain model in the thermodynamic equilibrium. Here we provide thealgorithm for computing the stresses as well as the elasticity tensor of this model. Given the time step ∆ t n ,and the value of ˜ S ∞ iso n and Q αn at each quadrature points, proceed through the following steps to computethe isochoric part of the second Piola-Kirchhoff stress and the isochoric part of the elasticity tensor.1. Calculate the deformation gradient and the strain measures based on the displacement U n +1 , F n +1 = I + ∇ X U n +1 , J n +1 = det ( F n +1 ) , C n +1 = F Tn +1 F n +1 , ˜ C n +1 = J − / n +1 C n +1 .
2. Calculate the fictitious second Piola-Kirchhoff stress˜ S ∞ iso n +1 = 2 (cid:18) ∂G ∞ iso ∂ ˜ C (cid:19) n +1 .
3. Calculate the fictitious elasticity tensor˜ C ∞ iso n +1 = 4 J − n +1 (cid:18) ∂ G ∞ iso ∂ ˜ C ∂ ˜ C (cid:19) n +1 .
4. For α = 1 , · · · , m , calculate Q αn +1 = β ∞ α exp ( ξ α ) ˜ S ∞ iso n +1 + exp ( ξ α ) (cid:16) exp ( ξ α ) Q αn − β ∞ α ˜ S ∞ iso n (cid:17) .
5. For α = 1 , · · · , m , calculate ˜ S α neq n +1 = J n +1 β ∞ α µ α ˜ C ∞ iso n +1 : Q αn +1 .
6. Calculate the fictitious stress ˜ S n +1 = ˜ S ∞ iso n +1 + m (cid:88) α =1 ˜ S α neq n +1 .
7. Calculate the projection tensor P n +1 = I − C − n +1 ⊗ C n +1 .
8. Calculate the isochoric part of the second Piola-Kirchhoff stress S iso n +1 = J − n +1 P n +1 : ˜ S n +1 .
9. Calculate the fictitious elasticity tensors for the non-equilibrium part,˜ C α neq n +1 = 4 β ∞ α J − n +1 µ α (cid:18) ∂ G ∞ iso ∂ ˜ C ∂ ˜ C ∂ ˜ C (cid:19) n +1 : Q αn +1 + δ α β ∞ α J n +1 µ α ˜ C ∞ iso n +1 : ˜ C ∞ iso n +1 .
10. Calculate the fictitious elasticity tensor˜ C n +1 = ˜ C ∞ iso n +1 + m (cid:88) α =1 ˜ C α neq n +1 .
11. Calculate the fourth-order tensor˜ P n +1 = C − n +1 (cid:12) C − n +1 − C − n +1 ⊗ C − n +1 .
12. Calculate the isochoric part of the elasticity tensor C iso n +1 = P n +1 : ˜ C n +1 : P Tn +1 + 23 Tr (cid:16) J − n +1 ˜ S n +1 (cid:17) ˜ P n +1 − (cid:0) C − n +1 ⊗ S iso n +1 + S iso n +1 ⊗ C − n +1 (cid:1) . eferences [1] J. Ferry, Viscoelastic properties of polymers, John Wiley & Sons, 1980.[2] J. Humphrey, Cardiovascular solid mechanics: cells, tissues, and organs, Springer Science & Business Media, 2013.[3] M. Shaw, W. MacKnight, Introduction to polymer viscoelasticity, John Wiley & Sons, 2018.[4] J. Ben´ıtez, F. Mont´ans, The mechanical behavior of skin: Structures and models for the finite element analysis, Computer& Structures 190 (2017) 75–107.[5] G. Holzapfel, T. Gasser, M. Stadler, A structural model for the viscoelastic behavior of arterial walls: Continuumformulation and finite element analysis, European Journal of Mechanics-A/Solids 21 (3) (2002) 441–463.[6] B. Coleman, M. Gurtin, Thermodynamics with internal state variables, Journal of Chemical Physics 47 (1967) 597–613.[7] M. Horstemeyer, D. Bammann, Historical review of internal state variable theory for inelasticity, International Journalof Plasticity 26 (2010) 1310–1334.[8] G. Maugin, The saga of internal variables of state in continuum thermo-mechanics (1893–2013), Mechanics ResearchCommunications 69 (2015) 79–86.[9] B. Coleman, W. Noll, Foundations of linear viscoelasticity, Reviews of modern physics 33 (1961) 239.[10] M. Gurtin, E. Sternberg, On the linear theory of viscoelasticity, Archive for Rational Mechanics and Analysis 11 (1962)291–356.[11] R. Christensen, A nonlinear theory of viscoelasticity for application to elastomers, Journal of Applied Mechanics 47(1980) 763.[12] C. Drapaca, S. Sivaloganathan, G. Tenti, Nonlinear constitutive laws in viscoelasticity, Mathematics and mechanics ofsolids 12 (2007) 475–501.[13] B. Coleman, Thermodynamics of materials with memory, Archive for Rational Mechanics and Analysis 17 (1964) 1–46.[14] M. Puso, J. Weiss, Finite element implementation of anisotropic quasilinear viscoelasticity, ASME Journal of Biome-chanical Engineering 120 (1998) 162–170.[15] Y. Fung, Biomechanics: mechanical properties of living tissues, Springer Science & Business Media, 2013.[16] Y. Yu, P. Perdikaris, G. Karniadakis, Fractional modeling of viscoelasticity in 3D cerebral arteries and aneurysms, Journalof Computational Physics 323 (2016) 219–242.[17] P. Perdikaris, G. Karniadakis, Fractional-order viscoelasticity in one-dimensional blood flow modelss, Annals of Biomed-ical Engineering 42 (2014) 1012–1023.[18] D. Craiem, F. Rojo, J. Atienza, R. Armentano, G. Guinea, Fractional-order viscoelasticity applied to describe uniaxialstress relaxation of human arteries, Physics in Medicine & Biology 53 (2008) 4543.[19] J. Humphrey, K. Rajagopal, A constrained mixture model for growth and remodeling of soft tissues, MathematicalModels and Methods in Applied Sciences 12 (2002) 407–430.[20] A. Valent´ın, J. Humphrey, G. Holzapfel, A finite element-based constrained mixture implementation for arterial growth,remodeling, and adaptation: Theory and numerical verification, International Journal for Numerical Methods in Biomed-ical Engineering 29 (2013) 822–849.[21] S. Park, Y. Kim, Fitting Prony-series viscoelastic models with power-law presmoothing, Journal of materials in civilengineering 13 (2001) 26–32.[22] X. Zeng, G. Scovazzi, N. Abboud, O. Colom´es, S. Rossi, A dynamic variational multiscale method for viscoelasticityusing linear tetrahedral elements, International Journal for Numerical Methods in Engineering 112 (2017) 1951–2003.[23] S. Fran¸cois, Un mod`ele visco´elastique non lin´eaire avec configuration interm´ediaire, J. M´ecanique 13 (1974) 679–713.[24] P. L. Tallec, C. rahler, Numerical models of steady rolling for non-linear viscoelastic structures in finite deformations,International Journal for Numerical Methods in Engineering 37 (1994) 1159–1186.[25] S. Reese, S. Govindjee, A theory of finite viscoelasticity and numerical aspects, International Journal of Solids andStructures 35 (1998) 3455–3482.[26] S. Reese, S. Govindjee, Theoretical and numerical aspects in the thermo-viscoelastic material behaviour of rubber-likepolymers, Mechanics of Time-Dependent Materials 1 (1998) 357–396.[27] D. Peri´c, W. Dettmer, A computational model for generalized inelastic materials at finite strains combining elastic,viscoelastic and plastic material behaviour, Engineering Computations (2003).[28] B. Nedjar, An anisotropic viscoelastic fibre–matrix model at finite strains: continuum formulation and computationalaspects, Computer Methods in Applied Mechanics and Engineering 196 (2007) 1745–1756.[29] H. Liu, G. Holzapfel, B. Skallerud, V. Prot, Anisotropic finite strain viscoelasticity: Constitutive modeling and finiteelement implementation, Journal of the Mechanics and Physics of Solids 124 (2019) 172–188.[30] W. Hong, X. Zhao, J. Zhou, Z. Suo, A theory of coupled diffusion and large deformation in polymeric gels, Journal ofMechanics and Physics of Solids 56 (2008) 1779–1793.[31] R. Ogden, Large deformation isotropic elasticity-on the correlation of theory and experiment for incompressible rubberlikesolids, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 326 (1972) 565–584.[32] J. Lubliner, A model of rubber viscoelasticity, Mechanics Research Communications 12 (1985) 93–99.[33] J. Simo, C. Miehe, Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implemen-tation, Computer Methods in Applied Mechanics and Engineering 98 (1992) 41–104.[34] M. Latorre, F. Mont´ans, Anisotropic finite strain viscoelasticity based on the sidoroff multiplicative decomposition andlogarithmic strains, Computational Mechanics 56 (2015) 503–531.[35] M. Latorre, F. Mont´ans, Fully anisotropic finite strain viscoelasticity based on a reverse multiplicative decompositionand logarithmic strains, Computers & Structures 163 (2016) 56–70.
36] J. Simo, On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects,Computer Methods in Applied Mechanics and Engineering 60 (1987) 153–173.[37] G. Holzapfel, On large strain viscoelasticity: continuum formulation and finite element applications to elastomeric struc-tures, International Journal for Numerical Methods in Engineering 39 (1996) 3903–3926.[38] G. Holzapfel, T. Gasser, A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, compu-tational aspects and applications, Computer Methods in Applied Mechanics and Engineering 190 (2001) 4379–4403.[39] T. Gasser, C. Forsell, The numerical implementation of invariant-based viscoelastic formulations at finite strains. Ananisotropic model for the passive myocardium, Computer Methods in Applied Mechanics and Engineering 200 (2011)3637–3645.[40] S. Govindjee, J. Simo, Mullins’ effect and the strain amplitude dependence of the storage modulus, International journalof solids and structures 29 (1992) 1737–1751.[41] O. G¨ultekin, G. Sommer, G. Holzapfel, An orthotropic viscoelastic model for the passive myocardium: continuum basisand numerical treatment, Computer methods in biomechanics and biomedical engineering 19 (2016) 1647–1664.[42] J. Simo, T. Hughes, Computational Inelasticity, Springer Science & Business Media, 2006.[43] R. Taylor, K. Pister, G. Goudreau, Thermomechanical analysis of viscoelastic solids, International Journal for NumericalMethods in Engineering 2 (1970) 45–59.[44] G. Holzapfel, Nonlinear Solid Mechanics: A Continuum Approach for Engineering, John Wiley & Sons, 2000.[45] M. Kaliske, H. Rothert, Formulation and implementation of three-dimensional viscoelasticity at small and finite strains,Computational Mechanics 19 (1997) 228–239.[46] G. Holzapfel, J. Simo, A new viscoelastic constitutive model for continuous media at finite thermomechanical changes,International Journal of Solids and Structures 33 (1996) 3019–3034.[47] E. Pe˜na, J. Pe˜na, M. Doblar´e, On modelling nonlinear viscoelastic effects in ligaments, Journal of Biomechanics 41 (2008)2659–2666.[48] J. Pe˜na, M. Mart´ınez, E. Pe˜na, A formulation to model the nonlinear viscoelastic properties of the vascular tissue, ActaMechanica 217 (2011) 63–74.[49] P. Haupt, Continuum mechanics and thoery of materials, Springer Science & Business Media, 2013.[50] S. Govindjee, T. Potter, J. Wilkening, Dynamic stability of spinning viscoelastic cylinders at finite deformation, Inter-national journal of solids and structures 51 (2014) 3589–3603.[51] B. Coleman, W. Noll, The thermodynamics of elastic materials with heat conduction and viscosity, Archive for RationalMechanics and Analysis 13 (1963) 167–178.[52] A. Green, P. Naghdi, A general theory of an elastic-plastic continuum, Archive for Rational Mechanics and Analysis 18(1965) 251–281.[53] A. Green, P. Naghdi, Some remarks on elastic-plastic deformation at finite strain, International Journal of EngineeringScience 9 (1971) 1219–1229.[54] F. Auricchio, L. B. da Veiga, C. Lovadina, A. Reali, R. Taylor, P. Wriggers, Approximation of incompressible largedeformation elastic problems: some unresolved issues, Computational Mechanics 52 (2013) 1153–1167.[55] J. Liu, A. Marsden, A unified continuum and variational multiscale formulation for fluids, solids, and fluid-structureinteraction, Computer Methods in Applied Mechanics and Engineering 337 (2018) 549–597.[56] G. Hauke, T. Hughes, A unified approach to compressible and incompressible flows, Computer Methods in AppliedMechanics and Engineering 113 (1994) 389–395.[57] G. Hauke, T. Hughes, A comparative study of different sets of variables for solving compressible and incompressible flows,Computer Methods in Applied Mechanics and Engineering 153 (1998) 1–44.[58] J. Liu, A. Marsden, A robust and efficient iterative method for hyper-elastodynamics with nested block preconditioning,Journal of Computational Physics 383 (2019) 72–93.[59] J. Liu, W. Yang, M. Dong, A. Marsden, The nested block preconditioning technique for the incompressible Navier-Stokesequations with emphasis on hemodynamic simulations, Computer Methods in Applied Mechanics and Engineering 367(2020) 113122.[60] J. Liu, A. Marsden, Z. Tao, An energy-stable mixed formulation for isogeometric analysis of incompressible hyperelasto-dynamics, International Journal for Numerical Methods in Engineering 120 (2019) 937–963.[61] J. Liu, W. Yang, I. Lan, A. Marsden, Fluid-structure interaction modeling of blood flow in the pulmonary arteries usingthe unified continuum and variational multiscale formulation, Mechanics Research Communications 107 (2020) 103556.[62] A. Buffa, C. de Falco, G. Sangalli, Isogeometric analysis: Stable elements for the 2d stokes equation, International Journalfor Numerical Methods in Fluids (2011) 1407–1422.[63] B. Hosseini, M. M¨oller, S. Turek, Isogeometric analysis of the Navier-Stokes equations with Taylor-Hood B-spline ele-ments, Applied Mathematics and Computation (2015) 264–281.[64] J. Liu, I. Lan, O. Tikenogullari, A. Marsden, A note on the accuracy of the generalized- α scheme for the incompressiblenavier-stokes equations, International Journal for Numerical Methods in Engineering (2021) 638–651.[65] T. Hughes, J. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and meshrefinement, Computer Methods in Applied Mechanics and Engineering 194 (2005) 4135–4195.[66] J. Cottrell, T. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons,2009.[67] J. Evans, Y. Bazilevs, I. Babuska, T. Hughes, n -widths, sup-infs and optimality ratios for the k -version of the isogeometricfinite element method, Computer Methods in Applied Mechanics and Engineering 198 (2009) 1726–1741.[68] S. Lipton, J. Evans, Y. Bazilevs, T. Elguedj, T. Hughes, Robustness of isogeometric structural discretizations undersevere mesh distortion, Computer Methods in Applied Mechanics and Engineering 199 (2010) 357–373.
69] T. Elguedj, Y. Bazilevs, V. Calo, T. Hughes, ¯ B and ¯ F projection methods for nearly incompressible linear and non-linearelasticity and plasticity using higher-order NURBS elements, Computer Methods in Applied Mechanics and Engineering197 (2008) 2732–2762.[70] J. Evans, T. Hughes, Isogeometric divergence-conforming B-splines for the Darcy-Stokes-Brinkman equations, Mathe-matical Models and Methods in Applied Sciences 23 (2013) 671–741.[71] B. Cockburn, F. Li, C. Shu, Locally divergence-free discontinuous galerkin methods for the maxwell equations, Journalof Computational Physics (2004) 588–610.[72] P. Flory, Thermodynamic relations for high elastic materials, Transactions of the Faraday Society 57 (1961) 829–838.[73] J. Simo, R. Taylor, K. Pister, Variational and projection methods for the volume constraint in finite deformation elasto-plasticity, Computer Methods in Applied Mechanics and Engineering 51 (1985) 177–208.[74] J. Marsden, T. Hughes, Mathematical foundations of elasticity, Dover Publications Inc. New York, 1993.[75] G. Scovazzi, T. Hughes, Lecture notes on continuum mechanics on arbitrary moving domains, Tech. Rep. SAND-2007-6312P, Sandia National Laboratories (2007).[76] D. Schroeder, An introduction to thermal physics, Addison Wesley Longman, San Francisco, 1999.[77] C. Truesdell, W. Noll, The Non-Linear Field Theories of Mechanics, Springer, 1965.[78] J. Criscione, J. Humphrey, A. Douglas, W. Hunter, An invariant basis for natural strain which yields orthogonal stressresponse terms in isotropic hyperelasticity, Journal of the Mechanics and Physics of Solids 48 (2000) 2445–2465.[79] C. Sansour, On the physical assumptions underlying the volumetric-isochoric split and the case of anisotropy, EuropeanJournal of Mechanics-A/Solids 27 (2008) 28–39.[80] J. Liu, C. Landis, H. Gomez, T. Hughes, Liquid-Vapor Phase Transition: Thermomechanical Theory, Entropy StableNumerical Formulation, and Boiling Simulations, Computer Methods in Applied Mechanics and Engineering 297 (2015)476–553.[81] J. Sengers, How Fluids Unmix: Discoveries by the School of Van der Waals and Kamerlingh Onnes, Edita-the PublishingHouse of the Royal, 2003.[82] K. Rajagopal, A. Srinivasa, An implicit thermomechanical theory based on a Gibbs potential formulation for describingthe response of thermoviscoelastic solids, International Journal of Engineering Science 70 (2013) 15–28.[83] K. Surana, Y. Mendoza, J. Reddy, Constitutive theories for thermoelastic solids in lagrangian description using gibbspotential, Acta Mechanica 224 (2013) 1019–1044.[84] J. Ball, Convexity conditions and existence theorems in nonlinear elasticity, Archive for rational mechanics and Analysi63 (1976) 337–403.[85] J. Helfenstein, M. Jabareen, E. Mazza, S. Govindjee, On non-physical response in models for fiber-reinforced hyperelasticmaterials, International Journal of Solids and Structures 47 (2010) 2056–2061.[86] O. G¨ultekin, H. Dal, G. Holzapfel, On the quasi-incompressible finite element analysis of anisotropic hyperelastic mate-rials, Computational Mechanics 63 (2019) 443–453.[87] Y. Fung, A First Course in Continuum Mechanics: For Physical and Biological Engineers and Scientists, Prentice Hall,Englewood Cliffs, 1994.[88] M. Gurtin, An introduction to continuum mechanics, Academic press, 1982.[89] A. Tobolsky, I. Prettyman, J. Dillon, Stress relaxation of natural and synthetic rubber stocks, Journal of Applied Physics15 (1944) 380–395.[90] O. G¨ultekin, G. Sommer, G. Holzapfel, An orthotropic viscoelastic model for the passive myocardium: continuum basisand numerical treatment, Computer methods in biomechanics and biomedical engineering 19 (2016) 1647–1664.[91] K. Valanis, Irreversible Thermodynamics of Continuous Media, Internal Variable Theory, Springer-Verlag, Wien, 1972.[92] E. Lee, Elastic-plastic deformation at finite strains, Journal of Applied Mechanics 36 (1969) 1–6.[93] X. Meng, T. Laursen, Energy consistent algorithms for dynamic finite deformation plasticity, Computer Methods inApplied Mechanics and Engineering 191 (2002) 1639–1675.[94] B. Eidel, C. Kuhn, Order reduction in computational inelasticity: Why it happens and how to overcome it-The ODE-caseof viscoelasticity, International Journal for Numerical Methods in Engineering 87 (2011) 1046–1073.[95] J. Bonet, Large strain viscoelastic constitutive models, International Journal of Solids and Structures 38 (2001) 2953–2968.[96] H. Liu, G. Holzapfel, B. Skallerud, V. Prot, Anisotropic finite strain viscoelasticity: Constitutive modeling and finiteelement implementation, Journal of the Mechanics and Physics of Solids 124 (2019) 172–188.[97] C. Miehe, J. Keck, Superimposed finite elastic-viscoelastic-plastoelastic stress response with damage in filled rubberypolymers. Experiments, modelling and algorithmic implementation, Journal of the Mechanics and Physics of Solids 48(2000) 323–365.[98] C. Miehe, S. G¨oktepe, A micro-macro approach to rubber-like materials. Part II: The micro-sphere model of finite rubberviscoelasticity, Journal of the Mechanics and Physics of Solids 53 (2005) 2231–2258.[99] J. Simo, N. Tarnow, K. Wong, Exact energy-momentum conserving algorithms and symmetric schemes for nonlineardynamics, Computer Methods in Applied Mechanics and Engineering 100 (1992) 63–116.[100] I. Romero, An analysis of stress formula for energy-momentum methods in nonlinear elastodynamics, ComputationalMechanics 50 (2012) 603–610.[101] M. Kr¨uger, M. Groß, P. Betsch, An energy-entropy-consistent time stepping scheme for nonlinear thermo-viscoelasticcontinua, ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift f¨ur Angewandte Mathematik und Mechanik96 (2016) 141–178.[102] E. Hinton, J. Campbell, Local and global smoothing of discontinuous finite element functions using a least squaresmethod, International Journal for Numerical Methods in Engineering 8 (1974) 461–480. α method, Journal of applied mechanics 60 (1993) 371–375.[105] K. Jansen, C. Whiting, G. Hulbert, A generalized- α method for integrating the filtered Navier-Stokes equations with astabilized finite element method, Computer Methods in Applied Mechanics and Engineering 190 (2000) 305–319.[106] C. Kadapa, W. Dettmer, D. Peri´c, On the advantages of using the first-order generalised-alpha scheme for structuraldynamic problems, Computers & Structures 193 (2017) 226–238.[107] H. Hilber, T. Hughes, Collocation, dissipation and ‘overshoot’ for time integration schemes in structural dynamics,Earthquake Engineering & Structural Dynamics 6 (1978) 99–117.[108] G. Scovazzi, B. Carnes, X. Zeng, S. Rossi, A simple, stable, and accurate linear tetrahedral finite element for transient,nearly, and fully incompressible solid dynamics: a dynamic variational multiscale approach, International Journal forNumerical Methods in Engineering 106 (2016) 799–839.[109] R. Ortigosa, M. Franke, A. Janz, A. Gil, P. Betsch, An energy-momentum time integration scheme based on a convexmulti-variable framework for non-linear electro-elastodynamics, Computer Methods in Applied Mechanics and Engineer-ing 339 (2018) 1–35.[110] P. Betsch, A. Janz, An energy-momentum consistent method for transient simulations with mixed finite elements devel-oped in the framework of geometrically exact shells, International Journal for Numerical Methods in Engineering 108(2016) 423–455.[111] D. Kuhl, M. Crisfield, Energy-conserving and decaying algorithms in non-linear structural dynamics, International Journalfor Numerical Methods in Engineering 45 (1999) 569–599.[112] E. Fancello, J. Ponthot, L. Stainier, A variational formulation of constitutive models and updates in non-linear finiteviscoelasticity, International Journal for Numerical Methods in Engineering 65 (11) (2006) 1831–1864.method for integrating the filtered Navier-Stokes equations with astabilized finite element method, Computer Methods in Applied Mechanics and Engineering 190 (2000) 305–319.[106] C. Kadapa, W. Dettmer, D. Peri´c, On the advantages of using the first-order generalised-alpha scheme for structuraldynamic problems, Computers & Structures 193 (2017) 226–238.[107] H. Hilber, T. Hughes, Collocation, dissipation and ‘overshoot’ for time integration schemes in structural dynamics,Earthquake Engineering & Structural Dynamics 6 (1978) 99–117.[108] G. Scovazzi, B. Carnes, X. Zeng, S. Rossi, A simple, stable, and accurate linear tetrahedral finite element for transient,nearly, and fully incompressible solid dynamics: a dynamic variational multiscale approach, International Journal forNumerical Methods in Engineering 106 (2016) 799–839.[109] R. Ortigosa, M. Franke, A. Janz, A. Gil, P. Betsch, An energy-momentum time integration scheme based on a convexmulti-variable framework for non-linear electro-elastodynamics, Computer Methods in Applied Mechanics and Engineer-ing 339 (2018) 1–35.[110] P. Betsch, A. Janz, An energy-momentum consistent method for transient simulations with mixed finite elements devel-oped in the framework of geometrically exact shells, International Journal for Numerical Methods in Engineering 108(2016) 423–455.[111] D. Kuhl, M. Crisfield, Energy-conserving and decaying algorithms in non-linear structural dynamics, International Journalfor Numerical Methods in Engineering 45 (1999) 569–599.[112] E. Fancello, J. Ponthot, L. Stainier, A variational formulation of constitutive models and updates in non-linear finiteviscoelasticity, International Journal for Numerical Methods in Engineering 65 (11) (2006) 1831–1864.