Physical response of hyperelastic models for composite materials and soft tissues

Abstract

A hyperelastic model must not only characterize the mechanical response of a composite material such as soft tissue, but also ensure numerical stability by a feasible set of material parameters. Apart from the well-known ill-conditioning problem caused by the incompressibility constraint, the paper indicates another ill-conditioning occurring in any general fibre-reinforced material model for tubular organs when unbalance between the fibre strain energy and the matrix strain energy becomes too large. Specifically, although the Holzapfel model is polyconvex, this problem can be observed as an unphysical behaviour in a physiological deformation range of a tissue such as arterial wall and intestine by thickening in the thickness direction associated with a volume growth of a specimen in a tension test. Particularly, the same problem for a polyconvex modified Fung-type model with the matrix characterized by the neo-Hookean model has been discussed for the first time. By investigating the influence of the shear modulus in these two models, we not only prove the cause of the ill-conditioning but also propose a solution to control the unbalance in the strain energy. The numerical results show significant enhancement of the model stability in overcoming the unphysical deformation.

Background

Composite materials such as biological soft tissues can be modelled using a phenomenological approach such as an exponential function by Fung et al. [1] and later fully completed in a more general form by Humphrey [2]. This classical Fung-type model has been widely used in biomechanics such as rabbit skin (Tong and Fung [3]), rabbit arteries (Chuong and Fung [4]) and a right coronary artery (Takamizawa [5]). On the contrary, a micro/structural approach models the microstructural constituents of tissues such as modelling fibres of connective tissue as long sinusoidal beams by Comninou and Yannas [6] and flat collagenous fibres (Lanir [7]). Particularly, Holzapfel et al. [8] postulated a strain–energy function (called the Holzapfel model) for arteries based on microstructures of the tissues by considering each layer of the artery as a composite material whose non-collagenous tissues, e.g. elastin, are considered as isotropic matrix substance. The collagen, which is assumed to be composed of two families of collagenous fibres, plays the role of a reinforcing material, and is assumed to be perfectly embedded in the isotropic matrix. The fibres and the matrix are characterized by different strain–energy functions and the material coefficients have clear physical meaning.

Besides the capability of representing the material behaviour of biological soft tissue accurately, constitutive equations must also ensure numerical stability in computer simulation by satisfying the Legendre–Hadamard or the ellipticity condition equivalent to a convexity or polyconvexity condition (Marsden et al. [9]).

Although the Holzapfel model is polyconvex (Schröder et al. [10]), a numerical instability known as unphysical behaviour can occur for an improper set of material coefficients as described in tension tests by Gasser et al. [11]. Considering the tests, the structural arrangement of the stiffened fibres in combination with the soft ground matrix activates a load-carrying mechanism in which the collagen fibres need to rotate towards the loading direction until they are able to carry significant load. Consequently, there is an unrealistic thickening of the sample strip in tension tests and it violates the incompressibility condition with a volume growth. To prevent this non-physical behaviour, Gasser et al. [11] took into account an effect of the dispersion of the collagen fibres whose orientations are not uniform. Admittedly, this proposal helps the model come closer to the real behaviour. In fact, by investigating terms of the strain–energy function, it is obvious that the incorporated dispersion is able to control the non-physical thickening by decreasing the difference between the orthotropic energy and the isotropic energy.

A set of material constants, which is identified by curve fitting and satisfies the goodness of fit, is not unique in general. For example, by fixing the shear modulus, Helfenstein et al. [12] identified two sets of material constants of the Holzapfel model with and without a volumetric–isochoric split. This split was also known as another source of the non-physical response (Ehlers and Eipper [13]). Shear moduli differing by a factor of 10 can be observed for two different arterial layers (Holzapfel et al. [8]). Obviously, a too soft ground matrix characterized by the small shear modulus of the neo-Hookean model may cause the constitutive matrix of the overall model to be ill-conditioned because the entries of this matrix are summations of both the small stiffness of the ground matrix and the large stiffness of the fibres, leading to a high condition number. Thus, this composite material model shows strongly directional behaviour and would behave non-physically in a direction along which the strain–energy dissipation is minimal, e.g. in directions orthogonal to the fibre orientations (Duong et al. [14]). In tubular organs this is observed in the radial direction.

In this paper, we show that the above thickening effect is concerned with the ill-conditioning problem and no recourse to the “fibre rotation”, as given by Gasser et al. [11], is needed. Observations of the ill-conditioning for nonlinear elastic models for certain classes of soft biological tissue have been mentioned earlier, see, e.g. Sun and Sacks [15]. Besides, the formulation of biological models is not similar, and hence mechanisms causing the ill-conditioning of each model are different. Herein, we would prove that the unrealistic phenomenon of the Holzapfel model is resulting from the ill-conditioned constitutive matrix caused by large differences of the magnitudes of the shear terms in the elasticity tensor. Moreover, the modified Fung-type model, abbreviated as MFH model (Holzapfel [16]), has been proven to be affected by this ill-conditioning problem in the physiological deformation range for the first time (Duong [17]). However, no indication for relation between this numerical instability and the mentioned “fibre rotation” or the source of the thickening effect has been provided in literature [11–13]. In addition, no suggestion for treating the ill-conditioning issue as well as for obtaining a lower condition number has been proposed because solving the condition number equation through matrix manipulation leads to a very complicated formulation. Thus, in this paper the numerical condition number of the elasticity tensor (hereafter it is shortly named the condition number), whose value changes as a function of material parameters and deformation, is also investigated. Numerical results show the significant influence of the shear moduli in the two considered models on the unphysical behaviour. Consequently, this suggests how to treat the unphysical mechanical response of the material law by controlling the unbalance of the energy with the highest possible shear modulus, leading to significant improvements in the numerical results.

Constitutive Equations

Models for soft biological tissues

