Multiplicative, Non-Newtonian Viscoelasticity Models for Rubber Materials and Brain Tissues: Numerical Treatment and Comparative Studies

In many aspects, elastomers and soft biological tissues exhibit similar mechanical properties such as a pronounced nonlinear stress–strain relation and a viscoelastic response to external loads. Consequently, many models use the same rheological framework and material functions to capture their behavior. The viscosity function is thereby often assumed to be constant and the corresponding free energy function follows that one of the long-term equilibrium response. This work questions this assumption and presents a detailed study on non-Newtonian viscosity functions for elastomers and brain tissues. The viscosity functions are paired with several commonly used free energy functions and fitted to two different types of elastomers and brain tissues in cyclic and relaxation experiments, respectively. Having identified suitable viscosity and free energy functions for the different materials, numerical aspects of viscoelasticity are addressed. From the multiplicative decomposition of the deformation gradient and ensuring a non-negative dissipation rate, four equivalent viscoelasticity formulations are derived that employ different internal variables. Using an implicit exponential map as time integration scheme, the numerical behavior of these four formulations are compared among each other and numerically robust candidates are identified. The fitting results demonstrate that non-Newtonian viscosity functions significantly enhance the fitting quality. It is shown that the choice of a viscosity function is even more important than the choice of a free energy function and the classical neo-Hooke approach is often a sufficient choice. Furthermore, the numerical investigations suggest the superiority of two of the four viscoelasticity formulations, especially when complex finite element simulations are to be conducted.


Introduction
An essential property of rubber materials is their rate and time dependent stress response which can be observed in relaxation, creep as well as cyclic experiments. This characteristic is highly non-linear and even more pronounced for filled rubbers. In material science and computational mechanics, such a behavior is commonly referred to as viscoelasticity and its modeling is still object of current research, see for instance [1][2][3][4]. Due to these hysteretic properties, the scope of engineering applications typically comprises damping, isolating and absorbing components especially for vehicles, aseismic structures or low-noise machineries.
Soft biological tissues exhibit qualitatively very similar properties to rubber materials. For example, soft tissues are nearly incompressible, can undergo large deformations and show viscoelastic behavior. Moreover, they exhibit preconditioning effects which lead to stress softening and permanent set, commonly called Mullins effect for elastomers. There are also crucial differences, e.g., many tissues exhibit different behavior in tension and compression what is called tension-compression asymmetry. Furthermore, tendons or blood vessels are reinforced by fibres in preferred directions and, hence, show a strong anisotropy. In contrast, some tissues, mostly non-load-bearing like brain and fat tissue, were shown to behave nearly isotropic, cf. Budday et al. [5].
The present paper compares the capability and properties of existing viscoelasticity models based on experimental data of rubbers and soft tissues. The similarities in the mechanical response of these materials suggest to reproduce their behavior using the same material models [6][7][8].
Specifically brain tissues are considered here due to the relevance of their viscoelastic behavior in simulations of head injuries like concussion and the existence of suitable experimental data sets for several types of brain tissues, e.g., Budday et al. [9].
For the comparison, a modeling framework is employed that aims to capture specifically viscoelastic material properties at a minimum level of model complexity. In terms of rheological models, it can be represented by a Maxwell element and a spring connected in parallel, see Fig. 1, which is often referred to as standard solid or Zener model. The Maxwell element itself is a series connection of a spring and a dashpot which yields the viscoelastic overstress, whereas the spring in parallel generates the long-term stress response. Commonly, two approaches exist to model the Maxwell element. The first is a convolution approach, i.e., the stress differential equation from the small strain theory is adopted. The second approach is a multiplicative decomposition of the deformation gradient into an elastic and an inelastic part, representing the deformation of the spring and the dashpot, respectively. As mentioned by Govindjee & Reese [10] and to the best of the authors' knowledge, the convolution approach has not been shown to satisfy the second law of thermodynamics in a general way. Therefore, the decomposition of the deformation gradient is employed here. 1 The standard solid in Fig. 1 requires four response functions to yield an executable model. The long-term stress response is comprised of a deviatoric and a hydrostatic contribution which are derived from the decoupled free energy functions Ψ eq and Ψ vol , respectively. On the other hand, the Maxwell element is defined via the free energy function Ψ neq as well as the viscosity function . Since the hydrostatic stress of rubbers exhibits only negligible viscoelastic properties [14] and the same is often assumed for brain tissue [8], Ψ neq contributes to the deviatoric stress only. Thus, the total deviatoric stress consists of an equilibrium (index eq ) and a non-equilibrium (index neq ) part. If the chosen free energy function Ψ neq of the Maxwell element leads to a constant shear modulus (i.e., neo-Hooke or Mooney-Rivlin model), then its behavior is referred to as Hookean viscoelasticity in the present paper. Moreover, a constant viscosity produces a so-called Newtonian viscoelasticity. 2 Since such simple approaches are insufficient to describe the complex behavior of rubbers or brain tissues, many authors proposed more sophisticated material functions.
In this work, such advanced approaches are compiled and the influence of the viscosity function and the non-equilibrium free energy function on the fitting quality is investigated. For a wide applicability, two distinct rubber materials (carbon black filled ethylene propylene diene rubber (EPDM) & natural rubber (NR)) as well as brain tissues (cortex (C) & corpus callosum (CC)) under relaxation as well as cyclic loading are considered. In addition, since the derivation of the viscoelastic evolution equation yields four equivalent formulations with different internal variables, also numerical tests are conducted. Favorable formulations in terms of efficiency and robustness are identified.
The manuscript is organized as follows. Firstly, in Sect. 2, the theory of multiplicative viscoelasticity is explained in detail and four equivalent formulations of the evolution equation are presented. Moreover, the viscosity functions described in the literature are summarized. Thereafter, in Sects. 3 and 4, the experimental data and the fitting procedure used herein are presented, before the applicability of multiplicative viscoelasticity to rubber materials and brain tissues is studied in Sects. 5 and 6. Finally, in Sect. 7, the numerical behavior of different formulations of viscoelastic evolution equations and update schemes is examined. The work closes with a summary and outlook in Sect. 8.

Fig. 1
Standard solid (also known as Zener model) being the rheological representation of the discussed viscoelasticity models 1 Note that alternative approaches which make use of an additive decomposition of the rate of deformation tensor [11] or of the stress power [12] lead to the identical constitutive equations as the multiplicative split of the deformation gradient in case of a Maxwell element [, 12, 13].

Basic Kinematics
To describe the deformation of the Maxwell element depicted in Fig. 1, the multiplicative decomposition of the deformation gradient into an inelastic (index i , dashpot) and a subsequent elastic (index e , spring) part is employed, see for instance Haupt [15]. This approach is often referred to as Sidoroff decomposition and has the advantage that the intermediate configuration between the elastic and inelastic deformation can be interpreted as the stress-free state which is instantaneously obtained after removal of external loads. Moreover, it is weakly invariant under an isochoric change of the reference configuration (w-invariant), see Shutov & Ihlemann [16] and Shutov [17]. That is, a transformation of the reference configuration can be counterbalanced by an appropriate transformation of the initial conditions. W-invariance is desirable for example for the simulation of multi-stage processes [18] or prestressed biological tissues [19].
The deformation gradients' polar decompositions read where R denote rotation and U stretch tensors. The corresponding spatial velocity gradients are defined as which can each be additively decomposed into a symmetric and a skew-symmetric part known as rate of deformation tensor d and spin tensor w , respectively, The symmetric and skew-symmetric parts of a tensor are defined as sym(X) = 1∕2 X + X T and skew(X) = 1∕2 X − X T . After some algebra, see for example Korelc & Wriggers [20], the identities are obtained. Further kinematic quantities, used in the following discussion, are the left ( b ) and right ( C ) Cauchy-Green tensors and their first principal invariants The deformation measures introduced in Eqs. (1), (9) and (10) can be multiplicatively split into an unimodular part X = det −1∕3 (X) X and a spherical part det 1∕3 (X) I with I being the identity tensor. This decomposition separates the isochoric, volume-preserving distortion from the volumetric dilation. The first isochoric principal invariants of the Cauchy-Green tensors are denoted by and the volume change by Furthermore, the right Cauchy-Green tensors in Eq. (10) and the rate of deformation tensors in Eq. (4) are linked by the identities

