0
Research Papers

# Numerical Approximation of Elasticity Tensor Associated With Green-Naghdi RateOPEN ACCESS

[+] Author and Article Information
Haofei Liu

Department of Mechanics,
Tianjin University,
Tianjin 300072, China

Wei Sun

Tissue Mechanics Laboratory,
The Wallace H. Coulter Department of Biomedical Engineering,
Georgia Institute of Technology,
Technology Enterprise Park,
Room 206, 387 Technology Circle,
Atlanta, GA 30313-2412
e-mail: wei.sun@bme.gatech.edu

Manuscript received July 26, 2016; final manuscript received May 3, 2017; published online June 16, 2017. Assoc. Editor: Hai-Chao Han.

J Biomech Eng 139(8), 081007 (Jun 16, 2017) (8 pages) Paper No: BIO-16-1317; doi: 10.1115/1.4036829 History: Received July 26, 2016; Revised May 03, 2017

## Abstract

Objective stress rates are often used in commercial finite element (FE) programs. However, deriving a consistent tangent modulus tensor (also known as elasticity tensor or material Jacobian) associated with the objective stress rates is challenging when complex material models are utilized. In this paper, an approximation method for the tangent modulus tensor associated with the Green-Naghdi rate of the Kirchhoff stress is employed to simplify the evaluation process. The effectiveness of the approach is demonstrated through the implementation of two user-defined fiber-reinforced hyperelastic material models. Comparisons between the approximation method and the closed-form analytical method demonstrate that the former can simplify the material Jacobian evaluation with satisfactory accuracy while retaining its computational efficiency. Moreover, since the approximation method is independent of material models, it can facilitate the implementation of complex material models in FE analysis using shell/membrane elements in abaqus.

<>

## Introduction

With an ever increasing use of finite element (FE) analysis to study biomechanics involved in biomedical applications, many constitutive models have been developed to simulate a variety of biomaterial behaviors. These models not only have nonlinear, hyperelastic properties, undergoing finite deformation, such as the Holzapfel–Gasser–Ogden (HGO) model [1], but can also incorporate microstructural [2], viscoelastic [3], growth and remodeling [4,5], active contraction [6], or fatigue damage characteristics [7]. Consequently, these models vary in functional forms and sometimes involve time series and complex integrations. The implementation of such experimentally derived, user-defined material models in commercial FE packages such as abaqus (Simulia, Providence, RI) can be challenging, which may pose a barrier for the widespread use of these material models.

To implement a user-defined material model, the user must explicitly specify, in the user-defined material subroutine (UMAT), the Cauchy stress and the tangent modulus tensor (also known as elasticity tensor or material Jacobian) associated with the material model. Correct derivation and coding of these two quantities in the user material subroutine are important to obtain converged simulation results. While the Cauchy stress derived from a constitutive model is relatively easy to obtain, the tangent modulus tensor is often problematic especially when multidecomposition of the constitutive model into volumetric and distortional parts is required. Tangent modulus tensor is the power-conjugate between a specific form of the objective stress rate and the deformation rate tensor. Objective stress rates, referring to the intrinsic requirement for the constitutive law to be frame-independent, can be expressed in different forms, including the Jaumann rate and Green-Naghdi rate of Kirchhoff or Cauchy stresses in commercial FE packages, such as abaqus.

Shell/membrane elements are often preferred in finite element analyses of thin membrane structures such as heart valves [712]. Shell and membrane elements possess significant advantages over 3D solid elements in simulating thin structures due to the simplification in topology and reduced demand for computational resources. Unfortunately, the tangent modulus tensor based on the Green-Naghdi rate of Kirchhoff stress, which is used for abaqus shell/membrane elements [13], is rather complex to evaluate analytically. Details on the analytical calculation of the elasticity tensor associated with Green-Naghdi rate of Kirchhoff stress can be found in Simo and Hughes [14], Prot et al. [15], and also Sec. 2.4 of this paper.

Since the tangent moduli serve as an iterative operator in solving boundary value problems using Newton's method, the exact analytical solution is preferred to provide quadratic rate of convergence [16], but not mandatory. This critical property of the tangent moduli in computational mechanics gives rise to an alternative method—approximation approaches of the tangent modulus tensor, which can be easy to implement without substantially compromising the convergence rate. Miehe [17] reported an approximation to the elasticity tensor which was based on the convective rate of Kirchhoff stress. Sun et al. [18] presented a method that approximated the elasticity tensor associated with the Jaumann rate of Kirchhoff stress tensor, which is used in 3D solid elements in abaqus. Liu and Sun [19] conducted an accuracy analysis and computational efficiency study of the method proposed in Ref. [18] and confirmed that the approximation was able to achieve reasonably good accuracy and comparable computational efficiency to the closed-form solution. Recently, Tanaka and Fujikawa [20] suggested an approximation method for tangent moduli based on Green-Naghdi rate of Kirchhoff stress and implemented the method using hyperelastic material models. However, to rigorously evaluate the accuracy and efficiency of the approximation method, a side-by-side comparison of the approximation and the closed-form solution is needed. In this paper, we present a detailed description of both methods for the evaluation of elasticity tensor associated with the Green-Naghdi rate of Kirchhoff stress. By comparing the approximation results with the closed-form solution using a popular hyperelastic material model in the biomechanics community, the accuracy and computational efficiency of the approximation method are evaluated. Finally, a more complex material model, involving the integration of a fiber distribution function, is incorporated to demonstrate the effectiveness of the method.

## Continuum Mechanics Framework

###### Kinematics.

