Anisotropic Constitutive Modeling of Compressible Biological Tissue

Fuzhang Zhao

APD Optima Study, Lake Forest, USA.

DOI: 10.4236/apm.2022.125027                

Abstract

The anisotropic continuum stored energy density (ACSED) functional is applied for accurate constitutive modeling of biological tissues and finite element implementation without the isochoric—volumetric split, the anisotropic—isotropic split, or the anisotropic invariant split. Related stress and elasticity tensors in the reference and current configurations are worked out. A new kinematic model is derived based on the tangent Poisson’s ratio as a cubic polynomial function of stretch. The ACSED model, along with the kinematic model, accurately fits uniaxial extension test data for compressible human skin, bovine articular cartilage, and human aorta samples.

Keywords

ACSED Functional, Biological Tissue, Kinematic Model, Nonlinear Elastic Deformation, Tangent Poisson’s Ratio

1. Introduction

Biological tissues, composed of collagen fibers and extracellular matrices, can be treated as nonlinear elastic fiber-reinforced composites. Significant progress in histological structures, mechanical characterizations, anisotropic constitutive models, and finite element implementations has been made for biological tissues. Finite element commercial packages play an important role in fundamentally solving practical boundary value problems. Essentially, we can better understand relations among histological compositions, anisotropic growth, mechanical properties, and performed functions of biological tissues. Continuum nonlinear solid mechanics was documented by Holzapfel (2000) [1], Continuum biomechanics was summarized by Humphrey (2003) [2], the theory of mixed finite element methods has been overviewed by Wriggers (2009) [3], anisotropic constitutive models have been reviewed by Chagnon, Rebouah, and Favier (2015) [4], and multi-scale structural modeling and mechanobiology have been reviewed by Lanir (2017) [5]. Despite extensive progress in the past few decades, many challenges in the constitutive modeling of biological tissues and its finite element implementation remain.

As a mathematical convenience, the three major assumptions of existing anisotropic constitutive models implemented into finite element commercial packages are the isochoric—volumetric split, the anisotropic—isotropic split, and the anisotropic invariant split. First, for the incompressibility condition, the multiplicative decomposition of the deformation gradient into isochoric and volumetric parts is enforced, establishing the isochoric—volumetric split. Second, existing anisotropic stored energy density functionals are generally divided into a deformation part for fiber and a deformation part for matrix, constructing the anisotropic—isotropic split. Third, irreducible invariants are directly applied and widely accepted, popularizing the anisotropic invariant split.

With the isochoric—volumetric split, along with the anisotropic—isotropic split, anisotropic stored energy functionals should be divided into four parts. However, the invalidity of the isochoric—volumetric split for anisotropic nonlinear elastic deformations has been irrefutably demonstrated by Cowin (1990) [6], Sansour (2008) [7], Ní Annaidh et al. (2013) [8], Wang and Liu (2018) [9], Jin and Stanciulescu (2019) [10], and others. Thus, most anisotropic stored energy functionals are generally divided into three parts: the isochoric-anisotropic deformation for fiber, the isochoric-isotropic deformation for matrix, and the volumetric-isotropic deformation for matrix. If the invalid volumetric-anisotropic deformation part for fiber is removed, the invalid isochoric-anisotropic deformation part for fiber should also be removed but be kept there, making the current implemented theory self-contradictory. The volumetric-isotropic deformation part for matrix is constructed as a function of volume ratio only, setting up a partial or incomplete or truncated formulation since pressure is a function of distance, area, and volume in general. Furthermore, the invalidity of the isochoric—volumetric split for isotropic nonlinear elastic deformations is demonstrated by Zhao (2020) [11]. With the anisotropic—isotropic split, the rule of a compatible deformation between collagen fibers and extracellular matrices is violated and other dilemmas have been addressed by Zhao (2019) [12]. The theory of representations for tensor polynomials has been extensively developed, as described by Spencer (1971) [13]. Tensor function representations have been developed in complete and irreducible invariant forms. The anisotropic constitutive framework is generally built from direct applications of irreducible invariants, which may be too general to reveal specific forms of constitutive models for nonlinear elastic deformations as indicated by Zheng (1994) [14].

In finite element implementations, the perfect incompressibility causes a numerical singularity. In displacement-based finite element methods, a locking occurs as Poisson’s ratio gets closer to 0.5. As a remedy, the reduced integration scheme is used. Spurious frequencies, however, are introduced due to inexact numerical integration of element stiffness matrices as emphasized by Bathe (1996) [15]. To overcome the dilemma, a mixed finite element formulation with independent fields of displacement and stress as unknowns was coined by Pian (1964) [16]. The incompressibility condition is weakened as the near incompressibility condition and the locking is alleviated at the cost of additional pressure degrees of freedom. The pressure field is determined by equilibrium equations and boundary conditions rather than constitutive and kinematic models. Different mixed finite element formulations for nearly incompressible nonlinear elastic deformations can be found in the articles by Oden and Key (1970) [17], by Häggblad and Sundberg (1983) [18], and by Zdunek and Rachowicz (2017) [19] among others. Furthermore, mixed finite element formulations are not stable in all displacement and pressure fields and spurious pressure oscillations are found by Legrain, Moësa, and Huerta (2008) [20].

In practical applications of finite element methods for biological tissues, small changes in the Poisson’s ratio near 0.5 result in significant overestimation or underestimation in stress which is probably due to spurious pressure oscillations. For small stretches, stresses can be predicted by more than or less than an order of magnitude, causing unphysical results as demonstrated by Ní Annaidh et al. (2013) [8]. Unphysical prediction of a sphere deformed into a larger sphere instead of an expected ellipsoid under hydrostatic tension was found by Vergori et al. (2013) [21], undesired volume growth in uniaxial extension was discovered by Helfenstein et al. (2010) [22], and a significant unphysical auxetic response for non-auxetic materials was also predicted by Horgan and Murphy (2020) [23]. Some improvement efforts within the three major splits have been made by Limbert and Taylor (2002) [24], Nolan et al. (2014) [25], Pierrat et al. (2016) [26], Cai et al. (2016) [27], Murphy and Rogerson (2018) [28], and Fereidoonnezhad, O’Connor, and McGarry (2020) [29] among others.

Thus, it is desirable to formulate a general constitutive model for biological tissues without the three major splits. The ACSED functional for ith fiber family, ${\Psi }_i$, based on the energy balance between the stored energy density and stress work done equivalents in continuum mechanics under isothermal processes is covariantly formulated by Zhao (2018) [30]

$P_i:F=S_i:C={\tau }_i:I={\sigma }_i:JI=2{\Psi }_i\left(I_{1,i},I_{2,i},I_{3,i}\right),$ (1)