Thermodynamics
The local, isothermal Clausius-Planck inequality per unit reference volume in terms of the Helmholtz free energy is given by with X ∶ Y = tr X ⋅ Y T , see for instance Haupt [15]. Basically, this inequality demands that for thermodynamic consistency the dissipation rate D m must be non-negative at any deformation state. D m is equal to the mechanical stress power ∶ d less the change of the free energy Ψ . Herein, denotes the symmetric Kirchhoff stress.
The total free energy of the standard solid in Fig. 1 reads The scalar state variable Ψ 0,max is needed for modeling the Mullins effect of rubber materials in Sect. 2.6.1, see also Ricker et al. [21] for details, and is omitted for brain tissue

Viscoelastic Evolution Equation
Material models that account for inelastic material behavior typically introduce internal state variables and associated evolution equations. For a physically plausible behavior, the evolution equations should always ensure a non-negative dissipation rate, cf. Eq. (15). In what follows, four equivalent viscoelasticity formulations with different internal variables are derived which satisfy inequality (27) by prescribing constitutive equations for D i and W i . They are consecutively referred to as formulation A to D.
Interpreting Eq. (3) 3 as a differential equation the evolution of the inelastic deformation gradient F i is determined by prescribing an inelastic velocity gradient l i or equivalently L i , cf. Eq. (6) 2 . To guarantee a nonnegative D m,neq in Eq. (27), one may define so that where ‖X‖ = √ X ∶ X ≥ 0 denotes the Frobenius norm. The viscosity > 0 is associated with the dashpot of the Maxwell element and will be discussed in Sect. 2.5. The choice W i = 0 is made for convenience rather than physically motivated. In contrast to non-affine, plastic deformations of metals where the spin may be related to the orientation of crystal slip systems, the micromechanical interpretation of W i for the viscoelastic distortion of filled rubbers or biological tissue remains unclear. Note that, in general, the constitutive equation (29) does not lead to a rotation-free elastic or inelastic deformation, viz., R e ≠ I and R i ≠ I . A detailed discussion on the rotation of the intermediate configuration is given for instance by Boyce et al. [23] or Dafalias [24].
Using Eqs. (6) 2 and (19) 1 , assumption (29) leads to the evolution equation of the inelastic deformation gradient In the present manuscript, Ψ neq is assumed to be an isotropic function of C e , see Sect. 1, and hence Ψ neq ∕ C e and C e commute. Thus, l i is symmetric implying l i = d i and w i = 0 . Furthermore, Ψ neq is assumed to be an isochoric function, i.e., Ψ neq is a function of the unimodular part C e . Then, C e ⋅ Ψ neq ∕ C e and, hence, l i are traceless. This implies a volume-preserving viscoelastic flow, viz., det F i = const. as well as a deviatoric non-equilibrium Kirchhoff stress, viz., tr neq = 0.
Alternatively, one can write Eq. (31) in terms of L i and neq reading As it will be shown in Sect. 2.3, contrary to formulation A, the corresponding time integrating scheme uses the elastic deformation gradient F e as the internal variable.
Since l i = d i holds true for isotropic Ψ neq , Eqs. (14) 3 and (23) 2 can be employed to rewrite the definition of l i in Eq. (31) as Moreover, to obtain a simplified stress calculation rule, the identity is applied leading to the 2nd Piola-Kirchhoff stress see for instance Haupt [15]. This is particularly helpful in case of invariant-based strain energy functions with Ψ neq C e = Ψ neq C ⋅ C −1 i yielding Here, C i is used as the internal variable and neither C e nor the unsymmetric F i , F e are present in the constitutive (32) formulation B: equations. Moreover, the stress definition in Eq. (36) shows that Ψ neq is in general not an isotropic function of C and, hence, S neq and C do not commute. That is, the model response exhibits anisotropic properties where C −1 i acts similarly to the structural tensor of transverse isotropy.
Finally, another representation for isotropic Ψ neq in terms of b e results from the identity Eq. (23) 1 applied to Eq. (33). It reads with the Lie-derivative The coaxiality and commutativity of neq and b e follow from the isotropy of the strain energy function. An overview of all four formulations and their properties is given in Table 1.

Time Integration
The implicit exponential map is a first order time integration scheme for the differential equations of the form where F (n) i denotes the inelastic deformation gradient at the beginning of the current time increment t (n) . For the sake of simplicity, the index (n + 1) referring to quantities at the end of the current increment is omitted: t = t (n+1) . The algebraic equation (39) has to be solved for F i in each increment. To prove the volume-preservation of the exponential map, the determinant is applied on both-hand sides of Eq. (39). Employing the identities det(X ⋅ Y) = det(X) det(Y) and exp(tr(X)) = det(exp(X)) , it can be shown for isochoric Ψ neq with tr T neq ⋅ C e = 0 , cf. Sect.
In the present manuscript, the initial condition F (0) i = I is prescribed so that det F i = 1 ∀t.
A p p l y i n g E q s . ( 6 ) a n d ( 1 ) a s w e l l as exp A −1 ⋅ X ⋅ A = A −1 ⋅ exp(X) ⋅ A , the integration scheme (39) can be equivalently rewritten as (37) formulation D: formulation A: to be solved for F e . This can be considered as an update scheme for for mulation B, Eq. (32). Wit h 1∕det(X) = det X −1 , it is proven to be volume-preserving i n t h e s a m e way a s fo r m u l a t i o n A , v i z . , e . An update scheme for formulation C in terms of C i can be derived by applying the exponential map to Eq. (36), viz., with the initial condition C (0) i = I . Shutov & Kreißig [28] showed that the exponential map is symmetry-preserving for differential equations of the form Ẋ (t) = f (X, t) ⋅ X(t) with det(X(t = 0)) = 1 if the following conditions are fulfilled: tr(f (X, t)) = 0 and (f (X, t)) k ⋅ X(t) is symmetric for all k = 1, 2, 3, … . Here, the first condition is satisfied for isochoric Ψ neq since tr T neq ⋅ C e = tr C ⋅ S neq = 0 . Besides, this property also proves volume-preservation of the update scheme (41). The second condition is satisfied for isotropic Ψ neq (w.r.t. C e ) where C e ⋅ T neq = T neq ⋅ C e y i e l d s C ⋅ S neq ⋅ C i = C i ⋅ S neq ⋅ C a n d , h e n c e , Note that in comparison to Eqs. (39) and (40), the argument of the exponential is not symmetric in general resulting in higher computational costs. Furthermore, inverting and forwarding update scheme (41) to the current configuration provides a time integration scheme for formulation D, see evolution equation (37), reading with the internal variable b e . Since the operations X −1 and F ⋅ X ⋅ F T preserve symmetry, formulation D is a symmetrypreserving update scheme, too. Moreover, its volumep r e s e r v a t i o n . The implicit exponential maps of all four formulations, viz. Eqs. (39)- (42) are iteratively solved applying a Newton-Raphson scheme with a line search algorithm such that descending Newton steps are ensured. The state variable at the beginning of each increment serves as the initial guess for the iteration of the updated state variable. The iteration is (40) formulation B: Table 1 Overview of formulations for multiplicative viscoelasticity with W i = 0 , exponential mapping and isotropic free energy function (even though Simo [25] and Korelc & Wriggers [20] deal with elastoplasticity rather than viscoealsticity, they use the same kinematics) Local residual with exponen- Use in literature [20] [ 15,26] [ 25,27] assumed to be converged if the norm of the residual as well as the change of the unknown within the current Newtonstep are lower than 10 −8 .