Deformation of a body is described by mapping $χ$ from a material point $X$ at its undeformed (referential) configuration $Ω0$ to its position $x=χ(X,t)$ at deformed (current) configuration $Ω$ at time t. The deformation can be described by the deformation gradient F, defined as $F=(∂x/∂X)$, with Jacobian $J=det(F)$. $C=FTF$ and are the right and left Cauchy–Green tensors, respectively. The Green–Lagrange strain tensor in material description is expressed as $E=(1/2)(C−I)$.

The deformation gradient can also be expressed as $F=RU=VR$ via Polar decomposition, where R is the rotation tensor, representing local rotation of a referential orientation. U and V denote the right and left stretch tensors, respectively. Due to the requirement of objectivity (i.e., frame indifference), the mechanical responses of materials are frequently addressed in corotational systems. A corotational system is defined so that the coordinate system is embedded in the material and rotates with it. Hence, the deformation gradient tensor in the corotational system can be expressed as: $F̂=RTF$, where the hats “ $̂$ ” are used to denote parameters in the corotational frame throughout the paper. It should be noted that the deformation gradient tensor F is a two-point tensor; therefore, only the first part of it describing the spatial coordinates is rotated back to counter the local rotation.

The rate of deformation can be defined as velocity (v) gradient: $L=grad(v)=F˙F−1$. It can be further decomposed into a symmetric and asymmetric part Display Formula

(1)$L=D+W, D=sym(L)=12(L+LT), W=skew(L)=12(L−LT)$

where D and W are known as the rate of deformation tensor and the rate of rotation tensor, respectively. A spin tensor $Ω:=R˙RT$ stands for the spin of the current orientation.

###### Hyperelastic Stress Response.

Among various measures of stress, three of them are employed in this paper. They are (1) the second Piola-Kirchhoff stress, S, defined in the reference configuration, (2) the Kirchhoff stress, $τ$, which is a push-forward of the second Piola-Kirchhoff stress, and (3) the Cauchy stress $σ$. These three measures of stress are expressed as [21] Display Formula

(2)$S=2∂W(C)∂C, τ=FSFT, σ=J−1τ$

where W denotes a scalar valued strain-energy function postulated to exist for hyperelastic materials.

###### Tangent Modulus Tensor.

Implementation of a user-defined material model in abaqus/umat usually involves the evaluation of (a) Cauchy stress $σ$ and (b) the tangent modulus tensor $C$ that is derived from a rate-form constitutive equation. There are various forms of elasticity tensors with respect to different forms of objective rate of stress. The following objective rates, namely, the Jaumann rate of the Kirchhoff stress $τ∇J$, the convective rate of the Kirchhoff stress $τ∇c$, and the Green-Naghdi rate of the Kirchhoff stress $τ∇G$ are frequently used in the literature and adopted in abaqusDisplay Formula

(3)$τ∇J=τ˙−Wτ−τWT=cτJ:D$
Display Formula
(4)$Lv(τ)=τ∇c=τ˙−Lτ−τLT=cτc:D$
Display Formula
(5)$τ∇G=τ˙−Ωτ−τΩT=cτG:D$

where the notation $∇$ in the superscript denotes objective derivative, while the notation “$⋅$” stands for material time derivative, $cτJ$ (or $cτc$, $cτG$) denotes the fourth-order elasticity tensor, with the superscripts representing the form of stress and stress rate utilized in the equation, i.e., $τJ$ indicates the Jaumann rate of the Kirchhoff stress, $τc$ is the convective rate of the Kirchhoff stress, $τG$ is the Green-Naghdi rate of the Kirchhoff stress, and $Lv(*)$ is the Lie time derivative of *. For 3D continuum elements, abaqus uses the elasticity tensor associated with the Jaumann rate of the Kirchhoff stress. For structural elements, abaqus uses the one associated with the Green-Naghdi rate of the Kirchhoff stress $cτG$. It is also worth noting that we represent the elasticity tensors in the current configuration with the lowercase $c$, whereas the elasticity tensor in the referential configuration is denoted using the uppercase $C$.

In Secs. 2.4 and 2.5, the elasticity tensor $cτG$ is calculated with both analytical and approximation approaches.

###### Analytical Method.

The elasticity tensor in the referential configuration is expressed as [14] Display Formula

(6)$C=∂S(C)∂C=4∂2W(C)∂C∂C$

The elasticity tensor $cτc$ defined in the current configuration is then obtainable after a push-forward operation [15,21] Display Formula

(7)$cτc=χ*(C), (cτc)abcd=FaAFbBFcCFdD(C)ABCD$

By comparing Eq. (3) with Eq. (4) and using the split $L=D+W$ in Eq. (1), the tangent moduli for the Jaumann rate of the Kirchhoff stress can be shown as Display Formula

(8)

Finally, the elasticity tensor for the Green-Naghdi rate of the Kirchhoff stress $CτG$ is expressible as Display Formula

(9)$cτG=cτJ+cspin$

where $Cspin$ is obtained by Display Formula

(10)$cspin:D=(W−Ω)τ−τ(W−Ω)$

Upon adopting the definition and notation in Simo and Hughes [14], we have Display Formula

(11)$W−Ω=Λ:D$

where Display Formula

(12)$Λabcd=c1(Vacδbd+Vadδbc−Vbdδac−Vbcδad)−c2(Bacδbd+Badδbc−Bbdδac−Bbcδad)+c3(BacVbd+BadVbc−VacBbd−VadBbc)$

with Display Formula

(13)$c1=IV22(IVIIV−IIIV), c2=IV2(IVIIV−IIIV), c3=12(IVIIV−IIIV)$

and $IV:=tr(V)$, $IIV:=(1/2)(tr(V)2−tr(V⋅V))$, and $IIIV:=det(V)$ are the invariants of the left stretch tensor V, respectively.

$cτG$ can thus be represented in its index form as Display Formula

