On the quasi-incompressible finite element analysis of anisotropic hyperelastic materials

Quasi-incompressible behavior is a desired feature in several constitutive models within the finite elasticity of solids, such as rubber-like materials and some fiber-reinforced soft biological tissues. The Q1P0 finite element formulation, derived from the three-field Hu–Washizu variational principle, has hitherto been exploited along with the augmented Lagrangian method to enforce incompressibility. This formulation typically uses the unimodular deformation gradient. However, contributions by Sansour (Eur J Mech A Solids 27:28–39, 2007) and Helfenstein et al. (Int J Solids Struct 47:2056–2061, 2010) conspicuously demonstrate an alternative concept for analyzing fiber reinforced solids, namely the use of the (unsplit) deformation gradient for the anisotropic contribution, and these authors elaborate on their proposals with analytical evidence. The present study handles the alternative concept from a purely numerical point of view, and addresses systematic comparisons with respect to the classical treatment of the Q1P0 element and its coalescence with the augmented Lagrangian method by means of representative numerical examples. The results corroborate the new concept, show its numerical efficiency and reveal a direct physical interpretation of the fiber stretches.


Introduction
That the constitutive models for finite elasticity of rubber-like materials and some fiber-reinforced soft biological tissues suffer from ill-conditioning of the global stiffness matrix, referred to as the locking phenomenon, is a well-known issue. Locking mainly arises when the standard displacement formulations are used, but is not directly related to a physical response, such as bending and quasi-incompressible elasticity. In fact, for such problems first-order shape functions (bior tri-linear interpolations) used to approximate the displacement field over a finite element exhibit convergence issues, B Gerhard A. Holzapfel holzapfel@tugraz.at 1 see, e.g., Hughes [19], Zienkiewicz and Taylor [40] and Wriggers [38]. One way to avoid this problem is to use a mixed variational formulation that hinges on the Hu-Washizu principle, where the master field appears together with additional subsidiary conditions, as first introduced by Nagtegaal et al. [21] and discussed in Brezzi and Fortin [4] for small strains. It was later extended to finite strain problems by Simo et al. [32] and is also well documented in the literature, by, e.g., Miehe [20] and Wriggers [38]. This approach falls into the category of the three-field Hu-Washizu principle with the constitutive function augmented by a volumetric constraint via the scalar conjugate pairs, i.e. pressure-dilatation. One of the relevant finite element formulations, known as the Q1P0 element, is robust for most problems in solid mechanics; however, it may cause numerical instabilities for the pressure for some specific loading and boundary conditions, see Wriggers [38].
Other methods used to circumvent the locking phenomenon include h-and p-refinement strategies (Düster et al. [5]), reduced integration schemes with stabilization techniques (Belytschko et al. [2,3]; Reese [24]) and enhanced strain formulations based on the Hu-Washizu principle (Simo and coworkers [28][29][30]). The above-stated methods, albeit effective in averting the locking issue, are connected with other problems, e.g., increased computational costs, artificial stabilizing parameters and numerical instabilities upon irregular mesh distortion. Therefore, they fall into disfavor among, e.g., the biomechanics community when a mere finite hyperelastic analysis of the material is sought. In addition, the more recent revelations by Schröder et al. [27] and Wriggers et al. [39] highlight the mixed variational principles for the treatment of the inextensibility limit in fiber-reinforced materials and soft biological tissues, which are useful in the presence of extremely stiff fibers, i.e. very high stiffness associated with fibers may lead to low convergence rates in the primary fields, e.g., the displacements. Since its introduction by Flory [7] the multiplicative decomposition of the deformation gradient has gained popularity in the context of the variational Hu-Washizu principle. This split generates nonphysical results in the case of simple tension and compression problems for isotropic hyperelastic models capturing rubber-like materials, which are known to be quasi-incompressible; for more details see Ehlers and Eipper [6]. By the same token, nonphysical responses may occur in the numerical analysis of anisotropic materials unless the analysis is based on a quasi-incompressible material formulation. In this respect, the multiplicative split of the deformation gradient leads to a twofold issue for fiber-reinforced hyperelastic materials: (i) Stresses in the fibers are expected to be one-dimensional since fibers are assumed to behave like one-dimensional springs, see Sansour [26]. However, the use of isochoric anisotropic invariants automatically yields a projection tensor that generates the stress components perpendicular to the alignment of the fibers, thereby violating the basic assumption. (ii) For the anisotropic contribution the use of isochoric anisotropic invariants leads to a 'competition' between the anisotropic part and the volumetric part of the free energies in the process of energy minimization. As a matter of fact, if the solid undergoes volumetric deformations, a much lower strain energy is stored in the system in comparison with that resulting from the fibers undergoing deformation for the same amount of global stretch, say λ, as illustrated in Fig. 1. As a consequence, the system favors the volumetric part and tends to generate spurious spherical deformations accompanied by a volume growth at relatively small stretches. Such a disparity is discernable in a typical numerical uniaxial extension test.
Several studies, e.g., Helfenstein et al. [14], Annaidh et al. [1] and Nolan et al. [22], have reported the erroneous analysis results of fiber-reinforced anisotropic material models for soft biological tissues (Weiss et al. [37], Holzapfel et al. [16] and Rubin and Bodner [25]) when they are mistakenly used in the compressible domain; e.g., a sphere reinforced with one family of fibers would be deformed into a sphere with a smaller size upon hydrostatic pressure instead of taking on an ellipsoidal shape. One remedy for (ii) is to implement the computationally (rather) expensive augmented Lagrangian method to bring the analysis towards the incompressibility limit, see Glowinski and Le Tallec [8,9] and Simo and Taylor [31] among others. Alternatively, the multiplicative decomposition of the deformation gradient can be avoided for the anisotropic part, as suggested by Sansour [26] and Helfenstein et al. [14], which solves the issue on the constitutive level without using any ad hoc algorithm.
The emphasis of the present article is placed upon the comparison of two remedies namely the augmented Lagrangian method and the use of an unsplit deformation gradient F for the anisotropic contribution; the consequences thereof are elucidated with simple examples in the context of the Q1P0 finite element formulation. In the authors' opinion, Nonlinear deformation of an anisotropic solid with the reference configuration B ∈ R 3 and the spatial configuration S ∈ R 3 . The bijective deformation map is ϕ : B × R + → S, which transforms a material point X ∈ B at time t 0 onto a spatial point x = ϕ(X, t) ∈ S at time t, with the deformation gradient F. The anisotropic micro-structure of the material point X is rendered by two families of fibers with unit vectors M and M . Likewise, the anisotropic micro-structure of the spatial point x is described by m and m , as the spatial counterparts of M and M , respectively a systematic comparison of the above-mentioned concepts is particularly relevant for highlighting the issues within the biomechanics community. We further emphasize that for the sake of brevity the classical Q1P0 element, its coalescence with the augmented Lagrangian method and the Q1P0 element without the multiplicative split in the anisotropic contribution are hereinafter denoted by Q1P0, Q1P0+AL (Q1P0 element along with Augmented Lagrangian method), and Q1P0+WAS (Q1P0 element Without Anisotropic Split), respectively. Section 2 summarizes in brief the background on the constitutive modeling of fibrous (soft) tissues where the collagen fibers are embedded in an otherwise isotropic matrix material. Subsequently, Sect. 3 documents simple yet representative boundary-value problems which demonstrate volume changes, isotropic and anisotropic energy contributions, and the associated Cauchy stresses for the considered material under uniaxial extension and extension-inflationtorsion tests. Finally, Sect. 4 concludes the paper with a brief summary and some remarks.