For the sake of comparison and investigation of the ill-conditioning, in this section two well-known material models as well as their modifications are recalled; the Fung model [1] and the Holzapfel model [8].

General Fung model

The strain–energy function of the classical Fung model is

$$W_{\text{F}} ({\mathbf{E}},J) = C[\exp (c_{1} E_{1}^{2} + c_{2} E_{2}^{2} + c_{3} E_{3}^{2} + 2c_{4} E_{1} E_{2} + 2c_{5} E_{2} E_{3} + 2c_{6} E_{3} E_{1} ) - 1] + 2U(J)$$

(1)

where $C,c_{i} \,(i = 1, \ldots ,6)$ are material parameters. C has the dimension of a modulus; all c _i are dimensionless. $E_{1,2,3}$ are the principal Green–Lagrange strains, the volumetric part $U(J) = \frac{1}{2}p(J - 1)$; $J = \det \,{\mathbf{F}}$ where F is the deformation gradient, for incompressible materials J = 1. In the finite element analysis p plays the role of a penalty factor which is typically expressed by the bulk modulus $k = \frac{2(1 + v)}{3(1 - 2v)}\mu$, with $\mu$ equivalent to the small strain shear modulus and $v = 0.5$ expresses the incompressibility similarly to the small strain limit.

Holzapfel model

The Fung model does not need any knowledge about the fibre orientation in the tissue. In contrast, in the Holzapfel model the isotropic part $W_{\text{iso}}$ presents the non-collagenous tissues, and the anisotropic part $W_{\text{ani}}$ describes the collagen fibres under an angle $\alpha$ typically with respect to the hoop direction of a hollow organ, see Fig. 1:

$$W_{\text{H}} (I_{1} ,I_{4} ,I_{6} ,J) = W_{\text{iso}} (I_{1} ) + W_{\text{ani}} (I_{4} ,I_{6} ) + U(J)$$

(2)

The isotropic term is the neo-Hookean model $W_{\text{iso}} (I_{1} ) = \frac{\mu }{2}(I_{1} - 3)$ and the anisotropic term is

$$W_{\text{ani}} (I_{4} ,I_{6} ) = \frac{{k_{1H} }}{{2k_{2H} }}\left\{ {\exp [k_{2H} (I_{4} - 1)^{2} ] - 1} \right\} + \frac{{k_{1H} }}{{2k_{2H} }}\left\{ {\exp [k_{2H} (I_{6} - 1)^{2} ] - 1} \right\}$$

(3)

with $k_{1H}$ a material parameter, having the dimension of a modulus, $k_{2H}$ a dimensionless material parameter; invariants of the right Cauchy–Green tensor ${\mathbf{C}}:I_{1} = {\text{tr}}{\mathbf{C}}$, $I_{4} = {\varvec{a}}_{0} \cdot {\mathbf{C}} {\varvec{a}}_{0}$, $I_{6} = \varvec{g}_{0} \cdot {\mathbf{C}} \varvec{g}_{0}$, $\varvec{a}_{0}$ and $\varvec{g}_{0}$ unit direction vectors of the two fibre families: $\varvec{a}_{0} = [0\,\cos \alpha \,\sin \alpha ]^{T}$, $\varvec{g}_{0} = [0\,\cos \alpha \, - \sin \alpha ]T$.

Fig. 1

Fibre structure of organs. a Human ureter and a strip for tension tests (Dorrell et al. [34]). b Small intestine (Gabella [35])

Some modified material laws

If the Fung potential is utilized to characterize the fibres in 3D cases, it can be formulated by combining with an isotropic model like the neo-Hookean model [18, 19].

$$W_{\text{MF}} ({\mathbf{E}},J) = \frac{\mu }{2}(I_{1} - 3) + W_{\text{F}} ({\mathbf{E}},J).$$

(4)

In addition, an isotropic Fung-type potential considered by Nguyen et al. [20] and Duong et al. [21] is formulated as

$$W_{\text{MF}} ({\mathbf{E}},J) = \frac{\mu }{2}(I_{1} - 3) + C[\exp \{ Q_{\text{iso}} ({\mathbf{E}})\} - 1],$$

(5)

where $Q_{\text{iso}} ({\mathbf{E}}) = c_{1} (E_{1}^{2} + E_{2}^{2} + E_{3}^{2} ) + 2c_{2} (E_{1} E_{2} + E_{2} E_{3} + E_{3} E_{1} )$.

If the parameters of a 2D Fung potential are chosen to represent fibres in a fibre plane, then the material law in (4) becomes the MFH (Holzapfel [16]) written as

$$W_{\text{MFH}} ({\mathbf{E}},J) = W_{\text{aniso}} = W_{\text{iso}} (I_{1} ) + 2U(J),$$

(6)

$$W_{\text{aniso}} = C[\exp \{ Q_{\text{fibre}} ({\mathbf{E}})\} - 1];\;W_{\text{iso}} (I_{1} ) = \mu (I_{1} - 3),$$

(7)

$$Q_{\text{fibre}} ({\mathbf{E}}) = c_{1} E_{2}^{2} + c_{2} E_{3}^{2} + c_{3} E_{2} E_{3} ,$$

(8)

where c ₁, c ₂ and c ₃ are dimensionless parameters so that the angle between two fibre families along the hoop direction is defined by $\tan (2\alpha ) = \frac{{2c_{3} }}{{c_{1} - c_{2} }}$, see Fig. 1. The component $W_{\text{aniso}}$ contributes mechanical properties only in the fibre plane O23 as an exponential function, whereas $W_{\text{iso}} (I_{1} ) = \mu (I_{1} - 3)$ characterizes the material response in all directions as a polynomial term. This would lead to the same ill-conditioning problem as described above. Considering a human ureter or a small intestine, see Fig. 1, the microstructure of these tubular organs may include layers of both one fibre and two fibre families. For this application, a slight modification of the Holzapfel model by Nguyen et al. [22] can be used. In general, the mechanical properties of the two fibre families might be different $(k_{2p} \ne k_{3p} ),$ and the orthotropic part (3) could be slightly modified as