(14)$(cτG)abcd=(cτJ)abcd+Λakcdτkb−τakΛkbcd$

It is noted from Eqs. (6)(8), (12), and (14) that while $cτJ$ possesses full symmetry (not only major symmetry but also minor symmetries), $cτG$, on the other hand, loses its major symmetry while retaining minor symmetries. It is also noted that the UMAT implementation of the Green-Naghdi rate of the Kirchhoff stress $τ∇G$ in abaqus requires all the equations from Eqs. (6)(14). The implementation procedure is thus complex, lengthy, and prone to coding errors, especially when the constitutive equations have complex functional forms.

###### Approximation Method.

Alternatively, Sun et al. [18] and Tanaka and Fujikawa [20] showed that numerical approximation of a rate form of elasticity tensors can make such evaluation much simpler to implement and the approach is independent of material models.

It is known that rotations taking place in the deformation are compensated in the corotational configuration, and the material stress rate in the configuration is objective. Therefore, the Green-Naghdi rate of the Kirchhoff stress, in the corotational configuration, is defined as Display Formula

(15)$τ̂∇G=τ̂˙−Ω̂τ̂−τ̂Ω̂T=ĉτG:D̂$

where $τ̂=RTτR$.

Due to the fact that the rotation R goes to identity in the corotational configuration, the spin tensor $Ω=R˙RT$ goes to zero, thus, Eq. (15) reduces to Display Formula

(16)$τ̂∇G=τ̂˙=ĉτG:D̂$

Equation (16) can then be linearized as Display Formula

(17)$Δτ̂=ĉτG:ΔD̂$

Following Miehe [17], Sun et al. [18], Liu and Sun [19], the perturbation to the rate of deformation can be further written as Display Formula

(18)$ΔD̂(îĵ)=12(ΔL̂+ΔL̂T)=12(ΔF̂(îĵ)F̂−1+(ΔF̂(îĵ)F̂−1)T)$

with the perturbation to the deformation gradient $ΔF̂$ defined in the following form: Display Formula

(19)$ΔF̂(îĵ)=ε2[(êî⊗êĵ)F̂+(êĵ⊗êî)F̂]$

where $ε$ is the perturbation parameter, and ${êî}î=1,2,3$ represents the basis vectors in the corotational frame.

Therefore, Display Formula

(20)$ΔD̂(îĵ)=ε2(êî⊗êĵ+êĵ⊗êî)$

The increment of the Kirchhoff stress $τ̂$ is taken as Display Formula

(21)$Δτ̂≅τ̂(F̂+ΔF̂(ij))−τ̂(F̂)$

Combining Eqs. (17), (20), and (21) and utilizing the (minor) symmetric feature of $cτG$, the forward Euler approximation of $cτG$ in the corotational frame has its index form Display Formula

(22)$(ĉτG)âb̂îĵ≅1ε[τ̂(F̂+ΔF̂(îĵ))âb̂−τ̂(F̂)âb̂]$

The tangent moduli for the Green-Naghdi rate of the Kirchhoff stress in the current configuration can be recovered via Display Formula

(23)$(ĉτG)abcd=RaâRbb̂RcĉRdd̂(ĉτG)âb̂ĉd̂$

By converting the Green-Naghdi rates (obtained from either analytical or approximate evaluation) to the form that is accepted by abaqus, we get the material Jacobian Display Formula

(24)$cMJ=1JcτG$

## Numerical Examples

In this section, we demonstrate the implementation of the approximation method with two constitutive models: (1) a modified version of the Holzapfel–Gasser–Ogden (HGO) model [1]. The HGO model (and its derivatives) is a fiber structural-based hyperelastic model that is commonly used to describe soft tissue material properties in biomedical applications [2224]. (2) A fiber structural-based model that employs a continuous fiber distribution function [25]. This model [25] has an integral function for which the analytical form of the tangent moduli is difficult to obtain. We use these two examples to demonstrate the effectiveness of the approximation method in abaqus.

###### Numeric Example 1—Modified HGO Model.

The HGO model [1] is modified as Display Formula

(25)$W(C,M,N,J)=c10(exp[c01(I1−3)]−1)+k12k2[exp(k2ε12)−1]+k12k2[exp(k2ε22)−1]−p(J−1)$

where

Display Formula

(26)$ε1=κ(I1−3)+(1−3κ)(I4−1) and ε2=κ(I1−3)+(1−3κ)(I6−1)$

$c01 and c10$ are the material parameters characterizing extracellular matrix and are the material parameters for embedded collagen fibers. $κ$, ranging from 0 to 1/3, determines the dispersion of fibers; the invariants $I1:=tr(C)$, $I4:=C:(M⊗M)$, and $I6:=C:(N⊗N)$, where M$(M=[ cos α sin α0])$ and N $(N=[ cos β sin β0])$ represent the referential fiber orientations, with $α$ and $β$ representing the angles between the orientations of the fibers in the undeformed state and the material axis; $J=1$ for incompressible material; and p is an indeterminate Lagrangian multiplier. The analytical calculations of stress tensor S, elasticity tensor $C$, and the determination of p for the modified HGO model are shown in the Appendix.

###### Finite Element Model Using the Modified HGO Model.

In this study, material parameters of the modified HGO model were obtained by means of data fitting to match the previously published planar biaxial mechanical testing [26] (resultant parameters: $c01$ = 14.87, $c10$ = 1.16 kPa, $k1$ = 4.48 kPa, $k2$ = 62.20, $κ$ = 0, $α$ = 32.64 deg, and $β$ = 56.02 deg). The modified HGO model was implemented for abaqus membrane elements via UMAT. The effectiveness of the implementation was examined by simulating the biaxial experiment under experimental loading condition. Detailed descriptions of the biaxial testing on soft tissue can be found in Sacks [27] and Sacks et al. [28].