Remark 1
Preliminary studies for the numerical benchmark tests in Sect. 7 showed that two aspects in the definition of the residual are advantageous for a fast convergence of the Newton-Raphson method: rearranging the update formulae such that the inverse of the exponential is present (i.e., a negative sign in the exponential's argument) and the inverse of the unknown is avoided. Appropriate definitions of the residual (denoted by R NR ) are given in Table 1.

Viscosity Functions
To fulfill the Clausius-Planck inequality, it is sufficient to define a positive viscosity , see Sect. 2.3. In the following, a list of viscosity functions presented in the literature is compiled 3 . Only models which have been proposed for rubbers or soft tissues are considered, i.e., approaches from fluid rheology or metal creep are beyond the scope of this work (although some viscosity functions in the list are borrowed from these fields of research). Equivalent representations of the employed stress and strain invariants for all four formulations A-D are given in Table 2. Material parameters are denoted by p, , , , , .
• The simplest, but very common Newtonian viscoelasticity with = const. is usually not the best choice to reproduce real material behavior. This approach serves as a reference herein. To improve the parameter identification procedure, the viscosity is fitted on a logarithmic scale, viz., This definition with the fitting parameter p leads to a fast, reproducible fitting and is applied to all viscosity functions below. Note that the units will be omitted in the following for the sake of readability. • A classical approach which is part of several approaches below is the Norton-type viscosity [29] with the power law Herein, ‖X‖ = √ X ∶ X denotes the Frobenius norm 4 . The sign of the stress exponent is chosen such that the viscosity decreases at large overstresses for positive -values. This effect is known as shear thinning or pseudoplasticity. A numerical drawback of this viscosity function is the behavior Hence, in the exponential map (Sect. 2.4), the term Δt∕ can readily produce overflows. To overcome such problems, either the parameter range of must be restricted or alternative time integration schemes must be considered. However, this issue is beyond the scope of this manuscript. • Haupt & Lion [30] as well as Plagge et al. [2] defined an exponentially stress dependent relaxation time in the context of convolution models for elastomers. It was also used earlier to describe the creep of metal crystals [31]. Herein, their approach is applied to a non-constant viscosity leading to • The hyperbolic sine power law by Garofalo [32] (43) = 10 p MPa s . Table 2 Equivalent stress and strain measures (note that for some conversions it was made use of the symmetric and traceless properties stemming from the choice of an isotropic and isochoric free energy function) Trace of Cauchy-Green tensor Note that for convenience some authors prescribe the fluidity = 1∕ (i.e., the reciprocal of the viscosity) or the effective creep rate ̇= ) rather than the viscosity itself. combines the behavior of a power law for small arguments and an exponential law for large arguments. Hurtado et al. [33] proposed its application to rubber materials. • Lion [26] proposed the viscosity function He stated that "for physical reasons, the material function [...] may depend on the arguments T neq and e i " is the Almansi strain. However, the interpretation of the intermediate stress T neq remains difficult since it is an artificial quantity obtained by a pull-back of the overstress neq (similar to the 2nd Piola-Kirchhoff stress S neq ). Moreover, a physically deviatoric stress does not lead to a traceless T neq . Hence, its Frobenius norm is not equivalent to the von Mises invariant (see footnote 3) such that the meaning of [6] derived from physical considerations the four parameter function which is multiplicatively comprised of the Nortontype viscosity, see Eq. (44), and a strain dependent contribution. The sign of the strain exponent is chosen so that the viscosity increases at large inelastic deformations for a positive . The material parameter is introduced due to "some difficulty in the numerical implementation" [6] and fixed to 0.01 for rubber materials. However, nearly zero viscosity values are still obtained for tr b i → 3 and > 1 such that numerical issues as discussed for the Norton-type viscosity can be aggravated. • The strain-rate dependent Carreau model is originally used to describe the non-Newtonian flow of fluids. It is adapted to brain tissue by Bilston et al. [34]. Note that the inelastic D i is employed here rather than the total rate of deformation tensor d as by Bilston et al. [34]. This choice is motivated by Fig. 1 and the assumption that the viscosity of the dashpot should only be affected by the corresponding inelastic deformation and not by the total deformation of the Maxwell element. 5 The parameter describes the viscosity at large deformation rates relative to small rates, viz., also originates from the fluid rheology. Hrapko et al. [35] employed the Ellis model for the simulation of brain tissue. • Prevost et al. [8] modelled the viscoelastic behavior of brain tissue with a modification of Eq. (48) [6]. The strain exponent is fixed to = 2 yielding after some rearrangements Contrary to Bergström & Boyce [6] and Prevost et al. [8], the parameter is herein not fixed to a pre-specified value. • Hur tado et al. [33] proposed the general form � and presented the formulation called "power law strain hardening model". Note that i is a monotonically increasing variable since ̇i ≥ 0 . Hence, the viscosity always tends to infinity (for 0 < < 1 ) or to zero (for < 0 ). The physical motivation for this irreversible, accumulating behavior induced by inelastic strain is not specified by Hurtado et al. [33]. • Kumar & Lopez-Pamies [36] proposed a six parameter approach for rubbers which is an extension of the Ellis model to a strain dependency. • Recently, Dal et al. [3] considered the relaxation kinetics of a single polymer chain resulting in Contrary to the other strain dependent viscosity functions, this approach employs an elastic invariant rather than an inelastic one. Note the similarity of the last term to a neo-Hooke free energy, cf. Table 4. The viscosity functions are summarized in Table 3. Approaches requiring additional numerical effort (for instance, numerical integration stemming from polymer chain statistics [37]) are beyond the scope of this comparison. Furthermore, since only decoupled, purely isochoric approaches are considered, viscosities depending on the hydrostatic presssure 6 p [39] or on the volume change J are not discussed. Summing up, hitherto many stress and strain dependent viscosity functions exist but there is no common opinion on the exact stress and strain measures to be used, i.e., intermediate or Kirchhoff stress and inelastic, elastic, total or accumulated strain. However, it can be noticed that in the source papers an increasing stress is always assumed to lead to a decreasing viscosity which can be motivated by shear thinning. Vice versa, increasing inelastic strain is associated with an increasing viscosity.
Hyperbolic sine power law Garofalo [32] Strain hardening power law Hurtado et al. [33] Kumar & Lopez-Pamies [36] 10 p 6 As shown by for instance by Champagne et al. [38], rubbers can approach the glass transition also due to a large compressive load. Thus, a pressure dependent viscosity can be a reasonable modeling approach to consider the concomitant change of the viscoelastic properties in the glassy state.  (1) is chosen with the material parameters m ∈ (0, ∞) and r ∈ [0, 1) , see Appendix 2 for the construction of H. The basic strain energy function is chosen to be as proposed by Plagge et al. [2]. It is called simplified extended tube model and provides a good fitting to the virgin curve of the tested materials, see Sect. 3, at low number of parameters, see also Ricker & Wriggers [43]. G 0 > 0 is the initial shear modulus and n > 0 is the polymer chain extensibility defining a pole at Ī 1 − 3 → n . The volumetric response is determined by the assumption of perfect incompressibility, see Sect. 4.
For the Maxwell element, all viscosity functions from the previous section 2.5 as well as the non-equilibrium free energy functions (a)-(e) in Table 4 are fitted to the experimental data to identify suitable approaches. The neo-Hooke potential is taken as reference free energy function. The detailed fitting procedure and the results are presented in Sect. 5. Compared to (a), the free energy function (b) introduces one additional parameter to generate an upturn of the viscoelastic overstress, i.e., an increasing shear modulus. In contrast, free energy function (c) with three parameters has a decreasing modulus (if B > 0 ). Approach (d) is identical to the basic free energy of the equilibrium spring and exhibits a pole. Approach (e) provides a Mullins-type damage for the non-equilibrium spring which is constructed in the same manner as the equilibrium free energy in Eq. (55) with H as in Eq. (58). Note that free energy function (e) is used only for rubber materials and (f) only for brain tissues due to the choice for the equilibrium free energy functions. (a) Neo-Hooke (55) and (58)