$$W_{\text{ani}} (I_{4} ,I_{6} ) = \frac{{k_{1p} }}{2}\{ \exp [k_{2p} (I_{4} - 1)^{2} ] - 1\} + \frac{{k_{1p} }}{2}\{ \exp [k_{3p} (I_{6} - 1)^{2} ] - 1\}$$

(9)

where $k_{1p}$ has the dimension of a modulus, k _2p and k _3p are dimensionless;

The orthotropic part (9) can re-use the material parameters of the original Holzapfel orthotropic part (3) in the published literature using the relations:

$$k_{1p} = k_{1H} /k_{2H} ,\,k_{2p} = k_{2H} ,\,k_{3p} = k_{2H}$$

(10)

When modelling a structure of only one fibre family, the second family in the model (9) is deactivated simply by imposing $k_{3p} = 0$.

Gasser et al. [11] proposed another modified Holzapfel model (called the Holzapfel–Gasser model) obtained by reformulating the orthotropic terms (3) as

$$W_{\text{ani}} (I_{1} ,I_{4} ,I_{6} ) = \sum\limits_{i = 4,6} {\frac{{k_{1H} }}{{2k_{2H} }}[\exp \{ k_{2H} [k_{G} I_{1} + (1 - 3k_{G} )I_{i} - 1]^{2} \} - 1]}$$

(11)

This model has incorporated a so-called fibre dispersion factor $k_{\text{G}} = [0,1/3]$. For $k_{\text{G}} = 0,$ the energy function (11) becomes the orthotropic part (3) of the original Holzapfel model. For $k_{\text{G}} = 1/3,$ the energy function (11) represents an isotropic material.

Non-physical behaviour

To obtain the unique solution of a boundary-value problem, the strain–energy function is required to satisfy the quasi-convexity condition. Ball [23] proposed a polyconvexity concept implying quasi-convexity and strong ellipticity. This guarantees material stability. For isotropic materials, the neo-Hookean model satisfies the concept of polyconvexity (Schröder and Neff [25]). Similarly, the Holzapfel model is also a polyconvex function with positive material parameters (Schröder and Neff [25]). In addition, the MFH model is also polyconvex (Balzani [24]). However, polyconvexity is not sufficient to ensure physical responses of the Holzapfel model and of the MFH as described in the following.

The Holzapfel model and the MFH model both suffer from the ill-conditioning problem even in the physiological deformation range of tissues with improper sets of material constants. This leads to a volume growth and hence the incompressibility constraint is violated. The unrealistic response of the orthotropic model and the transversely isotropic model in uniaxial tension tests are depicted in Fig. 2. A volumetric–isochoric split decomposes the strain energy into a volumetric part and a deviatoric term and might cause the non-physical behaviour of materials not restricted to (near) incompressibility (Ehlers and Eipper [13]). An alternative solution was also suggested by Sansour [26] and afterwards realized by Helfenstein et al. [12], in which they have adopted the original invariants of the right Cauchy–Green tensor for fibres, whereas the incompressibility was imposed on the isotropic strain energy. They also stated that the augmented Lagrange method is an alternative approach. However, we demonstrate that the augmented Lagrange method can just mitigate the ill-conditioning problem. Therefore, the paper presents a different source of the ill-conditioning through several tension tests employing the Holzapfel and the MFH.

Fig. 2

Non-physical responses of the Holzapfel model (the directions $r,\theta$ and z are denoted 1, 2 and 3, respectively)

Holzapfel model with two fibre families and one fibre family

We perform uniaxial tension tests on a rectangular adventitial strip exhibiting two fibre families, see Fig. 1. Gasser et al. [11]) used the Holzapfel model to conduct these tests showing the thickening by employing pseudo-material constants, which are thus unable to describe the mechanics of the adventitial strip. In contrast, for our numerical tests, we utilized the set of material constants obtained from curve fitting to experimental data. Therefore, the model using these material coefficients characterizes the mechanics of the adventitia. In addition, the tests were carried out in the physiological deformation range of the adventitial specimen, but more importantly, the ill-conditioning was persistently observed. To this end, the material parameters of the Holzapfel model for the case of two fibre families (Fig. 2a) were identified by curve fitting for two uniaxial tensile tests along the axial and hoop directions of an arterial wall (see Fig. 3) in which the data of an arterial wall were reproduced from Holzapfel et al. [27]); $\mu = 22.722$ (kPa), $k_{H1} = 94.3205$ (kPa), $k_{H2} = 240.5662$, α = 49.3552°. Moreover, the parameters for a transversely isotropic model are identical to those in the case of two fibre families, except for the fibre angle, see Fig. 2b.

Fig. 3

Curves fitted to experimental data of an adventitial layer (Holzapfel et al. [27]) using the Holzapfel and the MFH models

The simulation of a half of the specimen in Fig. 2 with Poisson’s ratio v = 0.49996 shows clearly the thickening of the strip at the middle, and this effect becomes more severe if the load continuously increases. Specifically, the thickness in Or direction (O1, displacement 1) orthogonal to the fibre plane $O\theta z$ (O23) for the two fibre families increases which is inconsistent with the physically expected transverse contraction. Interestingly, for the case of one fibre family, the thickness increases in all directions in a plane orthogonal to the fibre direction O3, see Fig. 2. Although this set of material parameters ensures the polyconvexity of the Holzapfel model, the numerical results are, however, incorrect. The MFH model is subjected to this problem in the same manner.

Ill-conditioning and condition number of constitutive matrix

Considering the ${\mathbb{C}}$, the ill-conditioning is due to the large condition number of the constitutive matrix. Since the elasticity tensor is symmetric and positive definite, its 2-norm condition is defined as