The FE model was then set up in accordance with the specimen (glutaraldehyde-treated bovine pericardium (GLBP)) with the dimension of 25 mm × 25 mm × 0.4 mm. Four evenly spaced nodal forces, with 5 mm apart from each other and 2.5 mm inside the specimen edge (Fig. 1(a)), were imposed on each side. Each nodal force was 2.5 N, which produced 1 MPa Lagrangian stress on each edge. The geometry domain (Fig. 1(a)) was discretized with 20 × 20 eight-node membrane elements (M3D8). Especially, the 4 × 4 elements (5 mm × 5 mm) in the central region (enclosed by the frame in Fig. 1(a)) were selected to represent the square region enclosed by the markers in the biaxial testing. Within this region, the averaged strain and stress from the numerical simulation were compared with those taken from the biaxial test. In the material model implementation, elasticity tensors were calculated with both approximate and analytical methods described in Secs. 2.4 and 2.5. The simulations were performed on a Windows 64-bit desktop computer equipped with Intel Core i7-3770 processor (3.40 GHz) and 16 GB RAM. The predicted stress contour at end loading is shown in Fig. 1(b).

Figure 1(c) shows that the stress–strain relations obtained from experiment and numerical simulation are in good agreement, highlighting the effectiveness of the material model in capturing the key features of the soft tissue property.

###### Accuracy Analysis of Approximation.

The accuracy of the approximation can be quantified using the relative error of the approximated tangent. The relative error was defined as follows [29]: Display Formula

(27)$Er=(∑a,b,c,d(cabcdMJ−cabcdMJ,apm)2)1/2/(∑a,b,c,d(cabcdMJ)2)1/2$

where the superscript “apm” represents the approximated values.

The relative error obtained using the approximation method was plotted against a range of perturbation parameter values of $ε$ in Fig. 2(a). As shown in Tanaka and Fujikawa [20], the approximation errors are of the same order as $ε$, which suggests a preference for a small value of $ε$. However, when $ε$ becomes too small, the evaluation is prone to numerical truncation errors. There exists an optimal value for $ε$. It is shown in Fig. 2(a) that the approximation achieved its optimal accuracy when $ε$ was selected to be $εo=$ 1.0 × 10−8.

###### Computational Efficiency Analysis.

The total iteration number and the CPU time used to achieve convergence were plotted against varying values of perturbation parameter, as shown in Figs. 2(b) and 2(c). For comparison purposes, UHYPER, the abaqus subroutine officially recommended for user-defined hyperelastic material, is utilized as a benchmark for performance. In general, the computational efficiency follows the approximation accuracy of the tangent, that is, the more accurate the approximation is, the less iteration number (the shorter CPU time) is taken to converge the simulation. To take a closer look at the computational efficiency, the residual norms of the largest force at each iteration are listed in Table 1 for all the tested schemes. A typical load step was selected to do the comparison. It is shown that the approximation achieves quadratic rate of convergence when $εo$ is selected for perturbation. It is observed that the performance of approximation is similar to that of UHYPER in terms of CPU time and convergence rate when optimal value is chosen for perturbation.

It is worth noting that the CPU time per iteration of an implicit FE simulation would be dominated by the factoring time of stiffness matrix when the simulation scale becomes sufficiently large [29]. Therefore, the closed-form solution should in general have better computational efficiency for large-scale simulation, due to the intrinsic slower convergence of the approximation method. However, given the severe stress concentration taking place in this simulation and that the approximation takes almost the same number of iterations as the UHYPER, the overall performance of the approximation is similar to that of the closed-form solution. Combining its feature of simple to use and material model independency, the approximation method remains attractive.

###### Numeric Example 2—The Material Model With Continuous Fiber Distribution Function.

Sun et al. [25] proposed a structural model for arterial tissue. It takes the form as Display Formula

(28)$W(C,M,N,J)=(1−df)d0(I1−3)+∫−π2π2df(R1(θ)Wf1(I4)+R2(θ)Wf2(I6))dθ−p(J−1)$

where Display Formula

(29)$Wf1(I4)=d1(ed2(I4−1)2−1), I4>1Wf2(I6)=d1(ed2(I6−1)2−1), I6>1$

$df$ is the fiber volume fraction; $d0$, $d1$, and $d2$ are the material constants; and $R1(θ)$ and $R2(θ)$ are the fiber distribution functions which are normal distributions with standard deviation (SD) and mean angle between fiber and circumferential direction, $Θ¯$. Since this model has an integral function, it is relatively complicated to obtain the analytical form of the tangent moduli. Thus, the approximation method is used.

###### Finite Element Model Using the Distributed Fiber Model.

Following Sun et al. [25], we construct an FE model of rat carotid artery segment based on the geometry described by Zulliger et al. [30]. The vessel is fixed at both ends in the axial direction and is inflated up to 25 kPa with pressure acting on the inner surface. The computational model is comprised of 1056 M3D4 elements. The material parameters are determined so as to match the experimental data (Figs. 3(c) and 3(d)). The baseline group of parameters are Figures 3(a) and 3(b), show the initial meshes and deformed meshes after inflation.

Figures 3(c) and 3(d) show the change of the outermost radius as the vessel is pressurized in experiment [30] and numerical simulations, with varying $Θ¯$ and SD. It demonstrates that when fibers were oriented more in the circumferential direction, with $Θ¯$ from 45 deg to 15 deg, the artery is less compliant, whereas as the fiber dispersion increases, SD from 1 to 45, the vessel becomes stiffer, which are in agreement with Sun et al. [25].

###### Computational Efficiency Analysis.