Material Model for Brain Tissues
Brain tissue exhibits stress softening similar to the Mullins effect of rubber materials as outlined in Sect. 1. Hence, the same modeling approach as in Sect. 2.6.2 could be used. However, the considered experimental data, see Sect. 3, include only one amplitude level for each loading mode which in preliminary studies turned out to be insufficient to provide reasonable material parameters. Therefore, a purely hyperelastic model for the equilibrium response is employed and, hence, preconditioning effects are not accounted for.
The first load cycle of the experimental data is omitted for the parameter identification and only considered as loading history.
A crucial difference between the behavior of rubbers and brain tissues is the pronounced asymmetric tensioncompression stress response of the latter [, 7,9]. That is, the material stiffness under uniaxial compression is much larger than under uniaxial tension. This behavior manifests itself when fitting a one term Ogden model to tension and compression data. In case of rubber materials, > 0 and c 10 > c 01 can typically be observed whereas for soft tissue the opposite behavior is present ( < 0 , c 10 < c 01 ). Furthermore, brain tissues show a pronounced strain stiffening, see also Budday et al. [5], leading to large absolute values of the Ogden parameter (typically, = −20 … − 40 ). It should be noted that such large exponents are numerically undesirable. However, the Mooney-Rivlin model (and most other free energy functions, especially Ī 1 -Ī 2 -based approaches) fail to capture properly these two effects. Therefore, the Ogden model in Eq. (60) is commonly used and also employed here.
The description of the volumetric response bases again on the incompressibility assumption. The viscosity functions to be investigated are identical to those for rubber materials, cf. Table 3. Furthermore, non-equilibrium free energy functions (a)-(d) and (f) from Table 4 are considered.

Material Model for Numerical Tests
For numerical tests under realistic conditions, a model setup which will be identified to be well-suited for rubber materials is chosen, see Sect. 5. The equilibrium free energy is given by Eqs. (55), (58) and (59). For the sake of simplicity, is chosen for the volumetric part leading to a linear p-Jrelation, i.e., to a constant bulk modulus K = 2 Ψ vol ∕ J 2 . As viscosity function, the stress and strain dependent function according to Eq. (51) is considered [8]. Moreover, a neo-Hooke model is chosen for the non-equilibrium stress contribution, cf.

Experiments
To find the best viscoelasticity model for rubber materials, cyclic, uniaxial tensile tests including relaxation phases and several amplitudes were conducted with two industrial compounds, see Figs. 2 and 3. The first compound is a sulfur crosslinked, carbon black filled EPDM for sealing applications which is loaded up to 155 % strain at a strain rate of 10 %∕s . Whereas the second one is a sulfur crosslinked, carbon black filled NR for vibration isolation with a maximum load of 210 % strain at 1 %∕s . Note that the latter shows a widening of the hysteresis loop at > 2 which is probably due to strain-induced crystallization. For brain tissue, experimental data are taken from Budday et al. [9] where the average material response from multiple brain tissue samples was presented. That paper is chosen as data source since it presents an extensive data base including both cyclic and relaxation experiments in multiple deformation modes. The cyclic tests are comprised of three displacement-driven, sinusoidal load cycles. The amplitudes are = 0.2 , = 0.9 and = 1.1 under simple shear, uniaxial compression and uniaxial tension, respectively. See Fig. 4 for the stress-stretch plots and the applied frequencies. For the relaxation tests, the tissue is loaded with a rate of 100 mm∕min to the above specified amplitudes and afterwards the displacement is kept constant for a duration of 300 s , see Fig. 5. Relaxation data are only provided for compression and shear in the source paper and consequently relaxation under tension is not considered herein. Furthermore, data for specimens from multiple brain regions are available. In the following, the viscoelastic models are compared using data from the cortex and the corpus callosum as they exhibit the stiffest and softest material responses, respectively, from all tested brain regions.

Parameter Fitting
To find the best parameter set for each model and material, a least square problem is formulated and solved by a Trust-Region algorithm. For the rubber materials with uniaxial tension data, the residual is defined in terms of the relative error between the model's and the experimental 1 st Piola-Kirchhoff stress in loading direction. That leads to the cost function m denotes the number of considered experimental observations, i.e., load steps. The model's 1 st Piola-Kirchhoff stress P mod,i is obtained from the Kirchhoff stress mod divided by the measured stretch exp,i (formulations A, B, D) or from the 2nd Piola-Kirchhoff stress S mod,i multiplied by exp,i (formulation C). The stress is computed from the experimental deformation gradient F exp,i , the state variable i−1 at the beginning of the increment and the material parameters p j , j = 1 … n . must be replaced by F i , F e , C i or b e depending on the considered formulation.
For the brain tissues with multiple test data and deformation modes, the cost function reads where F k are the contributions from the shear, compression and tension experiments. Here, a normalized error is employed in the residual for each experiment rather than the relative error, viz., since the cyclic brain data include stress values close or equal to zero where measurement noise would lead to arbitrarily large relative errors. The normalization w.r.t. the maximum stress of each experiment avoids a bias towards experiments with large stress values (typically compression experiments for brain tissue, cf. Sect. 3).
For the fitting procedure, perfect incompressibility is assumed. Consequently, the components of the experimental deformation gradient are constructed from det F exp,i = 1 as Moreover, the hydrostatic pressure p in the model is not computed from the volumetric free energy but from boundary conditions, i.e., the lateral directions are assumed to be stress-free. The Jacobian for the optimization problem is obtained as

Comparison of Viscelasticity Models for Rubber Materials
To find the best viscoelasticity model for the tested rubber materials, a two-step procedure is applied. First, the nonequilibrium free energy function is defined as a neo-Hooke potential Ψ neq = G neq ∕2 Ī 1,e − 3 and is combined with all viscosity functions given in Table 3, see Sect. 5.1 for the results. Then, the constant viscosity is chosen and combined with the non-equilibrium free energy functions (a)-(e) in Table 4, cf. Sect. 5.2 for the outcome. The equilibrium free energy function Ψ eq given by Eqs. (55), (58) and (59) is kept constant during this study. Since preliminary studies showed that the model framework with just one Maxwell element is insufficient to capture accurately both viscoelastic phenomena, i.e., hysteresis loops and stress relaxation, this procedure is done twice for each material. On the one hand, the models are fitted only to the cyclic data, i.e., the (65) tension/compression: non-transparent lines in Figs. 2 and 3. On the other hand, only the relaxation data are considered, i.e., the transparent lines. This manuscript sticks to a single Maxwell element because multiple Maxwell elements would lead to an unmanageable number of combinations of free energy and viscosity functions with hardly foreseeable interdependencies. Moreover, separated fittings to cyclic and relaxation data allow to identify appropriate response functions for each viscoelastic effect and to interpret the role of the material parameters. The results can then be employed to construct more sophisticated material models with two or even more Maxwell elements capturing complex load scenarios. The range of feasible parameter values for the fitting procedure is chosen as generous as possible without violating the second law of thermodynamics, see Sect. 2.2. That is, > 0 must be satisfied to ensure a non-negative dissipation rate of the dashpot. For instance, the stress and strain exponents and in the viscosity functions are allowed to be positive as well as negative. In some cases, that is more generous than originally proposed by the authors, e.g., Bergström & Boyce [47] limited ∈ [0, 1] for viscosity function 6 because of the micro-mechanically motivated background based on the reptational dynamics of polymer chains. Moreover, very large exponents can cause numerically impracticable model behavior. However, finding appropriate parameter limits for each viscosity function and free energy function is not reasonable for a large number of models. In addition, the parameters can be very different for rubbers and soft tissues and promising models can readily be overlooked due to too strong constraints. Thus, for the sake of a fair comparison, only strictly necessary parameter bounds are applied, see Table 3.