where $P_i$, $S_i$, ${\tau }_i$, and ${\sigma }_i$ are the first Piola-Kirchhoff stress tensor, second Piola-Kirchhoff stress tensor, Kirchhoff stress tensor, and Cauchy stress tensor for ith fiber family, respectively. The corresponding conjugates $F$, $C$, $I$, and $JI$ are the deformation gradient tensor, the right Cauchy—Green tensor, unit tensor, and product of volume ratio and unit tensor, respectively. For $i=1,2,\cdots ,n$, the index n is the total number of fiber families. The anisotropic invariants or principal invariant components (PIC) of the right Cauchy—Green tensor $C=F^{\text{T}}F$, $I_{\mathrm{1,}i}$, $I_{\mathrm{2,}i}$, and $I_{\mathrm{3,}i}$, are defined by

$I_{1,i}=\text{tr}\left(C{A}_{0,i}\right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{I}_{2,i}=\frac{1}{2}\left[I_1{I}_{1,i}-\text{tr}\left(C^2{A}_{0,i}\right)\right],\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{I}_{3,i}=\det C,$ (2)

where $I_1=\text{tr}C$ is the first principal invariant for isotropic hyperelastic materials and $I_{3,i}=I_3$ is noticed. Biological tissues are usually modeled as isotropic, transversely isotropic, and orthotropic hyperelastic materials. For the orthotropic symmetry as an example, the structure tensors $A_{0,i}$ are defined as

$A_{0,i}=a_{0,i}\otimes a_{0,i},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}i=1,2,3,$ (3)

and $a_{0,i}$ is a direction unit vector for fiber i in reference configuration denoted by the subscript 0 and the matrix operators, “tr”, and “det”, denote trace and determinant operations, respectively.

The first-order derivatives of PICs with respect to $C$ are given by

$\frac{\partial I_{1,i}}{\partial C}=A_{0,i},$ (4)

$\frac{\partial I_{2,i}}{\partial C}=\frac{1}{2}\left(I_1{A}_{0,i}+I_{1,i}I-C{A}_{0,i}-A_{0,i}C\right),$ (5)

$\frac{\partial I_3}{\partial C}=I_3{C}^{-1}.$ (6)

The second-order derivatives of PICs with respect to $C$ are further given by

$\frac{{\partial }^2{I}_{1,i}}{\partial C^2}=\mathbb{O},$ (7)

$\frac{{\partial }^2{I}_{\mathrm{2,}i}}{\partial C^2}=\frac{1}{2}\left(A_{\mathrm{0,}i}\otimes I+I\otimes A_{\mathrm{0,}i}-A_{\mathrm{0,}i}\odot I-I\odot A_{\mathrm{0,}i}\right)\mathrm{,}$ (8)

$\frac{{\partial }^2{I}_3}{\partial C^2}=I_3\left(C^{-1}\otimes C^{-1}-C^{-1}\odot C^{-1}\right)\mathrm{,}$ (9)

where $\mathbb{O}$ is the fourth-order null tensor.

The second Piola—Kirchhoff stress tensor and the Cauchy stress tensor fot ith fiber family desired in finite element implementations can be generally derived based on (1), (4), (5), and (6). The second Piola—Kirchhoff stress tensor in reference configuration turns out to be

$\begin{array}{c}{S}_i=2{\displaystyle \underset{j=1}{\overset{3}{\sum }}}\frac{\partial {\Psi }_i}{\partial I_{j,i}}\frac{\partial I_{j,i}}{\partial C}\\ =2\frac{\partial {\Psi }_i}{\partial I_{1,i}}{A}_{0,i}+\frac{\partial {\Psi }_i}{\partial I_{2,i}}\left(I_1{A}_{0,i}+I_{1,i}I-C{A}_{0,i}-A_{0,i}C\right)+2I_3\frac{\partial {\Psi }_i}{\partial I_3}{C}^{-1},\end{array}$ (10)

and the Cauchy stress tensor, ${\sigma }_i$, in current configuration can be converted by the push-forward operation of the second Piola—Kirchhoff stress tensor (10)₂ and readily obtained as

$\begin{array}{c}{\sigma }_i=\frac{1}{J}F{S}_i{F}^{\text{T}}\\ =2\frac{\partial {\Psi }_i}{J\partial I_{\mathrm{1,}i}}{A}_i+\frac{\partial {\Psi }_i}{J\partial I_{\mathrm{2,}i}}\left(I_1{A}_i+I_{\mathrm{1,}i}B-B{A}_i-A_{i}B\right)+2J\frac{\partial {\Psi }_i}{\partial I_3}I\mathrm{,}\end{array}$ (11)

where $A_i=F{A}_{\mathrm{0,}i}{F}^{\text{T}}$ is the structure tensor in current configuration, $B=F{F}^{\text{T}}$ is the left Cauchy—Green tensor, and $J=\sqrt{I_3}$ is the volume ratio or Jacobian.

Substituting (10)₂ or (11)₂ into (1) yields the partial differential equation for ${\Psi }_i$

${\Psi }_i=I_{1,i}\frac{\partial {\Psi }_i}{\partial I_{1,i}}+2I_{2,i}\frac{\partial {\Psi }_i}{\partial I_{2,i}}+3I_3\frac{\partial {\Psi }_i}{\partial I_3}.$ (12)

With the Lie group method, the general solution normalized by $I_{\mathrm{1,}i}$ reads

${\Psi }_i/I_{\mathrm{1,}i}=f\left(I_{\mathrm{2,}i}^{\frac{1}{2}}/I_{\mathrm{1,}i}\right)+g\left(I_3^{\frac{1}{3}}/I_{\mathrm{1,}i}\right)\mathrm{,}$ (13)

where functions f and g can be fixed by curvatures of deformations based on differential geometry. For normal, shear, and ellipsoidal stretches, the ACSED functional has been obtained by Zhao (2018) [30]

${\Psi }_i=c_{1,i}{I}_{1,i}+c_{2,i}\sqrt{I_{2,i}}+c_{3,i}\frac{I_{1,i}^{3c_{4,i}+1}}{I_3^{c_{4,i}}},$ (14)

where the constitutive parameters for ith fiber family, $c_{\mathrm{1,}i}$, $c_{\mathrm{2,}i}$, $c_{\mathrm{3,}i}$, and $c_{\mathrm{4,}i}$, are determined by experimental tests.

For biological tissues with simple kinematic relations, constant Poisson’s ratios ${\nu }_{1j}$ are used by Blatz and Ko (1962) [31],

${\nu }_{1j}=-\frac{\ln {\lambda }_j}{\ln {\lambda }_1}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{or}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\lambda }_j={\lambda }_1^{-{\nu }_{1j}},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}j=2,3.$ (15)

For biological tissues with complex kinematic relations, stretch dependent or tangent Poisson’s ratios are used by Fortes and Nogueira (1989) [32]

${\nu }_{1j}\left({\lambda }_1\right)=-\frac{\text{d}\left(\ln {\lambda }_j\right)}{\text{d}\left(\ln {\lambda }_1\right)}=-\frac{{\lambda }_1\text{d}{\lambda }_j}{{\lambda }_j\text{d}{\lambda }_1},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}j=2,3,$ (16)

where “ln” denotes natural logarithmic function.

The major objectives are to apply the ACSED functional, to develop the kinematic model based on tangent Poisson’s ratio, to derive the stress and elasticity tensors in reference and current configurations, and to fit uniaxial extension and kinematic test data for compressible human skin, bovine articular cartilage, and human aorta samples.

