- Email: [email protected]
A closed form solution for the uniaxial tension test of biological soft tissues
ARTICLE IN PRESS International Journal of Non-Linear Mechanics 45 (2010) 535–541
Contents lists available at ScienceDirect
International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm
A closed form solution for the uniaxial tension test of biological soft tissues F. Peyraut a,, Z.-Q. Feng b, N. Labed a, C. Renaud b a b
Laboratoire M3M, Universite´ de Technologie de Belfort-Montbe´liard, 90010 Belfort, France ´ nerge´tique d’E ´ vry, Universite´ d’E ´vry-Val d’Essonne, 40 rue du Pelvoux, 91020 E ´vry, France ´canique et d’E Laboratoire de Me
a r t i c l e in fo
abstract
Article history: Received 25 August 2009 Received in revised form 11 February 2010 Accepted 19 February 2010
This paper concerns the modeling of biological soft tissues in the framework of anisotropic hyperelasticity. A closed form solution is proposed in the special case where uniaxial unconstrained tension loading is applied to a strip modeled by Holzapfel–Gasser–Ogden’s (HGO) strain energy. The classical Cardano’s formula is used to calculate the solution. This solution is investigated for different values of b which represents the angle between the collagen fibers and the circumferential direction. This study allows for the understanding of the reason why a one-to-one correspondence does not exist between the principal stretch and the fourth invariant of the right Cauchy–Green deformation tensor. It is demonstrated that the relationship becomes non-bijective when b is greater than a critical angle of 54:733 . The HGO model is implemented into an in-house finite element program and numerical results are in good agreement with theoretical ones. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Anisotropic hyperelasticity Biological soft tissues HGO model Theory and modeling
1. Introduction To investigate internal deformation and stress of biological soft tissues such as ligaments, tendons or arterial walls, anisotropic hyperelastic constitutive laws are often used in the framework of finite element analysis [1–3]. The most used strain-energy functions take a power law form [4] or present an exponential behavior [5,6]. More recently, Balzani et al. [7] 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 been also recently developed to capture the inter-fiber angle change by adding to the strain energy the contribution of the fiber-matrix shear interaction [8]. 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 [9–11]. It is usually assumed that anisotropy is due to the collagen fibers behavior [12], 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 [1,7]. Each anisotropic density refers to a preferred direction of the material. For example, in Holzapfel–Gasser–Ogden’s constitutive law [6,13], two transversely isotropic energies with two distinct preferred directions corresponding to two superposed collagen fiber families.
Corresponding author. Tel.: + 33 3 84 58 31 91; fax: + 33 3 84 58 31 46.
E-mail address: [email protected] (F. Peyraut). 0020-7462/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijnonlinmec.2010.02.003
In this paper, a theoretical study related to a uniaxial unconstrained tension loading providing homogeneous deformation is considered with Holzapfel–Gasser–Ogden’s constitutive law. The nonlinear equilibrium equations are reduced into a cubic polynomial form by introducing an appropriate cylindrical coordinate system. The polynomial equation is solved thanks to Cardano’s formula and the analytical solution is successfully compared to numerical one obtained thanks to the university finite element program FER [14]. The results are investigated for different values of the material parameter b representing the non-deformed initial angle between the collagen fibers and the circumferential direction.
2. HGO hyperelastic model Most energy densities used to model transversely isotropic and orthotropic soft tissues take a power law form [4] or present an exponential behavior [5,6]. It is usually assumed that anisotropy is due to the collagen fibers behavior [12] while the ground substance, or matrix, behaves in an isotropic manner. The energy densities modeling transversely isotropic and orthotropic soft tissues are thus separated into isotropic and anisotropic parts [1] W ¼ Wiso þ
n X
a Wani
ð1Þ
a¼1
Each anisotropic density Waani 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,
ARTICLE IN PRESS 536
F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
a2
β β
a1
Fig. 1. Angle b for circumferentially oriented strip.
Holzapfel–Gasser–Ogden’s constitutive law [6,12] superposes two transversely isotropic energies with two distinct preferred directions a1 and a2 corresponding to two fiber families: 8 9 8 9 > > = < c > = ð2Þ a1 ¼ s , a2 ¼ s with c ¼ cosðbÞ, s ¼ sinðbÞ > > > : > ; : ; 0 0
Fig. 2. Cylindrical coordinate system.
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 (Fig. 1). b is predicted to be 433 39 in the work of Balzani et al. [7]. In our work, we consider b as a model parameter varying from 01 to 901. The strain measure adopted here is the Green–Lagrangian strain tensor 1 E ¼ ðCIÞ 2
a2
t
ð3Þ
@x @u ¼ Iþ , @X @X
J ¼ detðFÞ 4 0
I2 ¼ trðcof ðCÞÞ,
t Fig. 3. Unconstrained uniaxial tension test.
I3 ¼ detðCÞ,
J4 ¼ trðCMÞ,
J5 ¼ trðC2 MÞ ð5Þ
where cof(C) denotes the co-factor matrix of C and M is the so-called structural tensor representing the transverse-isotropy group and referring to a preferred direction a of the material M¼aa
ð6Þ
It is noted that (5) and (6) give J4 ¼ trðFT Fa aÞ ¼ JFaJ2
t
ð4Þ
X,x and u represent respectively the reference and the current positions and the displacement vector of a material point. According to Zhang–Rychlewski’s theorem [15], 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 I1, I2 and I3 of the right Cauchy–Green deformation tensor and two additional mixed invariants J4 and J5 [9–11] I1 ¼ trðCÞ,
x2
x1
where I is the second order unit tensor, C¼F F is the right Cauchy–Green deformation tensor and F the transformation gradient defined by
a1
t
x3
T
F¼
stress tensor r are given by S¼2
@W pC1 , @C
r ¼ J1 FSFT
ð8Þ
To uncouple the deviatoric part to the dilatational part of the response, the volume preserving part F^ ¼ J1=3 F of the deformation is introduced [1]. The modified invariants related to T C^ ¼ F^ F^ ¼ J 2=3 C are expressed from (5) by I^ 1 ¼ I1 J 2=3 ,
I^2 ¼ I2 J4=3 ,
a J^ 4 ¼ J4a J 2=3 ,
a J^ 5 ¼ J5a J 4=3
ð9Þ
The exponential type HGO density adopted in this work uses these modified invariants as follows: 8 a a k 2 > < 1 ½ek2 ðJ^ 4 1Þ 1 if J^ 4 Z 1 a ð10Þ Wani ¼ 2k2 > a :0 if J^ o 1 4
ð7Þ
The double brackets represent the usual Euclidian norm. The square root of J4 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 y represents the deformed angle between the collagen fibers and the circumferential direction (Fig. 2). In what follows, J4 and y will be used to solve the equilibrium equations in an analytical manner. In the case of hyperelastic materials, there exists an elastic potential function W which is a scalar function of the strain tensors. Generally, soft biological tissues are assumed to be incompressible. The incompressibility is characterized in terms of principal stretches l1 , l2 and l3 by J ¼ detðFÞ ¼ l1 l2 l3 ¼ 1 related to an extra pressure p. The second Piola–Kirchhoff stress tensor S and the corresponding Cauchy
This energy density is case sensitive with respect to Ja4 because the case of J4a o 1 represents the shortening of the fibers which is assumed to generate no stress. The proof of convexity of (10) with respect to F is given in [4,7]. The non-collagenous matrix of the media is modeled by the neo-Hookean isotropic density Wiso ¼ c1 ðI^1 3Þ
ð11Þ
It is noted that the volumetric-isochoric split of the above HGO model does only hold for (quasi) incompressible 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. [16] where a simple compressible anisotropic analytical model is developed. In order to conduct analytical solutions, we consider a uniaxial unconstrained tension test as depicted in Fig. 3.
ARTICLE IN PRESS F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
This loading is usually applied to biological soft tissues for identifying parameters of the material model. In this test, the parameters are set to c1 ¼10.2069 kPa, k1 ¼0.0017 kPa and k2 ¼882.847. These values are proposed in [7] to fit numerical models with experimental data. As the test leads to homogeneous deformations, the equilibrium equations are automatically satisfied and the boundary conditions related to the free traction surfaces are
s22 ¼ s33 ¼ 0
ð12Þ
The principal stretches l1 , l2 and l3 will be calculated in next sections by using an analytical approach and compared to finite element results. The incompressibility constraint is included in the finite element model by using a penalty term of the form k(J 1)2/2 with k¼ 105.
537
direction (Fig. 2) 1=2
cosðyÞ ¼ l1 cJ 4
,
1=2
sinðyÞ ¼ l2 sJ 4
ð19Þ
If (19) is reported in (17), the principal stretches l1 and l2 are replaced by two new variables y and J4. Eq. (17) becomes then a cubic polynomial equation: 3=2
s2 cJ 4 ¼A cosðyÞcos3 ðyÞ ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2k1 s2 2 1þ ðJ4 1Þek2 ðJ4 1Þ c1
ð20Þ
To calculate the roots of (20) with Cardano’s formula, three different cases have to be discussed with respect to the sign of the discriminant D
D ¼ A2 =41=27
ð21Þ
Case a: 3. Theoretical study
2
D 4 0 3 A 4 pffiffiffi
ð22Þ
3 3
By using (1), (8)–(11), it is easily demonstrated that 0 B J4 Z 1 : r ¼ p@
80 2 l > > > : 1
1
1
1 1
þ 4k1 ðJ^ 4 1Þe
k2 ðJ^ 4 1Þ2
0
C I^1 B C @ A 3 2
l22
1 1
l2 1 l2
80 2 2 > > > :
1
l22 s2 0
0 1 C J^ 4 B C @ A 3
1
19 > > = C A > > 1 ;
19 > > = C A > > 1 ;
ð13Þ 0 B J4 o 1 : r ¼ p@
80 2 l > > > : 1
1
1
1 1
l22
0 1 C I^1 B C @ A 3 2
l2 1 l2
1
19 > > = C A > > 1 ;
ð14Þ By considering s33 ¼ 0 to express the pressure p in terms of the principal stretch and taking into account the incompressibility condition J¼1, we obtain 8 < s11 ¼ 2½c1 ðl2 l2 l2 Þ þ 2k1 l2 c2 ðJ4 1Þek2 ðJ4 1Þ2 1 1 2 1 ð15Þ J4 Z 1 : 2 2 2 k2 ðJ4 1Þ2 : s22 ¼ 2½c1 ðl22 l2 l Þ þ 2k l ¼0 1 2 s ðJ4 1Þe 1 2 2
2 2
J4 o 1 : s11 ¼ 2c1 ðl1 l1 l2 Þ,
2 s22 ¼ 2c1 ðl22 l2 1 l2 Þ ¼ 0
ð16Þ
To determine the principal stretches in an analytical manner, we will demonstrate that the equilibrium equations can be transformed into a third degree polynomial form. The polynomial equation can then be algebraically solved by using Cardano’s formula. As the stresses expressions (15) and (16) depend on the value of J4, we logically start by considering two different cases.
It results from the definition (20) of A A rs2 c
ð23Þ
We consider the function of y representing the left-hand side of (20) f ðyÞ ¼ cosðyÞcos3 ðyÞ ¼ cosðyÞsin2 ðyÞ,
y A ½0,903
ð24Þ
It ispeasy ffiffiffi to demonstrate that the maximum value of f is equal to 2=3 3 and is reached if y ¼ bc with (Fig. 4) pffiffiffi bc ¼ Arccosð 3=3Þ 54:733 ð25Þ It then results from (23) that pffiffiffi A r2=3 3
ð26Þ
This result is contradictory with (22). Consequently, Case a will never occur. Case b: 2
D ¼ 0 3 A ¼ pffiffiffi
ð27Þ
3 3
In this case, Cardano’s formula gives one single and one double real root pffiffiffi pffiffiffi cosðy1 Þ ¼ 2 3=3, cosðy2 Þ ¼ 3=3 ð28Þ According to (24), the roots must be positive. The single root cosðy1 Þ is thus not acceptable. However, the double root cosðy2 Þ is acceptable as it is included in the range [0, 1]. The corresponding
0.4
2/3√3
0.3
The free traction surfaces condition s22 ¼ 0 applied to (15) gives 2
2 2
2
2
c1 ðl2 l1 l2 Þ þ2k1 l2 s2 ðJ4 1Þek2 ðJ4 1Þ ¼ 0
¼1
0.2
ð17Þ
0.1
Besides, it results from (2) and (5) that 1=2 1=2 ðl1 cJ 4 Þ2 þ ðl2 sJ 4 Þ2
f (θ)
3.1. Case J4 Z1
ð18Þ
This suggests that a cylindrical coordinate system could be appropriate to solve the non-linear equation (17). This cylindrical coordinate system is defined by setting J4½ as the radial coordinate and by introducing the polar angle y representing the deformed angle between the collagen fibers and the circumferential
0
βc 0
10
20
30
40 50 θ (°)
60
Fig. 4. Left-hand side of (20).
70
80
90
ARTICLE IN PRESS 538
F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
stretches are deduced from (19) pffiffiffi pffiffiffi pffiffiffi l1 ¼ J41=2 =ð 3cÞ, l2 ¼ 2J41=2 =ð 3sÞ
ð29Þ
As the maximum value of f is reached if y ¼ bc (Fig. 4), (27) gives A ¼ f ðbc Þ Z f ðbÞ
single root defined by cosðy1 Þ. If D is negative and J4 ¼1, there exists two distinct single roots defined by cosðy1 Þ and cosðy3 Þ. 3.2. Case J4 o1
ð30Þ
Besides, it results from the definitions (20) and (24) of A and f that If J4 41 : A o f ðbÞ
If J4 ¼ 1 : A ¼ f ðbÞ
ð31Þ
The conditions (30) and (31) are together consistent if and only if b ¼ bc and J4 ¼1. It is then deduced from (29) that Case b gives thus a double root corresponding to the evident solution l1 ¼ l2 ¼ 1. Case c: 2
D o0 3 Ao pffiffiffi
ð32Þ
3 3
In this case, according to Cardano’s formula, there exist three distinct real roots " pffiffiffi # pffiffiffi 2 3 a þ2ði1Þp 3 3A cos cosðyi Þ ¼ , a ¼ Arccos , i ¼ 1,2,3 3 3 2 ð33Þ The locations on the trigonometric circle of the three angles used in (33) are depicted in Fig. 5. To determine the acceptable roots, we first notice that (20) implies that A is positive. It then results from (32) and (33) that pffiffiffi 3 ocosðy1 Þ o1, cosðy2 Þ o 0 0 o cosðy3 Þ o ð34Þ 3 It is noted that the second root cosðy2 Þ is negative and must be excluded. Among the two positive roots cosðy1 Þ and cosðy3 Þ, we propose to demonstrate that the third is unacceptable except if J4 is equal to 1. To prove it, we remind that a tensile loading test is considered. It then results from (19) and (33) that pffiffiffi a þ4p ð35Þ ðc 3Þ Z1 l1 ¼ 2J41=2 cos 3 After some algebra, it can be demonstrated that (35) leads to J4 r 1. As the case of J4 Z1 is considered, we have J4 ¼1. The root cosðy3 Þ is thus excluded except if J4 ¼1. In conclusion, if D is negative and J4 is strictly greater than 1, there exists only one
This case is easier to study than the previous one because the orthotropic behavior is not taken into account and the model is reduced to the neo-Hookean energy density according to Eq. (10). By using the stress expression (16) in conjunction with the free traction surface loading s22 ¼ 0, the neo-Hookean solution is directly obtained pffiffiffiffiffi ð36Þ l2 ¼ 1= l1 To study the dependence of l1 and l2 on J4 in the same way as for the case J4 Z 1, it is needed to solve again a cubic polynomial equation. This equation is deduced from (18) and (36) 3
gðl1 Þ ¼ l1 J4 c2 l1 þ t 2 ¼ 0
with t ¼ s=c
ð37Þ
This function admits a local maximum and a local minimum as shown in Fig. 6. The local maximum and minimum points are qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi defined respectively by ð J4 =3c2 ,t 2 þ 4J43 =27c6 Þ and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð J4 =3c2 ,t 2 4J43 =27c6 Þ. To determine the real roots of g, Eq. (37) is solved analytically by using again Cardano’s formula. Three different cases occur depending on the sign of the discriminant D1
D1 ¼ t4 =4J43 =27c6
ð38Þ
Case a:
D1 40 3 t2 4
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4J43 =ð27c6 Þ
ð39Þ
In this case, Cardano’s formula indicates that one root is a real number while the two others are conjugate complex numbers. Besides, it is reminded that the condition l1 Z 1 must be satisfied because a tensile load is considered. As g increases from 1 to qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t2 þ 4J43 =ð27c6 Þ 40 for l1 varying from 1 to J4 =ð3c2 Þ, the real root is negative and should not be selected. Finally, using definitions (24) and (25), the condition (39) can be simplified to J4 o½f ðbÞ=f ðbc Þ2=3
ð40Þ
Case b:
sin
D1 ¼ 0 3 t2 ¼
α/3
ð41Þ
2.4
1.2
π/6
(α+2π)/3
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4J43 =ð27c6 Þ
cos
π/6
g
π/6 0
−1.2
−2.4 −1.8 (α+4π)/3 Fig. 5. Location of a=3, ða þ 2pÞ=3 and ða þ 4pÞ=3 on the trigonometric circle.
−1.2
−0.6
0 λ1
0.6
Fig. 6. Graph of g ðb ¼ 503 ,J4 ¼ 0:5Þ.
1.2
1.8
ARTICLE IN PRESS F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
In this case, Cardano’s formula gives a single and a double real root Single root : l1 ¼ 3s2 =J4 ,
Double root : l1 ¼ 3s2 =2J4
ð42Þ
The single root is negative and has to be excluded. By using (41), it is noted that the double root corresponds to the minimum of the g function. To make this root acceptable, the condition l1 Z1 must hold such that pffiffiffi pffiffiffi ðt4 =4Þ1=3 Z1 3 t Z 2 3 b Z bc ¼ Arccosð 3=3Þ 54:733 ð43Þ The double root defined by (42) is thus solution of Eq. (37) if and only if the phenomenological angle b is greater than the critical angle bc defined by (43). It is worth noting that this critical angle has also been found by Guo et al. with a composites-based hyperelastic constitutive model [17]. Moreover, it results from (19), (36), (41), (42) that cosy ¼ cosbc ,
siny ¼ sinbc ) y ¼ bc
ð45Þ
Case c:
D1 o 0 3 t2 o
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4J43 =ð27c6 Þ
ð46Þ
In this case, Cardano’s formula provides three distinct real roots expressed trigonometrically rffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffi g J4 J4 g þ 2p cos cos , , l1a ¼ 2 l ¼ 2 1b 2 2 3 3 3c 3c rffiffiffiffiffiffiffiffi J4 g þ 4p cos ð47Þ l1c ¼ 2 2 3 3c qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
g ¼ Arccos½t2 27c6 =ð4J43 Þ
In the same way, it is concluded that l1c o1. Consequently, the two positive roots are unacceptable. Case c2: pffiffiffi b Z bc ¼ Arccosð 3=3Þ ð53Þ As the tangent function increases in the interval ½0, p=2, (46) and (53) yield to qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð54Þ t 2 Z 2 ) J4 =ð3c2 Þ 4 1 To conclude that l1c 41, we consider Eq. (54) and we apply again the intermediate value theorem in the same way as in Case c1. It results next from (49) that l1a 4 1. Consequently, the two positive roots l1a and l1c defined by (47) are acceptable. Finally, by reasoning in the same manner as for Case a, the condition (46) can be simplified to J4 4½f ðbÞ=f ðbc Þ2=3
ð55Þ
ð44Þ
It is thus remarkable that the deformed angle y is equal to the critical angle bc in the case of a double root occurrence corresponding to an initial angle b greater than bc . Finally, it is noted that condition (41) can be simplified in the same way as for Case a: J4 ¼ ½f ðbÞ=f ðbc Þ2=3
539
ð48Þ
To determine the acceptable roots, we first deduce from (46)–(48) that qffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 o l1c o J4 =ð3c2 Þ o l1a o J4 =c2 , 2 J4 =ð3c2 Þ o l1b o J4 =c2 ð49Þ The second root l1b is negative and cannot be selected. To determine whether or not the positive roots l1a and l1c are greater than one, according to the condition of a tension test, two cases are considered Case c1: pffiffiffi b o bc ¼ Arccosð 3=3Þ ð50Þ Because the cosine function decreases in the interval ½0, p=2, (50) yields to qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð51Þ 3c2 41 ) 0 o J4 =ð3c2 Þ o1 Eq. (51) means that the minimum of g is reached between 0 and 1. We then apply the intermediate value theorem to g and we account for condition (46) sffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffi! 4J43 J4 1J4 2 ¼t g o0, gð1Þ ¼ 2 4 0 2 3c 27c6 c rffiffiffiffiffiffiffiffi !, J4 ) (l1a A ð52Þ ,1 gðl1a Þ ¼ 0 3c2
4. Analytical and numerical results The analytical formula applied to the uniaxial tension test determines the evolution of l2 and y vs. J4 as shown in Fig. 7. If b o bc (e.g., b ¼ 203 ), there is no solution corresponding to J4 o 1. The biological tissue can thus not be shortened in the fiber direction. On the contrary, if b 4 bc (e.g., b ¼ 703 ), it is observed that the biological tissue can be shortened in the fiber direction. However, the correspondence between l2 and J4 is not one-to-one because of the existence of two distinct solutions. These two solutions are given by (47) and (48). It is finally noted that the curve related to b ¼ 703 starts with J4 equal to a minimum value defined by (45). This minimum value corresponds to the double root defined by (42). Eq. (45) shows that the minimum value of J4 only depends on the initial angle b. As the function f decreases if b 4 bc (Fig. 4), this minimum value decreases if b increases. This trend is well observed in Fig. 7 (left). The minimum values of J4 are respectively equal to 0.98, 0.85 and 0.58 for b ¼ 603 ,703 and 803 . The occurrence of a double root with J4 o 1 corresponds to a remarkable result given by (44), i.e., the deformed angle y is equal to the critical angle bc when the minimum value of J4 given by (45) is reached as shown in Fig. 7 (right). Fig. 8 shows the variation of the cylindrical coordinates (J4, y) vs. l1 for two values of b (203 and 703 ). As a tension test is considered in this work, l1 is greater than 1. In the case of b ¼ 703 , it is observed that J4 is first lower than 1 (the biological tissue is first shortened in the fiber direction) and decreases to a minimum value. This case has been previously discussed in Section 3.2 (Case c2) and can only occur if b Z bc . The minimum value of J4 is equal to 0.85 according to (45). This minimum is reached if l1 is approximately equal to 1.56 (Fig. 8 a) and corresponds to a deformed angle y equal to the critical angle bc (Fig. 8 b). After this minimum value, J4 increases and becomes greater than 1 if l1 is approximately greater than 2.29. This value indicates the transition from a purely isotropic to an anisotropic behavior of biological pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tissues. It can be analytically calculated by using l1 ¼ ð 1þ 4t 2 1Þ=2. This formula is easily obtained by solving (37) together with J4 ¼1. It is worth noting that if b o bc (e.g., b ¼ 203 Þ,J4 increases with l1 (Fig. 8 a) so that the biological tissue is always stretched in the fiber direction ðJ4 Z1Þ. The critical angle bc plays obviously a key role for the understanding of the model behavior. However, it should be noted that the minimum value of J4 and the correspondence between the solution and J4 were only discussed in the case of J4 o1. As the orthotropic part of the HGO energy density is assumed to generate no stress if J4 o 1, the critical angle concerns only the neo-Hookean part of the model. The results related to bc given in this paper could thus be extended to any anisotropic
ARTICLE IN PRESS 540
F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
90
1
β=60° β=70° β=80°
80 0.8
70 60
0.6
βc=54.73°
θ
λ2
50
0.4 0.2 0 0.5
0.6
40 30
β=20° β=60° β=70° β=80°
20 10 0.7
0.8
0.9 J4
1
1.1
1.2
0 0.5
1.3
0.6
0.7
0.8
0.9 J4
1
1.1
1.2
1.3
Fig. 7. l2 and y vs. J4 for different values of b.
70
1.3
60
1.2 1.1 θ
J4
40
1
1
1.5
2
2.5
30 20
β=20° β=70°
0.9 0.8
β=20° β=70°
50
10 3
0
3.5
1
1.5
2
λ1
2.5
3
3.5
λ1
Fig. 8. Evolution of the cylindrical coordinates (J4, y) for two values of b.
1.4 1.2
λ2, λ3
1
1
1
0.8 3
0.6
3
Fig. 10. Finite element modeling: initial and deformed shapes (left: l1 o 2:29, right: l1 4 2:29).
0.4 0.2 0
2
2
λ2−analytical λ3−analytical λ2−numerical λ3−numerical
isotropic neo-Hookean model and the relationship between l1 and l2 is defined by the neo-Hookean solution (36).
1
1.5
2 λ1
2.5
3
Fig. 9. Comparison between numerical and analytical results ðb ¼ 703 Þ.
model superposing any orthotropic density with the isotropic neo-Hookean density. The HGO model has been implemented into an in-house finite element code FER [14]. Fig. 9 shows a good agreement between the numerical results obtained with FER and the analytical calculation presented in this work with b ¼ 703 . Beyond a threshold value of l1 2:29, it is observed that l3 increases while l2 decreases. In other words, the section shortens in the direction 2 and swells in the direction 3 as shown in Fig. 10. As outlined in [12], this trend typically represents an anisotropic behavior. Before this transition, the HGO model reduces to the
5. Conclusions In this paper, a theoretical study has been performed with Holzapfel–Gasser–Ogden’s (HGO) constitutive law and applied to a unconstrained tension test. The free traction loading condition related to this test is expressed as a non-linear algebraic equation to give a relationship between the principal stretches l1 and l2 . To change this equation into a cubic polynomial form, we propose to use a cylindrical coordinate system by setting the square root of the mixed invariant J4 as the radial coordinate. The deformed angle y, between the collagen fibers and the circumferential direction, is used as the polar angle. A closed form solution has been obtained thanks to Cardano’s formula. A critical angle ðbc 54:733 Þ has been analytically determined.
ARTICLE IN PRESS F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
To evaluate the dependence of the principal stretches on J4, the analytical solutions were studied for different values of the initial angle b between the collagen fibers orientation and the circumferential direction. It appears that the correspondence between the solution and J4 is not one to one if b is greater than bc while uniqueness holds if b is lower than bc . We have also established that the biological tissue cannot be shortened in the fiber direction more than a minimum value of J4 given by an analytical formula. We have finally demonstrated that this minimum value corresponds to a double root. In this case, it is remarkable that the deformed angle y is equal to the critical angle bc . However, it should be noted that the minimum value of J4 and the correspondence between the solution and J4 were only discussed in the case of J4 o1. As the orthotropic part of the HGO energy density is assumed to generate no stress if J4 o 1, the critical angle bc only concerns the neo-Hookean part of the model. The results related to bc given in this paper could thus be extended to other anisotropic hyperelasticity models superimposing any orthotropic density with the isotropic neo-Hookean density. The closed form solutions are also used to validate the implementation of the HGO model into the finite element code FER. A fair agreement has been found between theoretical and numerical solutions. If b is greater than bc , the results underline the transition between purely isotropic and orthotropic behavior. In particular, it is observed that the strip section shortens in one direction and swells in another direction beyond a threshold value of l1 corresponding to J4 ¼1. References [1] J.A. Weiss, B.N. Maker, S. Govindjee, Finite element implementation of incompressible, transversely isotropic hyperelasticity, Comput. Meth. Appl. Mech. Eng. 135 (1996) 107–128.
541
[2] E.S. Almeida, R.L. Spilker, Finite element formulation for hyperelastic transversely isotropic biphasic soft tissues, Comput. Meth. Appl. Mech. Eng. 151 (1998) 513–538. ¨ [3] M. Ruter, E. Stein, Analysis, finite element computation and error estimation in transversely isotropic nearly incompressible finite elasticity, Comput. Meth. Appl. Mech. Eng. 190 (2000) 519–541. ¨ [4] J. Schroder, P. Neff, D. Balzani, A variational approach for materially stable anisotropic hyperelasticity, Int. J. Solids Struct. 42 (2005) 4352–4371. [5] Y.C. Fung, K. Fronek, P. Patitucci, Pseudoelasticity of arteries and the choice of its mathematical expression, Am. J. Physiol. 237 (1979) H620–H631. [6] G.A. Holzapfel, T.C. Gasser, R.W. Ogden, A new constitutive framework for arterial wall mechanics and a comparative study of material models, J. Elasticity 61 (2000) 1–48. ¨ [7] D. Balzani, P. Neff, J. Schroder, G.A. Holzapfel, A polyconvex framework for soft biological tissues. Adjustment to experimental data, Int. J. Solids Struct. 43 (2006) 6052–6070. [8] X.Q. Peng, Z.Y. Guo, B. Moran, An anisotropic hyperelastic constitutive model with fiber-matrix shear interaction for the human annulus fibrosus, J. Appl. Mech. Trans. ASME 73 (5) (2006) 815–824. [9] A.J.M. Spencer, Isotropic polynomial invariants and tensor functions, in: J.P. Boehler (Ed.), Applications of Tensor Functions in Solids Mechanics, CISM Course No. 292, Springer Verlag, Wien, 1987. [10] J.P. Boehler, Introduction to the invariant formulation of anisotropic constitutive equations, in: J.P. Boehler (Ed.), Applications of Tensor Functions in Solids Mechanics, CISM Course No. 292, Springer Verlag, Wien, 1987. [11] Q.S. Zheng, A.J.M. Spencer, Tensors which characterize anisotropies, Int. J. Eng. Sci. 31 (5) (1993) 679–693. [12] T.C. Gasser, R.W. Ogden, G.A. Holzapfel, Hyperelastic modelling of arterial layers with distributed collagen fibre orientations, J. R. Soc. Interface 3 (2006) 15–35. [13] G.A. Holzapfel, T.C. Gasser, A viscoelastic model for fiber-reinforced composites at finite strains: continuum basis, computational aspects and applications, Comput. Meth. Appl. Mech. Eng. 190 (2001) 4379–4403. [14] Z.Q. Feng /http://lmee. univ-evry.fr/ feng/FerSystem.htmlS, 2008. [15] J.M. Zhang, J. Rychlewski, Structural tensors for anisotropic solids, Arch. Mech. 42 (1990) 267–277. [16] Z.Y. Guo, F.C. Caner, X.Q. Peng, B. Moran, On constitutive modelling of porous neo-Hookean composites, J. Mech. Phys. Solids 56 (6) (2008) 2338–2357. [17] Z.Y. Guo, X.Q. Peng, B. Moran, A composites-based hyperelastic constitutive model for soft tissue with application to the human annulus fibrosus, J. Mech. Phys. Solids 54 (9) (2006) 1952–1971.
Contents lists available at ScienceDirect
International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm
A closed form solution for the uniaxial tension test of biological soft tissues F. Peyraut a,, Z.-Q. Feng b, N. Labed a, C. Renaud b a b
Laboratoire M3M, Universite´ de Technologie de Belfort-Montbe´liard, 90010 Belfort, France ´ nerge´tique d’E ´ vry, Universite´ d’E ´vry-Val d’Essonne, 40 rue du Pelvoux, 91020 E ´vry, France ´canique et d’E Laboratoire de Me
a r t i c l e in fo
abstract
Article history: Received 25 August 2009 Received in revised form 11 February 2010 Accepted 19 February 2010
This paper concerns the modeling of biological soft tissues in the framework of anisotropic hyperelasticity. A closed form solution is proposed in the special case where uniaxial unconstrained tension loading is applied to a strip modeled by Holzapfel–Gasser–Ogden’s (HGO) strain energy. The classical Cardano’s formula is used to calculate the solution. This solution is investigated for different values of b which represents the angle between the collagen fibers and the circumferential direction. This study allows for the understanding of the reason why a one-to-one correspondence does not exist between the principal stretch and the fourth invariant of the right Cauchy–Green deformation tensor. It is demonstrated that the relationship becomes non-bijective when b is greater than a critical angle of 54:733 . The HGO model is implemented into an in-house finite element program and numerical results are in good agreement with theoretical ones. & 2010 Elsevier Ltd. All rights reserved.
Keywords: Anisotropic hyperelasticity Biological soft tissues HGO model Theory and modeling
1. Introduction To investigate internal deformation and stress of biological soft tissues such as ligaments, tendons or arterial walls, anisotropic hyperelastic constitutive laws are often used in the framework of finite element analysis [1–3]. The most used strain-energy functions take a power law form [4] or present an exponential behavior [5,6]. More recently, Balzani et al. [7] 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 been also recently developed to capture the inter-fiber angle change by adding to the strain energy the contribution of the fiber-matrix shear interaction [8]. 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 [9–11]. It is usually assumed that anisotropy is due to the collagen fibers behavior [12], 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 [1,7]. Each anisotropic density refers to a preferred direction of the material. For example, in Holzapfel–Gasser–Ogden’s constitutive law [6,13], two transversely isotropic energies with two distinct preferred directions corresponding to two superposed collagen fiber families.
Corresponding author. Tel.: + 33 3 84 58 31 91; fax: + 33 3 84 58 31 46.
E-mail address: [email protected] (F. Peyraut). 0020-7462/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijnonlinmec.2010.02.003
In this paper, a theoretical study related to a uniaxial unconstrained tension loading providing homogeneous deformation is considered with Holzapfel–Gasser–Ogden’s constitutive law. The nonlinear equilibrium equations are reduced into a cubic polynomial form by introducing an appropriate cylindrical coordinate system. The polynomial equation is solved thanks to Cardano’s formula and the analytical solution is successfully compared to numerical one obtained thanks to the university finite element program FER [14]. The results are investigated for different values of the material parameter b representing the non-deformed initial angle between the collagen fibers and the circumferential direction.
2. HGO hyperelastic model Most energy densities used to model transversely isotropic and orthotropic soft tissues take a power law form [4] or present an exponential behavior [5,6]. It is usually assumed that anisotropy is due to the collagen fibers behavior [12] while the ground substance, or matrix, behaves in an isotropic manner. The energy densities modeling transversely isotropic and orthotropic soft tissues are thus separated into isotropic and anisotropic parts [1] W ¼ Wiso þ
n X
a Wani
ð1Þ
a¼1
Each anisotropic density Waani 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,
ARTICLE IN PRESS 536
F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
a2
β β
a1
Fig. 1. Angle b for circumferentially oriented strip.
Holzapfel–Gasser–Ogden’s constitutive law [6,12] superposes two transversely isotropic energies with two distinct preferred directions a1 and a2 corresponding to two fiber families: 8 9 8 9 > >
Fig. 2. Cylindrical coordinate system.
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 (Fig. 1). b is predicted to be 433 39 in the work of Balzani et al. [7]. In our work, we consider b as a model parameter varying from 01 to 901. The strain measure adopted here is the Green–Lagrangian strain tensor 1 E ¼ ðCIÞ 2
a2
t
ð3Þ
@x @u ¼ Iþ , @X @X
J ¼ detðFÞ 4 0
I2 ¼ trðcof ðCÞÞ,
t Fig. 3. Unconstrained uniaxial tension test.
I3 ¼ detðCÞ,
J4 ¼ trðCMÞ,
J5 ¼ trðC2 MÞ ð5Þ
where cof(C) denotes the co-factor matrix of C and M is the so-called structural tensor representing the transverse-isotropy group and referring to a preferred direction a of the material M¼aa
ð6Þ
It is noted that (5) and (6) give J4 ¼ trðFT Fa aÞ ¼ JFaJ2
t
ð4Þ
X,x and u represent respectively the reference and the current positions and the displacement vector of a material point. According to Zhang–Rychlewski’s theorem [15], 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 I1, I2 and I3 of the right Cauchy–Green deformation tensor and two additional mixed invariants J4 and J5 [9–11] I1 ¼ trðCÞ,
x2
x1
where I is the second order unit tensor, C¼F F is the right Cauchy–Green deformation tensor and F the transformation gradient defined by
a1
t
x3
T
F¼
stress tensor r are given by S¼2
@W pC1 , @C
r ¼ J1 FSFT
ð8Þ
To uncouple the deviatoric part to the dilatational part of the response, the volume preserving part F^ ¼ J1=3 F of the deformation is introduced [1]. The modified invariants related to T C^ ¼ F^ F^ ¼ J 2=3 C are expressed from (5) by I^ 1 ¼ I1 J 2=3 ,
I^2 ¼ I2 J4=3 ,
a J^ 4 ¼ J4a J 2=3 ,
a J^ 5 ¼ J5a J 4=3
ð9Þ
The exponential type HGO density adopted in this work uses these modified invariants as follows: 8 a a k 2 > < 1 ½ek2 ðJ^ 4 1Þ 1 if J^ 4 Z 1 a ð10Þ Wani ¼ 2k2 > a :0 if J^ o 1 4
ð7Þ
The double brackets represent the usual Euclidian norm. The square root of J4 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 y represents the deformed angle between the collagen fibers and the circumferential direction (Fig. 2). In what follows, J4 and y will be used to solve the equilibrium equations in an analytical manner. In the case of hyperelastic materials, there exists an elastic potential function W which is a scalar function of the strain tensors. Generally, soft biological tissues are assumed to be incompressible. The incompressibility is characterized in terms of principal stretches l1 , l2 and l3 by J ¼ detðFÞ ¼ l1 l2 l3 ¼ 1 related to an extra pressure p. The second Piola–Kirchhoff stress tensor S and the corresponding Cauchy
This energy density is case sensitive with respect to Ja4 because the case of J4a o 1 represents the shortening of the fibers which is assumed to generate no stress. The proof of convexity of (10) with respect to F is given in [4,7]. The non-collagenous matrix of the media is modeled by the neo-Hookean isotropic density Wiso ¼ c1 ðI^1 3Þ
ð11Þ
It is noted that the volumetric-isochoric split of the above HGO model does only hold for (quasi) incompressible 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. [16] where a simple compressible anisotropic analytical model is developed. In order to conduct analytical solutions, we consider a uniaxial unconstrained tension test as depicted in Fig. 3.
ARTICLE IN PRESS F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
This loading is usually applied to biological soft tissues for identifying parameters of the material model. In this test, the parameters are set to c1 ¼10.2069 kPa, k1 ¼0.0017 kPa and k2 ¼882.847. These values are proposed in [7] to fit numerical models with experimental data. As the test leads to homogeneous deformations, the equilibrium equations are automatically satisfied and the boundary conditions related to the free traction surfaces are
s22 ¼ s33 ¼ 0
ð12Þ
The principal stretches l1 , l2 and l3 will be calculated in next sections by using an analytical approach and compared to finite element results. The incompressibility constraint is included in the finite element model by using a penalty term of the form k(J 1)2/2 with k¼ 105.
537
direction (Fig. 2) 1=2
cosðyÞ ¼ l1 cJ 4
,
1=2
sinðyÞ ¼ l2 sJ 4
ð19Þ
If (19) is reported in (17), the principal stretches l1 and l2 are replaced by two new variables y and J4. Eq. (17) becomes then a cubic polynomial equation: 3=2
s2 cJ 4 ¼A cosðyÞcos3 ðyÞ ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2k1 s2 2 1þ ðJ4 1Þek2 ðJ4 1Þ c1
ð20Þ
To calculate the roots of (20) with Cardano’s formula, three different cases have to be discussed with respect to the sign of the discriminant D
D ¼ A2 =41=27
ð21Þ
Case a: 3. Theoretical study
2
D 4 0 3 A 4 pffiffiffi
ð22Þ
3 3
By using (1), (8)–(11), it is easily demonstrated that 0 B J4 Z 1 : r ¼ p@
80 2 l > > > : 1
1
1
1 1
þ 4k1 ðJ^ 4 1Þe
k2 ðJ^ 4 1Þ2
0
C I^1 B C @ A 3 2
l22
1 1
l2 1 l2
80 2 2 > > > :
1
l22 s2 0
0 1 C J^ 4 B C @ A 3
1
19 > > = C A > > 1 ;
19 > > = C A > > 1 ;
ð13Þ 0 B J4 o 1 : r ¼ p@
80 2 l > > > : 1
1
1
1 1
l22
0 1 C I^1 B C @ A 3 2
l2 1 l2
1
19 > > = C A > > 1 ;
ð14Þ By considering s33 ¼ 0 to express the pressure p in terms of the principal stretch and taking into account the incompressibility condition J¼1, we obtain 8 < s11 ¼ 2½c1 ðl2 l2 l2 Þ þ 2k1 l2 c2 ðJ4 1Þek2 ðJ4 1Þ2 1 1 2 1 ð15Þ J4 Z 1 : 2 2 2 k2 ðJ4 1Þ2 : s22 ¼ 2½c1 ðl22 l2 l Þ þ 2k l ¼0 1 2 s ðJ4 1Þe 1 2 2
2 2
J4 o 1 : s11 ¼ 2c1 ðl1 l1 l2 Þ,
2 s22 ¼ 2c1 ðl22 l2 1 l2 Þ ¼ 0
ð16Þ
To determine the principal stretches in an analytical manner, we will demonstrate that the equilibrium equations can be transformed into a third degree polynomial form. The polynomial equation can then be algebraically solved by using Cardano’s formula. As the stresses expressions (15) and (16) depend on the value of J4, we logically start by considering two different cases.
It results from the definition (20) of A A rs2 c
ð23Þ
We consider the function of y representing the left-hand side of (20) f ðyÞ ¼ cosðyÞcos3 ðyÞ ¼ cosðyÞsin2 ðyÞ,
y A ½0,903
ð24Þ
It ispeasy ffiffiffi to demonstrate that the maximum value of f is equal to 2=3 3 and is reached if y ¼ bc with (Fig. 4) pffiffiffi bc ¼ Arccosð 3=3Þ 54:733 ð25Þ It then results from (23) that pffiffiffi A r2=3 3
ð26Þ
This result is contradictory with (22). Consequently, Case a will never occur. Case b: 2
D ¼ 0 3 A ¼ pffiffiffi
ð27Þ
3 3
In this case, Cardano’s formula gives one single and one double real root pffiffiffi pffiffiffi cosðy1 Þ ¼ 2 3=3, cosðy2 Þ ¼ 3=3 ð28Þ According to (24), the roots must be positive. The single root cosðy1 Þ is thus not acceptable. However, the double root cosðy2 Þ is acceptable as it is included in the range [0, 1]. The corresponding
0.4
2/3√3
0.3
The free traction surfaces condition s22 ¼ 0 applied to (15) gives 2
2 2
2
2
c1 ðl2 l1 l2 Þ þ2k1 l2 s2 ðJ4 1Þek2 ðJ4 1Þ ¼ 0
¼1
0.2
ð17Þ
0.1
Besides, it results from (2) and (5) that 1=2 1=2 ðl1 cJ 4 Þ2 þ ðl2 sJ 4 Þ2
f (θ)
3.1. Case J4 Z1
ð18Þ
This suggests that a cylindrical coordinate system could be appropriate to solve the non-linear equation (17). This cylindrical coordinate system is defined by setting J4½ as the radial coordinate and by introducing the polar angle y representing the deformed angle between the collagen fibers and the circumferential
0
βc 0
10
20
30
40 50 θ (°)
60
Fig. 4. Left-hand side of (20).
70
80
90
ARTICLE IN PRESS 538
F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
stretches are deduced from (19) pffiffiffi pffiffiffi pffiffiffi l1 ¼ J41=2 =ð 3cÞ, l2 ¼ 2J41=2 =ð 3sÞ
ð29Þ
As the maximum value of f is reached if y ¼ bc (Fig. 4), (27) gives A ¼ f ðbc Þ Z f ðbÞ
single root defined by cosðy1 Þ. If D is negative and J4 ¼1, there exists two distinct single roots defined by cosðy1 Þ and cosðy3 Þ. 3.2. Case J4 o1
ð30Þ
Besides, it results from the definitions (20) and (24) of A and f that If J4 41 : A o f ðbÞ
If J4 ¼ 1 : A ¼ f ðbÞ
ð31Þ
The conditions (30) and (31) are together consistent if and only if b ¼ bc and J4 ¼1. It is then deduced from (29) that Case b gives thus a double root corresponding to the evident solution l1 ¼ l2 ¼ 1. Case c: 2
D o0 3 Ao pffiffiffi
ð32Þ
3 3
In this case, according to Cardano’s formula, there exist three distinct real roots " pffiffiffi # pffiffiffi 2 3 a þ2ði1Þp 3 3A cos cosðyi Þ ¼ , a ¼ Arccos , i ¼ 1,2,3 3 3 2 ð33Þ The locations on the trigonometric circle of the three angles used in (33) are depicted in Fig. 5. To determine the acceptable roots, we first notice that (20) implies that A is positive. It then results from (32) and (33) that pffiffiffi 3 ocosðy1 Þ o1, cosðy2 Þ o 0 0 o cosðy3 Þ o ð34Þ 3 It is noted that the second root cosðy2 Þ is negative and must be excluded. Among the two positive roots cosðy1 Þ and cosðy3 Þ, we propose to demonstrate that the third is unacceptable except if J4 is equal to 1. To prove it, we remind that a tensile loading test is considered. It then results from (19) and (33) that pffiffiffi a þ4p ð35Þ ðc 3Þ Z1 l1 ¼ 2J41=2 cos 3 After some algebra, it can be demonstrated that (35) leads to J4 r 1. As the case of J4 Z1 is considered, we have J4 ¼1. The root cosðy3 Þ is thus excluded except if J4 ¼1. In conclusion, if D is negative and J4 is strictly greater than 1, there exists only one
This case is easier to study than the previous one because the orthotropic behavior is not taken into account and the model is reduced to the neo-Hookean energy density according to Eq. (10). By using the stress expression (16) in conjunction with the free traction surface loading s22 ¼ 0, the neo-Hookean solution is directly obtained pffiffiffiffiffi ð36Þ l2 ¼ 1= l1 To study the dependence of l1 and l2 on J4 in the same way as for the case J4 Z 1, it is needed to solve again a cubic polynomial equation. This equation is deduced from (18) and (36) 3
gðl1 Þ ¼ l1 J4 c2 l1 þ t 2 ¼ 0
with t ¼ s=c
ð37Þ
This function admits a local maximum and a local minimum as shown in Fig. 6. The local maximum and minimum points are qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi defined respectively by ð J4 =3c2 ,t 2 þ 4J43 =27c6 Þ and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð J4 =3c2 ,t 2 4J43 =27c6 Þ. To determine the real roots of g, Eq. (37) is solved analytically by using again Cardano’s formula. Three different cases occur depending on the sign of the discriminant D1
D1 ¼ t4 =4J43 =27c6
ð38Þ
Case a:
D1 40 3 t2 4
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4J43 =ð27c6 Þ
ð39Þ
In this case, Cardano’s formula indicates that one root is a real number while the two others are conjugate complex numbers. Besides, it is reminded that the condition l1 Z 1 must be satisfied because a tensile load is considered. As g increases from 1 to qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t2 þ 4J43 =ð27c6 Þ 40 for l1 varying from 1 to J4 =ð3c2 Þ, the real root is negative and should not be selected. Finally, using definitions (24) and (25), the condition (39) can be simplified to J4 o½f ðbÞ=f ðbc Þ2=3
ð40Þ
Case b:
sin
D1 ¼ 0 3 t2 ¼
α/3
ð41Þ
2.4
1.2
π/6
(α+2π)/3
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4J43 =ð27c6 Þ
cos
π/6
g
π/6 0
−1.2
−2.4 −1.8 (α+4π)/3 Fig. 5. Location of a=3, ða þ 2pÞ=3 and ða þ 4pÞ=3 on the trigonometric circle.
−1.2
−0.6
0 λ1
0.6
Fig. 6. Graph of g ðb ¼ 503 ,J4 ¼ 0:5Þ.
1.2
1.8
ARTICLE IN PRESS F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
In this case, Cardano’s formula gives a single and a double real root Single root : l1 ¼ 3s2 =J4 ,
Double root : l1 ¼ 3s2 =2J4
ð42Þ
The single root is negative and has to be excluded. By using (41), it is noted that the double root corresponds to the minimum of the g function. To make this root acceptable, the condition l1 Z1 must hold such that pffiffiffi pffiffiffi ðt4 =4Þ1=3 Z1 3 t Z 2 3 b Z bc ¼ Arccosð 3=3Þ 54:733 ð43Þ The double root defined by (42) is thus solution of Eq. (37) if and only if the phenomenological angle b is greater than the critical angle bc defined by (43). It is worth noting that this critical angle has also been found by Guo et al. with a composites-based hyperelastic constitutive model [17]. Moreover, it results from (19), (36), (41), (42) that cosy ¼ cosbc ,
siny ¼ sinbc ) y ¼ bc
ð45Þ
Case c:
D1 o 0 3 t2 o
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4J43 =ð27c6 Þ
ð46Þ
In this case, Cardano’s formula provides three distinct real roots expressed trigonometrically rffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffi g J4 J4 g þ 2p cos cos , , l1a ¼ 2 l ¼ 2 1b 2 2 3 3 3c 3c rffiffiffiffiffiffiffiffi J4 g þ 4p cos ð47Þ l1c ¼ 2 2 3 3c qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
g ¼ Arccos½t2 27c6 =ð4J43 Þ
In the same way, it is concluded that l1c o1. Consequently, the two positive roots are unacceptable. Case c2: pffiffiffi b Z bc ¼ Arccosð 3=3Þ ð53Þ As the tangent function increases in the interval ½0, p=2, (46) and (53) yield to qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð54Þ t 2 Z 2 ) J4 =ð3c2 Þ 4 1 To conclude that l1c 41, we consider Eq. (54) and we apply again the intermediate value theorem in the same way as in Case c1. It results next from (49) that l1a 4 1. Consequently, the two positive roots l1a and l1c defined by (47) are acceptable. Finally, by reasoning in the same manner as for Case a, the condition (46) can be simplified to J4 4½f ðbÞ=f ðbc Þ2=3
ð55Þ
ð44Þ
It is thus remarkable that the deformed angle y is equal to the critical angle bc in the case of a double root occurrence corresponding to an initial angle b greater than bc . Finally, it is noted that condition (41) can be simplified in the same way as for Case a: J4 ¼ ½f ðbÞ=f ðbc Þ2=3
539
ð48Þ
To determine the acceptable roots, we first deduce from (46)–(48) that qffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 o l1c o J4 =ð3c2 Þ o l1a o J4 =c2 , 2 J4 =ð3c2 Þ o l1b o J4 =c2 ð49Þ The second root l1b is negative and cannot be selected. To determine whether or not the positive roots l1a and l1c are greater than one, according to the condition of a tension test, two cases are considered Case c1: pffiffiffi b o bc ¼ Arccosð 3=3Þ ð50Þ Because the cosine function decreases in the interval ½0, p=2, (50) yields to qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð51Þ 3c2 41 ) 0 o J4 =ð3c2 Þ o1 Eq. (51) means that the minimum of g is reached between 0 and 1. We then apply the intermediate value theorem to g and we account for condition (46) sffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffi! 4J43 J4 1J4 2 ¼t g o0, gð1Þ ¼ 2 4 0 2 3c 27c6 c rffiffiffiffiffiffiffiffi !, J4 ) (l1a A ð52Þ ,1 gðl1a Þ ¼ 0 3c2
4. Analytical and numerical results The analytical formula applied to the uniaxial tension test determines the evolution of l2 and y vs. J4 as shown in Fig. 7. If b o bc (e.g., b ¼ 203 ), there is no solution corresponding to J4 o 1. The biological tissue can thus not be shortened in the fiber direction. On the contrary, if b 4 bc (e.g., b ¼ 703 ), it is observed that the biological tissue can be shortened in the fiber direction. However, the correspondence between l2 and J4 is not one-to-one because of the existence of two distinct solutions. These two solutions are given by (47) and (48). It is finally noted that the curve related to b ¼ 703 starts with J4 equal to a minimum value defined by (45). This minimum value corresponds to the double root defined by (42). Eq. (45) shows that the minimum value of J4 only depends on the initial angle b. As the function f decreases if b 4 bc (Fig. 4), this minimum value decreases if b increases. This trend is well observed in Fig. 7 (left). The minimum values of J4 are respectively equal to 0.98, 0.85 and 0.58 for b ¼ 603 ,703 and 803 . The occurrence of a double root with J4 o 1 corresponds to a remarkable result given by (44), i.e., the deformed angle y is equal to the critical angle bc when the minimum value of J4 given by (45) is reached as shown in Fig. 7 (right). Fig. 8 shows the variation of the cylindrical coordinates (J4, y) vs. l1 for two values of b (203 and 703 ). As a tension test is considered in this work, l1 is greater than 1. In the case of b ¼ 703 , it is observed that J4 is first lower than 1 (the biological tissue is first shortened in the fiber direction) and decreases to a minimum value. This case has been previously discussed in Section 3.2 (Case c2) and can only occur if b Z bc . The minimum value of J4 is equal to 0.85 according to (45). This minimum is reached if l1 is approximately equal to 1.56 (Fig. 8 a) and corresponds to a deformed angle y equal to the critical angle bc (Fig. 8 b). After this minimum value, J4 increases and becomes greater than 1 if l1 is approximately greater than 2.29. This value indicates the transition from a purely isotropic to an anisotropic behavior of biological pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tissues. It can be analytically calculated by using l1 ¼ ð 1þ 4t 2 1Þ=2. This formula is easily obtained by solving (37) together with J4 ¼1. It is worth noting that if b o bc (e.g., b ¼ 203 Þ,J4 increases with l1 (Fig. 8 a) so that the biological tissue is always stretched in the fiber direction ðJ4 Z1Þ. The critical angle bc plays obviously a key role for the understanding of the model behavior. However, it should be noted that the minimum value of J4 and the correspondence between the solution and J4 were only discussed in the case of J4 o1. As the orthotropic part of the HGO energy density is assumed to generate no stress if J4 o 1, the critical angle concerns only the neo-Hookean part of the model. The results related to bc given in this paper could thus be extended to any anisotropic
ARTICLE IN PRESS 540
F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
90
1
β=60° β=70° β=80°
80 0.8
70 60
0.6
βc=54.73°
θ
λ2
50
0.4 0.2 0 0.5
0.6
40 30
β=20° β=60° β=70° β=80°
20 10 0.7
0.8
0.9 J4
1
1.1
1.2
0 0.5
1.3
0.6
0.7
0.8
0.9 J4
1
1.1
1.2
1.3
Fig. 7. l2 and y vs. J4 for different values of b.
70
1.3
60
1.2 1.1 θ
J4
40
1
1
1.5
2
2.5
30 20
β=20° β=70°
0.9 0.8
β=20° β=70°
50
10 3
0
3.5
1
1.5
2
λ1
2.5
3
3.5
λ1
Fig. 8. Evolution of the cylindrical coordinates (J4, y) for two values of b.
1.4 1.2
λ2, λ3
1
1
1
0.8 3
0.6
3
Fig. 10. Finite element modeling: initial and deformed shapes (left: l1 o 2:29, right: l1 4 2:29).
0.4 0.2 0
2
2
λ2−analytical λ3−analytical λ2−numerical λ3−numerical
isotropic neo-Hookean model and the relationship between l1 and l2 is defined by the neo-Hookean solution (36).
1
1.5
2 λ1
2.5
3
Fig. 9. Comparison between numerical and analytical results ðb ¼ 703 Þ.
model superposing any orthotropic density with the isotropic neo-Hookean density. The HGO model has been implemented into an in-house finite element code FER [14]. Fig. 9 shows a good agreement between the numerical results obtained with FER and the analytical calculation presented in this work with b ¼ 703 . Beyond a threshold value of l1 2:29, it is observed that l3 increases while l2 decreases. In other words, the section shortens in the direction 2 and swells in the direction 3 as shown in Fig. 10. As outlined in [12], this trend typically represents an anisotropic behavior. Before this transition, the HGO model reduces to the
5. Conclusions In this paper, a theoretical study has been performed with Holzapfel–Gasser–Ogden’s (HGO) constitutive law and applied to a unconstrained tension test. The free traction loading condition related to this test is expressed as a non-linear algebraic equation to give a relationship between the principal stretches l1 and l2 . To change this equation into a cubic polynomial form, we propose to use a cylindrical coordinate system by setting the square root of the mixed invariant J4 as the radial coordinate. The deformed angle y, between the collagen fibers and the circumferential direction, is used as the polar angle. A closed form solution has been obtained thanks to Cardano’s formula. A critical angle ðbc 54:733 Þ has been analytically determined.
ARTICLE IN PRESS F. Peyraut et al. / International Journal of Non-Linear Mechanics 45 (2010) 535–541
To evaluate the dependence of the principal stretches on J4, the analytical solutions were studied for different values of the initial angle b between the collagen fibers orientation and the circumferential direction. It appears that the correspondence between the solution and J4 is not one to one if b is greater than bc while uniqueness holds if b is lower than bc . We have also established that the biological tissue cannot be shortened in the fiber direction more than a minimum value of J4 given by an analytical formula. We have finally demonstrated that this minimum value corresponds to a double root. In this case, it is remarkable that the deformed angle y is equal to the critical angle bc . However, it should be noted that the minimum value of J4 and the correspondence between the solution and J4 were only discussed in the case of J4 o1. As the orthotropic part of the HGO energy density is assumed to generate no stress if J4 o 1, the critical angle bc only concerns the neo-Hookean part of the model. The results related to bc given in this paper could thus be extended to other anisotropic hyperelasticity models superimposing any orthotropic density with the isotropic neo-Hookean density. The closed form solutions are also used to validate the implementation of the HGO model into the finite element code FER. A fair agreement has been found between theoretical and numerical solutions. If b is greater than bc , the results underline the transition between purely isotropic and orthotropic behavior. In particular, it is observed that the strip section shortens in one direction and swells in another direction beyond a threshold value of l1 corresponding to J4 ¼1. References [1] J.A. Weiss, B.N. Maker, S. Govindjee, Finite element implementation of incompressible, transversely isotropic hyperelasticity, Comput. Meth. Appl. Mech. Eng. 135 (1996) 107–128.
541
[2] E.S. Almeida, R.L. Spilker, Finite element formulation for hyperelastic transversely isotropic biphasic soft tissues, Comput. Meth. Appl. Mech. Eng. 151 (1998) 513–538. ¨ [3] M. Ruter, E. Stein, Analysis, finite element computation and error estimation in transversely isotropic nearly incompressible finite elasticity, Comput. Meth. Appl. Mech. Eng. 190 (2000) 519–541. ¨ [4] J. Schroder, P. Neff, D. Balzani, A variational approach for materially stable anisotropic hyperelasticity, Int. J. Solids Struct. 42 (2005) 4352–4371. [5] Y.C. Fung, K. Fronek, P. Patitucci, Pseudoelasticity of arteries and the choice of its mathematical expression, Am. J. Physiol. 237 (1979) H620–H631. [6] G.A. Holzapfel, T.C. Gasser, R.W. Ogden, A new constitutive framework for arterial wall mechanics and a comparative study of material models, J. Elasticity 61 (2000) 1–48. ¨ [7] D. Balzani, P. Neff, J. Schroder, G.A. Holzapfel, A polyconvex framework for soft biological tissues. Adjustment to experimental data, Int. J. Solids Struct. 43 (2006) 6052–6070. [8] X.Q. Peng, Z.Y. Guo, B. Moran, An anisotropic hyperelastic constitutive model with fiber-matrix shear interaction for the human annulus fibrosus, J. Appl. Mech. Trans. ASME 73 (5) (2006) 815–824. [9] A.J.M. Spencer, Isotropic polynomial invariants and tensor functions, in: J.P. Boehler (Ed.), Applications of Tensor Functions in Solids Mechanics, CISM Course No. 292, Springer Verlag, Wien, 1987. [10] J.P. Boehler, Introduction to the invariant formulation of anisotropic constitutive equations, in: J.P. Boehler (Ed.), Applications of Tensor Functions in Solids Mechanics, CISM Course No. 292, Springer Verlag, Wien, 1987. [11] Q.S. Zheng, A.J.M. Spencer, Tensors which characterize anisotropies, Int. J. Eng. Sci. 31 (5) (1993) 679–693. [12] T.C. Gasser, R.W. Ogden, G.A. Holzapfel, Hyperelastic modelling of arterial layers with distributed collagen fibre orientations, J. R. Soc. Interface 3 (2006) 15–35. [13] G.A. Holzapfel, T.C. Gasser, A viscoelastic model for fiber-reinforced composites at finite strains: continuum basis, computational aspects and applications, Comput. Meth. Appl. Mech. Eng. 190 (2001) 4379–4403. [14] Z.Q. Feng /http://lmee. univ-evry.fr/ feng/FerSystem.htmlS, 2008. [15] J.M. Zhang, J. Rychlewski, Structural tensors for anisotropic solids, Arch. Mech. 42 (1990) 267–277. [16] Z.Y. Guo, F.C. Caner, X.Q. Peng, B. Moran, On constitutive modelling of porous neo-Hookean composites, J. Mech. Phys. Solids 56 (6) (2008) 2338–2357. [17] Z.Y. Guo, X.Q. Peng, B. Moran, A composites-based hyperelastic constitutive model for soft tissue with application to the human annulus fibrosus, J. Mech. Phys. Solids 54 (9) (2006) 1952–1971.