Comparison of Viscosity Functions
Seeking for the best viscosity function, each approach from Table 3 is combined with a neo-Hooke non-equilibrium free energy for the Maxwell element. The equilibrium free energy is described in Sect. 2.6.1. The fitting results for cyclic and relaxation data of both rubber materials are visualized in Fig. 6.
The reference model 1 in this comparison employs a single, constant viscosity. This assumption clearly leads to the poorest outcome for both rubber materials and loading modes. Its stress response is depicted in Figs. 25 and 26 in Appendix 5 and reveals an insufficient reproduction of the material behavior. In contrast, viscosity function 11 by Kumar & Lopez-Pamies [36] ranks first in nearly all cases.
However, this approach requires six parameters whereas all others need four parameters at most. As a consequence, some parameters cannot be uniquely identified with the given data, especially for the NR compound. That is, a lot of variation in the optimized parameters can be observed between different initial guesses although the cost function value is nearly identical. This problem becomes also apparent in high parameter correlations, especially between the stress exponent and the stress scaling factor (as well as between and ).
Promising alternatives are the viscosity functions 5 by Lion [26], 6 by Bergström & Boyce [6] and 9 by Prevost et al. [8] which require just two or three parameters. All these top ranked models have in common a stress dependency as well as a strain dependency on Ī 1,i . In contrast, approach 12 by Dal et al. [3] employing Ī 1,e in a neo-Hooke-like contribution (cf. Sect. 2.5) provides similar results as the purely stress dependent viscosity functions 2-4 and 8. Their stress dependency can notably improve the calibration on the relaxation data but not on the cyclic data. The dependency on the accumulated viscous strain i of viscosity function 10 appears to be helpful particularly for the fitting to the relaxation behavior of the EPDM compound. Interestingly, only for these data, the parameter was optimized to a negative value, i.e., the viscosity increases ( → ∞ for ̇i → ∞ , cf. Sect. 2.5) whereas the viscosity tends to zero for the EPDM cyclic data and both NR data sets. Furthermore, the strain rate dependency of the Carreau model 7 does not seem to be advantageous for rubber materials.
Viscosity function 4 by Lion [26] with only two parameters differs from the other approaches in terms of the employed stress and strain measure ‖T‖ and ‖b −1 i ‖ , cf. Table 3. Its convincing performance for the NR compound and for the cyclic behavior of the EPDM compound reveals that these are suitable measures to capture the viscoelastic behavior despite the missing physical interpretability of T neq , cf. Sect. 2.5. However, the model is not among the most promising ones for the EPDM relaxation data. A reason might be the fixed ratio between the stress and strain measure, viz., ‖T‖ ∕‖b −1 i ‖ 3 . This is in contrast to the other top-performing viscosity functions 6, 9 and 11, where the stress and strain dependency can be adjusted separately. These three models show a comparatively weak strain dependency for the EPDM relaxation data, e.g., model 9 has a small -and a large -value, see Table 5.
In addition to the non-Newtonian approaches with one Maxwell-element, a parallel connection of 13 Maxwell elements with constant shear moduli and viscosities is considered in the comparison, see approach 13 named relaxation time spectrum in Fig. 6. Employing the implicit but iterative-free time integration scheme by Shutov et al. [48] yields a numerically efficient implementation, see Appendix 4 for details. For the NR relaxation data, the relaxation time spectrum achieves a small improvement compared to a single Maxwell element whereas the EPDM relaxation data benefit considerably from multiple Maxwell elements. The fittings to the cyclic data do not show a notable improvement due to the constant strain rate in the   can lead to kinks at turning points in the stress-strain plot. = 6 in Table 5 is barely acceptable. 4. All models struggle to predict large changes in stress after a relaxation phase, see for instance the time ranges 77 … 82 min and 83 … 88 min in Fig. 10a. 5. For certain parameter sets, the models can be highly sensitive to noise in the input strain data, e.g., see the second and third relaxation phase in Fig. 8a. This is due to the equilibrium spring rather than the Maxwell element. More specific, the steep slope at the upper turning points in Fig. 8b leads to the sensitive behavior and is caused by the large damage parameter m. 6. A visual inspection of the stress-strain curves of the calibrated model is always recommended. That is, the stress response in other deformation modes (e.g., simple shear, compression, biaxial tension), with different experimental protocols (e.g., cyclic, relaxation, creep), beyond the experimental strain range and with different strain rates should be checked for plausibility. Such predictability is of particular importance for three-dimensional simulations with arbitrary deformation states. For instance, in Fig. 8b, the large overstresses at the lower turning points and the steep stress upturn at the upper turning points of the re-and unloading curves do not reflect the real material behavior.

Comparison of Non-equilibrium Free Energy Functions
After the identification of promising viscosity functions in the previous section, a study on suitable free energy functions for the non-equilibrium spring is conducted. The model setup is comprised of the equilibrium free energy according to Sect. 2.6.1 and a dashpot with a constant viscosity. The non-equilibrium free energy functions (a)-(e) given in Table 4 are investigated. The neo-Hooke model (a) with a constant shear modulus serves again as a reference. The results are presented in Fig. 11. It can be noted that the non-Hookean viscoelasticity does not lead to an improved model calibration for the cyclic data of both materials. Moreover, the EPDM relaxation data benefit only from approach (c) by Yeoh [45]. In contrast, the relaxation behavior of the NR compound gains from all non-Hookean free energy functions. Concluding, the non-Newtonian approaches with a constant shear modulus (see previous section) are in general much more advantageous than a non-Hookean Maxwell element with a constant viscosity.
Free energy function (c) ranks first for the relaxation data of both materials and reduces the RMSE by 23 % and 15 % compared to the neo-Hooke spring. Its unique feature is its decreasing shear modulus. However, the visual inspection of the stress relaxation curves revealed a slightly concave shape for some relaxation phases, e.g., between 34 min and 38 min in Fig. 12. This is contrary

Remark 2
The parameter fittings for the rubber materials were rerun with the viscoelastic models employing the reverse order in the multiplicative split F = F i ⋅ F e , see Appendix 3. The optimized parameters and the goodness of fit did not notably change.

Remark 3
The study on suitable non-equilibrium free energy functions was repeated with viscosity function 9 by Prevost et al. [8] (i.e., a non-Newtonian viscoelasticity) rather than a constant viscosity. In this case, the improvements in comparison to the neo-Hooke potential (also combined with viscosity 9) were even smaller (3 % at most).
Remark 4 Furthermore, the model calibrations were redone with an absolute error in the cost function instead of a relative error. In this case, the stress and strain exponents in the viscosity functions tended to larger and hence numerically undesirable values, see also statement 3 in Sect. 5.1. Moreover, the cyclic data slightly benefited from a non-constant shear modulus when employing an absolute error, contrary to Fig. 11. However, the general tendencies and the ranking orders remained largely unchanged.