This paper is organized as follows. In Section 2, the stress—stretch relation and the kinematic relation in uniaxial extension mode are derived for constitutive modeling biological tissues. Related stress and elasticity tensors are derived for finite element implementation. In Section 3, the ACSED constitutive model, along with the new kinematic model, is applied to fit uniaxial extension test data for different compressible biological tissues. In Section 4, modeling of biological tissues, the isochoric—volumetric split, and the predictive ACSED model are discussed. In Section 5, conclusions are drawn and the desired future work is provided.

2. ACSED Constitutive Model

2.1. Stress—Stretch Relation and Kinematic Relation

For ith fiber family, $P_i$ is related to $F$ by

$P_i^{\text{T}}=\frac{\partial {\Psi }_i}{\partial F}=\frac{\partial {\Psi }_i}{\partial I_{\mathrm{1,}i}}\frac{\partial I_{\mathrm{1,}i}}{\partial F}+\frac{\partial {\Psi }_i}{\partial I_{\mathrm{2,}i}}\frac{\partial I_{\mathrm{2,}i}}{\partial F}+\frac{\partial {\Psi }_i}{\partial I_3}\frac{\partial I_3}{\partial F}\mathrm{.}$ (17)

With all three PICs generally treated as variables, the first-order derivatives of the ACSED functional are worked out as

$\frac{\partial {\Psi }_i}{\partial I_{1,i}}=c_{1,i}+c_{3,i}\left(3c_{4,i}+1\right)\frac{I_{1,i}^{3c_{4,i}}}{I_3^{c_{4,i}}},$ (18)

$\frac{\partial {\Psi }_i}{\partial I_{2,i}}=\frac{c_{2,i}}{2\sqrt{I_{2,i}}},$ (19)

$\frac{\partial {\Psi }_i}{\partial I_3}=-c_{3,i}{c}_{4,i}\frac{I_{1,i}^{3c_{4,i}+1}}{I_3^{c_{4,i}+1}},$ (20)

and the second-order derivatives of the ACSED functional are further obtained as

$\frac{{\partial }^2{\Psi }_i}{\partial I_{1,i}^2}=3c_{3,i}{c}_{4,i}\left(3c_{4,i}+1\right)\frac{I_{1,i}^{3c_{4,i}-1}}{I_3^{c_{4,i}}},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{{\partial }^2{\Psi }_i}{\partial I_{2,i}^2}=-\frac{c_{2,i}}{4\sqrt{I_{2,i}^3}},$ (21)

$\frac{{\partial }^2{\Psi }_i}{\partial I_3^2}=c_{3,i}{c}_{4,i}\left(c_{4,i}+1\right)\frac{I_{1,i}^{3c_{4,i}+1}}{I_3^{c_{4,i}+2}},$ (22)

$\frac{{\partial }^2{\Psi }_i}{\partial I_{1,i}\partial I_{2,i}}=\frac{{\partial }^2{\Psi }_i}{\partial I_{2,i}\partial I_3}=0,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{{\partial }^2{\Psi }_i}{\partial I_{1,i}\partial I_3}=-c_{3,i}{c}_{4,i}\left(3c_{4,i}+1\right)\frac{I_{1,i}^{3c_{4,i}}}{I_3^{c_{4,i}+1}}.$ (23)

Mechanical characterizations of anisotropic hyperelastic materials have been accurately and effectively conducted by uniaxial extension tests. Thus, the ACSED constitutive model in uniaxial extension mode is derived from the Equations (2), (17), (18), (19), and (20).

For intima, media, and advantitia layers of arteries, the structure tensors for two families of preferred fiber directions with angles $\pm \theta$ are given respectively by

${{\rm A}}_{\mathrm{0,1}}\left(\theta \right)=\left(\begin{array}{ccc}{\cos }^2\theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& \cos \theta \sin \theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& 0\\ \cos \theta \sin \theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& {\sin }^2\theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& 0\\ \mathrm{0\; <mtext>\; \hspace{0.17em}\; </mtext>}& \mathrm{0\; <mtext>\; \hspace{0.17em}\; </mtext>}& 0\end{array}\right)$ (24)

and

$A_{\mathrm{0,2}}\left(-\theta \right)=\left(\begin{array}{ccc}{\cos }^2\theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& -\cos \theta \sin \theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& 0\\ -\cos \theta \sin \theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& {\sin }^2\theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& 0\\ \mathrm{0\; <mtext>\; \hspace{0.17em}\; </mtext>}& \mathrm{0\; <mtext>\; \hspace{0.17em}\; </mtext>}& 0\end{array}\right)$ (25)

where $\pm \theta$ measure in-plane symmetric angles between fiber direction and principal stretch direction. In uniaxial extension tests of arteries, the related tensors are $C=\text{diag}\left[{\lambda }_1^2\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>}{\lambda }_2^2\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>}{\lambda }_3^2\right]$, $C^2=\text{diag}\left[{\lambda }_1^4\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>}{\lambda }_2^4\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>}{\lambda }_3^4\right]$. The related PICs are obtained as

$I_{1,i}={\lambda }_1^2{\cos }^2\theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}+{\lambda }_2^2{\sin }^2\theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>},$ (26)

$I_{2,i}=\frac{1}{2}\left({\lambda }_1^2{\lambda }_2^2+{\lambda }_2^2{\lambda }_3^2{\sin }^2\theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}+{\lambda }_3^2{\lambda }_1^2{\cos }^2\theta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\right),$ (27)

$I_3={\lambda }_1^2{\lambda }_2^2{\lambda }_3^2,$ (28)

and their derivatives with respect to ${\lambda }_1$ are worked out as