$$\kappa = \frac{{\lambda_{\hbox{max} } ({\mathbb{C}})}}{{\lambda_{\hbox{min} } ({\mathbb{C}})}}$$

(12)

where $\lambda_{\hbox{max} } ({\mathbb{C}})$ and $\lambda_{\hbox{min} } ({\mathbb{C}})$ are the maximal and minimal eigenvalues of ${\mathbb{C}}$, respectively. A large condition number $\kappa \gg 1$ indicates an ill-conditioned matrix and low accuracy of the results of matrix manipulations. Thus, the larger the condition number is, the worse the solution of the equation system becomes. In general, it can be simply said that the numerical solution of the linear equation system might lose about k decimal digits of precision when its condition number is $\kappa = 10^{k}$ (Liu and Chang [28]). In addition, if the solution is obtained from an iterative process, then its numerical error is iteratively accumulated. Iterative numerical analyses of biological soft tissues based on nonlinear materials subjected to large deformations severely suffer from such cumulated errors. An additional important issue to be discussed is the condition number of the elasticity tensor. Since the Holzapfel model is an additive function of a linear polynomial isotropic term and an exponential anisotropic part, the unbalance between the two components increases exponentially with increase of deformation. When the unbalance attains some bound, the ill-conditioning of the elasticity tensor ${\mathbb{C}}$ may become effective. As it was shown by Sun and Sacks [15], even if the material model is built up from convex functions, the resulting model may be subjected to the ill-conditioning, which in the case of the Holzapfel model accrues from magnitude differences amongst components of the elasticity tensor as shown below. In such a case, the numerical analysis is unable to give an accurate solution. Thus, the elasticity tensor must be well-conditioned, apart from the convexity or the polyconvexity condition. The physical behaviour in relation to the ill-conditioning is illustrated and analysed in the following sections.

Results

Numerical tests: physical behaviour of submucosa layer

In this section, we indicate the relation between the ill-conditioning and the non-physical behaviour of the Holzapfel model and the MFH model. Obviously, for a 3D boundary-value problem, it is impossible to compute analytically the condition number because of complicated matrix manipulations. Therefore, the condition number is calculated numerically for our tension tests.

Equibiaxial test of submucosa

Material parameters of the Holzapfel model are obtained by curve fitting to the intestinal experimental data which were reproduced from Ciarletta et al. [29]. Biaxial deformation is prescribed by imposing two in-plane stretches in the interval 1.0–1.3: the longitudinal stretch $\lambda_{\text{z}}$ and the hoop stretch $\lambda_{\theta }$ (by the incompressibility constraint $\lambda_{\text{r}} = 1/(\lambda_{\theta } ,\lambda_{\text{z}} )$ through the specimen thickness). Then, the Cauchy stresses in the longitudinal, hoop and radial directions are functions of the defined stretches. Material parameters of the Holzapfel model were identified using a nonlinear algorithm which minimizes the function:

$$\frac{1}{2}\sum\limits_{i = 1}^{n} {(\sigma_{\theta }^{e} - \sigma_{\theta }^{p} )_{i}^{2} + \sum\limits_{i = 1}^{n} {(\sigma_{z}^{e} - \sigma_{z}^{p} )_{i}^{2} } }$$

(13)

where $\sigma_{\theta }^{e} ,\sigma_{z}^{e}$ are stresses from the experimental data of Ciarletta et al. [29], $\sigma_{\theta }^{p} ,\sigma_{z}^{p}$ are stresses calculated from the Holzapfel model and $i = 1, \ldots ,n$ are data points.

Furthermore, it is well known that the curve fitting may result in non-unique sets of constants and some of these sets can produce unphysical behaviour (Sun and Sacks [15]; Schmidt [30]). Let us choose the shear modulus given by Ciarletta et al. [29] as the reference, $\mu^* = 1.58\;{\text{kPa}}$. Apart from the set of material parameters by Ciarletta et al. [29], similar to the ones by Holzapfel et al. [8], two other material sets with the shear moduli varying by factors of 10 and 0.1 were produced with the corresponding values of goodness of fit. Note that during the curve fitting, the angle of the model in Ciarletta et al. [29] was fixed ($\alpha = 60^{ \circ }$) so that the two new fitted material sets of the Holzapfel model can characterize the mechanical properties close to those of the original material constants in Ciarletta et al. [29]. Obviously, Fig. 4 shows a calibration of three material sets for the fitted curves which are very close together in the “toe region”. The stress–stretch curves are drawn by the analytical relationships in the incompressibility condition of the material. However, to avoid the high nonlinearity of the incompressibility constraint in finite element analysis, we can solve the problem for nearly incompressible materials through the penalty method imposed by the volumetric strain–energy function in which Poisson’s ratio is set to be as close as possible to 0.5 (in literature $\nu = \;0.4996$ is commonly used). Therefore, we need to write the material law for finite element implementation in the nearly incompressible form.

Fig. 4

Equibiaxial test: fitted curves with $\mu = \mu^{*}$, $\mu = \mu^{*} /10$ and $\mu = 10\mu^{*}$

The dimension of the simulated specimen was as follows: $L_{\text{r}} L_{\theta } L_{\text{z}} = 0.2 \times 25 \times 25\,{\text{mm}}^{3}$. For the sake of symmetry, the model was simulated by only one single, mixed finite element because there is only the volumetric locking effect in this case. A penalty method was used to enforce the incompressibility condition of the model. To prevent the ill-conditioning of the stiffness matrix that may occur at high levels of incompressibility, i.e. $0.4996 < \nu < 0.5$, the Uzawa algorithm, an augmented Lagrangian method, was utilized (Taylor [31]). The displacements $d = 0.35L_{\theta }$ were identically prescribed in both $\theta$ and z directions in 350 equal time steps.The predicted Cauchy stresses