Comparison of Viscelasticity Models for Brain Tissues
The procedure for comparing the viscosity and the non-equilibrium free energy functions for brain tissue is similar to that for rubber in the previous section 5. That is, the viscosity functions from Table 3 are combined with the neo-Hooke potential and fitted to the experimental data. Subsequently, the non-equilibrium free energy functions (a)-(d) & (f) from Table 4 coupled with a constant viscosity are calibrated on the data. In both studies, the equilibrium stress is obtained from an Ogden model and parameter bounds are prescribed only to guarantee a positive dissipation rate, i.e., a positive viscosity. These bounds are presented in Table 3. Moreover, the cyclic data and the relaxation data are fitted separately as a single Maxwell element is not sufficient to capture both phenomena simultaneously, see also the model calibrations by Budday et al. [9].

Comparison of Viscosity Functions
The results obtained with a neo-Hooke non-equilibrium spring and different viscosity functions are depicted in Fig. 13. Some conclusions similar to those for the rubber materials can be drawn. For instance, a constant viscosity provides the worst results and is not sufficient to reasonably capture the viscoelastic effects, see Figs. 27, 28 and 29 in Appendix 6. Moreover, the Careau model 7 depending on the strain rate ‖ ‖ D i ‖ ‖ provides only little improvement. The best results are obtained with the overparameterized viscosity function 11 [36].
In comparison to rubber materials, there are also some significant differences. Most of the models are suitable to notably improve the fit to the relaxation data of both brain tissues. The numerous good fits are possibly due to relaxation experiments with a single strain amplitude. Accordingly, tests with more strain amplitudes could refine this benchmark. Looking for the most promising candidate for relaxation, one approach is to be highlighted. The Ellis model 8 provides sound results, see Fig. 16, and reproducible parameters with different initial guesses. The cost function is reduced by 95 % (corpus callosum) and 54 % (cortex) compared to a constant viscosity. Interestingly, the -p a r a m e t e r w h i c h d e f i n e s t h e r a t i o is fitted for both tissues to zero, being the lower parameter bound. Therefore, to reduce the number of fitting parameters, = 0 can be fixed such that a reduced version with only three parameters is obtained. Looking at the fits to cyclic data, only viscosity functions 6 [6] and 11 [36] perform well. This can be explained by the strain and stress exponents which are both Table 5 Parameters obtained from fitting with viscosity function 9 by Prevost et al. [8] and a neo-Hooke free energy for the non-equilibrium spring (to anonymize the experimental data, the parameters are modified to match the normalized stress data in Figs. 2 and 3) Compound, loading Parameters fitting parameters. In contrast, the other strain dependent viscosity functions 5, 9 and 12 prescribe a fixed strain exponent of 3, 2 and −1 , respectively. Thus, how a change in stress and strain leads to an increase or decrease of the viscosity, is not adjustable. This effect becomes even more apparent when the well performing viscosity function 6 is compared to 9. The latter is derived from the former, see Sect. 2.5. While model 6 fixes the strain scaling factor , model 9 keeps the strain exponent constant, leading to quite different results. Note that fitting both the exponent and scaling factor, like in model 11, typically results in a large correlation between these parameters. Contrary to approach 11, viscosity function 6 [6] generates reproducible, less correlated parameters and is hence considered as a more suitable function to model the cyclic data of brain tissues. It reduces the cost function by 21 % (cortex) and 26 % (corpus callosum). The model stresses are illustrated in Fig. 14 for the cortex and in Fig. 15 for the corpus callosum. Looking into details, two critical aspects should be mentioned. Firstly, for the corpus callosum tissue, the viscosity function 6 generates a disputable kink at P = 0 in the stress-stretch curves for the tension and compression data, see Fig. 15b. This is due to the strong tension-compression asymmetry and the very large overstress at = 1 in the compression test which is even larger than the maximum stress of the tensile test. However, none of the models can capture this material behavior reasonably. Secondly, viscosity function 6 yields an unexpected kink for the cortex tissue at the upper turning point during the preconditioning cycle, see Fig. 14a. Indeed, the preconditioning data are not considered for fitting, but they reveal the limited predictability.
Concluding, using a non-Newtonian viscosity function for modeling brain tissues significantly enhances the fitting quality both for relaxation and cyclic loading modes. The cyclic data are best approximated when the viscosity function holds a strain and stress exponent as fitting parameter.

Comparison of Non-equilibrium Free Energy Functions
To identify suitable non-Hookean viscoelasticity models for brain tissues, the non-equilibrium free energy functions (a)- Table 4 are combined with a constant viscosity. The results are depicted in Fig. 17.
For the relaxation data, the largest improvements in the cost function compared to the neo-Hooke free energy are 25 % (for the cortex tissue with free energy function (c)) and 28 % (for the corpus callosum tissue with Ogden model (f), cf. Fig. 18). In contrast, the fit to the cyclic data of the cortex and corpus callosum gains at most by 12 % and 4 % , respectively, using free energy function (f) and (b). These findings promote similar conclusions as for the rubber materials. First of all, the effect of a non-Newtonian approach is much larger than of a non-Hookean approach, cf. Sect. 6.1. Moreover, the effect of the latter is more pronounced for relaxation data than for cyclic data. Furthermore, the simplified extended tube model (d) is again not a reasonable choice for any data set.

Remark 5
The study was redone employing viscosity function 6 [6] for the cyclic data and 8 (Ellis model) for the relaxation data. In this case, the improvements in comparison to the neo-Hooke potential (also combined with viscosity function 6 and 8, respectively) were even smaller. See also remark 3.

Comparison of Viscoelasticity Formulations
To identify the numerically most robust and fastest viscoelasticity formulation, the model versions A-D were implemented in the updated Lagrange framework of the commercial finite element software MSC Marc via the HypEla2 subroutine. Details on the HypEla2 subroutine are given in Appendix 1. The programming was done with the Mathematica Add-On AceGen [20] allowing an automated code generation from the symbolic Wolfram language such that identical conditions and analytical material tangents are used for all formulations. The computation of the matrix exponential according to Hudobivnik & Korelc [49] is employed. For all numerical tests, the model definition from Sect. 2.6.3 was used.
In theory, formulations C and D are preferable in terms of numerical efficiency because their update schemes with the symmetric state variables C i and b e lead to a nonlinear system of only six equations, see discussion in Sect. 2.4. Moreover, in case of formulation C, the material tangent for an implicit finite element implementation can be computed more efficiently in terms of C instead of the unsymmetric F . However, the matrix exponential has to be applied to the unsymmetric argument Δt∕ C ⋅ S neq resulting in higher computational costs. See Table 1 for an overview of these properties. The following investigations show which formulations indeed perform well in terms of number of iterations and robustness.
Firstly, a cylindrical geometry is simulated, see Fig. 19a. The aspect ratio h∕d = 1∕4 is taken from Shutov et al. [48]. A cyclic, both-sided shear test with a strain rate of 10 %∕s and amplitudes of 60 % nominal strain according to Fig. 19b is simulated. Subsequently, a cyclic tensioncompression load is applied with a strain rate of 10 %∕s and +30 %/−10 % strain amplitudes. Making use of the symmetry, one half of the geometry is meshed with 540 H1P0 elements. Moreover, three temporal discretizations are tested with 160, 320 and 640 increments for the shear as well as tension/compression load with a simulation time of 120 s and 40 s , respectively. No automatic incrementation (cutbacks or adaptive time stepping) is allowed for the sake of equal conditions. The material parameters optimized for NR cyclic data are used, see Table 5.
It is observed that all formulations yield the same global response for each temporal discretization, proving their equivalence. Furthermore, the total number of global iterations (i.e., to find the nodal unknowns of the assembled system) is nearly identical for all formulations. In contrast, the number of local iterations (i.e., to find the updated state variable at the integration points) differs significantly, see Fig. 20. Formulation D needs the lowest number of iterations, closely followed by B. On the other hand, formulation A and C are far behind with  Table 4 combined with a constant viscosity Fig. 12 Stress response with a constant viscosity and the non-equilibrium free energy function (c) [45] fitted to the relaxation data of the EPDM compound almost 40 % and 50 % more iterations for the largest time step size. For the smallest time steps, these differences are less pronounced but still around 20 %.
The well performing formulations B and D have in common that they use elastic state variables F e and b e . In contrast, the choice of using either a deformation gradient (A & B) or an Cauchy-Green tensor (C & D) as the state variable is of lower importance. These observations are probably due to smaller changes in the elastic deformation than in the inelastic deformation, viz., the non-equilibrium spring is stiffer than the dashpot. Thus, the initial guess for an iteration in formulation B or D is more likely closer to the solution and, hence, less iterations are needed. Note that a generalization of these conclusions must be made carefully as the model behavior is material parameter dependent. However, since realistic material parameters are chosen, this outcome has a practical relevance.
For a second numerical study, a twisted beam with a quadratic cross-section a 2 and a length = 5a is simulated. An increasing amplitude up to 540 • with a constant frequency is applied, see Fig. 21. A mesh of 2 × 2 × 14 = 56 H1P0 elements and a time stepping with 160, 320 and 640 increments are considered. The material parameters identified for EPDM cyclic data are used, see Table 5.
The global response is depicted in Fig. 22 showing that the maximum applicable twist depends on the formulation and the time step size. This ranking is consistent with the previous results of the sheared cylinder. That is, a robust behavior is observed with formulations D and B and poor results are obtained with formulations A and C. Interestingly, the former two fail to converge near the turning points, i.e., when changing from loading to unloading whereas the latter two struggle near the unloaded state, i.e., when crossing zero load.

