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

Abstract

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.

1 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 (bi- or tri-linear interpolations) used to approximate the displacement field over a finite element exhibit convergence issues, 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].

Fig. 1

Schematic views of a typical uniaxial extension of a transversely isotropic solid with unit dimensions: a with the unimodular invariant $\bar{I}_4 = \bar{\lambda }^2$; b with the invariant $I_4 = \lambda ^2$. Qualitatively (a) reveals erroneous transverse deformations along the y- and z-directions, i.e. isotropic directions upon a uniaxial stretch $\lambda$ in the x-direction, which is in contrast to (b)

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 $\lambda$, 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.

Fig. 2

Nonlinear deformation of an anisotropic solid with the reference configuration ${\mathcal {B}}\in \mathbb {R}^3$ and the spatial configuration ${\mathcal {S}}\in \mathbb {R}^3$. The bijective deformation map is ${\varvec{\varphi }}: {\mathcal {B}}\times \mathbb {R}^+ \rightarrow {\mathcal {S}}$, which transforms a material point $\mathbf{X}\in \mathcal {B}$ at time $t_0$ onto a spatial point $\mathbf{x}= {\varvec{\varphi }}(\mathbf{X}, t) \in \mathcal {S}$ at time t, with the deformation gradient $\mathbf{F}$. The anisotropic micro-structure of the material point $\mathbf{X}$ is rendered by two families of fibers with unit vectors $\mathbf{M}$ and $\mathbf{M}^\prime$. Likewise, the anisotropic micro-structure of the spatial point $\mathbf{x}$ is described by $\mathbf{m}$ and $\mathbf{m}^\prime$, as the spatial counterparts of $\mathbf{M}$ and $\mathbf{M}^\prime$, respectively

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 $\mathbf{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, 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–inflation–torsion tests. Finally, Sect. 4 concludes the paper with a brief summary and some remarks.

2 Motion and deformation in an anisotropic continuum

Let ${\mathcal {B}} \subset \mathbb {R}^3$ be a solid body of interest in the reference configuration parametrized by the material point $\mathbf{X}\in {\mathcal {B}}$, while $\partial \mathcal {B}\subset \mathbb {R}^{2}$ denotes the boundary of the reference configuration $\mathcal {B}\subset \mathbb {R}^{3}$. The spatial configuration is denoted by $\mathcal {S} \subset {\mathbb {R}}^3$ with the spatial point $\mathbf{x}\in \mathcal {S}$, and its boundary $\partial \mathcal {S}\subset \mathbb {R}^{2}$. The bijective deformation map ${\varvec{\varphi }}_t(\mathbf{X}) : \mathcal {B}\rightarrow \mathcal {S}$, at time $t \in \mathbb {R}^{+}$, maps a material point $\mathbf{X}$ onto a spatial point $\mathbf{x}$, i.e. ${\varvec{\varphi }}_t(\mathbf{X}) : \mathbf{X}\mapsto \mathbf{x}$, see Fig. 2. Let $\mathbf{E}_I$ and $\mathbf{e}_i$ denote rectangular Cartesian base vectors in the reference and spatial configuration, respectively. The fundamental deformation measure, i.e. the deformation gradient, reads

$$\begin{aligned} \mathbf{F}= \nabla _\mathrm{\mathbf X} {\varvec{\varphi }}_t (\mathbf{X}), \end{aligned}$$

(1)

mapping the unit tangent of a reference point onto its counterpart in the spatial configuration. The deformation gradient $\mathbf{F}$, its cofactor $\mathrm{cof}\mathbf{F}= J\mathbf{F}^\mathrm{-T}$, and its Jacobian $J = \mathrm{det}\mathbf{F}$ relate the deformation of the infinitesimal line ($\mathrm{d}\mathbf{X}$ and $\mathrm{d}\mathbf{x}$), area ($\mathrm{d}\mathbf{A}$ and $\mathrm{d}\mathbf{a}$), and volume ($\mathrm{d}V$ and $\mathrm{d}v$) elements, i.e.

$$\begin{aligned} \mathrm{d}\mathbf{x}= \mathbf{F}\mathrm{d}\mathbf{X}, \qquad \mathrm{d}\mathbf{a}= \mathrm{cof}\mathbf{F}\mathrm{d}\mathbf{A}, \qquad \mathrm{d}v = J\mathrm{d}V. \end{aligned}$$

(2)

The deformations are non-penetrable for $J > 0$. Following Flory [7], one may decouple the deformation gradient $\mathbf{F}= \mathbf{F}_\mathrm{vol}\overline{\mathbf{F}}$ into its dilatational $\mathbf{F}_\mathrm{vol}= J^{1/3}\mathbf{I}$ and unimodular $\overline{\mathbf{F}} = J^{-1/3}\mathbf{F}$ parts. Subsequently, we introduce the right and left Cauchy–Green tensors

$$\begin{aligned} \displaystyle \mathbf{C}= \mathbf{F}^\mathrm{T}\mathbf{F}, \qquad \mathbf{b}= \mathbf{F}\mathbf{F}^\mathrm{T}, \end{aligned}$$

(3)

in the reference and spatial configuration, respectively. Their unimodular parts are defined as $\overline{\mathbf{C}} = J^{-2/3}\mathbf{C}$ and $\overline{\mathbf{b}} = J^{-2/3}\mathbf{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

$$\begin{aligned} \displaystyle I_1 = \mathrm{tr}(\mathbf{b}), \qquad I_2 = \frac{1}{2}\left[ I_1^2 - {\text {tr}}(\mathbf{b}^2) \right] , \qquad I_3 = \det \mathbf{b}\end{aligned}$$

(4)

reflecting the isotropic response of the solid and building up an integrity basis of $\mathbf{F}$, see Spencer [33]. Accordingly, the physically meaningful fourth and sixth invariants read

$$\begin{aligned} \displaystyle I_\mathrm{4} = \mathbf{m}\cdot \mathbf{m}, \qquad I_\mathrm{6} = \mathbf{m}^\prime \cdot \mathbf{m}^\prime , \end{aligned}$$

(5)

which capture the anisotropic response of the solid where $\mathbf{m}= \mathbf{F}\mathbf{M}$ and $\mathbf{m}^\prime = \mathbf{F}\mathbf{M}^\prime$ describe the spatial counterparts of the reference unit vectors $\mathbf{M}$ and $\mathbf{M}^\prime$, as shown in Fig. 2. The corresponding unimodular forms of the aforementioned invariants read $\bar{I}_i = J^{-2/3} I_i$, where $i \in \{1,2,4,6\}$. For the unimodular invariants describing anisotropy, the unimodular spatial vectors $\overline{\mathbf{m}} = \overline{\mathbf{F}} \mathbf{M}$ and $\overline{\mathbf{m}}^\prime = \overline{\mathbf{F}} \mathbf{M}^\prime$ are emphasized.

3 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 $\Psi = \hat{\Psi }(\ldots )$, with function-specific arguments. The Q1P0 and the $Q1P0 {+}{AL }$ basically rely on the same constitutive equations where the multiplicative decomposition of $\mathbf{F}= \mathbf{F}_\mathrm{vol}\overline{\mathbf{F}}$ holds for both isotropic and anisotropic contributions. Hence, we may express the volumetric U(J), the isotropic $\hat{\Psi }_\mathrm{iso}$ and the anisotropic $\hat{\Psi }_\mathrm{ani}$ parts as

$$\begin{aligned} \displaystyle \hat{\Psi }\left( \mathbf{F}, \mathbf{A}_{{\mathbf{{m}}}}, \mathbf{A}_{{\mathbf{{m}}}^\prime } \right) = U(J) + \hat{\Psi }_\mathrm{iso}\left( \overline{\mathbf{F}}\right) + \hat{\Psi }_\mathrm{ani}\left( \overline{\mathbf{F}}, \overline{\mathbf{A}}_{{\mathbf{{m}}}}, \overline{\mathbf{A}}_{{\mathbf{{m}}}^\prime }\right) , \end{aligned}$$

(6)

and the formulation for ${Q1P0 {+}{} WAS }$ reads

$$\begin{aligned} \displaystyle \hat{\Psi }(\mathbf{F}, \mathbf{A}_{{\mathbf{{m}}}}, \mathbf{A}_{{\mathbf{{m}}}^\prime } ) = U(J) + \hat{\Psi }_\mathrm{iso}(\overline{\mathbf{F}}) + \hat{\Psi }_\mathrm{ani}(\mathbf{F}, \mathbf{A}_{{\mathbf{{m}}}}, {\mathbf{A}}_{{\mathbf{{m}}}^\prime }), \end{aligned}$$

(7)

where we have introduced the Eulerian form of the structure tensors $\overline{\mathbf{A}}_{{\mathbf{{m}}}}$, $\overline{\mathbf{A}}_{{\mathbf{{m}}}^\prime }$ and $\mathbf{A}_{{\mathbf{{m}}}}$, $\mathbf{A}_{{\mathbf{{m}}}^\prime }$ as $\overline{\mathbf{A}}_{{\mathbf{{m}}}} = \overline{\mathbf{m}} \otimes \overline{\mathbf{m}}$, $\overline{\mathbf{A}}_{{\mathbf{{m}}}^\prime } = \overline{\mathbf{m}}^\prime \otimes \overline{\mathbf{m}}^\prime$ and $\mathbf{A}_{{\mathbf{{m}}}} = \mathbf{m}\otimes \mathbf{m}$, $\mathbf{A}_{{\mathbf{{m}}}^\prime } = \mathbf{m}^\prime \otimes \mathbf{m}^\prime$, 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

$$\begin{aligned} \displaystyle U(J)= \kappa (J - \text {ln}J - 1), \end{aligned}$$

(8)

whereas the isotropic part $\hat{\Psi }_\mathrm{iso}$ follows from the neo-Hookean material as (see, e.g., Ogden [23])

$$\begin{aligned} \hat{\Psi }_\mathrm{iso}\left( \overline{\mathbf{F}}\right) = \hat{\Psi }_\mathrm{iso}\left( \bar{I}_1\right) = \frac{\mu }{2} \left( \bar{I}_1 - 3\right) . \end{aligned}$$

(9)

In view of (6) the anisotropic part $\hat{\Psi }_\mathrm{ani}$ is a function of the unimodular invariants $\bar{I}_4$ and $\bar{I}_6$ such that

$$\begin{aligned} \hat{\Psi }_\mathrm{ani}(\overline{\mathbf{F}}, \overline{\mathbf{A}}_{{\mathbf{{m}}}}, \overline{\mathbf{A}}_{{\mathbf{{m}}}^\prime })= & {} \hat{\Psi }_\mathrm{ani}(\bar{I}_4, \bar{I}_6)\nonumber \\= & {} \frac{ { \bar{k}_1 }}{2 { \bar{k}_2 } } \sum _{i = 4,6} \left\{ \text {exp}\left[ { \bar{k}_2} \left( \bar{I}_i - 1\right) ^2\right] - 1\right\} ,\nonumber \\ \end{aligned}$$

(10)

as suggested by Holzapfel et al. [16]. However, for (7) $\hat{\Psi }_\mathrm{ani}$ becomes a function of the invariants $I_4$ and $I_6$, i.e.

$$\begin{aligned} \hat{\Psi }_\mathrm{ani}(\mathbf{F}, \mathbf{A}_{{\mathbf{{m}}}}, \mathbf{A}_{{\mathbf{{m}}}^\prime })= & {} \hat{\Psi }_\mathrm{ani}( I_4, I_6)\nonumber \\= & {} \frac{k_1}{2k_2} \sum _{i = 4,6} \left\{ \text {exp}\left[ k_2 \left( I_i - 1\right) ^2\right] - 1\right\} .\nonumber \\ \end{aligned}$$

(11)

In (8) $\kappa$ denotes the bulk modulus, whereas $\mu$ refers to the shear modulus in (9). The anisotropic material parameters $\bar{k}_1$, $\bar{k}_2$ in (10) and $k_1$, $k_2$ in (11) represent a stress-like material parameter and a dimensionless parameter, respectively. In general, $\bar{k}_1$, $\bar{k}_2$ and $k_1$, $k_2$ are different from one another.

3.1 Stress expressions

From the Coleman–Noll procedure we arrive at the definition of the Kirchhoff stress tensor ${\varvec{\tau }}$ according to

$$\begin{aligned} \displaystyle {\varvec{\tau }}= & {} { 2 \mathbf{F}(\partial _{{\mathbf{{C}}}} \hat{\Psi }) \mathbf{F}^\mathrm{T}} = {\varvec{\tau }}_\mathrm{vol} + {\varvec{\tau }}_\mathrm{iso} + {\varvec{\tau }}_\mathrm{ani}, \end{aligned}$$

(12)

where ${\varvec{\tau }}_\mathrm{vol}$, ${\varvec{\tau }}_\mathrm{iso}$, ${\varvec{\tau }}_\mathrm{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

(13)

Table 1 Elastic material parameters used for the analysis of an anisotropic material for Sects. 4.1–4.3

where $\hat{p} = J \mathrm{d} U/\mathrm{d} J = \kappa (J-1)$, and the projection tensor is defined as in which the symmetric fourth-order identity tensor $\mathbb {I}$ has the index form $(\mathbb {I})_{ij}^{kl} = \frac{1}{2}(\delta _i^k \delta _j^l + \delta _i^l \delta _j^k)$. 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 ${\varvec{\tau }}_\mathrm{ani}$ reads

(14)

where the deformation-dependent scalar coefficients are defined as $\bar{\psi }_4 = \partial _{\bar{I}_4} \hat{\Psi }_\mathrm{ani}(\bar{I}_4, \bar{I}_6)$ and $\bar{\psi }_6 = \partial _{\bar{I}_6} \hat{\Psi }_\mathrm{ani}(\bar{I}_4, \bar{I}_6)$. The counterpart of (14) for ${Q1P0 {+}{} WAS }$ is given by

$$\begin{aligned} \displaystyle {\varvec{\tau }}_\mathrm{ani} = { 2 \mathbf{F}[\partial _{{\mathbf{{C}}}} \hat{\Psi }_\mathrm{ani}(I_4, I_6)] \mathbf{F}^\mathrm{T}} = 2 \psi _4 \mathbf{m}\otimes \mathbf{m}+ 2\psi _6 \mathbf{m}^\prime \otimes \mathbf{m}^\prime , \end{aligned}$$

(15)

where the coefficients in this case read $\psi _4 = \partial _{I_4} \hat{\Psi }_\mathrm{ani}(I_4, I_6)$ and $\psi _6 = \partial _{I_6} \hat{\Psi }_\mathrm{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].

4 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.

4.1 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 $\mathbf{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^\mathrm{avg}$ (averaged over the 8 finite elements; subsequently the index $(\bullet )^\mathrm{avg}$ stands for average) with respect to stretch $\lambda _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^\mathrm{avg}$ with $\hat{\Psi }_\mathrm{ani}^\mathrm{avg}$. In fact, high values of $U^\mathrm{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.

Fig. 3

Uniaxial extension of a transversely isotropic material with a single family of fibers $\mathbf{M}$ aligned in the x-direction, the direction of the displacement-driven loading. The plots show the Jacobian $J^\mathrm{avg}$ averaged over 8 unstructured hexahedral elements, the average volumetric free energy $U^\mathrm{avg}$ (in kPa) and the average anisotropic free energy $\hat{\Psi }_\mathrm{ani}^\mathrm{avg}$ (in kPa) as a function of stretch $\lambda _x$. Results are obtained on the basis of: Q1P0 (solid curve with empty triangles), $Q1P0 {+}{} AL$ (solid curve with empty circles), and the approach avoiding the multiplicative decomposition of $\mathbf{F}$ for the anisotropic contribution $Q1P0 {+}{} WAS$ (red dashed curve). The comparison is shown for two different sets of material parameters, namely for case a referring to a and for case b referring to b, compare with Table 1. (Color figure online)

Fig. 4

Uniaxial extension of a transversely isotropic material with a single family of fibers $\mathbf{M}$ aligned in the x-direction, with displacement-driven loading applied in the y-direction. The plots show the Jacobian $J^\mathrm{avg}$ averaged over 8 unstructured hexahedral elements, the average volumetric free energy $U^\mathrm{avg}$ (in kPa) and the average isotropic free energy $\hat{\Psi }_\mathrm{iso}^\mathrm{avg}$ (in kPa) as functions of stretch $\lambda _y$. Results are obtained on the basis of: Q1P0 (solid curve with empty triangles), ${Q1P0 {+}{} AL }$ (solid curve with empty circles), and the approach avoiding the multiplicative decomposition of $\mathbf{F}$ for the anisotropic contribution $Q1P0 {+}{} WAS$ (red dashed curve). The comparison is shown for two different sets of material parameters, namely for case a referring to a and for case c referring to c, compare with Table 1. (Color figure online)

4.2 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 y-direction, 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^\mathrm{avg}$, the volumetric free energy $U^\mathrm{avg}$ and the isotropic free energy $\hat{\Psi }_\mathrm{iso}^\mathrm{avg}$, 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 $\lambda _y$. All treatments are able to yield, to a large extent, the relevant isotropic response.

4.3 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 $\mathbf{M}$ and $\mathbf{M}^\prime$, where $40^\circ$ is the angle between the fibers and the circumferential axis $\theta$, 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 $\hat{u}_z = 2$ mm; an inner pressure growing up to $\hat{p}_\mathrm{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 $\hat{\gamma } = 60^\circ$, 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 $\tilde{u}_x=\tilde{u}_y=\tilde{u}_z=0$. For the purpose of comparison we only examine case b, see Table 1.

Fig. 5

a Cylindrical tube with a hypothetical tissue and geometry with dimensions $H = 10$, $R = 8$, and $T = 2$ mm with related boundary and loading conditions $\tilde{u}_x = \tilde{u}_y = \tilde{u}_z = 0$ mm, $\hat{u}_z = 2$ mm, $\hat{p}_\mathrm{i} = 500$ mmHg and $\hat{\gamma } = 60^\circ$; b finite element mesh of the corresponding geometry; c visual representation of the first family of fibers $\mathbf{M}$ at each node; d visual representation of the second family of fibers $\mathbf{M}^\prime$ at each node

Fig. 6

Distributions of the radial $\sigma _{rr}$, the circumferential $\sigma _{\theta \theta }$ and the axial $\sigma _{zz}$ Cauchy stress components for aQ1P0; b$Q1P0 {+}{} WAS$; c$Q1P0 {+}{AL}$. (Color figure online)

Figure 6 reveals dramatic differences in the Cauchy stress responses illustrated along the cylindrical coordinates r, $\theta$ and z. The Q1P0 element formulation generates much lower stress components ($\sigma _{rr}$, $\sigma _{\theta \theta }$, $\sigma _{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 ${\texttt {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.

Table 2 Quadratic rate of convergence of the global residual for $Q1P0 {+}{} WAS$, Sect. 4.3

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.

5 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 $\mathbf{F}$ for the anisotropic part of the constitutive equations. The results corroborate the new concept namely $Q1P0 {+}{} WAS$.

Several anisotropic constitutive models presume quasi-incompressibility. Hence, there should be no discrepancy between the theoretical structure and the related numerical application. Having scrutinized the results in Figs. 3, 4, and 6 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 $\bar{I}_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.