Incorporating fiber recruitment in hyperelastic modeling of vascular tissues by means of kinematic average

Abstract

We present a framework for considering the gradual recruitment of collagen fibers in hyperelastic constitutive modeling. An effective stretch, which is a response variable representing the true stretch at the tissue-scale, is introduced. Properties of the effective stretch are discussed in detail. The effective stretch and strain invariants derived from it are used in selected hyperelastic constitutive models to describe the tissue response. This construction is investigated in conjunction with Holzapfel-Gasser-Ogden family strain energy functions. The ensuing models are validated against a large body of uniaxial and bi-axial stress–strain response data from human aortic aneurysm tissues. Both the descriptive and the predictive capabilities are examined. The former is evaluated by the quality of constitutive fitting, and the latter is assessed using finite element simulation. The models significantly improve the quality of fitting, and reproduce the experiment displacement, stress, and strain distributions with high fidelity in the finite element simulation.

Appendices

Appendix 1 Proof of properties of ${{\bar{\lambda }}}$

Properties 1, 2, 3, and 4 hinge on the derivative of ${{\bar{\lambda }}}$:

$$\begin{aligned} {{\bar{\lambda }}}^\prime = \frac{\alpha }{\alpha +\beta } F(\lambda ^{-1}, \alpha +1,\beta ,\text {upper}) + \frac{d{}}{d{\lambda }} \int _0^{\frac{1}{\lambda }} \frac{p^{\alpha -1}( 1- p)^{\beta -1} }{B(\alpha ,\beta )} \,dp +\lambda \frac{d{}}{d{\lambda }} \int _{\frac{1}{\lambda }}^1 \frac{p^{\alpha }( 1- p)^{\beta -1} }{B(\alpha ,\beta )} \,dp \end{aligned}$$

(4.10)

A straightforward computation shows that the second and third terms cancel out, leaving

$$\begin{aligned} {{\bar{\lambda }}}^\prime = \frac{\alpha }{\alpha +\beta } F(\lambda ^{-1};\alpha +1,\beta ,\text {upper}) \end{aligned}$$

(4.11)

Clearly,

$$\begin{aligned} {{\bar{\lambda }}}^\prime (1) = 0, \quad \lim _{\lambda \rightarrow \infty }{{\bar{\lambda }}}^\prime (\lambda ) = \frac{\alpha }{\alpha +\beta } \end{aligned}$$

(4.12)

Properties 3 and 4 are proved. Further, it is evident that

$$\begin{aligned} 0< {{\bar{\lambda }}}^\prime < 1 \quad \forall \lambda \in (1, \infty ) \end{aligned}$$

(4.13)

and

$$\begin{aligned} {{\bar{\lambda }}}(1) = 1 \end{aligned}$$

(4.14)

Properties 1 and 2 are the consequences of these results. First, ${{\bar{\lambda }}}^\prime >0$ indicates that ${{\bar{\lambda }}}$ is a strictly increasing function in $(1,\infty )$. Consequently,

$$\begin{aligned} {{\bar{\lambda }}}(\lambda ) > {{\bar{\lambda }}}(1) = 1 \quad \forall \lambda \in (1, \infty ) \end{aligned}$$

(4.15)

Similarly, ${{\bar{\lambda }}}^\prime < 1$ implies that ${{\bar{\lambda }}}(\lambda ) - \lambda$ is a strictly decreasing function in $(1,\infty )$, and thus

$$\begin{aligned} {{\bar{\lambda }}}(\lambda ) - \lambda < {{\bar{\lambda }}}(1) - 1 = 0 \quad \forall \lambda \in (1, \infty ) \end{aligned}$$

(4.16)

We thus conclude that

$$\begin{aligned} 1< {{\bar{\lambda }}}(\lambda ) < \lambda \quad \forall \lambda \in (1, \infty ) \end{aligned}$$

(4.17)

To show that ${{\bar{\lambda }}}(\mathbf{F})$ is a convex function of $\mathbf{F}$, first observe that ${{\bar{\lambda }}}^\prime = \int _{\frac{1}{\lambda }}^1 \frac{p^{\alpha }( 1- p)^{\beta -1} }{B(\alpha ,\beta )} \,dp$, and hence