$\frac{\partial I_{1,i}}{\partial {\lambda }_1}=2{\lambda }_1{\cos }^2\theta ,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{\partial I_{2,i}}{\partial {\lambda }_1}={\lambda }_1\left({\lambda }_2^2+{\lambda }_3^2{\cos }^2\theta \right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\frac{\partial I_3}{\partial {\lambda }_1}=2{\lambda }_1{\lambda }_2^2{\lambda }_3^2.$ (29)

For two fiber families with $\pm \theta$, substituting the derivatives (18), (19), and (20), PICs (26), (27), and (28) and derivatives of PICs (29) into (17) and summing up stresses yields the total nominal stress equation

$\begin{array}{l}{P}_{11}={\displaystyle \underset{i=1}{\overset{2}{\sum }}}{P}_{11,i}=4c_{1,i}{\lambda }_1{\cos }^2\theta +\frac{c_{2,i}{\lambda }_1({\lambda }_2^2+{\lambda }_3^2{\cos }^2\theta )}{\sqrt{\frac{1}{2}\left({\lambda }_1^2{\lambda }_2^2+{\lambda }_2^2{\lambda }_3^2{\sin }^2\theta +{\lambda }_3^2{\lambda }_1^2{\cos }^2\theta \right)}}\\ +4c_{3,i}\frac{{\left({\lambda }_1^2{\cos }^2\theta +{\lambda }_2^2{\sin }^2\theta \right)}^{3c_{4,i}}}{{\left({\lambda }_1{\lambda }_2{\lambda }_3\right)}^{2c_{4,i}}}\left[\left(2c_{4,i}+1\right){\lambda }_1{\cos }^2\theta -c_{4,i}\frac{{\lambda }_2^2}{{\lambda }_1}{\sin }^2\theta \right],i=1,2.\end{array}$ (30)

The corresponding stress normalization condition reads

$4c_{1,i}{\cos }^2\theta +c_{2,i}\left(1+{\cos }^2\theta \right)+4c_{3,i}\left[\left(2c_{4,i}+1\right){\cos }^2\theta -c_{4,i}{\sin }^2\theta \right]=0.$ (31)

For one fiber family with $\theta =0$, the constitutive model (30) combined with the stress normalization condition (31) is simplified as

$P_{11,1}=2c_{1,1}\left({\lambda }_1-1\right)+c_{2,1}\left(\sqrt{\frac{{\lambda }_2^2+{\lambda }_3^2}{2}}-1\right)+2c_{3,1}\left(2c_{4,1}+1\right)\left[\frac{{\lambda }_1^{4c_{4,1}+1}}{{\left({\lambda }_2{\lambda }_3\right)}^{2c_{4,1}}}-1\right].$ (32)

In the constitutive Equations (30) and (32), transverse stretches ${\lambda }_2$ and ${\lambda }_3$ as a function of longitudinal stretch ${\lambda }_1$ are needed. Thus, kinematic relations, ${\lambda }_2\left({\lambda }_1\right)$ and ${\lambda }_3\left({\lambda }_1\right)$, need to be modeled. Integrating the equation (16) yields the general kinematic model

${\lambda }_j=\exp \left[-{\displaystyle \int }\frac{{\nu }_{1j}\left({\lambda }_1\right)}{{\lambda }_1}\text{d}{\lambda }_1\right]\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}j=\mathrm{2,3,}$ (33)

where “exp” or “e” denotes an exponential function.

In order to relate transverse stretches ${\lambda }_2$ and ${\lambda }_3$ with longitudinal stretch ${\lambda }_1$, a cubic polynomial is assumed for tangent Poisson’s ratios as

${\nu }_{1j}\left({\lambda }_1\right)=k_{0j}+k_{1j}{\lambda }_1+k_{2j}{\lambda }_1^2+k_{3j}{\lambda }_1^3,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}j=2,3,$ (34)

and substituting (34) into (33), integrating, and applying the normalization condition of ${{\lambda }_j|}_{{\lambda }_1=1}=1$ yields the general kinematic model

${\lambda }_j={\lambda }_1^{-k_{0j}}{\text{e}}^{-\left[k_{1j}\left({\lambda }_1-1\right)+k_{2j}\left({\lambda }_1^2-1\right)/2+k_{3j}\left({\lambda }_1^3-1\right)/3\right]},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}j=2,3,$ (35)

where parameters $k_{0j}$, $k_{1j}$, $k_{2j}$, and $k_{3j}$ need to be determined by kinematic characterization in uniaxial extension tests.

2.2. Anisotropic Stress Tensors

For the ACSED functional with ith fiber family (14), the second Piola-Kirchhoff stress tensor, $S_i$, can be obtained as

$\begin{array}{c}{S}_i=2c_{1,i}{A}_{0,i}+c_{2,i}\frac{\left(I_1{A}_{0,i}+I_{1,i}I-C{A}_{0,i}-A_{0,i}C\right)}{2\sqrt{I_{2,i}}}\\ +2c_{3,i}\frac{I_{1,i}^{3c_{4,i}}\left[\left(3c_{4,i}+1\right)A_{0,i}-c_{4,i}{I}_{1,i}{C}^{-1}\right]}{I_3^{c_{4,i}}},\end{array}$ (36)

and the corresponding Cauchy stress tensor, ${\sigma }_i$, in current configuration can be converted by the push-forward operation (11)₁ of the second Piola-Kirchhoff stress tensor (36) as

$\begin{array}{c}{\sigma }_i=2c_{1,i}\frac{A_i}{J}+c_{2,i}\frac{\left(I_1{A}_i+I_{1,i}B-B{A}_i-A_{i}B\right)}{2\sqrt{I_{2,i}}J}\\ +2c_{3,i}\frac{I_{1,i}^{3c_{4,i}}\left[\left(3c_{4,i}+1\right)A_i-c_{4,i}{I}_{1,i}I\right]}{I_3^{c_{4,i}}J},\end{array}$ (37)

and the total second Piola-Kirchhoff and Cauchy stress tensors for n number of fiber families can respectively be summed up as

$S={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{S}_i,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{and}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\sigma ={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\sigma }_i.$ (38)

2.3. Anisotropic Elasticity Tensors

The fourth-order elasticity tensor in reference configuration can be readily derived by

${ℂ}_i=4{\displaystyle \underset{j=1}{\overset{3}{\sum }}}\left[{\displaystyle \underset{k=1}{\overset{3}{\sum }}}\frac{{\partial }^2{\Psi }_i}{\partial I_{j,i}\partial I_{k,i}}\left(\frac{\partial I_{j,i}}{\partial C}\otimes \frac{\partial I_{k,i}}{\partial C}\right)+\frac{\partial {\Psi }_i}{\partial I_{j,i}}\frac{{\partial }^2{I}_{j,i}}{\partial C^2}\right].$ (39)

Substituting the first-order and second-order derivatives of PICs (4) through (9), simplifying, and rearranging produces

$\begin{array}{c}{ℂ}_i={\Delta }_{\mathrm{1,}i}{A}_{\mathrm{0,}i}\otimes A_{\mathrm{0,}i}+{\Delta }_{\mathrm{2,}i}\left(A_{\mathrm{0,}i}\otimes I+I\otimes A_{\mathrm{0,}i}\right)+{\Delta }_{\mathrm{3,}i}\left(A_{\mathrm{0,}i}\otimes H_{\mathrm{0,}i}+H_{\mathrm{0,}i}\otimes A_{\mathrm{0,}i}\right)\\ +{\Delta }_{\mathrm{4,}i}\left(A_{\mathrm{0,}i}\otimes C^{-1}+C^{-1}\otimes A_{\mathrm{0,}i}\right)+{\Delta }_{\mathrm{5,}i}I\otimes I+{\Delta }_{\mathrm{6,}i}\left(I\otimes H_{\mathrm{0,}i}+H_{\mathrm{0,}i}\otimes I\right)\\ +{\Delta }_{\mathrm{7,}i}{H}_{\mathrm{0,}i}\otimes H_{\mathrm{0,}i}+{\Delta }_{\mathrm{8,}i}\left(I\otimes C^{-1}+C^{-1}\otimes I\right)+{\Delta }_{\mathrm{9,}i}\left(C^{-1}\otimes H_{\mathrm{0,}i}+H_{\mathrm{0,}i}\otimes C^{-1}\right)\\ +{\Delta }_{\mathrm{10,}i}{C}^{-1}\otimes C^{-1}+{\Delta }_{\mathrm{11,}i}\left(A_{\mathrm{0,}i}\odot I+I\odot A_{\mathrm{0,}i}\right)+{\Delta }_{\mathrm{12,}i}{C}^{-1}\odot C^{-1}\mathrm{,}\end{array}$ (40)

where $H_{\mathrm{0,}i}=0.5\left(C{A}_{\mathrm{0,}i}+A_{\mathrm{0,}i}C\right)$ is an interim second-order tensor for abbreviation and the twelve coefficients ${\Delta }_{\mathrm{1,}i}\mathrm{,}{\Delta }_{\mathrm{2,}i}\mathrm{,}\cdots \mathrm{,}{\Delta }_{\mathrm{12,}i}$ for the ACSED functional, using (21), (22), and (23), are given by

${\Delta }_{1,i}=12c_{3,i}{c}_{4,i}\left(3c_{4,i}+1\right)\frac{I_{1,i}^{3c_{4,i}-1}}{I_3^{c_{4,i}}}-\frac{c_{2,i}{I}_1^2}{I_{2,i}^{1.5}},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\Delta }_{2,i}=\frac{c_{2,i}}{I_{2,i}^{1.5}}-\frac{c_{2,i}{I}_1{I}_{1,i}}{4I_{2,i}^{1.5}},$ (41)