Motion and deformation in an anisotropic continuum
Let B ⊂ R 3 be a solid body of interest in the reference configuration parametrized by the material point X ∈ B, while ∂B ⊂ R 2 denotes the boundary of the reference configuration B ⊂ R 3 . The spatial configuration is denoted by S ⊂ R 3 with the spatial point x ∈ S, and its boundary ∂S ⊂ R 2 . The bijective deformation map ϕ t (X) : B → S, at time t ∈ R + , maps a material point X onto a spatial point x, i.e. Fig. 2. Let E I and e i denote rectangular Cartesian base vectors in the reference and spatial configuration, respectively. The fundamental deformation measure, i.e. the deformation gradient, reads mapping the unit tangent of a reference point onto its counterpart in the spatial configuration. The deformation gradient F, its cofactor cofF = J F −T , and its Jacobian J = detF relate the deformation of the infinitesimal line (dX and dx), area (dA and da), and volume (dV and dv) elements, i.e.
The deformations are non-penetrable for J > 0. Following Flory [7], one may decouple the deformation gradient F = F vol F into its dilatational F vol = J 1/3 I and unimodular F = J −1/3 F parts. Subsequently, we introduce the right and left Cauchy-Green tensors in the reference and spatial configuration, respectively. Their unimodular parts are defined as C = J −2/3 C and b = J −2/3 b. For an elaborate treatment see, e.g., Truesdell and Noll [36], Holzapfel [15] and Gurtin et al. [13]. Next, we exploit the representation theorem of invariants according to Spencer [34] and Holzapfel [15], and define the three irreducible isotropic invariants reflecting the isotropic response of the solid and building up an integrity basis of F, see Spencer [33]. Accordingly, the physically meaningful fourth and sixth invariants read