Remark 6
The numerical tests were additionally conducted with a constant viscosity instead of the viscosity function 9 by Prevost et al. [8]. This much simpler material model showed similar results for all formulations.

Remark 7
The numerical studies were also carried out with viscoelastic models employing the reverse order in the multiplicative split F = F i ⋅ F e , see Eqs. (84) and (85) in Appendix 3. The global reaction forces and torques changed only slightly but the numerical behavior was poor. That is, they failed to simulate the sheared cylinder in Fig. 19 even with a constant viscosity and the smallest time step size. Moreover, for the twisted beam in Fig. 21, the maximum applicable twist was less than for formulation C.

Remark 8
In addition, the update scheme by Shutov et al. [48] based on formulation C, see Eq. 91 in Appendix 4, was tested with viscosity function 9 and a constant shear modulus. The global responses for all three time step sizes in both numerical tests were quasi-identical to the exponential maps. The performance of the scheme was good, similar to formulation B and D. For the sheared cylinder, the number of local iterations were 3 … 6 % larger than for formulation D, cf.

Conclusion
The applicability of non-Newtonian and non-Hookean viscoelasticity to the simulation of rubber materials as well as brain tissue under large deformations was analyzed. For this purpose, the modeling framework of a standard solid based on the multiplicative decomposition of the deformation gradient into an inelastic and an elastic part was derived. As time integration scheme for the evolution equation, an exponential map was applied. On the one hand, the fitting quality of twelve viscosity functions as well as five free energy functions for the non-equilibrium spring were studied for two rubber compounds and two brain tissues. On the other hand, the numerical properties of four equivalent formulations of this framework were tested to identify the most robust and fastest implementation.
For the tested rubber materials, the studies revealed that the use of non-Newtonian viscosities is beneficial compared to the use of a constant viscosity. Moreover, these viscosity functions also outperform a generalized Maxwell element which is often used in the literature. The stress as well as strain dependent viscosity functions by Prevost et al. [8], Bergström & Boyce [6] and Lion [26] were identified as a good compromise between fitting capability and model complexity in terms of number of parameters. That is, they reproduce well both viscoelastic phenomenons, i.e., hysteresis loops and stress relaxation, and the optimized parameters are reproducible with different initial guesses. The study of different free energy functions for the non-equilibrium spring showed only small improvements compared to a neo-Hooke approach. In addition, some free energy functions provided physically doubtful behavior, i.e., a concave instead of a convex shaped stress-time curve during a relaxation phase.
Some conclusions for the rubber materials also apply for brain tissue. For instance, a non-Newtonian modeling has a larger impact on the fitting results than a non-Hookean approach and a strain dependent viscosity  Fitting results for brain tissues with non-Hookean viscoelasticity: free energy functions from Table 4 combined with a constant viscosity 1 3 function is important to capture particularly the cyclic behavior. For brain tissue, the Ellis model (with = 0 fixed) and the viscosity function by Bergström & Boyce [6] are the most promising candidates for relaxation and cyclic data, respectively.
The study on the numerical properties showed that the viscoelasticity formulations which employ the elastic deformation gradient or the elastic left Cauchy-Green tensor as internal state variable allow large mesh distortions and need few local iterations. The former of these two formulation has the advantage that it is applicable to anisotropic free energy functions but it needs slightly more iterations and uses an unsymmetric state variable.
For future works there are some interesting aspects which were beyond the scope of the present manuscript. On the one hand, experiments with multiple strain rates should be conducted to test the range of the models' validity and predictability. For rubber materials, experiments under pure shear or equibiaxial tension can be of interest to investigate the effect of Ī 2 -dependent free energy functions.  On the other hand, it can be worth to add plastic behavior to the model framework since a permanent set can be observed for rubber and soft tissue. In the current model framework, these effects are also treated by the viscoelastic modeling. Moreover, the framework should be extended by one or more additional Maxwell elements such that both the cyclic as well as relaxation behavior are captured  simultaneously. For this purpose, the findings of this paper are a useful basis for the design of a sophisticated model.

Implementation of Material Models in MSC Marc via the HypEla2 Subroutine
The commercial finite element software MSC Marc provides a total Lagrange (TL) as well as an updated Lagrange (UL) framework for large strain simulations. When implementing user-defined material models via the HypEla2 subroutine, the variable iupdat is passed in to indicate which formulation is chosen, see Table 6 for an overview of relevant input variables. In both cases, the deformation gradient and the temperature at the beginning and at the end of the current time increment are available. Moreover, the state variables at the beginning of the increment are provided. The user must return the stress at the end of the increment, the material tangent as well as the change of the state variables. In case of the TL framework, the 2nd Piola-Kirchhoff stress S and the TL tangent = 2 dS∕dC are required. For updated Lagrange simulations, the Cauchy stress = 1∕J and the consistent UL tangent with Jaumann correction denoted by ℂ must be returned, see Table 7 for an overview of required output variables. A detailed derivation of the tangents is given by Ihlemann [50] and summarized by Ji et al. [51], see their Eq. (40). Note that MSC Marc employs the Voigt notation for the stress and tangent with the storage order xx, yy, zz, xy, yz, zx. Special considerations are needed when nearly incompressible materials like rubbers or soft tissues are simulated since MSC Marc makes use of Herrmann elements to overcome volumetric locking. In the following, additional programming aspects within the HypEla2 subroutine are briefly explained in terms of the TL framework (see Table 7 for the corresponding UL formulation). Assuming a decoupled isochoric-volumetric stress response, a pressure-like primary unknown p * (with lowered shape function order) and a corresponding constraint equation g(J, , p * ) = 0 are introduced where J and are the volume change and absolute temperature, respectively. Rather than defining p * to be equal to the hydrostatic pressure p obtained from the constitutive equations, it is reasonable to demand that p * = p∕f with a correction function f (J, ) , see Landgraf [52]. Thus, the constraint may be written as