$$\begin{aligned} &\sigma _{{\text{r}}} = 0 \Rightarrow p = S_{{\text{r}}} \lambda _{{\text{r}}}^{2} , \\ & \sigma _{\theta } = S_{\theta } \lambda _{\theta }^{2} - p,\sigma _{{\text{z}}} = S_{{\text{z}}} \lambda _{{\text{z}}}^{2} - p \\ \end{aligned}$$

(14)

will be considered as the gold standard to check the numerical behaviour of the Holzapfel model in which S _r, $S_{\theta }$ and $S_{\text{z}}$ are the principal second Piola–Kirchhoff stresses. The structural tensors $\varvec{a}_{0} \otimes \varvec{a}_{0}$ and $\varvec{g}_{0} \otimes \varvec{g}_{0}$ define the principal orientations of the orthotropic material, indicating the stiffest and the mid-compliant orientations, respectively. Thus, the stress–stretch curves along $\theta$ and z directions are different. In the equibiaxial tension test, see Figs. 5, 6, for $\mu = \mu^*/10,$ there is a good agreement between the curve 2 and the predicted one for stretches smaller than a critical stretch $\lambda_{\text{c}}^{\mu^*/10} = 1.25$. This is the largest stretch at which the response of the model is still sufficiently close to the analytical stress. For larger stretches $\lambda > \lambda_{\text{c}}^{\mu^*/10}$, the thickness of the simulated specimen begins to increase and therefore the incompressibility condition is violated. In contrast, the neo-Hookean model is immune from this thickening, see curve 5 in Fig. 5. Thus, the unphysical deformation directly relates to the co-existence of the isotropic part and the anisotropic term in the strain–energy function. As a consequence, the unphysical deformation results in incorrect stresses; see the stress–stretch curves 3–6 in Fig. 6. However, employing a larger shear modulus $\mu$ ($\mu = \mu^*$ or $\mu = 10\mu^*$) can greatly improve the simulations by increasing the critical stretch. Moreover, the ratio $W_{\text{ani}} /W_{\text{iso}}$ between the exponential function $W_{\text{ani}}$ and the linear polynomial function $W_{\text{iso}}$ is approximately linear for small deformation, but increases exponentially at higher stretch levels, see Fig. 7. As seen in Fig. 5, this rapid growth for $\mu = \mu^*/10$ occurs at stretches larger than $\lambda_{c}^{\mu^*/10}$ and is accordingly associated with the increase of the condition number, see Fig. 8. The smallest shear modulus $\mu = \mu^*/10$ leads to the smallest critical stretch compared to the other cases, see Fig. 5.

Fig. 5

Equibiaxial test: comparison of displacements in thickness direction of the strip of an intestinal layer for three values of the shear modulus ($\mu^{*} = 1.58$ kPa; v = 0.4996)

Fig. 6

Equibiaxial test: effects of the shear modulus on stress–stretch curves; FEM results compared to analytical ones (prediction)

Fig. 7

Equibiaxial test: comparison of the predicted strain–energy ratios for the three cases of $\mu$

Fig. 8

Equibiaxial test: condition numbers and relative difference between condition numbers in the cases; $\mu^{*} /10$ vs. $\mu^{*}$, $\mu^{*}$ vs. $10\mu^{*}$ and $\mu^{*} /10$ vs. $10\mu^{*}$ ($\mu^{*} = 1.58$ kPa; v = 0.4996)

Uniaxial tension test of submucosa

A uniaxial tension test in the axial direction of the submucosa layer was investigated with the material parameters in Table 1. The specimen has the dimension of $L_{\text{r}} L_{\theta } L_{\text{z}} = 0.5 \times 3 \times 12\,{\text{mm}}^{3}$. Due to symmetry, only one-eighth of the specimen is modelled, see Fig. 9, by 2376 hexahedral mixed elements (Taylor [31]). An axially given displacement $0.15L_{\text{z}}$ was prescribed at one end of the specimen model. For the shear moduli $\mu = \mu^{*} /10$ and $\mu = \mu^{*}$ the thickness of the specimen increases for the first step of the simulation; hence, the critical stretches are very close to 1, see Fig. 9. This occurs because the axial tension increases the stiffness of the fibres faster than in the case of the equibiaxial tension since the fibres form a small angle (30°) with respect to the axial axis. For $\mu = 10\mu^{*}$, the thickness displacement of the specimen is physically decreased as shown in Fig. 9 since the stiffness difference between the fibres and the ground matrix is smaller and therefore the incompressibility condition is ensured. This phenomenon can be avoided or mitigated when the material parameters are suitably chosen. Specifically, the shear modulus should be identified as large as possible.

Table 1 Material parameters for the submucosa layer of large intestine (μ* = 1.58 kPa)

Fig. 9

Uniaxial tensile test in z direction. a Displacement in thickness direction r [mm] and b Cauchy stress [kPa] in z (3) direction: left side $\mu /\mu^{*} = 10$—without thickening; right side $\mu /\mu^{*} = 0.1$—thickening (with deformed and un-deformed shape)

Unrealistic response of the MFH model with two fibre families

To obtain the material constants for the MFH model (6), a curve fitting phase was performed on the adventitial data (Holzapfel et al. [27]), see Fig. 3. All coefficients are checked for the convexity conditions which are derived as $c_{3} ,c_{3} ,c_{3} > 0$, $c_{1} c_{2} - c_{3}^{2} > 0$. Specifically, the coefficients are identified as follows: $\mu^* = 22.722$, $C = 0.185$, $c_{1} = 180.082$, $c_{2} = 373.143$, $c_{3} = 245.052$. Similar to the investigation of the Holzapfel model in equibiaxial tensile test with the three values of the shear modulus, we also adopted $\mu = 10\mu^{*}$, $\mu = 0.1\mu^{*}$ and $\mu = 10\mu^{*}$ to simulate the radial stretch in uniaxial tensile tests as illustrated in Fig. 10b. Although the model is polyconvex, it is strongly subjected to the ill-conditioning problem in the physiological deformation range of the tissue as shown in Fig. 10 in the case $\mu = \mu^{*}$. The results show that using a higher value of the shear modulus leads to a better numerical stability. This finding is similar to the one in the case that the Holzapfel model is used.