${\Delta }_{3,i}=\frac{c_{2,i}{I}_1}{2I_{2,i}^{1.5}},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\Delta }_{4,i}=-4c_{3,i}{c}_{4,i}\left(3c_{4,i}+1\right)\frac{I_{1,i}^{3c_{4,i}}}{I_3^{c_{4,i}}},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\Delta }_{5,i}=-\frac{c_{2,i}{I}_{1,i}^2}{4I_{2,i}^{1.5}},$ (42)

${\Delta }_{6,i}=\frac{c_{2,i}{I}_{1,i}}{2I_{2,i}^{1.5}},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\Delta }_{7,i}=-\frac{c_{2,i}}{I_{2,i}^{1.5}},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\Delta }_{8,i}={\Delta }_{9,i}=0,$ (43)

${\Delta }_{10,i}=4c_{3,i}{c}_{4,i}^2\frac{I_{1,i}^{3c_{4,i}+1}}{I_3^{c_{4,i}}},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\Delta }_{11,i}=-\frac{c_{2,i}}{I_{2,i}^{0.5}},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\Delta }_{12,i}=4c_{3,i}{c}_{4,i}\frac{I_{1,i}^{3c_{4,i}+1}}{I_3^{c_{4,i}}}.$ (44)

The elasticity tensor in current configuration, ${\text{c}}_i^{|}$, can be converted from the elasticity tensor in reference configuration, ${ℂ}_i$, by the following push-forward operation of (40)

${\text{c}}_i^{|}=J^{-1}\left(F\odot F\right)\mathrm{<mo>\; :\; </mo>}{ℂ}_i\mathrm{<mo>\; :\; </mo>}\left(F^{\text{T}}\odot F^{\text{T}}\right)\mathrm{.}$ (45)

With the properties of double contraction operations between fourth-order tensors [11], the general coupled fourth-order elasticity tensor in current configuration (45) can be converted as

$\begin{array}{c}{\text{c}}_i^{|}={\delta }_{\mathrm{1,}i}{A}_i\otimes A_i+{\delta }_{\mathrm{2,}i}\left(A_i\otimes B+B\otimes A_i\right)+{\delta }_{\mathrm{3,}i}\left(A_i\otimes H_i+H_i\otimes A_i\right)\\ +{\delta }_{\mathrm{4,}i}\left(A_i\otimes I+I\otimes A_i\right)+{\delta }_{\mathrm{5,}i}B\otimes B+{\delta }_{\mathrm{6,}i}\left(B\otimes H_i+H_i\otimes B\right)\\ +{\delta }_{\mathrm{7,}i}{H}_i\otimes H_i+{\delta }_{\mathrm{8,}i}\left(I\otimes B+B\otimes I\right)+{\delta }_{\mathrm{9,}i}\left(I\otimes H_i+H_i\otimes I\right)\\ +{\delta }_{\mathrm{10,}i}I\otimes I+{\delta }_{\mathrm{11,}i}\left(A_i\odot B+B\odot A_i\right)+{\delta }_{\mathrm{12,}i}I\odot I\mathrm{,}\end{array}$ (46)

where $H_i=F{H}_{\mathrm{0,}i}{F}^{\text{T}}=0.5\left(B{A}_i+A_{i}B\right)$ is for abbreviation and the corresponding coefficients ${\delta }_{\mathrm{1,}i}\mathrm{,}{\delta }_{\mathrm{2,}i}\mathrm{,}\cdots \mathrm{,}{\delta }_{\mathrm{12,}i}$ for the ACSED functional in current configuration are given by

${\delta }_{j,i}={\Delta }_{j,i}{J}^{-1},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}j=1,2,\cdots ,12,$ (47)

and the total elasticity tensors for n number of fiber families can be combined as

$ℂ={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{ℂ}_i,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{or}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\text{c}}_{}^{|}={\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\text{c}}_i^{|}.$ (48)

The coupled fourth-order elasticity tensors in both reference and current configurations (40) and (46) are derived for the ACSED functional.

3. Application of ACSED Model

3.1. Modeling of Human Abdominal Skin

Uniaxial extension tests with a full kinematic characterization of compressible human abdominal skin samples have been conducted at the stretch rate of 0.001 s⁻¹ by Wahlsten et al. (2019) [33]. A self-developed graphics digitizer with MATLAB has been used to read out the averaged experimental data. The converted nominal stress—stretch data with kinematic test data have been selected for the ACSED model (32) along with the general kinematic model (35).

The kinematic model (35) fits the test data as

${\lambda }_2={\text{e}}^{\left[-805.158\ln {\lambda }_1+2183.826\left({\lambda }_1-1\right)-981.877\left({\lambda }_1^2-1\right)+194.660\left({\lambda }_1^3-1\right)\right]}\mathrm{,}$ (49)

${\lambda }_3={\text{e}}^{\left[587.242\ln {\lambda }_1-1593.726\left({\lambda }_1-1\right)+720.646\left({\lambda }_1^2-1\right)-145.065\left({\lambda }_1^3-1\right)\right]}\mathrm{.}$ (50)

The stretch dependent volume ratio can be readily obtained as

$J={\lambda }_1^{-216.916}{\text{e}}^{\left[590.100\left({\lambda }_1-1\right)-261.231\left({\lambda }_1^2-1\right)+49.595\left({\lambda }_1^3-1\right)\right]}\mathrm{.}$ (51)

The models for tangent Poisson’s ratios are obtained as

${\nu }_{12}\left({\lambda }_1\right)=805.158-2183.826{\lambda }_1+1963.753{\lambda }_1^2-583.980{\lambda }_1^3\mathrm{,}$ (52)