A particular form of the model by Holzapfel et al. [16]
In the following we assume the existence of a Helmholtz free-energy function, say =ˆ (. . .), with function-specific arguments. The Q1P0 and the Q1P0+AL basically rely on the same constitutive equations where the multiplicative decomposition of F = F vol F holds for both isotropic and anisotropic contributions. Hence, we may express the volumetric U (J ), the isotropicˆ iso and the anisotropicˆ ani parts aŝ and the formulation for Q1P0+WAS readŝ where we have introduced the Eulerian form of the structure tensors A m , A m and A m , A m as A m = m ⊗ m, A m = m ⊗m and A m = m⊗m, A m = m ⊗m , respectively. Note that in (7) the multiplicative decomposition is solely used upon the isotropic response, while it is completely excluded from the fiber response. The volumetric part can simply be chosen as whereas the isotropic partˆ iso follows from the neo-Hookean material as (see, e.g., Ogden [23]) In view of (6) the anisotropic partˆ ani is a function of the unimodular invariantsĪ 4 andĪ 6 such thatˆ as suggested by Holzapfel et al. [16]. However, for (7)ˆ ani becomes a function of the invariants I 4 and I 6 , i.e.
In (8) κ denotes the bulk modulus, whereas μ refers to the shear modulus in (9). The anisotropic material parametersk 1 , k 2 in (10) and k 1 , k 2 in (11) represent a stress-like material parameter and a dimensionless parameter, respectively. In general,k 1 ,k 2 and k 1 , k 2 are different from one another.

Stress expressions
From the Coleman-Noll procedure we arrive at the definition of the Kirchhoff stress tensor τ according to where τ vol , τ iso , τ ani represent the volumetric, isotropic and anisotropic terms associated with the Kirchhoff stress tensor.
With the chain rule the volumetric and isotropic parts yield the forms , and the projection tensor is defined as P = I − 1 3 (I ⊗ I) in which the symmetric fourth-order identity tensor I has the index form . It needs to be emphasized that a mean dilation approach in the context of the Q1P0 formulation results in an averaged uniform pressure field, and a dilatation field over a finite element. Therefore, all volumetric responses essentially emerge on an element level instead of a constitutive level. It is also worth mentioning that all three formulations (Q1P0, Q1P0+AL, Q1P0+WAS) assume the same expressions for their volumetric and isotropic stress responses, as stated in (13). The mere difference originates from their distinct anisotropic constitutive forms. In fact, in case of Q1P0 and Q1P0+AL the anisotropic Kirchhoff stress tensor τ ani reads Case a κ = 5000 kPa μ = 10 kPak 1 = k 1 = 50 kPak 2 = k 2 = 2.0 Case b κ = 5000 kPa μ = 10 kPak 1 = k 1 = 500 kPak 2 = k 2 = 2.0 Case c κ = 5000 kPa μ = 100 kPak 1 = k 1 = 50 kPak 2 = k 2 = 2.0 where the deformation-dependent scalar coefficients are defined asψ 4 = ∂Ī 4ˆ ani (Ī 4 ,Ī 6 ) andψ 6 = ∂Ī 6ˆ ani (Ī 4 ,Ī 6 ). The counterpart of (14) for Q1P0+WAS is given by where the coefficients in this case read ψ 4 = ∂ I 4ˆ ani (I 4 , I 6 ) and ψ 6 = ∂ I 6ˆ ani (I 4 , I 6 ). The solution of nonlinear problems necessitates the consideration of spatial elasticity tensors which can be readily derived, see, e.g., [12,15].