$$\begin{aligned} {{\bar{\lambda }}}^{''}= \frac{d{}}{d{\lambda }} \int _{\frac{1}{\lambda }}^1 \frac{p^{\alpha }( 1- p)^{\beta -1} }{B(\alpha ,\beta )} \,dp = \frac{\left( \frac{1}{\lambda }\right) ^{\alpha +2}\left( 1- \left( \frac{1}{\lambda }\right) \right) ^{\beta -1} }{B(\alpha ,\beta )} > 0 \quad \text {for}\quad \lambda \in (1,\infty ) \end{aligned}$$

(4.18)

By construction, ${{\bar{\lambda }}} = 1$ for $\lambda \le 1$, which is convex. Thus, ${{\bar{\lambda }}}$ is a convex function of $\lambda$ over $(0,\infty )$. It follows that

$$\begin{aligned} \left. {{\bar{\lambda }}}(\lambda _2) - {{\bar{\lambda }}}(\lambda _1) \ge \frac{d{{{\bar{\lambda }}}}}{d{\lambda }}\right|_{\lambda =\lambda _1} (\lambda _2 - \lambda _1) \quad \text {for}\ \lambda \in (0,\infty ) \end{aligned}$$

(4.19)

At the same time, $\lambda =\sqrt{\mathbf{F}\mathbf{N}\cdot \mathbf{F}\mathbf{N}}$ is a convex function of $\mathbf{F}$ for any fixed $\mathbf{N}$ (Hartmann and Neff 2003), so that

$$\begin{aligned} \left. \lambda (\mathbf{F}_2) - \lambda (\mathbf{F}_1) \ge \frac{d{\lambda (\mathbf{F})}}{d{\mathbf{F}}}\right|_{\mathbf{F}= \mathbf{F}_1}:(\mathbf{F}_2 -\mathbf{F}_1) \end{aligned}$$

(4.20)

Combining these two inequalities and noticing $\frac{d{{{\bar{\lambda }}}}}{d{\lambda }} \ge 0$, we conclude

$$\begin{aligned} \left. {{\bar{\lambda }}}(\mathbf{F}_2) - {{\bar{\lambda }}}(\mathbf{F}_1) \ge \frac{d{{{\bar{\lambda }}}}}{d{\lambda }} \frac{d{\lambda (\mathbf{F})}}{d{\mathbf{F}}}\right|_{\mathbf{F}= \mathbf{F}_1}: \left. (\mathbf{F}_2 -\mathbf{F}_1) = \frac{d{{{\bar{\lambda }}}(\mathbf{F})}}{d{\mathbf{F}}}\right|_{\mathbf{F}= \mathbf{F}_1}:(\mathbf{F}_2 -\mathbf{F}_1) \end{aligned}$$

(4.21)

which proves the convexity.

Appendix 2 Lebedev quadrature

Six-point quadrature

$$\begin{aligned} \frac{1}{4\pi } \int _{S^2} f\, da = \frac{1}{6} \sum _{i=1}^6 f(\mathbf{a}_i^{(1)}), \quad \mathbf{a}_i^{(1)}=\{\pm 1,0,0\}, \{0\pm 1,0\}, \{0,0,\pm 1\} \end{aligned}$$

(4.22)

14-point quadrature

$$\begin{aligned} \frac{1}{4\pi } \int _{S^2} f\, da= & {} \frac{1}{15} \sum _{i=1}^6 f(\mathbf{a}_i^{(1)}) + \frac{3}{40} \sum _{i=1}^8 f(\mathbf{a}_i^{(3)}) \end{aligned}$$

(4.23)

$$\begin{aligned} \mathbf{a}_i^{(3)}= & {} \{ \pm \frac{ \sqrt{3}}{3} , \pm \frac{\sqrt{3}}{3} , \pm \frac{\sqrt{3}}{3}\}, \{\mp \frac{\sqrt{3}}{3} \pm \frac{\sqrt{3}}{3},\pm \frac{\sqrt{3}}{3}\}, \{\pm \frac{\sqrt{3}}{3}, \mp \frac{\sqrt{3}}{3},\pm \frac{\sqrt{3}}{3}\}, \{ \pm \frac{\sqrt{3}}{3}, \pm \frac{\sqrt{3}}{3}, \mp \frac{\sqrt{3}}{3} \} \end{aligned}$$

(4.24)

Appendix 3 Response function of the GOH model with recruitment

The derivative of $\mathbf{S}_f$ in (3.17) to $\frac{\mathbf{C}}{2}$ yields the material tangent tensor