${\nu }_{13}\left({\lambda }_1\right)=-587.242+1593.726{\lambda }_1-1441.291{\lambda }_1^2+435.196{\lambda }_1^3\mathrm{.}$ (53)

The constitutive parameters have been determined by the trial-and-error-on-digit (TED) method and the linear least square (LLSQ) method combined by Zhao (2019) [12]. The ACSED model (32) accurately fits the nominal stress-stretch test data as shown in Figure 1(a). Kinematic models (49) and (50) have been used to plot Figure 1(b). As shown, the kinematic model with the tangent Poisson’s ratios accurately fits the kinematic test data. The tangent Poisson’s ratio models (52) and (53) have been used to plot Figure 1(c). As shown, the tangent Poisson’s ratio ${\nu }_{12}$ varies between 1.022 and 3.253 and the Poisson’s ratio ${\nu }_{13}$ varies from 0.388 to 1.787 as ${\lambda }_1$ increases from 1.000 to 1.200. In Figure 1(d), both test and model (51) show the volume reduction of about 30.6%.

Figure 1. Comparison between ACSED model and uniaxial extension test of human skin.

3.2. Modeling of Bovine Articular Cartilage

Uniaxial extension tests for bovine articular cartilage of humeral joint samples with full kinematic characterizations have been conducted at the stretch rate of 0.01 s⁻¹ by Woo et al. (1979) [34]. The converted nominal stress—stretch data and kinematic data at 0˚ middle zone have been selected to fit the ACSED model (32) and the kinematic model (35). The kinematic relation and the tangent Poisson’s ratio for lateral contraction read

${\lambda }_2={\lambda }_1^{-55.194}{\text{e}}^{\left[135.311\left({\lambda }_1-1\right)-54.151\left({\lambda }_1^2-1\right)+9.334\left({\lambda }_1^3-1\right)\right]}\mathrm{,}$ (54)

${\nu }_{12}\left({\lambda }_1\right)=55.194-135.311{\lambda }_1+108.302{\lambda }_1^2-28.002{\lambda }_1^3\mathrm{,}$ (55)

and ${\lambda }_3={\lambda }_1^{-{\nu }_{13}}$ with ${\nu }_{13}=2.241$ for theckness reduction.

The ACSED model (32) fits the nominal stress—stretch test data shown in Figure 2(a). The kinematic models (54) and (15)₂ with ${\nu }_{13}=2.241$ depict Figure 2(b). The tangent and constant Poisson’s ratio models (55) depict Figure 2(c). Both test and model show the volume reduction of about 43% as shown in Figure 2(d).

Figure 2. Comparison between ACSED model and uniaxial extension test of bovine articular cartilage.

3.3. Modeling of Human Abdominal Aorta

Uniaxial extension tests for circumferential and axial strips of anatomically separated intima, media, and advantitia layers from human abdominal aorta, along with the kinematic characterization of circumferential contraction—axial stretch relation, have been studied by Holzapfel (2006) [35]. The test results of media axial strips with two families of collagen fiber angles of $\pm \theta ={52.2}^{\circ }$ at a constant crosshead speed of 1 mm/min and 37˚C have been selected. A self-developed graphics digitizer with MATLAB has been used to read out experimental data and converted to axial nominal stress $P_1$ and axial stretch ${\lambda }_1$ and the kinematic relation between circumferential contraction ${\lambda }_2$ and axial stretch ${\lambda }_1$. The circumferential contraction and radial reduction are assumed to relate axial stretch by constant Poisson’s ratios in (15). With curve-fitting, the Poisson’s ratios are ${\nu }_{12}=0.359$ and ${\nu }_{13}=0.491$ for the volume ratio at maximum stretch $J_m=1.04$. Uniaxial extension test data for media strips with $J_m=1.04$ have been fitted by the anisotropic CSE model (30) with the stress normalization condition (31). The constitutive parameters have been solved by the TED-LLSQ method. The comparison between test data and model for stress—stretch relation and kinematic relations with $J_m=1.04$ is shown in Figure 3.

In the physiological pressure range from 50 to 200 mmHg, a relative volume gain for most arteries is 2% to 6% as studied by Yosibash et al. (2014) [36] and Yossef et al. (2017) [37]. Thus, the average value of relative volume increase 4% is used.

4. Discussion

4.1. Modeling of Biological Tissues

The constitutive parameters resolved from fitting uniaxial extension tests of compressible human abdominal skins, bovine articular cartilages, and human abdominal aortas are listed in Table 1.

Figure 3. Comparison between ACSED model and uniaxial extension test of human aorta.

Table 1. Constitutive parameters of ACSED models for biological tissues.

Kinematic relations are indispensable in constitutive modeling compressible biological tissues. Poisson’s ratio plays a crucial role in relating stretches without work done to the stretch with a work done. The historical development on Poisson’s ratio has been reviewed by Greaves (2013) [38]. Constant, tangent, as well as other definitions of Poisson’s ratios have been summarized by Smith, Wootton, and Evans (1999) [39]. Their applications are extensive and some of which are studied by Beatty and Stalnaker (1986) [40], Carniel et al. (2014) [41], and Wahlsten et al. (2019) [33]. In this study, the tangent Poisson’s ratio as a cubic polynomial function of stretch is assumed in (34). The general kinematic model in (35) has been successfully applied to biological tissues.

By and large, the ACSED functional (14), along with the general kinematic model (35), offers a new approach in constitutive modeling of compressible as well as incompressible biological tissues.

4.2. The Isochoric-Volumetric Split

The $c_{\mathrm{3,}i}$ term with $I_{\mathrm{1,}i}^{3c_{\mathrm{4,}i}+1}/I_3^{c_{\mathrm{4,}i}}$ in the ACSED functional (14) represents ellipsoidal stretch deformation rather than spherical deformation.

The $c_{\mathrm{3,}i}$ term with $I_{\mathrm{1,}i}^{3c_{\mathrm{4,}i}}/I_3^{c_{\mathrm{4,}i}+0.5}$ in the Cauchy stress tensor (37) makes the multiplicative decomposition of the deformation gradient incomplete. Thus, the Cauchy stress tensor cannot be decoupled into a deviatoric stress tensor and a hydrostatic stress tensor.

The assumption of the isochoric—volumetric split is equivalent to assume that hydrostatic stress should only depend on volume change as emphasized by Charrier et al. (1988) [42] and by Murphy and Rogerson (2018) [28]. Taking one-third of the trace of Cauchy stress tensor (11) yields the general hydrostatic stress

$\frac{1}{3}\text{tr}{\sigma }_i=\frac{2}{3J}\left(I_{\mathrm{1,}i}\frac{\partial {\Psi }_i}{\partial I_{\mathrm{1,}i}}+2I_{\mathrm{2,}i}\frac{\partial {\Psi }_i}{\partial I_{\mathrm{2,}i}}+3I_3\frac{\partial {\Psi }_i}{\partial I_3}\right)\mathrm{.}$ (56)