Representative numerical examples
We hereby touch upon the quasi-incompressible hyperelastic performances of the Q1P0, the Q1P0+AL, and the Q1P0+WAS formulations. Comparisons regarding the uniaxial extension tests are analyzed by considering three different sets of material parameters, as summarized in Table 1. The cases b and c provide stiffer constitutive responses than case a in the sense of anisotropy and isotropy, respectively. The final example presents a thick-walled cylindrical tube idealizing a single-layered hypothetical arterial tissue, which is extended, inflated and twisted simultaneously.

Numerical investigations of Q1P0, Q1P0+AL, and Q1P0+WAS along a fiber direction
In order to scrutinize the constitutive responses associated with Q1P0, Q1P0+AL, and Q1P0+WAS, a simple unit cube discretized with 8 unstructured hexahedral elements is adopted, see Fig. 3. The isotropic ground matrix is reinforced by a single family of fibers with orientation M aligned in the x-direction, which is also the loading direction. Hence, the unit cube undergoes uniaxial extension. In order to better discuss the effect on the volumetric response, comparisons are analyzed by using two different sets of anisotropic material parameters (case a and case b), as provided in Table 1. The corresponding results are depicted in Fig. 3a, b, respectively. As far as the Q1P0 study is concerned, both Fig. 3a, b indicate a tremendous increase of the Jacobian J avg (averaged over the 8 finite elements; subsequently the index (•) avg stands for average) with respect to stretch λ x for both sets of material parameters-solid curves with empty triangles in Fig. 3. It is evident that the anisotropic split creates a kind of 'competition' between the volumetric part and the anisotropic part through the minimization principle, thereby favoring an upsurge in the volumetric free energy, while diminishing the anisotropic contribution, which can be grasped by comparing U avg withˆ avg ani . In fact, high values of U avg are essentially the result of the minimization principle preferring an increase in the volumetric response over the anisotropic one. This predicament can be overcome by applying either Q1P0+AL (solid curves with empty circles in Fig. 3) or Q1P0+WAS (red dashed curves in Fig. 3). Note that for Q1P0+AL we prescribe a maximum number of 20 augmented Lagrangian iterations for each incremental step which can also be increased in order to achieve a better performance for the case b.

Numerical investigations of Q1P0, Q1P0+AL, and Q1P0+WAS along an isotropic direction
This benchmark is identical to that described in Sect. 4.1 in regard to geometry, structure and mesh data, see Fig. 4. Yet, the uniaxial extension is in this case applied in the ydirection, hence in an isotropic direction. As the comparison of anisotropic responses become trivial due to the setup of the problem, different isotropic responses are probed in line with two distinct sets of material parameters, cases a and c in Table 1. The corresponding results are depicted in Fig. 4a, c, respectively. In all two cases Q1P0 (solid curves with empty triangles) and Q1P0+WAS (red dashed curves) yield identical results in regard to the average values of the Jacobian J avg , the volumetric free energy U avg and the isotropic free energyˆ avg iso , as expected. The other treatment Q1P0+AL, however, does not create any growth in volume, thereby providing the most physical response in a rigorous sense, see Fig. 4. Nonetheless, the volume growth associated with Q1P0 and Q1P0+WAS is not pronounced and is minimal when compared with the enormous swelling for Q1P0 in the previous example, see Fig. 3. As a result, the differences in the isotropic response are practically unnoticeable even at relatively large stretches λ y . All treatments are able to yield, to a large extent, the relevant isotropic response.