Fig. 22
Results of the benchmark from Fig. 21: robustness of the viscoelastic models (the square, circle, cross and star mark the last increment before failure for each formulation; formulation C fails for the coarsest discretization at the first increment so that the cross is missing in (a)) For the sake of completeness, the treatment of isotropic thermal expansion is included here with the thermal 8 and mechanical volume change J and J m = J∕J . Moreover, denotes the initial bulk modulus. The presented definition of p * leads to the modified 2nd Piola-Kirchhoff stress with S iso and S vol = −f p * J C −1 being the isochoric and volumetric stress contribution. (70) The motivation for the correction function f stems from the additionally required linearizations 2 dg∕dC , dg∕dp * and dS∕dp * which have to be returned together with via the array cf. Table 7. f is needed to ensure a numerically more efficient, symmetric tangent matrix, viz., 2 dg∕dC = dS∕dp * . With Eqs. (69) and (71) the requirement g∕ J J C −1 = f J C −1 is obtained, yielding the differential equation g∕ J = f to be solved for the correction function f. Note that the solution depends on the chosen volumetric free energy Ψ vol J m . To (72) Herrmann formulation must be used   8 The thermal volume change may be defined as J = exp 3 ( Δ ) ≈ (1 + Δ ) 3 with the coefficient of linear thermal expansion , cf. Lu and Pister [53]. ensure that f = 1 and p = p * holds true in case of a constant bulk modulus with Ψ vol = K 0 ∕2 J m − 1 2 , cf. Eq. (61), the factor J 2 ∕K 0 was added in Eq. (69). Unfortunately, for some volumetric free energy functions, the solution for f leads to an indeterminate or non-continuously differentiable behavior for J, J → 1 requiring additional numerical treatment. On the other hand, Table 8 shows the correction functions of non-critical Ψ vol . The required output of HypEla2 subroutines for Herrmann elements is summarized in Table 7. In addition, conversions from the TL to the UL framework and vice versa are provided, see also Ihlemann [50]. The operators T 34 and sym 34 ( ) in the table denote the transposition and the symmetrization of the third and fourth indices, viz.
The operator ⊙ produces a fourth order tensor which exhibits minor symmetries allowing a representation in Voigt notation. Applying this operator to the deformation gradient and its inverse/transposed provides with a simplified calculation Appendix 2

Mullins-Type Damage
As outlined in Sect. 2.6.1, the equilibrium free energy function with discontinuous damage is obtained from the scalar function H Ψ 0 , Ψ 0,max where Ψ 0 describes the virgin load response and Ψ 0,max is its current maximum value. This section summarizes how to construct such a function and provides a small comparison of feasible approaches.
Ricker et al. [21] identified two generic classes used in the existing literature, namely H Ψ 0 ∕Ψ 0,max and H Ψ 0,max − Ψ 0 . The drawback of the former class is its behavior when going back to the unloaded configuration with Ψ 0 = 0 , for instance in a cyclic shear or tensile test with increasing amplitudes.
In this load scenario, the material stiffness is independent of the maximum load since H Ψ 0 ∕Ψ 0,max is then independent of Ψ 0,max . This behavior is contrary to experimental findings where the stiffness at the unloaded configuration progressively reduces for higher amplitudes, cf. Fig. 2. Therefore, the latter class of functions H Ψ 0,max − Ψ 0 is considered here. For readability, the definition x = Ψ 0,max − Ψ 0 holds for the following discussion. Constructing such a function, the following two requirements have to be fulfilled: For thermodynamic consistency, the requirement Finally, some requirements are added for a convenient handling: • the H-function employs two parameters m > 0 and 1 > r > 0 , of which the former scales the argument, viz., H(m x) and the latter defines the lower limit as lim x→∞ H(x) = 1 − r • the curvature does not change ( H �� (x) > 0 , i.e., strict convexity) • a closed-form representation of Ψ eq given by Eq. (55) in terms of elementary functions exists Such functions can be constructed for instance from the convex branch of a sigmoid function. A list of suitable H-functions is compiled in Table 9. They are plotted in Fig. 23a with r = 1 (leading to lim x→∞ H(x) = 0 ). Moreover, m is chosen such that H � (0) = 1 , i.e., all approaches show the same initial slope.
The results for the EPDM compound are summarized in Fig. 24a and sorted by the cost function value. Apparently, a good fitting result correlates with a steep slope of the H-function at x = 0 and a slow convergence for x → ∞ in Fig. 23b. An exception is the 1∕ 3 √ x-approach that exhibits the steepest initial slope and convergences slowly but does not rank among the top candidates. The fitted stress-stretch curves of the 1∕ √ x -, arcoth -and erf-approach and the cyclic data of the EPDM compound are depicted in Fig. 24b. The former two models provide a sound fit to the un-and reload curves, particularly near the turning points. At the upper turning point, they generate a steeper slope than the erf-approach, what is caused by the steeper initial slope of the H-function. At the lower turning point, the erf -function yields the same stress response for the third and fourth amplitude. This is due to its fast convergence, cf. Fig. 23b, and does not agree with the experimental observation. The other two approaches can qualitatively capture the real material behavior. Comparing the 1∕ √ x -to the arcoth-approach, the former shows a slightly better fit to the virgin curve whereas their remaining stress-stretch curves are nearly identical. However, an underestimation of the virgin curve is acceptable to some extent since the viscoelasticity model will generate some additional overstress.

Remark 9
The fitting was also conducted with the cyclic data of the NR compound. The ranking order was the same as in Figs. 24a.
Applying the exponential map to Eq. (82) yields an update scheme in terms of F i which can be rearranged to obtain an update formula in terms of F e corresponding to Eq. (83)  Table 9 for the cyclic data of the EPDM compound: (a) ranking and (b) stress-stretch plot of the two top-ranked as well as the reference model ( erf)

Generalized Maxwell Element with Relaxation Time Spectrum
The viscoelastic behavior of rubber materials is oftentimes modeled using a generalized Maxwell element (i.e., several Maxwell elements connected in parallel) with a discrete spectrum of constant relaxation times. To avoid the high computational effort of several exponential mappings, the time stepping method by Shutov et al. [48] is employed. It is an implicit time integration scheme based on a backward Euler method applied to formulation C and is an iterativefree update formula in case of a constant relaxation time. Its derivation is recapped in the following. With the derivatives of Ī 1,e = tr C e = tr C ⋅C . and the evolution equation is obtained as where = ∕ 2 Ψ � neq denotes the relaxation time. The Euler backward method provides the implicit, symmetrypreserving time integration scheme in terms of C i However, in general, it is not a volume-preserving time stepping method. To overcome this drawback, Shutov et al. [48] applied unimodularization to the right-hand side of Eq. (90). With the property a X =X for any non-zero, real scalar a and non-singular tensor X , the update formula simplifies to which is volume-preserving and, in addition, iterative-free for , Ψ � neq = const. Further numerically efficient update schemes are presented by Shutov [13] (for a Mooney-Rivlin strain energy function at the cost of a tensor square root) and Landgraf et al. [63] (for general Ī 1 -dependent strain energy functions at the cost of a scalar equation to be iteratively solved).
For the comparison in Secs. 5.1 and 6.1, 13 Maxwell elements in parallel are considered with 26 material parameters: 13 relaxation times k and 13 shear moduli G neq,k = 2 Ψ � neq,k . To minimize the number of fitting parameters and their correlations, the relaxation times are spaced equidistantly on a logarithmic scale, viz., with the fitting parameters lg 1 and Δlg( ) . The distribution of the shear moduli is approximated by the normal distribution introducing four parameters a, b, , . Thus, in total, only six fitting parameters are needed instead of 24. Since only a small range of the full relaxation time spectrum plays a role in the given experiments, the function is further reduced by setting Δlg( ) = 1∕2 and = −6 and b = 0 such that (88) S neq = 2 Ψ � neq C −1 ⋅ dev C ⋅ C −1 i (89)

Appendix 6
Fitting Results of the Reference Model for Brain Tissue