Hydrostatic stress in (56) is a function of $I_{\mathrm{1,}i}$, $I_{\mathrm{2,}i}$, and $I_3$, indicating the assumption of hydrostatic stress only as a function of $I_3$ is false. Thus, the isochoric—volumetric split is invalid.

4.3. Predictive ACSED Model

The ACSED functional was formulated by three different grouping processes. First, irreducible invariants are grouped into three PICs (2). Second, the PICs are lumped into dimensionless PIC groups (13) based on energy balance through Lie group method by Zhao (2016) [43]. Third, the dimensionless PIC groups are conducted based on curvatures in elasticity tensors through differential geometry.

The ACSED model can accurately fit and predict test data. Fitting capability of ACSED model has been demonstrated here but predictive capability is yet to be conducted. However, both fitting and predicting capabilities for the isotropic continuum stored energy model formulated through the same three grouping processes have been thoroughly examined through Treloar’s test data in uniaxial extension, uniaxial compression, equibiaxial extension, biaxial extension, pure shear, and simple shear modes by Zhao (2021) [44].

Biological tissues exhibit complex mechanical behaviors. Complexities include highly nonlinear stress—stretch response, anisotropy, inhomogeneity, rate dependency, active and passive responses, pre-stress, residual stress, in vivo and ex vivo characterizations, and others. Both accurate and complete set of experimental tests in multiple deformation modes for anisotropic biological tissues are more difficult to conduct but are still desirable for examining the predictive capability of anisotropic constitutive models.

5. Conclusion

The ACSED functional was formulated through three grouping processes rather than the three major splits. The one-dimensional normal stretch, two-dimensional shear stretch, and three-dimensional ellipsoidal stretch deformations become the new metric of nonlinear elastic deformation. The stress and elasticity tensors in both reference and current configurations have been derived. With the cubic polynomial for a tangent Poisson’s ratio, the general kinematic model is established. The ACSED model, along with the kinematic model, has been accurately fitted to the uniaxial extension test data of compressible human abdominal skin, bovine articular cartilage, and human abdominal aorta samples. In the future, the predictive capability of the ACSED model should be examined for biological tissues in different deformation modes.

Acknowledgements