Fig. 10

a Unphysical response b Unrealistic numerical results of stretch in tension test of the MFH model for the adventitial layer (Holzapfel et al. [27])

Discussion

The equibiaxial tests are investigated since it is convenient to study the condition number of the elasticity tensors of a single finite element model. The stiffness of the submucosa similar to that of the arterial adventitia is essentially determined by the material parameters for collagen fibres, whereas the stiffness of the bulk matrix is characterized by the shear modulus $\mu$. Although this stiffness is much smaller than the one of collagen fibres, a small change of $\mu$ leads to a considerable change in the simulation results. This numerical effect is discussed with the influence of the condition number of the elasticity tensor, which is composed of the three contributing terms ${\mathbb{C}} = {\mathbb{C}}_{\text{iso}} + {\mathbb{C}}_{\text{ani}} + {\mathbb{C}}_{\text{vol}}$, see also the "Appendix". For all considered simulations a fixed value of v, say v = 0.5, was used. Thus, the difference in the condition numbers corresponding to the three considered cases with three different values of the shear modulus is unaffected by the magnitudes of the components of the volumetric part ${\mathbb{C}}_{\text{vol}}$. In the anisotropic part ${\mathbb{C}}_{\text{ani}}$ there are certain entries that strongly dominate the others due to the large mechanical contribution by the reinforcing fibre parameters. Specifically, since all of the fibres are so-called in-plane fibres in the plane ($\theta - z$), the in-plane shear term ${\mathbb{C}}_{{{\text{ani}},\theta z}}$ and all normal terms are much larger than the remaining out-of-plane shear terms, e.g. ${\mathbb{C}}_{{{\text{ani,r}}\theta }}$ and ${\mathbb{C}}_{{{\text{ani,}}\;{\text{r}}z}}$, see Fig. 11. On the contrary to the irregular entries of ${\mathbb{C}}_{\text{ani}}$, the corresponding magnitude difference between the in-plane shear terms and the out-of-plane shear terms of ${\mathbb{C}}_{\text{iso}}$ (by the neo-Hookean model) is small, see Fig. 11, curves 5 and 6. As a consequence, the three contributing terms ${\mathbb{C}}_{\text{iso}} + {\mathbb{C}}_{\text{ani}} + {\mathbb{C}}_{\text{vol}}$ in the Holzapfel model lead to the elasticity tensor ${\mathbb{C}}$ possessing highly different magnitudes of its components, see Fig. 11. At large deformation, the in-plane shear term ${\mathbb{C}}_{{\theta {\text{z}}}}$ can be even larger than the normal terms. Therefore, the ill-conditioning of the elasticity tensor ${\mathbb{C}}$ is caused by the magnitude differences between the shear term ${\mathbb{C}}_{{\theta {\text{z}}}}$ and the shear terms (${\mathbb{C}}_{{{\text{r}}\theta }}$,${\mathbb{C}}_{\text{rz}}$). Moreover, it is also observed in the "Appendix" that ${\mathbb{C}}_{\text{ani}}$ contains exponential terms [terms A and B in (A.3) and (A.5)]. Consequently, once the deformation is sufficiently large, these differences between the components of ${\mathbb{C}}$ are extremely large (increasing exponentially) and in turn lead to an ill-conditioned elasticity tensor. In contrast, when the neo-Hookean law was used the resulting condition number increases almost linearly with increasing stretch, similarly to the difference between the neo-Hookean normal terms and the shear terms shown in Fig. 11. The neo-Hookean elasticity tensor has a low condition number curve, see the dotted curve in Fig. 8 and consequently it produces the correct solution, see curve 5 in Fig. 5.

Fig. 11

Equibiaxial test: ratio between components of the elasticity tensor. The magnitude of $\mu$ strongly affects the difference between components of the elasticity tensor of the Holzapfel model, but has no effect on the neo-Hookean model

For the Holzapfel model, the difference amongst the components of ${\mathbb{C}}$ can be adjusted through the shear modulus $\mu$. If the magnitudes of the components of ${\mathbb{C}}_{\text{iso}}$ are more comparable to the corresponding ones of ${\mathbb{C}}_{\text{ani}}$, e.g. using a larger shear modulus $\mu$, then the relative difference amongst the components of ${\mathbb{C}}$ are reduced, see Fig. 11. This helps decrease the condition number of the elasticity tensor in the equibiaxial simulation and hence the critical stretch is significantly increased, see Fig. 8.

For the case of the smallest shear modulus $\mu = \mu^{*} /10$, the unphysical deformation begins from the smallest critical stretch $\lambda = \lambda_{\text{c}}^{\mu^*/10}$, see Fig. 5. When larger values of the shear modulus $\mu = \mu^{*}$ and $\mu = 10\mu^{*}$ were used, the critical stretches increases significantly, denoted by $\lambda_{\text{c}}^{\mu^* }$ and $\lambda_{\text{c}}^{10\mu^* }$, respectively. Consequently, we have $\lambda_{\text{c}}^{{10\mu^{*} }} > \lambda_{\text{c}}^{{\mu^{*} }} > \lambda_{\text{c}}^{{\mu^{*} /10}}$. Correspondingly, Fig. 8 shows the relative differences of the condition numbers for the three cases. These differences always start from the corresponding lower critical stretch, at which the unphysical deformation begins to take place. For example, if $\mu = 10\mu^{*}$ and $\mu = 100\mu^{*}$ the condition number is decreased by 55 and 75 %, respectively. These moderate reductions already ensure the physical behaviour of the material model. When the augmented Lagrange method was not utilized, the smallest critical stretch was decreased, $\lambda_{\text{c}}^{{\mu^{*} /10}} = 1.2$, since the other ill-conditioning is caused by the nearly incompressible constraint. Therefore, the augmented Lagrange method can only mitigate the ill-conditioning problem but cannot prevent this in the deformation range of interest as an alternative solution as suggested by Helfenstein et al. [12]. Thus, apart from the source of the thickening effect proposed by Helfenstein et al. [12], we demonstrated that it is caused by the large difference between the two energy terms. Figure 2 shows that the unphysical deformation is observed in a plane or in directions orthogonal to the fibre direction. This is also supported by the equibiaxial test with the unrealistic displacement perpendicular to the fibre plane.