Figure 3(e) shows the iteration number taken by the simulations for different perturbation parameters. When the perturbation parameter is in the range of 10−12 and 10−4, the iteration number required is unchanged. The effect of perturbation parameter becomes obvious only when the perturbation is chosen outside this range, and the iteration number increases substantially until no convergence is obtained. CPU time follows the same trend.

From these two numerical examples, we have demonstrated that the approximation method simplifies the evaluation of $cτG$. One only needs to calculate the stresses and their transformation between frames. In addition, a major advantage of the approximation method lies in its feature of material model independence, i.e., a UMAT code written for a particular material model can be easily adapted for another material model with minimum effort on coding—by simply updating the stress calculations of the new model. It enables the implementation of complex constitutive models in abaqus.

## Conclusion

This paper presents a detailed description of an approximation method to evaluate the tangent moduli associated with the Green-Naghdi rate of the Kirchhoff stress. By means of implementing a user-defined fiber-reinforced hyperelastic material model in abaqus/Standard via user-defined subroutine UMAT, the aforementioned method was successfully applied and its performance was examined. Upon a proper selection of the perturbation parameter, the approximation method achieved good accuracy and its computational efficiency was similar to UHYPER in terms of the total number of iterations. We believe the method presented here will facilitate the implementation of user-defined material models due to the fact that calculating and coding stress and rotation of frame are much easier than directly evaluating the tangent moduli based on Green-Naghdi rate of Kirchhoff stress itself, especially as material models are becoming more complicated in biomedical engineering applications.

## Acknowledgements

This project was funded in part by Grant Nos. NIH HL104080 and HL108240.

## Appendices

###### Appendix

The analytical derivation of the second P-K stress along with that of the Lagrangian multiplier p and the elasticity tensor can be presented as below:

The second Piola-Kirchhoff stress can be expressed as Display Formula

(A1)$S=2∂W(I1,I4,I6,J)∂C=2(W1I+W4M⊗M+W6N⊗N)−pC−1$

where Display Formula

(A2)$W1=∂W(I1,I4,I6,J)∂I1=c10c01(exp[c01(I1−3)]−1)+k1κ[ε1 exp(k2ε12)+ε2 exp(k2ε22)]W4=∂W(I1,I4,I6,J)∂I4=k1(1−3κ)ε1 exp(k2ε12)W6=∂W(I1,I4,I6,J)∂I6=k1(1−3κ)ε2 exp(k2ε22)$

According to the plane stress condition, the out of plane stress S33 is zero, so that we get $p=2W1C33$ and Display Formula

(A3)$S=2(W1I+W4M⊗M+W6N⊗N)−2W1C33C−1$

To calculate the elasticity tensor associated with Green-Naghdi rate, we start with the calculation of the referential elasticity tensor $C$Display Formula

(A4)$C=2∂S(C)∂C=4(W11I⊗∂I1∂C+W14I⊗∂I4∂C+W16I⊗∂I6∂C+W41M⊗M⊗∂I1∂C+W44M⊗M⊗∂I4∂C+W61N⊗N⊗∂I1∂C+W66N⊗N⊗∂I6∂C)−4C33C−1⊗(W11∂I1∂C+W14∂I4∂C+W16∂I6∂C)−4W1C−1⊗∂C33∂C−4W1C33∂C−1∂C$

where Display Formula

(A5)$W11=∂2W(I1,I4,I6,J)∂I12=c10c012 exp[c01(I1−3)]+k1κ2[(1+2k2ε12)exp(k2ε12)+(1+2k2ε22)exp(k2ε22)]W44=∂2W(I1,I4,I6,J)∂I42=k1(1−3κ)2(1+2k2ε12)exp(k2ε12)W66=∂2W(I1,I4,I6,J)∂I62=k1(1−3κ)2(1+2k2ε22)exp(k2ε22)W14=∂2W(I1,I4,I6,J)∂I1∂I4=k1κ(1−3κ)(1+2k2ε12)exp(k2ε12)W16=∂2W(I1,I4,I6,J)∂I1∂I6=k1κ(1−3κ)(1+2k2ε22)exp(k2ε22)W41=W14W61=W16$

It should be noted that the incompressibility condition in planar biaxial experiment implies $C33=1/(C11C22−C122)$. This condition should apply in all the evaluations involved in Eq. (A4). It follows: Display Formula

(A6)$∂C33∂C11=−C22C332∂C33∂C22=−C11C332∂C33∂C12=∂C33∂C21=C12C332∂I1(C)∂C11=∂I1∂C11+∂I1∂C33∂C33∂C11=1−C22C332∂I1(C)∂C22=∂I1∂C22+∂I1∂C33∂C33∂C22=1−C11C332∂I1(C)∂C12=∂I1(C)∂C21=∂I1∂C33∂C33∂C12=C12C332∂I4(C,M)∂C11=∂I4∂C11+∂I4∂C33∂C33∂C11=M12∂I4(C,M)∂C22=∂I4∂C22+∂I4∂C33∂C33∂C22=M22∂I4(C,M)∂C12=∂I4(C,M)∂C21=∂I4∂C12+∂I4∂C33∂C33∂C12=M1M2∂I6(C,N)∂C11=∂I6∂C11+∂I6∂C33∂C33∂C11=N12∂I6(C,N)∂C22=∂I6∂C22+∂I6∂C33∂C33∂C22=N22∂I6(C,N)∂C12=∂I6(C,N)∂C21=∂I6∂C12+∂I6∂C33∂C33∂C12=N1N2(∂C−1∂C)abcd=−(C−1)am∂Cmn∂Ccd(C−1)nb$

## References