Extension-inflation-torsion test for Q1P0, Q1P0+AL, and Q1P0+WAS
The aim of this benchmark test is to compare the formulations Q1P0, Q1P0+AL, and Q1P0+WAS in regard to how much the mechanical responses differ from each other under extreme loading conditions. To this end we consider a cylindrical tube as a single-layered hypothetical arterial tissue with the geometry provided in Fig. 5a. The morphology renders anisotropy via two symmetric families of fibers M and M , where 40 • is the angle between the fibers and the circumferential axis θ , see Fig. 5c, d. The domain is discretized with 960 solid brick elements connected by 1320 nodes, see Fig. 5b. As for loading, monotonically increasing displacements are applied on the top plane in the z-direction up toû z = 2 mm; an inner pressure growing up top i = 500 mmHg is monotonically exerted on the inner layer of the tissue, while a monotonically increasing torsion is applied on the top plane reaching an angle of twistγ = 60 • , as depicted in Fig. 5a. Figure 5a also illustrates the boundary conditions, i.e. the displacements on the bottom plane are constrained along all three directionũ x =ũ y =ũ z = 0. For the purpose of comparison we only examine case b, see Table 1. Figure 6 reveals dramatic differences in the Cauchy stress responses illustrated along the cylindrical coordinates r , θ and z. The Q1P0 element formulation generates much lower stress components (σ rr , σ θθ , σ zz ) since the anisotropic contribution, likewise in Fig. 3, is disfavored during the energy minimization process, see Fig. 6a, which leads to relatively low stress components. This example ultimately pinpoints how predicted stress values become spurious under a supra-physiological loading scenario an arterial wall might undergo when a classical Q1P0 finite element formulation is chosen. Q1P0+WAS and Q1P0+AL yield very similar responses for this numerical example, with  differences less than 5%, compare Fig. 6b with Fig. 6c. It should also be highlighted that the maximum number of augmented Lagrangian iterations is set to 50 with the augmenting factor augf = 5 (see in FEAP [35]) in order to enforce incompressibility for Q1P0+AL in a strict sense. The augmented Lagrangian method in its turn leads to a major drawback, i.e. the computational time required to simulate the Q1P0+AL is around 20 minutes, while only 1 minute is needed for the Q1P0+WAS, rendering Q1P0+WAS nearly 20 times faster for the problem considered.
In addition, we briefly report on the computational performance of Q1P0+WAS; the quadratic rate of convergence behavior at t = 0.25, 0.5, 0.75, and 1 s is summarized in Table 2. All simulations are carried out on a single 3.2 GHz Intel ® Core TM i7-3930K CPU on a 64bit Linux operating system with 32 GB RAM.

Summary and concluding remarks
Some important aspects in regard to physically relevant and computationally efficient analyses of fiber-reinforced materials were presented. Following an extensive literature overview, a brief theoretical synopsis of anisotropic hyperelasticity was provided. Subsequently, the numerical performance of the classical Q1P0 element with and without the augmented Lagrangian method was examined together with a rather new concept introduced by Sansour [26] and Helfenstein et al. [14], who proposed the use of the (unsplit) deformation gradient tensor F for the anisotropic part of the constitutive equations. The results corroborate the new concept namely Q1P0+WAS.
Several anisotropic constitutive models presume quasiincompressibility. Hence, there should be no discrepancy between the theoretical structure and the related numeri- it is palpably shown that the concept Q1P0+WAS is able to generate quasi-incompressible responses which are in good conformity with those of Q1P0+AL even under extreme loading cases, as in Fig. 6. The feature through which Q1P0+WAS gains the upper hand over Q1P0+AL is the numerical efficiency. In fact, given the size of the mesh domain in Sect. 4.3, Q1P0+WAS appears to be approximately 20 times as fast as Q1P0+AL. It should be underlined that we also performed numerical studies with Q1P0+WAS using I 1 instead ofĪ 1 for the isotropic part. For this case almost the same stress results as in Fig. 6 were obtained with the quadratic rate of convergence behavior retained. Another credible aspect of Q1P0+WAS is that the formulation leads to a physical interpretation of the fiber stretches through I 4 and I 6 . This allows exclusion of compressed fibers, a significant issue which may cause erroneous model considerations, as pointed out by Holzapfel and Ogden [17,18]. At this point, we highlight that no systematic experimental data are yet available to substantiate or rebut incompressibility/compressibility of rupturing tissues. Moreover, enforcing   incompressibility via an expedient method such as the augmented Lagrangian method often generates convergence issues with such highly inelastic constitutive responses. For recent investigations of the rupture behavior of aortic tissues we refer to Gültekin et al. [10][11][12] and references therein. For soft biological tissues, it is of utmost importance to carry out fast yet reliable numerical analyses which enable a precise validation of experimental data obtained in vivo or ex vivo. This will better inform and guide medical monitoring and rupture risk assessment of diseases such as aortic dissection, atherosclerosis and aneurysms.