Moreover, the model of Weiss et al. [32] is constructed by an isotropic strain–energy term as a polynomial function in terms of the invariants $I_{1}$ and $I_{2}$, whereas the anisotropic term is characterized as an exponential function for one embedded fibre family. Therefore, this model would have a large difference in the mechanical contribution between the isotropic energy and the anisotropic energy, leading to the same instability problem as reported in Helfenstein et al. [12]. The thickening effect was also similarly observed for pure shear tests (Duong [17]). More importantly, the same is found for simulation of tension tests using of the MFH model. To this end, the influence of the shear modulus was also investigated for the MFH model as shown in Fig. 10. The model shows the bad effect in the physiological deformation range of the tension test with $\mu = \mu^{*}$ and $\mu = 0.1\mu^{*}$. The higher the shear modulus ($\mu = 10\mu^{*}$) is chosen, the more stable the MFH model becomes. Thus, the critical stretch increases with the increase of the shear modulus which should be chosen as large as possible for numerical stability. For the case of an incompressible Fung model in-plane stress, Sun and Sacks [15] imposed both the upper bound for the condition number and the convexity for the model to achieve numerical stability. However, this approach is limited to a specific case of plane stress. As shown in Fig. 8, the condition number is a function of deformation and the material parameters. Thus, preventing the ill-conditioning by fixing a predicted upper bound for the condition number is not always plausible.

Figure 5 shows an equibiaxial test simulated by the neo-Hookean model with only a linear polynomial function. The corresponding condition number curve of the neo-Hookean elasticity tensor is regular in Fig. 8. Moreover, Helfenstein et al. [12] have even utilized the real measure of fibre stretch ($I_{4}$ and $I_{6}$) instead of ($\overline{{I_{4} }}$ and $\overline{{I_{6} }}$) for the model by Rubin and Bodner [33], but there is of course negligible effect observed. This is due to the fact that there are only exponential terms in the strain–energy formulation. Thus, the instability problem originates from the unbalance of strain energy. Moreover, a model composed of the neo-Hookean and an isotropic version of the Fung-type potential (Nguyen et al. [20]; Duong et al. [21]) [see (5)] is also immune from the ill-conditioning even at very high stretch up to 2.0 due to the absence of an anisotropic term. A model proposed by Peng et al. [37] has both the fibre and the matrix terms described by polynomial functions and a shear interaction is modelled by a product of a weak exponential and a polynomial function in terms of all principal invariants. Therefore, the difference between the anisotropic and isotropic strain energy is equilibrated, leading to a stable model. In addition, Gasser et al. [11] have modified the structure tensor characterizing dispersion of fibre orientations to overcome the unphysical deformation. As discussed before, the fibre dispersion helps reduce the large difference between the isotropic strain energy and the anisotropic strain energy and therefore ensure stability of their model. However, the fibre dispersion is not always easily obtained due to limits of experiments

Conclusions

Generally, the larger difference in mechanical contributions apparently gives rise to the larger differences amongst the components of the material tangent stiffness matrix and therefore leads to the ill-conditioned stiffness matrix. The similar ill-conditioning is known from the simulation of composite structures comprising two main constituents with a great difference in mechanical stiffness. Typical examples of this are the organic-coated metals. Throughout the investigation of the ill-conditioning problem, the magnitude of $\mu$ strongly affects the simulation accuracy of the Holzapfel model but has no influence on the accuracy of the neo-Hookean model. Therefore, the thickening effect does not originate from the “fibre rotation” as discussed by Gasser et al. [11]. Moreover, the shear modulus also significantly influences the numerical stability of the MFH model but has no effect on any isotropic Fung-type model (Duong [17]). It accounts for the ill-conditioning problem, which results from the large condition numbers of the elasticity tensors of the Holzapfel model and the MFH model. Since the shear modulus $\mu$ of soft tissue usually varies in a bounded set, the shear modulus $\mu$ should be chosen closer to the physically upper bound (its natural value). Alternatively, Shore hardness indentation tests could be used to measure the shear modulus of the matrix directly (Duong [17]). These measured values of the shear modulus help ensure the numerical stability of the material law (Duong [17]).

Appendix

To implement the material laws for nearly incompressible materials into finite element schemes, we need to employ the penalty method using the volumetric strain–energy function as presented above. To this end, we split the deformation gradient F into the isochoric part $\overline{{\mathbf{F}}} = J^{ - 1/3} {\mathbf{F}}$ and the volumetric part $J^{1/3}$. Thus, $\overline{{\mathbf{F}}}$ is the modified deformation gradient tensor associated with the modified invariants $\overline{I}_{a}$ $(a = 1,4,6).$ The second Piola–Kirchhoff stress for the Holzapfel model is calculated as

$${\mathbf {S}} = 2\frac{\partial W}{\partial {\mathbf C}} = 2\sum\limits_{a = 1,4,6} {\frac{\partial W}{{\partial \overline{{I_{a} }} }}} \left( {\frac{{\partial \overline{I}_{a} }}{{\partial I_{a} }}\frac{{\partial I_{a} }}{\partial {\mathbf C}} + \frac{{\partial \overline{I}_{a} }}{{\partial I_{3} }}\frac{{\partial I_{3} }}{\partial {\mathbf C}}} \right) = {\mathbf {S}}_{\text{iso}} + {\mathbf {S}}_{\text{vol}} + {\mathbf {S}}_{\text{ani}}$$

