This article was downloaded by: [Tulane University] On: 29 September 2014, At: 12:56 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK
Computer Methods in Biomechanics and Biomedical Engineering Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/gcmb20
Anisotropic hyperelastic behavior of soft biological tissues a
a
Z.-W. Chen , P. Joli & Z.-Q. Feng
ab
a
Laboratoire de Mécanique et d'Énergétique d'Évry, Université d'Évry-Val d'Essonne, Évry, France b
School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu, P.R. China Published online: 15 Aug 2014.
To cite this article: Z.-W. Chen, P. Joli & Z.-Q. Feng (2014): Anisotropic hyperelastic behavior of soft biological tissues, Computer Methods in Biomechanics and Biomedical Engineering, DOI: 10.1080/10255842.2014.915082 To link to this article: http://dx.doi.org/10.1080/10255842.2014.915082
PLEASE SCROLL DOWN FOR ARTICLE Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content. This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http:// www.tandfonline.com/page/terms-and-conditions
Computer Methods in Biomechanics and Biomedical Engineering, 2014 http://dx.doi.org/10.1080/10255842.2014.915082
Anisotropic hyperelastic behavior of soft biological tissues Z.-W. Chena, P. Jolia and Z.-Q. Fenga,b* Laboratoire de Me´canique et d’E´nerge´tique d’E´vry, Universite´ d’E´vry-Val d’Essonne, E´vry, France; b School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu, P.R. China
a
Downloaded by [Tulane University] at 12:56 29 September 2014
(Received 13 January 2014; accepted 10 April 2014) Constitutive laws are fundamental to the studies of the mechanically dominated clinical interventions involving soft biological tissues which show a highly anisotropic hyperelastic mechanical properties. The purpose of this paper was to develop an improved constitutive law based on the Holzapfel –Gasser – Ogden’s model: to replace the isotropic part with Gent constitutive law so as to model the noncollagenous matrix of the media due to its generality and capability to reproduce the Neo-Hookean model. This model is implemented into an in-house finite element program. A uniaxial tension test is considered to study the influence of material parameter J m in Gent model and b which represents the angle between the collagen fibers and the circumferential direction. A simulation of an adventitial strip specimen under tension is performed to show the applicability of this constitutive law. Keywords: anisotropic hyperelasticity; Gent model; HGO model; soft biological tissues; uniaxial tension
1.
Introduction
Soft biological tissues such as ligaments, tendons, pelvic organs, or arterial walls are usually subjected to large deformations with negligible volume change and show a highly nonlinear anisotropic mechanical properties (Weiss et al. 1996; Almeida and Spilker 1998; Ru¨ter and Stein 2000; Pen˜a et al. 2006; Holzapfel and Ogden 2006). To determine the strain and stress in the soft biological tissues, Spencer (1954) first developed the framework of continuum mechanics by means of the definition of a strain energy density function expressed in terms of kinematic invariants to model the purely elastic response of these tissues. This approach was further adapted to finite element formulations of soft collagenous biological tissues (see e.g. Weiss et al. 1996; Pen˜a et al. 2007) for ligaments and (Holzapfel et al. 2000 for arteries). Several constitutive models have been proposed based on Spencer’s theoretical approach for fiber-reinforced composites (Spencer 1972, 1984). Most energy densities used to model transversely isotropic and orthotropic soft tissues take a power law form (Schro¨der et al. 2005) or present an exponential behavior (Fung et al. 1979; Holzapfel et al. 2000). Modeling of connective tissues often includes collagen fibers which play an important role in strengthening the tissue’s stiffness. These fiber are loaded in tension and buckle under compression (Holzapfel et al. 2004). Therefore, the inter-fiber arrangement mainly influences the mechanical behavior of connective tissues in tension. A theoretical study has been performed with Holzapfel –Gasser – Ogden’s (HGO) constitutive law and applied to an unconstrained tension test by Peyraut et al.
*Corresponding author. Email: [email protected] q 2014 Taylor & Francis
(2010). They performed a finite element analysis which has a fair agreement with theoretical solutions. Holzapfel et al. (2000) and Eberlein et al. (2001) proposed an anisotropic hyperelastic fiber-reinforced constitutive model for soft tissues. It is usually assumed that anisotropy is due to the collagen fibers behavior (Gasser et al. 2006), while the ground substance, or matrix, behaves in an isotropic manner, so the energy densities modeling transversely isotropic and orthotropic soft tissues are separated into isotropic and anisotropic parts (Weiss et al. 1996; Balzani et al. 2006). Each anisotropic density refers to a preferred direction of the material. For example, in the Holzapfel– Gasser – Ogden’s constitutive law (Holzapfel et al. 2000; Holzapfel and Gasser 2001), there were two transversely isotropic energies with two distinct preferred directions corresponding to two superposed collagen fiber families. In this paper, we propose using Gent constitutive law (Gent 1996) as the isotropic part to model the noncollagenous matrix of the media due to its generality and capability to reproduce the Neo-Hookean model, combining with Holzapfel –Gasser – Ogden’s constitutive law to represent the hyperelastic anisotropic behavior of soft biological tissues. Our intention is to present numerical investigations on the anisotropic hyperelastic models under large deformation. This particular model has been implemented into an in-house finite element code FER (Feng 2008). A uniaxial tension test is considered to study the influence of material parameter J m in the Gent model and b which represents the angle between the collagen fibers and the circumferential direction.
2
Z.-W. Chen et al. two fiber families: a 1 ¼ fcosðbÞ; sinðbÞ; 0}T a 2 ¼ fcosðbÞ; 2sinðbÞ; 0}T :
Figure 1.
Angle b for circumferentially oriented strip.
Downloaded by [Tulane University] at 12:56 29 September 2014
A simulatin of an adventitial strip specimen which is adopted in (Gasser et al. 2006) under tension is performed to show the applicability of this constitutive law.
2. Anisotropic hyperelastic constitutive law To investigate internal deformation and stress of soft biological tissues such as ligaments, tendons, or arterial walls, anisotropic hyperelastic constitutive laws are often used in the framework of finite element analysis (Weiss et al. 1996; Almeida and Spilker 1998; Ru¨ter and Stein 2000). Balzani et al. (2006) have proposed polyconvex strain energy functions combining an exponential form with a power law to take care of the tissues behavior in the low load domain. More realistic models have also been recently developed to capture the inter-fiber angle change by adding to the strain energy the contribution of the fiber – matrix shear interaction (Peng et al. 2006). In general, the anisotropy can be represented via the introduction of a so-called structural tensor, which allows a coordinate-invariant formulation on the constitutive equations (Spencer 1987; Boehler 1987; Zheng and Spencer 1993). 2.1 Anisotropic hyperelasticity Soft biological tissues are usually taken to be hyperelastic and often undergo large deformations. It is usually assumed that anisotropy is due to the collagen fibers behavior (Gasser et al. 2006), while the ground substance, or matrix, behaves in an isotropic manner, so the energy densities modeling transversely isotropic and orthotropic soft tissues are separated into isotropic and anisotropic parts (Weiss et al. 1996; Balzani et al. 2006). W ¼ W iso þ
n X
W aani
ð1Þ
ð2Þ
The phenomenological angle b represents the angle between the collagen fibers and the circumferential direction for a strip extracted, for example, from the media of a human abdominal aorta (Figure 1). b is predicted to be 43839 in the work of Balzani et al. (2006). In our work, we consider b as a model parameter varying from 08 to 908. In the continuation, the cosine and sine of this angle will be noted, respectively, c and s. In order to describe the geometrical transformation problems, the deformation gradient tensor is introduced by F¼
›x ›u ¼Iþ : ›X ›X
ð3Þ
X; x, and u represent, respectively, the reference and the current positions and the displacement vector of a material point. The constraint of incompressibility (isochoric deformation) is given by (Ogden 1984) J ¼ detðFÞ ¼ 1:
ð4Þ
Because of large displacements and rotations, Green – Lagrangian strain is adopted for the nonlinear relationships between strains and displacements. We note C the stretch tensor or the right Cauchy –Green strain tensor (C ¼ F T F). The Green –Lagrangian strain tensor E is defined by 1 E ¼ ðC 2 IÞ: 2
ð5Þ
According to the Zhang – Rychlewski’s theorem (Zhang and Rychlewski 1990), the condition of material symmetry is satisfied if structural tensors are additionally included in the strain energy density representation. Transversely isotropic densities can then be expressed with the three invariants I 1 ; I 2 , and I 3 of the right Cauchy – Green deformation tensor and two additional mixed invariants J 4 and J 5 (Spencer 1987; Boehler 1987; Zheng and Spencer 1993).
a¼1
Each anisotropic density W aani refers to a preferred direction of the material. The number of fiber families n is generally set to 1 to model tissues as ligaments or tendons while it is set to 2 to represent the behavior of arterial walls. For example, to model the embedded collagen fibers of soft biological arterial tissues, HGO’s constitutive law (Holzapfel et al. 2000; Gasser et al. 2006) superposes two transversely isotropic energies with two distinct preferred directions a 1 and a 2 corresponding to
I 1 ¼ trðCÞ; I 2 ¼ trðcof ðCÞÞ; I 3 ¼ detðCÞ; J 4 ¼ trðCMÞ; J 5 ¼ trðC 2 MÞ;
ð6Þ
where cof ðCÞ denotes the co-factor matrix of C, and M is the so-called structural tensor representing the transverseisotropy group and referring to a preferred direction a of the material: M ¼ a ^ a:
ð7Þ
Computer Methods in Biomechanics and Biomedical Engineering
3
The first strain invariant of C iso is defined by I^1 ¼ J 22=3 I 1 which replaces I 1 in WðI 1 Þ. ^ ¼ F^ T F^ ¼ J 22=3 C The modified invariants related to C are expressed from Equation (6) by a I^1 ¼ I 1 J 22=3 ; I^2 ¼ I 2 J 24=3 ; J^ 4 ¼ J a4 J 22=3 ; a J^ 5 ¼ J a5 J 24=3 :
ð13Þ
The exponential type HGO density adopted in this work uses these modified invariants as follows:
Downloaded by [Tulane University] at 12:56 29 September 2014
Figure 2.
Cylindrical coordinate system.
^ I^1 J^ a Þ þ W H ðJÞ; W ¼ Wð 4
ð14Þ
1 J2 2 1 2 lnJ ; W H ðJÞ ¼ d 2
ð15Þ
It is noted that Equations (6) and (7) give J 4 ¼ trðF Fa ^ aÞ ¼ kFak : 2
T
^ I^1 ; J^ a Þ ¼ W iso ðI^1 Þ þ Wð 4
ð8Þ
2 X
W ani ðJ^ 4 Þ; a
ð16Þ
a¼1
The double brackets represent the usual Euclidian norm. The square root of J 4 represents thus the stretch in the fiber direction. It can also be interpreted as the radial coordinate of Fa in a cylindrical coordinate system where the polar angle u represents the deformed angle between the collagen fibers and the circumferential direction (Figure 2). In the case of hyperelastic law, there exists an elastic potential function W (or strain energy density function) which is a scale function of one of the strain tensors, whose derivative with respect to a strain component determines the corresponding stress component. This can be expressed by S¼
›W ›W ¼2 ; ›E ›C
ð9Þ
where S is the second Piola– Kirchoff stress tensor. For numerical purpose, it proves useful to separate deformation in volumetric and isochoric parts by a multiplicative split of a deformation gradient as in (Flory 1961): F ¼ F iso F vol ;
ð10Þ
W aani
8 i k1 h k2 ðJ^a4 21Þ2 > < e 21 ¼ 2k2 > : 0
a if J^ 4 $ 1; a if J^ 4 , 1
:
ð17Þ
Generally, soft biological tissues are assumed to be incompressible. Equation (15), which was proposed by Horgan and Saccomandi (2003), represents a penalty term added to the finite element model to account for the incompressible behavior of the material. Here, d is the material incompressibility parameter. The initial bulk modulus K is defined as K ¼ 2=d. Some other models can be found in (Horgan and Saccomandi 2004, 2006). This a energy density is case-sensitive with respect to J^ 4 because a the case of J^ 4 , 1 represents the shortening of the fibers which is assumed to generate no stress. The proof of convexity of Equation (17) with respect to F is given in (Schro¨der et al. 2005; Balzani et al. 2006). The noncollagenous matrix of the media is modeled by the Gent isotropic density (Gent 1996; Feng et al. 2010): W iso ðI^1 Þ ¼ 2
I^1 2 3 m J m ln 1 2 ; Jm 2
ð18Þ
where F
vol
¼J
1=3
I;
F
iso
¼J
21=3
F:
ð11Þ
This decomposition is such that detðF iso Þ ¼ 1. It is easy to see that F and F iso have the same eigenvectors. The isochoric part of the right Cauchy – Green strain tensor C can then be defined as C iso ¼ J 22=3 C:
ð12Þ
where m is the shear modulus and J m is the constant limiting value for I^1 2 3. It is noted that for J m ! 1, we recover the well-known Neo-Hookean strain energy density: WðI 1 Þ ¼
m ðI 1 2 3Þ: 2
ð19Þ
It is noted that the volumetric –isochoric split of the above HGO model does only hold for (quasi)incompres-
4
Z.-W. Chen et al.
sible deformations. An extension to compressible deformations would require that the volumetric part of the strain energy function includes a dependency on the structural tensor. This is proved recently by Guo et al. (2008) where a simple compressible anisotropic analytical model is developed. We remind finally that by deriving W from Equation (9) and introducing the matrix of cofactors of C, Cof ðCÞ ¼ I 3 C 2T , it is conventionally obtained:
›W ›W ›W ›W 1 Iþ ðI 1 I 2 CÞ þ cof ðCÞ þ 1 M a ›I 1 ›I 2 ›I 3 ›J 4 ›W ›W ›W 2 1 1 2 2 þ 2 M a þ 1 ðCM a þ M a CÞ þ 2 ðCM a þ M a CÞ : ›J 4 ›J 5 ›J 5 ð20Þ
Downloaded by [Tulane University] at 12:56 29 September 2014
S¼2
In our particular case, Equation (20) is reduced to ›W ›W ›W a 1 ›W a 2 S¼2 Iþ cof ðCÞ þ 1 M þ 2 M : ð21Þ ›I 1 ›I 3 ›J 4 ›J 4 The derivatives of the density energy W with respect to the invariants are evaluated with equations (17,18):
›W 21=3 dW iso ¼ I3 ; ›I 1 d I^1
1 1 12 2 ; 2d J
ð22Þ
" # 2 ›W 1 ^ dW iso X a dW ani I1 J^ 4 ¼2 þ a 3I 3 ›I 3 d I^1 dJ^ 4 a¼1 þ
object-oriented programming language. The geometrically nonlinear analysis may be described by using the total or the updated Lagrangian formulations. The total Lagrangian formulation is derived with respect to the initial configuration. The updated Lagrangian formulation is derived with respect to the current configuration. In other words, the total Lagrangian formulation constructs the tangent stiffness matrix with respect to the initial configuration. On the other hand, the updated Lagrangian formulation constructs the tangent stiffness matrix with respect to the current configuration. The updated Lagrangian formulation is computationally effective (Zienkiewicz and Taylor 1991) because it does not include the initial displacement matrix. In the total Lagrangian formulation, the initial configuration remains constant. This simplifies the computation (Crisfield 1991). Therefore, the total Lagrangian formulation was selected in this work to describe nonlinear behavior. In the context of the finite element method and with Equations (3) and (5), the Green – Lagrangian strain can be formally written with linear and nonlinear terms in function of nodal displacements (Simo and Hughes 1998):
ð23Þ
E¼
1 B L þ B NL ðuÞ u; 2
where B L is the matrix which relates the linear strain term to the nodal displacements and B NL ðuÞ the matrix which relates the nonlinear strain term to the nodal displacements. From Equation (28), the incremental form of the strain – displacement relationship is
dE ¼ ðB L þ B NL ðuÞÞdu:
›W 21=3 dW ani ¼ I3 a ; ›J a4 dJ^ 4
ð24Þ
dW iso mJ m ; ¼ dI^1 2ðJ m 2 I^1 þ 3Þ
ð25Þ
ð28Þ
ð29Þ
Using the principle of virtual displacement, the virtual work equation is given as ð
dW ¼
SdEdV 2 F ext du;
ð30Þ
V0 a dW ani k2 ðJ^ 4 21Þ2 ^a : a ¼ k1 ðJ4 2 1Þe ^ dJ
ð26Þ
4
The Cauchy stress (or true stress) tensor s can be calculated from the second Piola – Kirchoff stress tensor S as follows: 1 s ¼ FSF T : J 2.2
ð27Þ
Finite element implementation
This model have been implemented and tested in the finite element code FER. This code is implemented using Cþ þ
where V 0 is the volume of the initial configuration and F ext the vector of external loads. Using Equations (9) and (29), the incremental form of the stress dS can be linked to the incremental form of the strain dE as follows:
dS ¼
›2 W ¼ D : dE ¼ D : ðB L þ B NL ðuÞÞdu; ›E 2
ð31Þ
where D denotes the constitutive tangent matrix. This fourth-order tensor is obtained from the derivative of S
Computer Methods in Biomechanics and Biomedical Engineering
5
90 Gent+HGO constitutive model HGO constitutive model Experiment
80 70 σ [kPa]
60 50 40 30 20 10 0
Downloaded by [Tulane University] at 12:56 29 September 2014
Figure 3.
þ cðM ^ MÞ; ð32Þ with 2mJ m J 24=3 ; ðJ m 2 I^1 þ 3Þ2
b¼
2mJ m J 22=3 ; 3ðJ m 2 I^1 þ 3Þ
c ¼ 4k1 ð1 þ 2k2 ðJ 4 2 1Þ Þe 2
k2 ðJ 4 21Þ2
ð33Þ
:
Substituting dE from Equation (29) into Equation (30) results in ð dW ¼ SðB L þ B NL ðuÞÞdudV 2 F ext du V0
¼ F int du 2 F ext du: The vector of internal forces is defined by ð SðB L þ B NL ðuÞÞdV: F int ¼
ð34Þ
ð35Þ
V0
Since du is arbitrary, a set of nonlinear equations can be obtained as F ext 2 F int ¼ 0
0.05
0.1
0.15 0.2 Δ l/l0
0.25
0.3
0.35
Figure 4. Cauchy stress sðkPaÞ versus strain Dl=l0 of the experiment and the constitutive model.
Unconstrained uniaxial tension test.
with respect to E in Equation (20): I1 I1 D ¼ a I 2 C 21 ^ I 2 C 21 ; 3 3 I1 ›C21 21 21 ; 2 b C ^ I þ ð1 þ ÞI ^ C 2 I 1 3 ›C 2 ›C21 J 2 C 21 ^ C 21 þ ðJ 2 2 1Þ þ d ›C
a¼
0
ð36Þ
This equation is strongly nonlinear, because of finite strains and large displacements of soft biological tissues.
A typical solution procedure for this type of nonlinear analysis is obtained by using the Newton – Raphson iterative procedure (Simo and Hughes 1998; Belytschko et al. 2000). It should be noted that all the models presented here have been implemented in the finite element code FER.
3.
Numerical results
To illustrate the results of the simulation using this model, we consider here two uniaxial tensile numerical tests. First, we consider a uniaxial unconstrained tension test as depicted in Figure 3. This loading is usually applied to soft biological tissues for identifying parameters of the material model. To fit this numerical model with experimental data in (Balzani et al. 2006), the material parameters for the best fit are chosen as m ¼ 17:41 kPa, J m ¼ 60:0 kPa, k1 ¼ 0:0017 kPa, k2 ¼ 882:847, and b ¼ 43:398. These last three values are proposed in (Balzani et al. 2006). Figure 4 illustrates the response of the constitutive model we proposed here (Gent þ HGO constitutive model), as well as the experiment data and response of the HGO constitutive model which are given by Balzani et al. (2006). As can be seen, the match is quite good, even better than HGO constitutive model using Neo-Hookean model, the exponential character can be presented better by using the Gent model. Whereas the Holzapfel– Gasser – Ogden’s constitutive law uses the Neo-Hookean model as the isotropic part, the Gent model has one more material parameter: J m , when it turns to infinity, we can get the Neo-Hookean strain energy density. So it is important to consider the influence of parameter J m . The other material parameters are the same as the last tension test except b ¼ 70:08 and several different values of J m were considered. Figure 5 illustrates the different
6
Z.-W. Chen et al. 24
1.6 Jm=3.0 Jm=5.0 Jm=6.0 Jm=7.0 Jm=10.0
σ [MPa]
16 12
Jm=10.0 Jm=20.0 Jm=50.0 Jm=100.0
1.2 σ [MPa]
20
0.8
8 0.4 4 0
0 0
0.4
0.8
1.2
1.6
2
0
0.4
0.8
Δ l/l0
1.6
2
10 ≤ Jm ≤ 100
Influence of J m on the mechanical response.
et al. (2010) for different values of the initial angle b between the collagen fibers orientation and the circumferential direction. This theoretical study was performed with the HGO constitutive law. A critical angle (bc < 54:738) has been analytically determined. It appears that the correspondence between the solution and J 4 is not one to one if b is greater than bc while uniqueness holds if b is lower than bc . A relationship between l2 and J 4 for different values of b is shown by Figure 7. Three different initial angles 408, 608, and 708 are considered in this test, the material parameters are chosen from (Holzapfel et al. 2000): m ¼ 6:0 kPa, J m ¼ 50:0 kPa, k1 ¼ 2:3632 kPa, and k2 ¼ 0:8393. The evolution of three principal stretches for different values of b is shown in Figure 8. As can be seen from Figure 7, if b , bc (e.g. b ¼ 208), during the whole tension process, there is no
mechanical responses due to the different J m . The influence of J m on the three principal stretches is shown in Figure 6. It seems that when J m , 6, the strength of material increases as J m grows, but the material has not been stretched to the length we expected. When J m . 6, the strength decreases with an increase in J m , especially when 6 , J m , 10, it decreases a lot while when J m . 50, and it stays stable and almost not influenced by J m . It is the same situation when it comes to deformation as shown in Figure 6 (right), when J m . 10, the deformation stays the same. We should notice that the influence of J m on the mechanical behavior of material models is also dependent on the other parameters such as m and b. To evaluate the dependence of the principal stretches on J 4 , the analytical solutions were studied by Peyraut
1.5
1.5
1
1 λ2, λ3
λ2, λ3
Downloaded by [Tulane University] at 12:56 29 September 2014
1 < Jm ≤ 10
Figure 5.
1.2 Δ l/l0
0.5
0.5 Jm=3.0, λ2 Jm=3.0, λ3 Jm=5.0, λ2
Jm=10.0, λ2 Jm=10.0, λ3 Jm=50.0, λ2
Jm=5.0, λ3 Jm=10.0, λ2 Jm=10.0, λ3
0
0 1
1.4
1.8
2.2
2.6
λ1 1 < Jm ≤ 10
Figure 6.
Jm=50.0, λ3 Jm=100.0, λ2 Jm=100.0, λ3
Influence of J m on the principal stretches.
3
1
1.4
1.8
2.2 λ1
10 ≤ Jm ≤ 100
2.6
3
Computer Methods in Biomechanics and Biomedical Engineering
7
1 0.8
λ2
0.6 0.4 β=20° β=60° β=70° β=80°
0.2 0 0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
J4
way that J 4 , 1, so J 4 is always $ 1, which means that the biological tissue is always stretched in the fiber direction and the anisotropic part works from the beginning, which explains why in Figure 8 (left) the two principle stretches l2 and l3 are not the same from the first place. The biological tissue can thus not be shortened in the fiber direction. On the contrary, if b . bc (e.g. b ¼ 708), J 4 can be less than 1 during some time, which means that the HGO model reduces to the isotropic Gent model. As can be seen in Figure 7, this time grows as the initial degree increases, so that it takes more time for b ¼ 708 to transfer into anisotropic behavior than b ¼ 608 as shown in Figure 8 (right). The biological tissue can be shortened in the fiber direction. At last, a numerical model of adventitial strip adopted in (Gasser et al. 2006) is considered to perform a uniaxial tensile test from the axial direction. The material parameters are m ¼ 15:28 kPa, J m ¼ 50:0 kPa, k1 ¼ 996:6 kPa, and k2 ¼ 524:6, and the structure parameter is
Figure 9. time.
Initial configurations and von-Mises stress at last
b ¼ 49:988. Adventitial strip model of length L ¼ 10.0 mm, width W ¼ 3.0 mm, and thickness T ¼ 0.5 mm is considered in this work. This model consists of 3200 hexahedral elements as shown in Figure 9 (left).
6
1.25 β=40°, λ2 β=40°, λ3
5
β=60°, λ2 β=60°, λ3 β=70°, λ2 β=70°, λ3
1
λ2, λ3
4 λ2, λ3
Downloaded by [Tulane University] at 12:56 29 September 2014
Figure 7. l2 versus J 4 for different values of b (Peyraut et al. 2010).
3
0.75
2 0.5 1 0
1
1.4
1.8
2.2
2.6
3
0.25
λ1
Figure 8.
Evolution of three principal stretches for different values of b.
1
1.4
1.8
2.2 λ1
2.6
3
8
Z.-W. Chen et al.
Downloaded by [Tulane University] at 12:56 29 September 2014
Figure 10.
Displacement in direction z.
Figure 9 (right) represents the deformed configuration and computed von-Mises stress of the specimen. According to the structural assumptions in (Holzapfel et al. 2000) and (Zulliger et al. 2004), the collagen is embedded as two families of fibers that are symmetrically disposed relative to the tensile (axial) direction. The embedded collagen fibers need to be rotated towards the direction of the loading until they are able to carry significant load. This leads to an increase in thickness in the middle of specimen as shown in Figure 9 (right). Due to the incompressibility character, the width of the specimen decreases. This numerical result is very similar to that observed in (Gasser et al. 2006). The displacement in the thickness direction (z) is shown in Figure 10. We can see clearly that the displacement in direction z is totally symmetrical.
4.
Conclusion
The development of constitutive laws is fundamental and a prerequisite for understanding the highly nonlinear deformation mechanisms and stress distributions of soft biological tissues under different loading conditions. In this work, the improved Gent þ HGO anisotropic hyperelastic law has been proposed and numerically implemented in the context of uniaxial tension tests. The exponential character of the stress – strain behavior of the considered tissue is modeled by the Gent þ HGO constitutive law and the curve for the uniaxial unconstrained tension test fits the experimental data better than the pure HGO model. An investigation is made on the effect of material parameter J m on the mechanical response and principal stretches. It turns out that when J m is small, the influence is obvious while with an increase in J m , the strain energy density gets more and more like Neo-Hookean model and the influence is paltry. Three different initial angles have been considered in this context; the result has a good agreement with the analytical solutions which was studied by Peyraut et al. (2010). A numerical simulation of adventitial strip has
been performed, in which the numerical result is very similar to that observed in (Gasser et al. 2006).
Funding This work was supported by a French ANR Grant [No. ANR-09SYSC-008] and the Natural Science Foundation of China [No. 11372260].
References Almeida ES, Spilker RL. 1998. Finite element formulation for hyperelastic transversely isotropic biphasic soft tissues. Comput Meth Appl Mech Eng. 151:513 – 538. Balzani D, Neff P, Schro¨der J, Holzapfel GA. 2006. A polyconvex framework for soft biological tissues. Adjustment to experimental data. Int J Solids Struct. 43:6052 – 6070. Belytschko T, Liu WK, Moran B. 2000. Nonlinear finite elements for continua and structures. Chichester: Wiley. Boehler JP. 1987. Introduction to the invariant formulation of anisotropic constitutive equations. In: Boehler JP, editor. Applications of tensor functions in solid mechanics, CISM Course No. 292. Vienna: Springer Verlag. Crisfield MA. 1991. Non-linear finite element analysis of solid and structures. Chichester: Wiley. Eberlein R, Holzapfel GA, Schulze-Bauer CAJ. 2001. An anisotropic model for annulus tissue and enhanced finite element analyses of intact lumbar disc bodies. Comput Methods Biomech Biomed Eng. 4:209– 230. Feng ZQ. 2008. A finite element research system. Available from: http://lmee.univ-evry.fr/, feng/FerSystem.html Feng ZQ, Renaud C, Cros JM, Zhang HW, Guan ZQ. 2010. A finite element finite-strain formulation for modeling colliding blocks of Gent materials. Int J Solids Struct. 47:2215– 2222. Flory PJ. 1961. Thermodynamic relations for high elastic materials. Trans Faraday Soc. 57:829– 838. Fung YC, Fronek K, Patitucci P. 1979. Pseudoelasticity of arteries and the choice of its mathematical expression. Am J Physiol. 237:H620 – H631. Gasser TC, Ogden RW, Holzapfel GA. 2006. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J R Soc Interface. 3:15 – 35. Gent AN. 1996. A new constitutive relation for rubber. Rubber Chem Technol. 69:59– 61.
Downloaded by [Tulane University] at 12:56 29 September 2014
Computer Methods in Biomechanics and Biomedical Engineering Guo ZY, Caner F, Peng XQ, Moran B. 2008. On constitutive modelling of porous neo-Hookean composites. J Mech Phys Solids. 56:2338– 2357. Holzapfel GA, Gasser TC. 2001. A viscoelastic model for fiberreinforced composites at finite strains: continuum basis, computational aspects and applications. Comput Meth Appl Mech Eng. 190:4379– 4403. Holzapfel GA, Gasser TC, Ogden RW. 2000. A new constitutive framework for arterial wall mechanics and a comparative study of material models. J Elast. 61:1 – 48. Holzapfel GA, Gasser TC, Ogden RW. 2004. Comparison of a multilayer structural model for arterial walls with a Fungtype model, and issues of material stability. J Biomech Eng. 126:264 – 275. Holzapfel GA, Ogden RW. 2006. Mechanics of biological tissue. Berlin: Springer. Horgan CO, Saccomandi G. 2003. Finite thermoelasticity with limiting chain extensibility. J Mech Phys Solids. 51:1127 – 1146. Horgan CO, Saccomandi G. 2004. Constitutive models for compressible nonlinearly elastic materials with limiting chain extensibility. J Elast. 77:123– 138. Horgan CO, Saccomandi G. 2006. Phenomenological hyperelastic strain-stiffening constitutive models for rubber. Rubber Chem Technol. 79:152– 169. Ogden RW. 1984. Non-linear elastic deformations. Chichester: Ellis Horwood. Pen˜a E, del Palomar AP, Calvo B, Martı´nez MA, Doblare´ M. 2007. Computational modelling of diarthrodial joints. Physiol, pathological and post-surgery simulations. Arch Comput Meth Eng. 14(1):1 – 54. Pen˜a E, Martı´nez MA, Calvo B, Doblare´ M. 2006. On the numerical treatment of initial strains in soft biological tissues. Int J Numer Meth Eng. 68:836 – 860. Peng XQ, Guo ZY, Moran B. 2006. An anisotropic hyperelastic constitutive model with fiber-matrix shear interaction for the human annulus fibrosus. J Appl Mech. 73:815 – 824.
9
Peyraut F, Feng ZQ, Labed N, Renaud C. 2010. A closed form solution for the uniaxial tension test of biological soft tissues. Int J Nonlinear Mech. 45:535 – 541. Ru¨ter M, Stein E. 2000. Analysis, finite element computation and error estimation in transversely isotropic nearly incompressible finite elasticity. Comput Meth Appl Mech Eng. 190:519 – 541. Schro¨der J, Neff P, Balzani D. 2005. A variational approach for materially stable anisotropic hyperelasticity. Int J Solids Struct. 42:4352 –4371. Simo JC, Hughes TJR. 1998. Computational inelasticity. New York: Springer-Verlag. Spencer AJM. 1954. Theory of invariants. Continuum Physics. New York: Academic Press; p. 239– 253. Spencer AJM. 1972. Deformations of fiber-reinforced materials. New York: Oxford University Press. Spencer AJM. 1984. Continuum theory of the mechanics of fiberreinforced composites. Vienna: Springer-Verlag. Spencer AJM. 1987. Isotropic polynomial invariants and tensor functions. In: Boehler JP, editor. Applications of tensor functions in solids mechanics, CISM Course No. 282. Vienna: Springer Verlag. Weiss JA, Maker BN, Govindjee S. 1996. Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput Meth Appl Mech Eng. 135:107 – 128. Zhang JM, Rychlewski J. 1990. Structural tensors for anisotropic solids. Arch Mech. 42:267 – 277. Zheng QS, Spencer AJM. 1993. Tensors which characterize anisotropies. Int J Eng. Sci. 31(5):679– 693. Zienkiewicz OC, Taylor RL. 1991. The finite element method. Maidenhead Berkshire: McGraw-Hill. Zulliger MA, Fridez P, Hayashi K, Stergiopulos N. 2004. A strain energy function for arteries accounting for wall composition and structure. J Biomech. 37:989 – 1000.
Computer Methods in Biomechanics and Biomedical Engineering Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/gcmb20
Anisotropic hyperelastic behavior of soft biological tissues a
a
Z.-W. Chen , P. Joli & Z.-Q. Feng
ab
a
Laboratoire de Mécanique et d'Énergétique d'Évry, Université d'Évry-Val d'Essonne, Évry, France b
School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu, P.R. China Published online: 15 Aug 2014.
To cite this article: Z.-W. Chen, P. Joli & Z.-Q. Feng (2014): Anisotropic hyperelastic behavior of soft biological tissues, Computer Methods in Biomechanics and Biomedical Engineering, DOI: 10.1080/10255842.2014.915082 To link to this article: http://dx.doi.org/10.1080/10255842.2014.915082
PLEASE SCROLL DOWN FOR ARTICLE Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content. This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http:// www.tandfonline.com/page/terms-and-conditions
Computer Methods in Biomechanics and Biomedical Engineering, 2014 http://dx.doi.org/10.1080/10255842.2014.915082
Anisotropic hyperelastic behavior of soft biological tissues Z.-W. Chena, P. Jolia and Z.-Q. Fenga,b* Laboratoire de Me´canique et d’E´nerge´tique d’E´vry, Universite´ d’E´vry-Val d’Essonne, E´vry, France; b School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu, P.R. China
a
Downloaded by [Tulane University] at 12:56 29 September 2014
(Received 13 January 2014; accepted 10 April 2014) Constitutive laws are fundamental to the studies of the mechanically dominated clinical interventions involving soft biological tissues which show a highly anisotropic hyperelastic mechanical properties. The purpose of this paper was to develop an improved constitutive law based on the Holzapfel –Gasser – Ogden’s model: to replace the isotropic part with Gent constitutive law so as to model the noncollagenous matrix of the media due to its generality and capability to reproduce the Neo-Hookean model. This model is implemented into an in-house finite element program. A uniaxial tension test is considered to study the influence of material parameter J m in Gent model and b which represents the angle between the collagen fibers and the circumferential direction. A simulation of an adventitial strip specimen under tension is performed to show the applicability of this constitutive law. Keywords: anisotropic hyperelasticity; Gent model; HGO model; soft biological tissues; uniaxial tension
1.
Introduction
Soft biological tissues such as ligaments, tendons, pelvic organs, or arterial walls are usually subjected to large deformations with negligible volume change and show a highly nonlinear anisotropic mechanical properties (Weiss et al. 1996; Almeida and Spilker 1998; Ru¨ter and Stein 2000; Pen˜a et al. 2006; Holzapfel and Ogden 2006). To determine the strain and stress in the soft biological tissues, Spencer (1954) first developed the framework of continuum mechanics by means of the definition of a strain energy density function expressed in terms of kinematic invariants to model the purely elastic response of these tissues. This approach was further adapted to finite element formulations of soft collagenous biological tissues (see e.g. Weiss et al. 1996; Pen˜a et al. 2007) for ligaments and (Holzapfel et al. 2000 for arteries). Several constitutive models have been proposed based on Spencer’s theoretical approach for fiber-reinforced composites (Spencer 1972, 1984). Most energy densities used to model transversely isotropic and orthotropic soft tissues take a power law form (Schro¨der et al. 2005) or present an exponential behavior (Fung et al. 1979; Holzapfel et al. 2000). Modeling of connective tissues often includes collagen fibers which play an important role in strengthening the tissue’s stiffness. These fiber are loaded in tension and buckle under compression (Holzapfel et al. 2004). Therefore, the inter-fiber arrangement mainly influences the mechanical behavior of connective tissues in tension. A theoretical study has been performed with Holzapfel –Gasser – Ogden’s (HGO) constitutive law and applied to an unconstrained tension test by Peyraut et al.
*Corresponding author. Email: [email protected] q 2014 Taylor & Francis
(2010). They performed a finite element analysis which has a fair agreement with theoretical solutions. Holzapfel et al. (2000) and Eberlein et al. (2001) proposed an anisotropic hyperelastic fiber-reinforced constitutive model for soft tissues. It is usually assumed that anisotropy is due to the collagen fibers behavior (Gasser et al. 2006), while the ground substance, or matrix, behaves in an isotropic manner, so the energy densities modeling transversely isotropic and orthotropic soft tissues are separated into isotropic and anisotropic parts (Weiss et al. 1996; Balzani et al. 2006). Each anisotropic density refers to a preferred direction of the material. For example, in the Holzapfel– Gasser – Ogden’s constitutive law (Holzapfel et al. 2000; Holzapfel and Gasser 2001), there were two transversely isotropic energies with two distinct preferred directions corresponding to two superposed collagen fiber families. In this paper, we propose using Gent constitutive law (Gent 1996) as the isotropic part to model the noncollagenous matrix of the media due to its generality and capability to reproduce the Neo-Hookean model, combining with Holzapfel –Gasser – Ogden’s constitutive law to represent the hyperelastic anisotropic behavior of soft biological tissues. Our intention is to present numerical investigations on the anisotropic hyperelastic models under large deformation. This particular model has been implemented into an in-house finite element code FER (Feng 2008). A uniaxial tension test is considered to study the influence of material parameter J m in the Gent model and b which represents the angle between the collagen fibers and the circumferential direction.
2
Z.-W. Chen et al. two fiber families: a 1 ¼ fcosðbÞ; sinðbÞ; 0}T a 2 ¼ fcosðbÞ; 2sinðbÞ; 0}T :
Figure 1.
Angle b for circumferentially oriented strip.
Downloaded by [Tulane University] at 12:56 29 September 2014
A simulatin of an adventitial strip specimen which is adopted in (Gasser et al. 2006) under tension is performed to show the applicability of this constitutive law.
2. Anisotropic hyperelastic constitutive law To investigate internal deformation and stress of soft biological tissues such as ligaments, tendons, or arterial walls, anisotropic hyperelastic constitutive laws are often used in the framework of finite element analysis (Weiss et al. 1996; Almeida and Spilker 1998; Ru¨ter and Stein 2000). Balzani et al. (2006) have proposed polyconvex strain energy functions combining an exponential form with a power law to take care of the tissues behavior in the low load domain. More realistic models have also been recently developed to capture the inter-fiber angle change by adding to the strain energy the contribution of the fiber – matrix shear interaction (Peng et al. 2006). In general, the anisotropy can be represented via the introduction of a so-called structural tensor, which allows a coordinate-invariant formulation on the constitutive equations (Spencer 1987; Boehler 1987; Zheng and Spencer 1993). 2.1 Anisotropic hyperelasticity Soft biological tissues are usually taken to be hyperelastic and often undergo large deformations. It is usually assumed that anisotropy is due to the collagen fibers behavior (Gasser et al. 2006), while the ground substance, or matrix, behaves in an isotropic manner, so the energy densities modeling transversely isotropic and orthotropic soft tissues are separated into isotropic and anisotropic parts (Weiss et al. 1996; Balzani et al. 2006). W ¼ W iso þ
n X
W aani
ð1Þ
ð2Þ
The phenomenological angle b represents the angle between the collagen fibers and the circumferential direction for a strip extracted, for example, from the media of a human abdominal aorta (Figure 1). b is predicted to be 43839 in the work of Balzani et al. (2006). In our work, we consider b as a model parameter varying from 08 to 908. In the continuation, the cosine and sine of this angle will be noted, respectively, c and s. In order to describe the geometrical transformation problems, the deformation gradient tensor is introduced by F¼
›x ›u ¼Iþ : ›X ›X
ð3Þ
X; x, and u represent, respectively, the reference and the current positions and the displacement vector of a material point. The constraint of incompressibility (isochoric deformation) is given by (Ogden 1984) J ¼ detðFÞ ¼ 1:
ð4Þ
Because of large displacements and rotations, Green – Lagrangian strain is adopted for the nonlinear relationships between strains and displacements. We note C the stretch tensor or the right Cauchy –Green strain tensor (C ¼ F T F). The Green –Lagrangian strain tensor E is defined by 1 E ¼ ðC 2 IÞ: 2
ð5Þ
According to the Zhang – Rychlewski’s theorem (Zhang and Rychlewski 1990), the condition of material symmetry is satisfied if structural tensors are additionally included in the strain energy density representation. Transversely isotropic densities can then be expressed with the three invariants I 1 ; I 2 , and I 3 of the right Cauchy – Green deformation tensor and two additional mixed invariants J 4 and J 5 (Spencer 1987; Boehler 1987; Zheng and Spencer 1993).
a¼1
Each anisotropic density W aani refers to a preferred direction of the material. The number of fiber families n is generally set to 1 to model tissues as ligaments or tendons while it is set to 2 to represent the behavior of arterial walls. For example, to model the embedded collagen fibers of soft biological arterial tissues, HGO’s constitutive law (Holzapfel et al. 2000; Gasser et al. 2006) superposes two transversely isotropic energies with two distinct preferred directions a 1 and a 2 corresponding to
I 1 ¼ trðCÞ; I 2 ¼ trðcof ðCÞÞ; I 3 ¼ detðCÞ; J 4 ¼ trðCMÞ; J 5 ¼ trðC 2 MÞ;
ð6Þ
where cof ðCÞ denotes the co-factor matrix of C, and M is the so-called structural tensor representing the transverseisotropy group and referring to a preferred direction a of the material: M ¼ a ^ a:
ð7Þ
Computer Methods in Biomechanics and Biomedical Engineering
3
The first strain invariant of C iso is defined by I^1 ¼ J 22=3 I 1 which replaces I 1 in WðI 1 Þ. ^ ¼ F^ T F^ ¼ J 22=3 C The modified invariants related to C are expressed from Equation (6) by a I^1 ¼ I 1 J 22=3 ; I^2 ¼ I 2 J 24=3 ; J^ 4 ¼ J a4 J 22=3 ; a J^ 5 ¼ J a5 J 24=3 :
ð13Þ
The exponential type HGO density adopted in this work uses these modified invariants as follows:
Downloaded by [Tulane University] at 12:56 29 September 2014
Figure 2.
Cylindrical coordinate system.
^ I^1 J^ a Þ þ W H ðJÞ; W ¼ Wð 4
ð14Þ
1 J2 2 1 2 lnJ ; W H ðJÞ ¼ d 2
ð15Þ
It is noted that Equations (6) and (7) give J 4 ¼ trðF Fa ^ aÞ ¼ kFak : 2
T
^ I^1 ; J^ a Þ ¼ W iso ðI^1 Þ þ Wð 4
ð8Þ
2 X
W ani ðJ^ 4 Þ; a
ð16Þ
a¼1
The double brackets represent the usual Euclidian norm. The square root of J 4 represents thus the stretch in the fiber direction. It can also be interpreted as the radial coordinate of Fa in a cylindrical coordinate system where the polar angle u represents the deformed angle between the collagen fibers and the circumferential direction (Figure 2). In the case of hyperelastic law, there exists an elastic potential function W (or strain energy density function) which is a scale function of one of the strain tensors, whose derivative with respect to a strain component determines the corresponding stress component. This can be expressed by S¼
›W ›W ¼2 ; ›E ›C
ð9Þ
where S is the second Piola– Kirchoff stress tensor. For numerical purpose, it proves useful to separate deformation in volumetric and isochoric parts by a multiplicative split of a deformation gradient as in (Flory 1961): F ¼ F iso F vol ;
ð10Þ
W aani
8 i k1 h k2 ðJ^a4 21Þ2 > < e 21 ¼ 2k2 > : 0
a if J^ 4 $ 1; a if J^ 4 , 1
:
ð17Þ
Generally, soft biological tissues are assumed to be incompressible. Equation (15), which was proposed by Horgan and Saccomandi (2003), represents a penalty term added to the finite element model to account for the incompressible behavior of the material. Here, d is the material incompressibility parameter. The initial bulk modulus K is defined as K ¼ 2=d. Some other models can be found in (Horgan and Saccomandi 2004, 2006). This a energy density is case-sensitive with respect to J^ 4 because a the case of J^ 4 , 1 represents the shortening of the fibers which is assumed to generate no stress. The proof of convexity of Equation (17) with respect to F is given in (Schro¨der et al. 2005; Balzani et al. 2006). The noncollagenous matrix of the media is modeled by the Gent isotropic density (Gent 1996; Feng et al. 2010): W iso ðI^1 Þ ¼ 2
I^1 2 3 m J m ln 1 2 ; Jm 2
ð18Þ
where F
vol
¼J
1=3
I;
F
iso
¼J
21=3
F:
ð11Þ
This decomposition is such that detðF iso Þ ¼ 1. It is easy to see that F and F iso have the same eigenvectors. The isochoric part of the right Cauchy – Green strain tensor C can then be defined as C iso ¼ J 22=3 C:
ð12Þ
where m is the shear modulus and J m is the constant limiting value for I^1 2 3. It is noted that for J m ! 1, we recover the well-known Neo-Hookean strain energy density: WðI 1 Þ ¼
m ðI 1 2 3Þ: 2
ð19Þ
It is noted that the volumetric –isochoric split of the above HGO model does only hold for (quasi)incompres-
4
Z.-W. Chen et al.
sible deformations. An extension to compressible deformations would require that the volumetric part of the strain energy function includes a dependency on the structural tensor. This is proved recently by Guo et al. (2008) where a simple compressible anisotropic analytical model is developed. We remind finally that by deriving W from Equation (9) and introducing the matrix of cofactors of C, Cof ðCÞ ¼ I 3 C 2T , it is conventionally obtained:
›W ›W ›W ›W 1 Iþ ðI 1 I 2 CÞ þ cof ðCÞ þ 1 M a ›I 1 ›I 2 ›I 3 ›J 4 ›W ›W ›W 2 1 1 2 2 þ 2 M a þ 1 ðCM a þ M a CÞ þ 2 ðCM a þ M a CÞ : ›J 4 ›J 5 ›J 5 ð20Þ
Downloaded by [Tulane University] at 12:56 29 September 2014
S¼2
In our particular case, Equation (20) is reduced to ›W ›W ›W a 1 ›W a 2 S¼2 Iþ cof ðCÞ þ 1 M þ 2 M : ð21Þ ›I 1 ›I 3 ›J 4 ›J 4 The derivatives of the density energy W with respect to the invariants are evaluated with equations (17,18):
›W 21=3 dW iso ¼ I3 ; ›I 1 d I^1
1 1 12 2 ; 2d J
ð22Þ
" # 2 ›W 1 ^ dW iso X a dW ani I1 J^ 4 ¼2 þ a 3I 3 ›I 3 d I^1 dJ^ 4 a¼1 þ
object-oriented programming language. The geometrically nonlinear analysis may be described by using the total or the updated Lagrangian formulations. The total Lagrangian formulation is derived with respect to the initial configuration. The updated Lagrangian formulation is derived with respect to the current configuration. In other words, the total Lagrangian formulation constructs the tangent stiffness matrix with respect to the initial configuration. On the other hand, the updated Lagrangian formulation constructs the tangent stiffness matrix with respect to the current configuration. The updated Lagrangian formulation is computationally effective (Zienkiewicz and Taylor 1991) because it does not include the initial displacement matrix. In the total Lagrangian formulation, the initial configuration remains constant. This simplifies the computation (Crisfield 1991). Therefore, the total Lagrangian formulation was selected in this work to describe nonlinear behavior. In the context of the finite element method and with Equations (3) and (5), the Green – Lagrangian strain can be formally written with linear and nonlinear terms in function of nodal displacements (Simo and Hughes 1998):
ð23Þ
E¼
1 B L þ B NL ðuÞ u; 2
where B L is the matrix which relates the linear strain term to the nodal displacements and B NL ðuÞ the matrix which relates the nonlinear strain term to the nodal displacements. From Equation (28), the incremental form of the strain – displacement relationship is
dE ¼ ðB L þ B NL ðuÞÞdu:
›W 21=3 dW ani ¼ I3 a ; ›J a4 dJ^ 4
ð24Þ
dW iso mJ m ; ¼ dI^1 2ðJ m 2 I^1 þ 3Þ
ð25Þ
ð28Þ
ð29Þ
Using the principle of virtual displacement, the virtual work equation is given as ð
dW ¼
SdEdV 2 F ext du;
ð30Þ
V0 a dW ani k2 ðJ^ 4 21Þ2 ^a : a ¼ k1 ðJ4 2 1Þe ^ dJ
ð26Þ
4
The Cauchy stress (or true stress) tensor s can be calculated from the second Piola – Kirchoff stress tensor S as follows: 1 s ¼ FSF T : J 2.2
ð27Þ
Finite element implementation
This model have been implemented and tested in the finite element code FER. This code is implemented using Cþ þ
where V 0 is the volume of the initial configuration and F ext the vector of external loads. Using Equations (9) and (29), the incremental form of the stress dS can be linked to the incremental form of the strain dE as follows:
dS ¼
›2 W ¼ D : dE ¼ D : ðB L þ B NL ðuÞÞdu; ›E 2
ð31Þ
where D denotes the constitutive tangent matrix. This fourth-order tensor is obtained from the derivative of S
Computer Methods in Biomechanics and Biomedical Engineering
5
90 Gent+HGO constitutive model HGO constitutive model Experiment
80 70 σ [kPa]
60 50 40 30 20 10 0
Downloaded by [Tulane University] at 12:56 29 September 2014
Figure 3.
þ cðM ^ MÞ; ð32Þ with 2mJ m J 24=3 ; ðJ m 2 I^1 þ 3Þ2
b¼
2mJ m J 22=3 ; 3ðJ m 2 I^1 þ 3Þ
c ¼ 4k1 ð1 þ 2k2 ðJ 4 2 1Þ Þe 2
k2 ðJ 4 21Þ2
ð33Þ
:
Substituting dE from Equation (29) into Equation (30) results in ð dW ¼ SðB L þ B NL ðuÞÞdudV 2 F ext du V0
¼ F int du 2 F ext du: The vector of internal forces is defined by ð SðB L þ B NL ðuÞÞdV: F int ¼
ð34Þ
ð35Þ
V0
Since du is arbitrary, a set of nonlinear equations can be obtained as F ext 2 F int ¼ 0
0.05
0.1
0.15 0.2 Δ l/l0
0.25
0.3
0.35
Figure 4. Cauchy stress sðkPaÞ versus strain Dl=l0 of the experiment and the constitutive model.
Unconstrained uniaxial tension test.
with respect to E in Equation (20): I1 I1 D ¼ a I 2 C 21 ^ I 2 C 21 ; 3 3 I1 ›C21 21 21 ; 2 b C ^ I þ ð1 þ ÞI ^ C 2 I 1 3 ›C 2 ›C21 J 2 C 21 ^ C 21 þ ðJ 2 2 1Þ þ d ›C
a¼
0
ð36Þ
This equation is strongly nonlinear, because of finite strains and large displacements of soft biological tissues.
A typical solution procedure for this type of nonlinear analysis is obtained by using the Newton – Raphson iterative procedure (Simo and Hughes 1998; Belytschko et al. 2000). It should be noted that all the models presented here have been implemented in the finite element code FER.
3.
Numerical results
To illustrate the results of the simulation using this model, we consider here two uniaxial tensile numerical tests. First, we consider a uniaxial unconstrained tension test as depicted in Figure 3. This loading is usually applied to soft biological tissues for identifying parameters of the material model. To fit this numerical model with experimental data in (Balzani et al. 2006), the material parameters for the best fit are chosen as m ¼ 17:41 kPa, J m ¼ 60:0 kPa, k1 ¼ 0:0017 kPa, k2 ¼ 882:847, and b ¼ 43:398. These last three values are proposed in (Balzani et al. 2006). Figure 4 illustrates the response of the constitutive model we proposed here (Gent þ HGO constitutive model), as well as the experiment data and response of the HGO constitutive model which are given by Balzani et al. (2006). As can be seen, the match is quite good, even better than HGO constitutive model using Neo-Hookean model, the exponential character can be presented better by using the Gent model. Whereas the Holzapfel– Gasser – Ogden’s constitutive law uses the Neo-Hookean model as the isotropic part, the Gent model has one more material parameter: J m , when it turns to infinity, we can get the Neo-Hookean strain energy density. So it is important to consider the influence of parameter J m . The other material parameters are the same as the last tension test except b ¼ 70:08 and several different values of J m were considered. Figure 5 illustrates the different
6
Z.-W. Chen et al. 24
1.6 Jm=3.0 Jm=5.0 Jm=6.0 Jm=7.0 Jm=10.0
σ [MPa]
16 12
Jm=10.0 Jm=20.0 Jm=50.0 Jm=100.0
1.2 σ [MPa]
20
0.8
8 0.4 4 0
0 0
0.4
0.8
1.2
1.6
2
0
0.4
0.8
Δ l/l0
1.6
2
10 ≤ Jm ≤ 100
Influence of J m on the mechanical response.
et al. (2010) for different values of the initial angle b between the collagen fibers orientation and the circumferential direction. This theoretical study was performed with the HGO constitutive law. A critical angle (bc < 54:738) has been analytically determined. It appears that the correspondence between the solution and J 4 is not one to one if b is greater than bc while uniqueness holds if b is lower than bc . A relationship between l2 and J 4 for different values of b is shown by Figure 7. Three different initial angles 408, 608, and 708 are considered in this test, the material parameters are chosen from (Holzapfel et al. 2000): m ¼ 6:0 kPa, J m ¼ 50:0 kPa, k1 ¼ 2:3632 kPa, and k2 ¼ 0:8393. The evolution of three principal stretches for different values of b is shown in Figure 8. As can be seen from Figure 7, if b , bc (e.g. b ¼ 208), during the whole tension process, there is no
mechanical responses due to the different J m . The influence of J m on the three principal stretches is shown in Figure 6. It seems that when J m , 6, the strength of material increases as J m grows, but the material has not been stretched to the length we expected. When J m . 6, the strength decreases with an increase in J m , especially when 6 , J m , 10, it decreases a lot while when J m . 50, and it stays stable and almost not influenced by J m . It is the same situation when it comes to deformation as shown in Figure 6 (right), when J m . 10, the deformation stays the same. We should notice that the influence of J m on the mechanical behavior of material models is also dependent on the other parameters such as m and b. To evaluate the dependence of the principal stretches on J 4 , the analytical solutions were studied by Peyraut
1.5
1.5
1
1 λ2, λ3
λ2, λ3
Downloaded by [Tulane University] at 12:56 29 September 2014
1 < Jm ≤ 10
Figure 5.
1.2 Δ l/l0
0.5
0.5 Jm=3.0, λ2 Jm=3.0, λ3 Jm=5.0, λ2
Jm=10.0, λ2 Jm=10.0, λ3 Jm=50.0, λ2
Jm=5.0, λ3 Jm=10.0, λ2 Jm=10.0, λ3
0
0 1
1.4
1.8
2.2
2.6
λ1 1 < Jm ≤ 10
Figure 6.
Jm=50.0, λ3 Jm=100.0, λ2 Jm=100.0, λ3
Influence of J m on the principal stretches.
3
1
1.4
1.8
2.2 λ1
10 ≤ Jm ≤ 100
2.6
3
Computer Methods in Biomechanics and Biomedical Engineering
7
1 0.8
λ2
0.6 0.4 β=20° β=60° β=70° β=80°
0.2 0 0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
J4
way that J 4 , 1, so J 4 is always $ 1, which means that the biological tissue is always stretched in the fiber direction and the anisotropic part works from the beginning, which explains why in Figure 8 (left) the two principle stretches l2 and l3 are not the same from the first place. The biological tissue can thus not be shortened in the fiber direction. On the contrary, if b . bc (e.g. b ¼ 708), J 4 can be less than 1 during some time, which means that the HGO model reduces to the isotropic Gent model. As can be seen in Figure 7, this time grows as the initial degree increases, so that it takes more time for b ¼ 708 to transfer into anisotropic behavior than b ¼ 608 as shown in Figure 8 (right). The biological tissue can be shortened in the fiber direction. At last, a numerical model of adventitial strip adopted in (Gasser et al. 2006) is considered to perform a uniaxial tensile test from the axial direction. The material parameters are m ¼ 15:28 kPa, J m ¼ 50:0 kPa, k1 ¼ 996:6 kPa, and k2 ¼ 524:6, and the structure parameter is
Figure 9. time.
Initial configurations and von-Mises stress at last
b ¼ 49:988. Adventitial strip model of length L ¼ 10.0 mm, width W ¼ 3.0 mm, and thickness T ¼ 0.5 mm is considered in this work. This model consists of 3200 hexahedral elements as shown in Figure 9 (left).
6
1.25 β=40°, λ2 β=40°, λ3
5
β=60°, λ2 β=60°, λ3 β=70°, λ2 β=70°, λ3
1
λ2, λ3
4 λ2, λ3
Downloaded by [Tulane University] at 12:56 29 September 2014
Figure 7. l2 versus J 4 for different values of b (Peyraut et al. 2010).
3
0.75
2 0.5 1 0
1
1.4
1.8
2.2
2.6
3
0.25
λ1
Figure 8.
Evolution of three principal stretches for different values of b.
1
1.4
1.8
2.2 λ1
2.6
3
8
Z.-W. Chen et al.
Downloaded by [Tulane University] at 12:56 29 September 2014
Figure 10.
Displacement in direction z.
Figure 9 (right) represents the deformed configuration and computed von-Mises stress of the specimen. According to the structural assumptions in (Holzapfel et al. 2000) and (Zulliger et al. 2004), the collagen is embedded as two families of fibers that are symmetrically disposed relative to the tensile (axial) direction. The embedded collagen fibers need to be rotated towards the direction of the loading until they are able to carry significant load. This leads to an increase in thickness in the middle of specimen as shown in Figure 9 (right). Due to the incompressibility character, the width of the specimen decreases. This numerical result is very similar to that observed in (Gasser et al. 2006). The displacement in the thickness direction (z) is shown in Figure 10. We can see clearly that the displacement in direction z is totally symmetrical.
4.
Conclusion
The development of constitutive laws is fundamental and a prerequisite for understanding the highly nonlinear deformation mechanisms and stress distributions of soft biological tissues under different loading conditions. In this work, the improved Gent þ HGO anisotropic hyperelastic law has been proposed and numerically implemented in the context of uniaxial tension tests. The exponential character of the stress – strain behavior of the considered tissue is modeled by the Gent þ HGO constitutive law and the curve for the uniaxial unconstrained tension test fits the experimental data better than the pure HGO model. An investigation is made on the effect of material parameter J m on the mechanical response and principal stretches. It turns out that when J m is small, the influence is obvious while with an increase in J m , the strain energy density gets more and more like Neo-Hookean model and the influence is paltry. Three different initial angles have been considered in this context; the result has a good agreement with the analytical solutions which was studied by Peyraut et al. (2010). A numerical simulation of adventitial strip has
been performed, in which the numerical result is very similar to that observed in (Gasser et al. 2006).
Funding This work was supported by a French ANR Grant [No. ANR-09SYSC-008] and the Natural Science Foundation of China [No. 11372260].
References Almeida ES, Spilker RL. 1998. Finite element formulation for hyperelastic transversely isotropic biphasic soft tissues. Comput Meth Appl Mech Eng. 151:513 – 538. Balzani D, Neff P, Schro¨der J, Holzapfel GA. 2006. A polyconvex framework for soft biological tissues. Adjustment to experimental data. Int J Solids Struct. 43:6052 – 6070. Belytschko T, Liu WK, Moran B. 2000. Nonlinear finite elements for continua and structures. Chichester: Wiley. Boehler JP. 1987. Introduction to the invariant formulation of anisotropic constitutive equations. In: Boehler JP, editor. Applications of tensor functions in solid mechanics, CISM Course No. 292. Vienna: Springer Verlag. Crisfield MA. 1991. Non-linear finite element analysis of solid and structures. Chichester: Wiley. Eberlein R, Holzapfel GA, Schulze-Bauer CAJ. 2001. An anisotropic model for annulus tissue and enhanced finite element analyses of intact lumbar disc bodies. Comput Methods Biomech Biomed Eng. 4:209– 230. Feng ZQ. 2008. A finite element research system. Available from: http://lmee.univ-evry.fr/, feng/FerSystem.html Feng ZQ, Renaud C, Cros JM, Zhang HW, Guan ZQ. 2010. A finite element finite-strain formulation for modeling colliding blocks of Gent materials. Int J Solids Struct. 47:2215– 2222. Flory PJ. 1961. Thermodynamic relations for high elastic materials. Trans Faraday Soc. 57:829– 838. Fung YC, Fronek K, Patitucci P. 1979. Pseudoelasticity of arteries and the choice of its mathematical expression. Am J Physiol. 237:H620 – H631. Gasser TC, Ogden RW, Holzapfel GA. 2006. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J R Soc Interface. 3:15 – 35. Gent AN. 1996. A new constitutive relation for rubber. Rubber Chem Technol. 69:59– 61.
Downloaded by [Tulane University] at 12:56 29 September 2014
Computer Methods in Biomechanics and Biomedical Engineering Guo ZY, Caner F, Peng XQ, Moran B. 2008. On constitutive modelling of porous neo-Hookean composites. J Mech Phys Solids. 56:2338– 2357. Holzapfel GA, Gasser TC. 2001. A viscoelastic model for fiberreinforced composites at finite strains: continuum basis, computational aspects and applications. Comput Meth Appl Mech Eng. 190:4379– 4403. Holzapfel GA, Gasser TC, Ogden RW. 2000. A new constitutive framework for arterial wall mechanics and a comparative study of material models. J Elast. 61:1 – 48. Holzapfel GA, Gasser TC, Ogden RW. 2004. Comparison of a multilayer structural model for arterial walls with a Fungtype model, and issues of material stability. J Biomech Eng. 126:264 – 275. Holzapfel GA, Ogden RW. 2006. Mechanics of biological tissue. Berlin: Springer. Horgan CO, Saccomandi G. 2003. Finite thermoelasticity with limiting chain extensibility. J Mech Phys Solids. 51:1127 – 1146. Horgan CO, Saccomandi G. 2004. Constitutive models for compressible nonlinearly elastic materials with limiting chain extensibility. J Elast. 77:123– 138. Horgan CO, Saccomandi G. 2006. Phenomenological hyperelastic strain-stiffening constitutive models for rubber. Rubber Chem Technol. 79:152– 169. Ogden RW. 1984. Non-linear elastic deformations. Chichester: Ellis Horwood. Pen˜a E, del Palomar AP, Calvo B, Martı´nez MA, Doblare´ M. 2007. Computational modelling of diarthrodial joints. Physiol, pathological and post-surgery simulations. Arch Comput Meth Eng. 14(1):1 – 54. Pen˜a E, Martı´nez MA, Calvo B, Doblare´ M. 2006. On the numerical treatment of initial strains in soft biological tissues. Int J Numer Meth Eng. 68:836 – 860. Peng XQ, Guo ZY, Moran B. 2006. An anisotropic hyperelastic constitutive model with fiber-matrix shear interaction for the human annulus fibrosus. J Appl Mech. 73:815 – 824.
9
Peyraut F, Feng ZQ, Labed N, Renaud C. 2010. A closed form solution for the uniaxial tension test of biological soft tissues. Int J Nonlinear Mech. 45:535 – 541. Ru¨ter M, Stein E. 2000. Analysis, finite element computation and error estimation in transversely isotropic nearly incompressible finite elasticity. Comput Meth Appl Mech Eng. 190:519 – 541. Schro¨der J, Neff P, Balzani D. 2005. A variational approach for materially stable anisotropic hyperelasticity. Int J Solids Struct. 42:4352 –4371. Simo JC, Hughes TJR. 1998. Computational inelasticity. New York: Springer-Verlag. Spencer AJM. 1954. Theory of invariants. Continuum Physics. New York: Academic Press; p. 239– 253. Spencer AJM. 1972. Deformations of fiber-reinforced materials. New York: Oxford University Press. Spencer AJM. 1984. Continuum theory of the mechanics of fiberreinforced composites. Vienna: Springer-Verlag. Spencer AJM. 1987. Isotropic polynomial invariants and tensor functions. In: Boehler JP, editor. Applications of tensor functions in solids mechanics, CISM Course No. 282. Vienna: Springer Verlag. Weiss JA, Maker BN, Govindjee S. 1996. Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput Meth Appl Mech Eng. 135:107 – 128. Zhang JM, Rychlewski J. 1990. Structural tensors for anisotropic solids. Arch Mech. 42:267 – 277. Zheng QS, Spencer AJM. 1993. Tensors which characterize anisotropies. Int J Eng. Sci. 31(5):679– 693. Zienkiewicz OC, Taylor RL. 1991. The finite element method. Maidenhead Berkshire: McGraw-Hill. Zulliger MA, Fridez P, Hayashi K, Stergiopulos N. 2004. A strain energy function for arteries accounting for wall composition and structure. J Biomech. 37:989 – 1000.