Gasser, T. C. , Ogden, R. W. , and Holzapfel, G. A. , 2006, “ Hyperelastic Modelling of Arterial Layers With Distributed Collagen Fibre Orientations,” J. R. Soc. Interface, 3(6), pp. 15–35. [PubMed]
Fan, R. , and Sacks, M. S. , 2014, “ Simulation of Planar Soft Tissues Using a Structural Constitutive Model: Finite Element Implementation and Validation,” J. Biomech., 47(9), pp. 2043–2054. [PubMed]
Holzapfel, G. A. , and Gasser, T. C. , 2001, “ A Viscoelastic Model for Fiber-Reinforced Composites at Finite Strains: Continuum Basis, Computational Aspects and Applications,” Comput. Methods Appl. Mech. Eng., 190(34), pp. 4379–4403.
Genet, M. , Rausch, M. K. , Lee, L. C. , Choy, S. , Zhao, X. , Kassab, G. S. , Kozerke, S. , Guccione, J. M. , and Kuhl, E. , 2015, “ Heterogeneous Growth-Induced Prestrain in the Heart,” J. Biomech., 48(10), pp. 2080–2089. [PubMed]
Lubarda, V. A. , and Hoger, A. , 2002, “ On the Mechanics of Solids With a Growing Mass,” Int. J. Solids Struct., 39(18), pp. 4627–4664.
Nordsletten, D. , McCormick, M. , Kilner, P. J. , Hunter, P. , Kay, D. , and Smith, N. P. , 2011, “ Fluid–Solid Coupling for the Investigation of Diastolic and Systolic Human Left Ventricular Function,” Int. J. Numer. Methods Biomed. Eng., 27(7), pp. 1017–1039.
Martin, C. , and Sun, W. , 2013, “ Modeling of Long-Term Fatigue Damage of Soft Tissue With Stress Softening and Permanent Set Effects,” Biomech. Model. Mechanobiol., 12(4), pp. 645–655. [PubMed]
Haj-Ali, R. , Dasi, L. P. , Kim, H. S. , Choi, J. , Leo, H. , and Yoganathan, A. P. , 2008, “ Structural Simulations of Prosthetic Tri-Leaflet Aortic Heart Valves,” J. Biomech., 41(7), pp. 1510–1519. [PubMed]
Kim, H. , Lu, J. , Sacks, M. S. , and Chandran, K. B. , 2008, “ Dynamic Simulation of Bioprosthetic Heart Valves Using a Stress Resultant Shell Model,” Ann. Biomed. Eng., 36(2), pp. 262–275. [PubMed]
Shi, Y. , Yao, J. , Young, J. M. , Fee, J. A. , Perucchio, R. , and Taber, L. A. , 2014, “ Bending and Twisting the Embryonic Heart: A Computational Model for C-Looping Based on Realistic Geometry,” Front. Physiol., 5, p. 297. [PubMed]
Sun, W. , Li, K. , and Sirois, E. , 2010, “ Simulated Elliptical Bioprosthetic Valve Deformation: Implications for Asymmetric Transcatheter Valve Deployment,” J. Biomech., 43(16), pp. 3085–3090. [PubMed]
Venkatasubramaniam, A. , Fagan, M. , Mehta, T. , Mylankal, K. , Ray, B. , Kuhan, G. , Chetter, I. C. , and McCollum, P. T. , 2004, “ A Comparative Study of Aortic Wall Stress Using Finite Element Analysis for Ruptured and Non-Ruptured Abdominal Aortic Aneurysms,” Eur. J. Vasc. Endovasc. Surg., 28(2), pp. 168–176. [PubMed]
ABAQUS, 2011, “ ABAQUS/Standard: User's Manual,” Dassault Systèmes Simulia, Johnston, RI.
Simo, J. , and Hughes, T. J. R. , 1998, Computational Inelasticity, Springer, New York.
Prot, V. , Skallerud, B. , and Holzapfel, G. , 2007, “ Transversely Isotropic Membrane Shells With Application to Mitral Valve Mechanics. Constitutive Modelling and Finite Element Implementation,” Int. J. Numer. Methods Eng., 71(8), pp. 987–1008.
ABAQUS, 2011, “ ABAQUS Theory Manual, ABAQUS 611 Documentation,” Dassault Systèmes Simulia, Johnston, RI.
Miehe, C. , 1996, “ Numerical Computation of Algorithmic (Consistent) Tangent Moduli in Large-Strain Computational Inelasticity,” Comput. Methods Appl. Mech. Eng., 134(3–4), pp. 223–240.
Sun, W. , Chaikof, E. L. , and Levenston, M. E. , 2008, “ Numerical Approximation of Tangent Moduli for Finite Element Implementations of Nonlinear Hyperelastic Material Models,” ASME J. Biomech. Eng., 130(6), p. 061003.
Liu, H. , and Sun, W. , 2015, “ Computational Efficiency of Numerical Approximations of Tangent Moduli for Finite Element Implementation of a Fiber-Reinforced Hyperelastic Material Model,” Comput. Methods Biomech. Biomed. Eng., 19(11), pp. 1–10.
Tanaka, M. , and Fujikawa, M. , 2011, “ Numerical Approximation of Consistent Tangent Moduli Using Complex-Step Derivative and Its Application to Finite Deformation Problems,” Trans. Jpn. Soc. Mech. Eng., Ser. A, 77(773), pp. 27–38.
Holzapfel, G. A. , 2000, Nonlinear Solid Mechanics: A Continuum Approach for Engineering, Wiley, Chichester, UK.
Davis, F. M. , Luo, Y. , Avril, S. , Duprey, A. , and Lu, J. , 2015, “ Pointwise Characterization of the Elastic Properties of Planar Soft Tissues: Application to Ascending Thoracic Aneurysms,” Biomech. Model. Mechanobiol., 14(5), pp. 967–978. [PubMed]
Schiavone, A. , and Zhao, L. G. , 2015, “ A Study of Balloon Type, System Constraint and Artery Constitutive Model Used in Finite Element Simulation of Stent Deployment,” Mech. Adv. Mater. Modern Processes, 1(1), pp. 1–15.
Smoljkić, M. , Vander Sloten, J. , Segers, P. , and Famaey, N. , 2015, “ Non-Invasive, Energy-Based Assessment of Patient-Specific Material Properties of Arterial Tissue,” Biomech. Model. Mechanobiol., 14(5), pp. 1045–1056. [PubMed]
Sun, W. , Chaikof, E. L. , and Levenston, M. E. , 2008, “ Development and Finite Element Implementation of a Nearly Incompressible Structural Constitutive Model for Artery Substitute Design,” ASME Paper No. SBC2008-193164.
Sun, W. , and Sacks, M. S. , 2005, “ Finite Element Implementation of a Generalized Fung-Elastic Constitutive Model for Planar Soft Tissues,” Biomech. Model. Mechanobiol., 4(2), pp. 190–199. [PubMed]
Sacks, M. , 1999, “ A Method for Planar Biaxial Mechanical Testing That Includes In-Plane Shear,” ASME J. Biomech. Eng., 121(5), pp. 551–555.
Sacks, M. S. , Smith, D. B. , and Hiester, E. D. , 1997, “ A Small Angle Light Scattering Device for Planar Connective Tissue Microstructural Analysis,” Ann. Biomed. Eng., 25(4), pp. 678–689. [PubMed]
Tanaka, M. , Fujikawa, M. , Balzani, D. , and Schröder, J. , 2014, “ Robust Numerical Calculation of Tangent Moduli at Finite Strains Based on Complex-Step Derivative Approximation and Its Application to Localization Analysis,” Comput. Methods Appl. Mech. Eng., 269, pp. 454–470.
Zulliger, M. A. , Fridez, P. , Hayashi, K. , and Stergiopulos, N. , 2004, “ A Strain-Energy Function for Arteries Accounting for Wall Composition and Structure,”J. Biomech., 37(7), pp. 989–1000. [PubMed]
View article in PDF format.