(A.1)

$\overline{I}_{a}$ $(a = 1,4,6)$ are the modified invariants when the volumetric–isochoric split is applied to the strain–energy function. The two terms $S_{\text{iso}}$ and $S_{\text{vol}}$ are derived from the isotropic energy and the volumetric energy, respectively. They are available in textbooks (Holzapfel [36]). The anisotropic part of the second Piola–Kirchhoff stress is

$$\begin{aligned} {\mathbf {S}}_{\text{ani}} &=& 2\frac{{\partial W_{\text{ani}} }}{\partial {\mathbf C}} = 2\left[ {\frac{{\partial W_{\text{ani}} }}{{\partial \overline{I}_{4} }}\frac{{\partial \overline{I}_{4} }}{{\partial I_{4} }}\frac{{\partial I_{4} }}{\partial {\mathbf C}} + \frac{{\partial W_{\text{ani}} }}{{\partial \overline{I}_{6} }}\frac{{\partial \overline{I}_{6} }}{{\partial I_{6} }}\frac{{\partial I_{6} }}{\partial {\mathbf C}}} \right] \\&&+ 2 {\left[ {\frac{{\partial W_{\text{ani}} }}{{\partial \overline{I}_{4} }}\frac{{\partial \overline{I}_{4} }}{{\partial I_{3} }}\frac{{\partial I_{3} }}{\partial {\mathbf C}} + \frac{{\partial W_{\text{ani}} }}{{\partial \overline{I}_{6} }}\frac{{\partial \overline{I}_{6} }}{{\partial I_{3} }}\frac{{\partial I_{3} }}{\partial {\mathbf C}}} \right]}\end{aligned}$$

(A.2)

The strain–energy function $W_{\text{ani}}$, Eq. (9), is an exponential function and so are all its derivatives. The first derivatives with respect to the invariants are

$$\frac{{\partial W_{\text{ani}} }}{{\partial \overline{I}_{4} }} = k_{1p} k_{2p} (\overline{{I_{4} }} - 1)\underbrace {{\exp [k_{2p} (\overline{{I_{4} }} - 1)^{2} ]}}_{\text{A}},\,\,\,\frac{{\partial W_{\text{ani}} }}{{\partial \overline{I}_{6} }} = k_{1p} k_{3p} (\overline{{I_{6} }} - 1)\underbrace {{\exp [k_{3p} (\overline{{I_{6} }} - 1)^{2} ]}}_{\text{B}}$$

(A.3)

The constitutive matrix is calculated as a second derivative of the strain energy function with respect to the right Cauchy–Green strain as:

$${\mathbb{C}} = 4\frac{{\partial^{2} W}}{{\partial {\mathbf C}^{2} }} = {\mathbb{C}}_{\text{iso}} + {\mathbb{C}}_{\text{vol}} + {\mathbb{C}}_{\text{ani}}$$

(A.4)

The two terms ${\mathbb{C}}_{\text{iso}}$ and ${\mathbb{C}}_{\text{vol}}$ are derived from the isotropic energy and the volumetric energy, respectively. They are available in textbooks (Holzapfel [36]). The anisotropic term is the exponential function

$$\begin{aligned} \mathbb{C}_{{{\text{ani}}}} & = 2\frac{{\partial S_{{{\text{ani}}}} }}{{\partial {\mathbf C}}} = 4k_{{1p}} \frac{\partial }{{\partial {\mathbf C}}}\left[ {k_{{2p}} A(\bar{I}_{4} - 1)I_{3}^{{ - 1/3}} \left( {\user2{a}_{0} \otimes {\text{ }}\user2{a}_{0} - \frac{1}{3}\bar{I}_{4} {\mathbf C}^{{ - 1}} } \right)} \right. \\ & \left. {\quad + k_{{3p}} B(\bar{I}_{6} - 1)I_{3}^{{ - 1/3}} \left( {\user2{g}_{0} \otimes \user2{g}_{0} - \frac{1}{3}\bar{I}_{6} {\mathbf C}^{{ - 1}} } \right)} \right] \\ \end{aligned}$$

(A.5)

In (A.5) we need to use following derivatives:

$$\begin{aligned} &\frac{{\partial A}}{{\partial {\mathbf C}}} = \frac{{\partial A}}{{\partial \overline{{I_{4} }} }}\frac{{\partial \overline{{I_{4} }} }}{{\partial {\mathbf C}}} = 2k_{{2p}} (\overline{{I_{4} }} - 1)A\frac{{\partial \overline{{I_{4} }} }}{{\partial {\mathbf C}}},\,\frac{{\partial B}}{{\partial {\mathbf C}}} = \frac{{\partial B}}{{\partial \overline{{I_{6} }} }}\frac{{\partial \overline{{I_{6} }} }}{{\partial {\mathbf C}}} = 2k_{{3p}} (\overline{{I_{6} }} - 1)B\frac{{\partial \overline{{I_{6} }} }}{{\partial {\mathbf C}}},\,\\ &\frac{{\partial ^{2} W_{{{\text{ani}}}} }}{{\partial \overline{{I_{4} }} ^{2} }} = k_{{1p}} k_{{2p}} A[1 + 2k_{{2p}} (\overline{{I_{4} }} - 1)^{2} ],\,\frac{{\partial ^{2} W_{{ani}} }}{{\partial \overline{{I_{6}^{2} }} }} = k_{{1p}} k_{{3p}} B[1 + 2k_{{3p}} (\overline{{I_{6} }} - 1)^{2} ] \end{aligned}$$

(A.6)

These above quantities are suitable for the total Lagrange formulations. For the updated Lagrange formulations, a push-forward step is furthermore required.