$$\begin{aligned} \begin{aligned} {\mathbb {D}}_f&= \sum _{i=1}^2 \biggl [ 2\mu _2 \frac{d{ \left( e^{ \gamma ({\bar{I}}_k-1)^2 } ({\bar{I}}_k-1) \right) }}{d{{\bar{I}}_\kappa }}\mathbf{H}_i\otimes \mathbf{H}_i \\&+ 2\mu _2 e^{ \gamma ({\bar{I}}_k-1)^2 } ({\bar{I}}_k-1) \left( \kappa {\mathbb {H}} + (1-2\kappa ) \frac{ d^2{{\bar{J}}}_i }{d J^2_i } \mathbf{M}_i\otimes \mathbf{M}_i \otimes \mathbf{M}_i\otimes \mathbf{M}_i \right) \biggr ] \end{aligned} \end{aligned}$$

(4.25)

where

$$\begin{aligned} \mathbf{H}_i = \kappa \left( \sum _{I=1}^2 \frac{d{{({{\bar{\lambda }}}^2(\mathbf{E}'_I))}}}{d{ (\lambda ^2(\mathbf{E}'_I))}} \mathbf{E}'_I\otimes \mathbf{E}'_I \right) + (1-2\kappa ) \frac{d{{{\bar{J}}}_i}}{d{ J_i}} \mathbf{M}_i\otimes \mathbf{M}_i, \ {\mathbb {H}} = \sum _{I=1}^2 \frac{d^2({{\bar{\lambda }}}^2(\mathbf{E}'_I))}{ d(\lambda ^2(\mathbf{E}'_I))^2} \mathbf{E}'_I\otimes \mathbf{E}'_I\otimes \mathbf{E}'_I\otimes \mathbf{E}'_I \end{aligned}$$

(4.26)

Response function from the alternative approximation of ${{\bar{I}}}_\kappa$. The stress function arising from the representation (3.18) reads

$$\begin{aligned} \mathbf{S}_f = \sum _{i=1}^2 \mu _{2} e^{ \gamma ({\bar{I}}_{ki}-1)^2 } ({\bar{I}}_{ki}-1) \left( \kappa _i {\bar{\mathbf{I}}} + (1-3\kappa _i) \frac{d{\bar{J_i}}}{d{ J_i}} \mathbf{M}_i\otimes \mathbf{M}_i \right) \end{aligned}$$

(4.27)

where ${\bar{\mathbf{I}}}:= \frac{d{{\bar{I}}_1}}{d{\mathbf{C}}} = \sum _{I=1}^3 \frac{f(\lambda _I) f'(\lambda _I)}{\lambda _I} \mathbf{N}_I\otimes \mathbf{N}_I$, and $\{\mathbf{N}_I\}_{I=1}^3$ are the principal vectors. The material tangent is

$$\begin{aligned} \begin{aligned} {\mathbb {D}}_f&= \sum _{i=1}^2 \biggl [ 2\mu _2 \frac{d{ \left( e^{ \gamma ({\bar{I}}_k-1)^2 } ({\bar{I}}_k-1) \right) }}{d{{\bar{I}}_\kappa }}\mathbf{H}_i\otimes \mathbf{H}_i \\&+ 2\mu _2 e^{ \gamma ({\bar{I}}_k-1)^2 } ({\bar{I}}_k-1) \left( \kappa \frac{d{{\bar{\mathbf{I}}}}}{d{\mathbf{C}}} + (1-3\kappa ) \frac{ d^2{{\bar{J}}}_i }{d J^2_i } \mathbf{M}_i\otimes \mathbf{M}_i \otimes \mathbf{M}_i\otimes \mathbf{M}_i \right) \biggr ] \end{aligned} \end{aligned}$$

(4.28)

where $\mathbf{H}_i = \kappa {\bar{\mathbf{I}}} + (1-3\kappa ) \frac{d{{{\bar{J}}}_i}}{d{ J_i}} \mathbf{M}_i\otimes \mathbf{M}_i$. The fourth-order tensor $\frac{d{{\bar{\mathbf{I}}}}}{d{\mathbf{C}}}$ can be obtained by a straightforward computation.

$$\begin{aligned} \begin{aligned} \frac{d{ {{\bar{\mathbf{I}}}}}}{d{\mathbf{C}}}&= \sum _{I=1}^3 \frac{1}{2\lambda _I} \frac{d{ }}{d{\lambda _I}}\left( \frac{f( {{\bar{\lambda }}}_I) f'({{\bar{\lambda }}}_I)}{\lambda _I} \right) \mathbf{N}_I\otimes \mathbf{N}_I \otimes \mathbf{N}_I\otimes \mathbf{N}_I \\&+ \sum _{I=1}^3 \sum _{J>I} \frac{1}{2( \lambda _I^2 - \lambda _J^2)} \left( \frac{f(\lambda _I)f'(\lambda _I)}{\lambda _I} - \frac{f(\lambda _J)f'(\lambda _J)}{\lambda _J} \right) \varvec{\Psi }_{IJ} \otimes \varvec{\Psi }_{IJ} \\ \end{aligned} \end{aligned}$$

(4.29)

where $\varvec{\Psi }_{IJ} = \mathbf{N}_I\otimes \mathbf{N}_J + \mathbf{N}_J\otimes \mathbf{N}_I$. In the limiting case of $\lambda _I \rightarrow \lambda _J$:

$$\begin{aligned} \frac{1}{\lambda _I^2 - \lambda _I^2} \left( \frac{{{\bar{\lambda }}}_I}{\lambda _I} \frac{d{ {{\bar{\lambda }}}_I }}{d{\lambda _I}} - \frac{{{\bar{\lambda }}}_J}{\lambda _J} \frac{d{ {{\bar{\lambda }}}_J }}{d{\lambda _J}} \right) = \frac{1}{2\lambda _I} \frac{d{ }}{d{\lambda _I}}\left( \frac{f( {{\bar{\lambda }}}_I) f'({{\bar{\lambda }}}_I)}{\lambda _I} \right) \end{aligned}$$

(4.30)

Volumetric-deviatoric splitting. In finite element implementation, a multiplicative volumetric-deviatoric splitting is applied to the deformation gradient. The isochoric factor ${\tilde{\mathbf{F}}} = J^{-\frac{1}{3}}\mathbf{F}$ is passed to the strain energy function. A volumetric term U(J) is added to energy function, giving

$$\begin{aligned} W(\mathbf{F}) = \frac{\mu _1}{2} ( {\text {tr}}{{\tilde{\mathbf{C}}}} -3 ) + W_f({{\tilde{\mathbf{C}}}}) + U(J), \quad {{\tilde{\mathbf{C}}}} = J^{-\frac{2}{3}}\mathbf{C}\end{aligned}$$

(4.31)

Let ${{\tilde{\varvec{\sigma }}}}$ and the $\tilde{{\mathbb {C}}}$ be the algorithm Cauchy stress and material tangent tensor arising from $W_f+W_g$, namely the response evaluating at $\mathbf{F}= {{\tilde{\mathbf{F}}}}$. The actual Cauchy stress and spatial material tensor are obtained through deviatoric projections (Simo et al. 1985) of the algorithm response, adding the contributions form U(J):

$$\begin{aligned} \begin{aligned} \mathbf{\sigma }&= {{\tilde{\varvec{\sigma }}}}_{\text{dev}} + U'(J) \varvec{1}, \quad \text {where} \ {{\tilde{\varvec{\sigma }}}}_{\text{dev}} = {\mathbb {I}}_{\text{dev}} {{\tilde{\varvec{\sigma }}}} \\ {\mathbb {C}}&= {\mathbb {I}}_{\text{dev}} \tilde{{\mathbb {C}}} {\mathbb {I}}_{\text{dev}} -\frac{2}{3}{{\tilde{\varvec{\sigma }}}}_{\text{dev}}\otimes \varvec{1}- \frac{2}{3} \varvec{1}\otimes {{\tilde{\varvec{\sigma }}}}_{\text{dev}} + \frac{2{\text {tr}}{{\tilde{\varvec{\sigma }}}}}{3} {\mathbb {I}}_{\text{dev}} -2 U'(J) {\mathbb {I}} + (JU'(J))' \varvec{1}\otimes \varvec{1}\end{aligned} \end{aligned}$$

(4.32)

where ${\mathbb {I}}$ and $\varvec{1}$ are the fourth- and second-order identity tensors, ${\mathbb {I}}_{\text{dev}} = {\mathbb {I}} -\frac{1}{3}\varvec{1}\otimes \varvec{1}$ is the deviatoric projector.