## References

Gasser, T. C. , Ogden, R. W. , and Holzapfel, G. A. , 2006, “ Hyperelastic Modelling of Arterial Layers With Distributed Collagen Fibre Orientations,” J. R. Soc. Interface, 3(6), pp. 15–35. [PubMed]
Fan, R. , and Sacks, M. S. , 2014, “ Simulation of Planar Soft Tissues Using a Structural Constitutive Model: Finite Element Implementation and Validation,” J. Biomech., 47(9), pp. 2043–2054. [PubMed]
Holzapfel, G. A. , and Gasser, T. C. , 2001, “ A Viscoelastic Model for Fiber-Reinforced Composites at Finite Strains: Continuum Basis, Computational Aspects and Applications,” Comput. Methods Appl. Mech. Eng., 190(34), pp. 4379–4403.
Genet, M. , Rausch, M. K. , Lee, L. C. , Choy, S. , Zhao, X. , Kassab, G. S. , Kozerke, S. , Guccione, J. M. , and Kuhl, E. , 2015, “ Heterogeneous Growth-Induced Prestrain in the Heart,” J. Biomech., 48(10), pp. 2080–2089. [PubMed]
Lubarda, V. A. , and Hoger, A. , 2002, “ On the Mechanics of Solids With a Growing Mass,” Int. J. Solids Struct., 39(18), pp. 4627–4664.
Nordsletten, D. , McCormick, M. , Kilner, P. J. , Hunter, P. , Kay, D. , and Smith, N. P. , 2011, “ Fluid–Solid Coupling for the Investigation of Diastolic and Systolic Human Left Ventricular Function,” Int. J. Numer. Methods Biomed. Eng., 27(7), pp. 1017–1039.
Martin, C. , and Sun, W. , 2013, “ Modeling of Long-Term Fatigue Damage of Soft Tissue With Stress Softening and Permanent Set Effects,” Biomech. Model. Mechanobiol., 12(4), pp. 645–655. [PubMed]
Haj-Ali, R. , Dasi, L. P. , Kim, H. S. , Choi, J. , Leo, H. , and Yoganathan, A. P. , 2008, “ Structural Simulations of Prosthetic Tri-Leaflet Aortic Heart Valves,” J. Biomech., 41(7), pp. 1510–1519. [PubMed]
Kim, H. , Lu, J. , Sacks, M. S. , and Chandran, K. B. , 2008, “ Dynamic Simulation of Bioprosthetic Heart Valves Using a Stress Resultant Shell Model,” Ann. Biomed. Eng., 36(2), pp. 262–275. [PubMed]
Shi, Y. , Yao, J. , Young, J. M. , Fee, J. A. , Perucchio, R. , and Taber, L. A. , 2014, “ Bending and Twisting the Embryonic Heart: A Computational Model for C-Looping Based on Realistic Geometry,” Front. Physiol., 5, p. 297. [PubMed]
Sun, W. , Li, K. , and Sirois, E. , 2010, “ Simulated Elliptical Bioprosthetic Valve Deformation: Implications for Asymmetric Transcatheter Valve Deployment,” J. Biomech., 43(16), pp. 3085–3090. [PubMed]
Venkatasubramaniam, A. , Fagan, M. , Mehta, T. , Mylankal, K. , Ray, B. , Kuhan, G. , Chetter, I. C. , and McCollum, P. T. , 2004, “ A Comparative Study of Aortic Wall Stress Using Finite Element Analysis for Ruptured and Non-Ruptured Abdominal Aortic Aneurysms,” Eur. J. Vasc. Endovasc. Surg., 28(2), pp. 168–176. [PubMed]
ABAQUS, 2011, “ ABAQUS/Standard: User's Manual,” Dassault Systèmes Simulia, Johnston, RI.
Simo, J. , and Hughes, T. J. R. , 1998, Computational Inelasticity, Springer, New York.
Prot, V. , Skallerud, B. , and Holzapfel, G. , 2007, “ Transversely Isotropic Membrane Shells With Application to Mitral Valve Mechanics. Constitutive Modelling and Finite Element Implementation,” Int. J. Numer. Methods Eng., 71(8), pp. 987–1008.
ABAQUS, 2011, “ ABAQUS Theory Manual, ABAQUS 611 Documentation,” Dassault Systèmes Simulia, Johnston, RI.
Miehe, C. , 1996, “ Numerical Computation of Algorithmic (Consistent) Tangent Moduli in Large-Strain Computational Inelasticity,” Comput. Methods Appl. Mech. Eng., 134(3–4), pp. 223–240.
Sun, W. , Chaikof, E. L. , and Levenston, M. E. , 2008, “ Numerical Approximation of Tangent Moduli for Finite Element Implementations of Nonlinear Hyperelastic Material Models,” ASME J. Biomech. Eng., 130(6), p. 061003.
Liu, H. , and Sun, W. , 2015, “ Computational Efficiency of Numerical Approximations of Tangent Moduli for Finite Element Implementation of a Fiber-Reinforced Hyperelastic Material Model,” Comput. Methods Biomech. Biomed. Eng., 19(11), pp. 1–10.
Tanaka, M. , and Fujikawa, M. , 2011, “ Numerical Approximation of Consistent Tangent Moduli Using Complex-Step Derivative and Its Application to Finite Deformation Problems,” Trans. Jpn. Soc. Mech. Eng., Ser. A, 77(773), pp. 27–38.
Holzapfel, G. A. , 2000, Nonlinear Solid Mechanics: A Continuum Approach for Engineering, Wiley, Chichester, UK.
Davis, F. M. , Luo, Y. , Avril, S. , Duprey, A. , and Lu, J. , 2015, “ Pointwise Characterization of the Elastic Properties of Planar Soft Tissues: Application to Ascending Thoracic Aneurysms,” Biomech. Model. Mechanobiol., 14(5), pp. 967–978. [PubMed]
Schiavone, A. , and Zhao, L. G. , 2015, “ A Study of Balloon Type, System Constraint and Artery Constitutive Model Used in Finite Element Simulation of Stent Deployment,” Mech. Adv. Mater. Modern Processes, 1(1), pp. 1–15.
Smoljkić, M. , Vander Sloten, J. , Segers, P. , and Famaey, N. , 2015, “ Non-Invasive, Energy-Based Assessment of Patient-Specific Material Properties of Arterial Tissue,” Biomech. Model. Mechanobiol., 14(5), pp. 1045–1056. [PubMed]
Sun, W. , Chaikof, E. L. , and Levenston, M. E. , 2008, “ Development and Finite Element Implementation of a Nearly Incompressible Structural Constitutive Model for Artery Substitute Design,” ASME Paper No. SBC2008-193164.
Sun, W. , and Sacks, M. S. , 2005, “ Finite Element Implementation of a Generalized Fung-Elastic Constitutive Model for Planar Soft Tissues,” Biomech. Model. Mechanobiol., 4(2), pp. 190–199. [PubMed]
Sacks, M. , 1999, “ A Method for Planar Biaxial Mechanical Testing That Includes In-Plane Shear,” ASME J. Biomech. Eng., 121(5), pp. 551–555.
Sacks, M. S. , Smith, D. B. , and Hiester, E. D. , 1997, “ A Small Angle Light Scattering Device for Planar Connective Tissue Microstructural Analysis,” Ann. Biomed. Eng., 25(4), pp. 678–689. [PubMed]
Tanaka, M. , Fujikawa, M. , Balzani, D. , and Schröder, J. , 2014, “ Robust Numerical Calculation of Tangent Moduli at Finite Strains Based on Complex-Step Derivative Approximation and Its Application to Localization Analysis,” Comput. Methods Appl. Mech. Eng., 269, pp. 454–470.
Zulliger, M. A. , Fridez, P. , Hayashi, K. , and Stergiopulos, N. , 2004, “ A Strain-Energy Function for Arteries Accounting for Wall Composition and Structure,”J. Biomech., 37(7), pp. 989–1000. [PubMed]