The author is thankful to his family, Jianming and Jiesi Zhao, for their encouragement, support, and help.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Holzapfel, G.A. (2000) Nonlinear Solid Mechanics: A Continuum Approach for Engineering. John Wiley & Sons, Chichester.
[2]	Humphrey, J.D. (2003) Continuum Biomechanics of Soft Biological Tissues. Proceedings of the Royal Society of London A, 459, 3-46. https://doi.org/10.1098/rspa.2002.1060
[3]	Wriggers, P. (2009) Mixed Finite Element Methods-Theory and Discretization. In: Carstensen, C. and Wriggers, P., Eds., Mixed Finite Element Technologies, CISM International Centre for Mechanical Sciences, Vol. 509, Springer, Wienna, 131-177. https://doi.org/10.1007/978-3-211-99094-0_5
[4]	Chagnon, G., Rebouah, M. and Favier, D. (2015) Hyperelastic Energy Densities for Soft Biological Tissues: A Review. Journal of Elasticity, 120, 129-160. https://doi.org/10.1007/s10659-014-9508-z
[5]	Lanir, Y. (2017) Multi-Scale Structural Modeling of Soft Tissues Mechanics and Mechanobiology. Journal of Elasticity, 129, 7-48. https://doi.org/10.1007/s10659-016-9607-0
[6]	Cowin, S.C. (1990) Deviatoric and Hydrostatic Mode Interaction in Hard and Soft Tissue. Journal of Biomechanics, 23, 11-14. https://doi.org/10.1016/0021-9290(90)90364-9
[7]	Sansour, C. (2008) On the Physical Assumptions Underlying the Volumetric-Isochoric Split and the Case of Anisotropy. European Journal of Mechanics A/Solids, 27, 28-39. https://doi.org/10.1016/j.euromechsol.2007.04.001
[8]	Ní Annaidh, A., Destrade, M., Gilchrist, M.D. and Murphy, J.G. (2013) Deficiencies in Numerical Models of Anisotropic Nonlinearly Elastic Materials. Biomechanics and Modeling in Mechanobiology, 12, 781-791. https://doi.org/10.1007/s10237-012-0442-3
[9]	Wang, M. and Liu, F. (2018) Compressible Hyperelastic Models for Soft Biological Tissue: A Review. Journal of Biomaterials and Tissue Engineering, 8, 1375-1389. https://doi.org/10.1166/jbt.2018.1888
[10]	Jin, T. and Stanciulescu, I. (2019) A Note on the Volumetric-Deviatoric Split on the Anisotropoc Constitutive Model for Fiber-Reinforced Materials. Biomedical Engineering International, 1, 16-24. https://doi.org/10.33263/BioMed12.016024
[11]	Zhao, F. (2020) Modeling and Implementing Compressible Isotropic Finite Deformation without the Isochoric-Volumetric Split. Journal of Advances in Applied Mathematics, 5, 57-70. https://doi.org/10.22606/jaam.2020.52002
[12]	Zhao, F. (2019) On Constitutive Modeling of Arteries. Journal of Advances in Applied Mathematics, 4, 54-68. https://doi.org/10.22606/jaam.2019.42003
[13]	Spencer, A.J.M. (1971) Theory of Invariants. In: Eringen, A.C., Ed., Continuum Physics, Vol. 1, Academic Press, New York, 239-353. https://doi.org/10.1016/B978-0-12-240801-4.50008-X
[14]	Zheng, Q.-S. (1994) Theory of Representations for Tensor Functions—A Unified Invariant Approach to Constitutive Equations. Applied Mechanics Reviews, 47, 545-587. https://doi.org/10.1115/1.3111066
[15]	Bathe, K.-J. (1996) Finite Element Procedures. Prentice-Hall, Inc., Englewood Cliffs.
[16]	Pian, T.H.H. (1964) Derivation of Element Stiffness Matrices by Assumed Stress Distributions. AIAA Journal, 2, 1333-1336. https://doi.org/10.2514/3.2546
[17]	Oden, J.T. and Key, J.E. (1970) Numerical Analysis of Finite Axisymmetrical Deformations of Incompressible Elastic Solids of Revolution. International Journal of Solids and Structures, 6, 497-518. https://doi.org/10.1016/0020-7683(70)90027-2
[18]	Häggblad, B. and Sundberg, J.A. (1983) Large Strain Solutions of Rubber Components. Computers and Structures, 17, 835-843. https://doi.org/10.1016/0045-7949(83)90097-4
[19]	Zdunek, A. and Rachowicz, W. (2017) A Mixed Higher Order FEM for Fully Coupled Compressible Transversely Isotropic Finite Hyperelasticity. Computers and Mathematics with Applications, 74, 1727-1750. https://doi.org/10.1016/j.camwa.2017.02.042
[20]	Legrain, G., Moësa, N. and Huerta, A. (2008) Stability of Incompressible Formulations Enriched with X-FEM. Computer Methods in Applied Mechanics and Engineering, 197, 1835-1849. https://doi.org/10.1016/j.cma.2007.08.032
[21]	Vergori, L., Destrade, M., McGarry, P. and Ogden, R.W. (2013) On Anisotropic Elasticity and Questions Concerning Its Finite Element Implementation. Computational Mechanics, 52, 1185-1197. https://doi.org/10.1007/s00466-013-0871-6
[22]	Helfenstein, J., Jabareen, M., Mazza, E. and Govindjee, S. (2010) On Non-Physical Response in Models for Fiber-Reinforced Hyperelastic Materials. International Journal of Solids and Structures, 47, 2056-2061. https://doi.org/10.1016/j.ijsolstr.2010.04.005
[23]	Horgan, C.O. and Murphy, J.G. (2020) Some Unexpected Predictions from Strongly Anisotropic Hyperelastic Constitutive Models of Soft Tissue. Mechanics of Soft Materials, 2, 9. https://doi.org/10.1007/s42558-020-00024-5
[24]	Limbert, G. and Taylor, M. (2002) On the Constitutive Modeling of Biological Soft Connective Tissues: A General Theoretical Framework and Explicit Forms of the Tensors of Elasticity for Strongly Anisotropic Continuum Fiber-Reinforced Composites at Finite Strain. International Journal of Solids and Structures, 39, 2343-2358. https://doi.org/10.1016/S0020-7683(02)00084-7
[25]	Nolan, D.R., Gower, A.L., Destrade, M., Ogden, R.W. and McGarry, J.P. (2014) A Robust Anisotropic Hyperelastic Formulation for the Modelling of Soft Tissue. Journal of the Mechanical Behavior of Biomedical Materials, 39, 48-60. https://doi.org/10.1016/j.jmbbm.2014.06.016
[26]	Pierrat, B., Murphy, J.G., MacManus, D.B. and Gilchrist, M.D. (2016) Finite Element Implementation of a New Model of Slight Compressibility for Transversely Isotropic Materials. Computer Methods in Biomechanics and Biomedical Engineering, 19, 745-758. https://doi.org/10.1080/10255842.2015.1061513
[27]	Cai, R., Holweck, F., Feng, Z.-Q. and Peyraut, F. (2016) A New Hyperelastic Model for Anisotropic Hyperelastic Materials with One Fiber Family. International Journal of Solids and Structures, 84, 1-16. https://doi.org/10.1016/j.ijsolstr.2015.11.008
[28]	Murphy, J.G. and Rogerson, G.A. (2018) Modelling Slight Compressibility for Hyperelastic Anisotropic Materials. Journal of Elasticity, 131, 171-181. https://doi.org/10.1007/s10659-017-9650-5
[29]	Fereidoonnezhad, B., O’Connor, C. and McGarry, J.P. (2020) A New Anisotropic Soft Tissue Model for Elimination of Unphysical Auxetic Behavior. Journal of Biomechanics, 111, Article ID: 110006. https://doi.org/10.1016/j.jbiomech.2020.110006
[30]	Zhao, F. (2018) Anisotropic Continuum Stored Energy Functional Solved by Lie Group and Differential Geometry. Advances in Pure Mathematics, 8, 631-651. https://doi.org/10.4236/apm.2018.87037
[31]	Blatz, P.J. and Ko, W.L. (1962) Application of Finite Elastic Theory to the Deformation of Rubbery Materials. Transactions of the Society of Rheology, 6, 223-251. https://doi.org/10.1122/1.548937
[32]	Fortes, M.A. and Nogueira, M.T. (1989) The Poisson Effect in Cork. Materials Science and Engineering: A, 122, 227-232. https://doi.org/10.1016/0921-5093(89)90634-5
[33]	Wahlsten, A., Pensalfini, M., Stracuzzi, A., Restivo, G., Hopf, R. and Mazza, E. (2019) On the Compressibility and Poroelasticity of Human and Murine Skin. Biomechanics and Modeling in Mechanobiology, 18, 1079-1093. https://doi.org/10.1007/s10237-019-01129-1
[34]	Woo, S.L.-Y., Lubock, P., Gomez, M.A., Jemmott, G.F., Kuei, S.C. and Akeson, W.H. (1979) Large Deformation Nonhomogeneous and Directional Properties of Articular Cartilage in Uniaxial Tension. Journal of Biomechanics, 12, 437-466. https://doi.org/10.1016/0021-9290(79)90028-9
[35]	Holzapfel, G.A. (2006) Determination of Material Models for Arterial Walls from Uniaxial Extension Tests and Histological Structure. Journal of Theoretical Biology, 238, 290-302. https://doi.org/10.1016/j.jtbi.2005.05.006
[36]	Yosibash, Z., Manor, I., Gilad, I. and Willentz, U. (2014) Experimental Evidence of the Compressibility of Arteries. Journal of the Mechanical Behavior of Biomedical Materials, 39, 339-354. https://doi.org/10.1016/j.jmbbm.2014.07.030
[37]	Yossef, O.E., Farajian, M., Gilad, I., Willenz, U., Gutman, N. and Yosibash, Z. (2017) Further Experimental Evidence of the Compressibility of Arteries. Journal of the Mechanical Behavior of Biomedical Materials, 65, 177-189. https://doi.org/10.1016/j.jmbbm.2016.08.013
[38]	Greaves, G.N. (2013) Poisson’s Ratio over Two Centuries: Challenging Hypotheses. Notes and Records of the Royal Society of London, 67, 37-58. https://doi.org/10.1098/rsnr.2012.0021
[39]	Smith, C.W., Wootton, R.J. and Evans, K.E. (1999) Interpretation of Experimental Data for Poisson’s Ratio of Highly Nonlinear Materials. Experimental Mechanics, 39, 356-362. https://doi.org/10.1007/BF02329817
[40]	Beatty, M.F. and Stalnaker, D.O. (1986) The Poisson Function of Finite Elasticity. Journal of Applied Mechanics, 53, 807-813. https://doi.org/10.1115/1.3171862
[41]	Carniel, E.L., Gramigna, V., Fontanella, C.G., Frigo, A., Stefanini, C., Rubini, A. and Natali, A.N. (2014) Characterization of the Anisotropic Mechanical Behaviour of Colonic Tissues: Experimental Activity and Constitutive Formulation. Experimental Physiology, 99, 759-771. https://doi.org/10.1113/expphysiol.2013.076091
[42]	Charrier, P., Dacorogna, B., Hanouzet, B. and Laborde, P. (1988) An Existence Theorem for Slightly Compressible Materials in Nonlinear Elasticity. SIAM Journal on Mathematical Analysis, 19, 70-85. https://doi.org/10.1137/0519005
[43]	Zhao, F. (2016) Continuum Constitutive Modeling for Isotropic Hyperelastic Materials. Advances in Pure Mathematics, 6, 571-582. https://doi.org/10.4236/apm.2016.69046
[44]	Zhao, F. (2021) Predictive Continuum Constitutive Modeling of Unfilled and Filled Rubbers. Journal of Advances in Applied Mathematics, 6, 144-159. https://doi.org/10.22606/jaam.2021.63002