## Figures

Fig. 1

FE meshes used in the biaxial test simulation: the arrows denote the loading directions with their round roots showing the sites of loading, the box enclosed region in the center represents the area delineated by the four markers in the biaxial test (a), Von Mises stress contour in the deformed configuration (b), and curves of Cauchy stress versus Green strain resulted from experiment (dashed line) and numerical simulation (solid line) (c)

Fig. 2

Relative error of the Jacobian approximation at different logarithmic perturbations (a), the CPU time (b), and the iteration number (c) at different logarithmic perturbations using themodified HGO model in the biaxial testing simulation. Dashed lines (—) represent the level achieved using UHYPER subroutine.

Fig. 3

FE meshes used in the vessel inflation simulation, before (a) and after (b) deformation; loading pressure versus the outermost radius of the deformed vessel, observed in experiment (dotted line) and numerical simulations. Material models with various mean fiber directions while the other parameters are kept the same as baseline values (c); material models with various standard deviations while the other parameters are kept the same as baseline values (d); and the iteration number at different logarithmic perturbations using the distributed fiber model in the rat carotid artery simulation (e).

## Tables

Table 1 Norm of the largest force residual at each iteration (taken at the load step #32 of the simulation)

## Errata

Some tools below are only available to our subscribers or users with an online account.

### Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related Proceedings Articles
Related eBook Content
Topic Collections