- Email: [email protected]
On the tension–compression switch in soft fibrous solids
European Journal of Mechanics A/Solids 49 (2015) 561e569
Contents lists available at ScienceDirect
European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol
On the tensionecompression switch in soft fibrous solids Gerhard A. Holzapfel a, *, Ray W. Ogden b a b
Institute of Biomechanics, Graz University of Technology, Kronesgasse 5-I, 8010 Graz, Austria School of Mathematics and Statistics, University of Glasgow University Gardens, Glasgow G12 8QW, Scotland, UK
a r t i c l e i n f o
a b s t r a c t
Article history: Received 20 May 2014 Accepted 18 September 2014 Available online
Many soft fibrous solids consist of an isotropic matrix with strong embedded fibers. In continuum models of such materials it is often considered that the fibers do not support compression and that their contribution to the overall mechanical response should be neglected when they are under compression. A tensionecompression ‘switch’ is used to turn off the contribution of the fibers in the model, but it is apparent that an inappropriate switch is often adopted and the predictions of the mechanical response can then be significantly different for different switching criteria. On the basis of a simple illustration this paper draws attention to the consequences of using an inappropriate switch, and identifies possible solutions. One of the possible solutions for rectifying the problem is provided in some detail by introducing a modified fiber distribution model and applying it to both plane strain and three-dimensional transversely isotropic fiber distributions. © 2014 Elsevier Masson SAS. All rights reserved.
Keywords: Fiber distribution Fibrous solid Constitutive modeling
1. Introduction Many biological tissues and engineering materials are composites that consist of fibers embedded in a matrix, which in many cases can be treated as an isotropic material. In the application of continuum mechanics to such structures it is often assumed that the fibers do not support compression. The implementation of this condition requires the use of a tensionecompression ‘switch’ that turns off the (anisotropic) mechanical contribution of the fibers when they undergo a compressive strain. However, there is some confusion in the literature and in finite element software concerning the nature of this switch. For example, an incorrect switch is adopted for the frequently-used dispersed fiber model from Gasser et al. (2006) in the Abaqus finite element software (Dassault mes Simulia Corp, 2013) which does not exclude fibers under Syste compression. One aim of this paper, therefore, is to highlight and explain this problem and then to offer possible modifications that circumvent the problem. In the paper by Gasser et al. (2006) we developed a model for describing the elastic properties of soft fibrous materials with dispersed fiber orientations. This model was based on the notion of a generalized structure tensor H defined, in slightly different notation from that used in Gasser et al. (2006), by
* Corresponding author. E-mail address: [email protected] (G.A. Holzapfel). http://dx.doi.org/10.1016/j.euromechsol.2014.09.005 0997-7538/© 2014 Elsevier Masson SAS. All rights reserved.
H¼
1 4p
Z rðQ; FÞNðQ; FÞ5NðQ; FÞdU;
(1)
U
where Q and F are spherical polar angles, U is the unit sphere defined by U ¼ S[0, 2p], with Q2S ¼ [0, p], r(Q, F) is an orientation density function normalized according to
1 4p
Z rðQ; FÞsinQdQdF ¼ 1;
(2)
U
and N(Q, F) is a unit vector representing the direction of an arbitrary fiber within a distribution of fibers with diverse orientations and weight r(Q, F). When specialized to a transversely isotropic distribution with r independent of F, H can be reduced to the form
H ¼ kI þ ð1 3kÞM5M;
(3)
where I is the identity tensor, the unit vector M is the mean direction of the distribution, k2[0, 1/3] is a dispersion parameter defined by
k¼
1 4
Z S
rðQÞsin3 QdQ;
(4)
562
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
and Z the normalization of r reduces to
1 2
and
rðQÞsinQdQ ¼ 1:
(5)
S
Note, in particular, the limiting case k ¼ 1/3, for which r(Q)≡1. The other limiting case, k ¼ 0, will be discussed in detail in Section 2.2.2. In principle the upper limit of k could be 1/2, but as shown in Holzapfel and Ogden (2010) the range (1/3,1/2] of values of k can lead to unphysical results. In its incompressible form,1 the model in Gasser et al (2006) has a strain-energy (or free-energy) function J given by
JðC; Hi Þ ¼ Jg ðCÞ þ
X
Jf i ðC; Hi Þ;
(6)
i¼1;2
where C ¼ FTF is the right CauchyeGreen tensor, F is the deformation gradient with the incompressibility constraint detF ¼ 1, Jg is the strain energy of the isotropic groundmatrix, Jfi is the free energy associated with the i th family of fibers (i ¼ 1,2), Hi is a generalized structure tensor of the form (3), with k taken to be the same for each family of fibers and with M replaced by a0i, i ¼ 1,2, which define the mean fiber directions for families i ¼ 1,2 in the reference configuration. Note that k represents a weighting parameter between isotropy (in the limit k ¼ 1/3) and transverse isotropy (for k ¼ 0). More specifically, Jg and Jfi were taken to have the forms
Jg ¼
1 mðI 3Þ; 2 1
Jf i ¼
k1 2k2
exp k2 Ei2 1 ;
(7)
where the constant m denotes the neo-Hookean parameter, k1 is a parameter with dimension of stress, k2 is a dimensionless parameter, I1 ¼ tr C is the first principal invariant of C and
JðC; H1 ; H2 Þ ¼ UðJÞ þ
i 1 k Xh 2 m I1 3 þ 1 exp k2 Ei 1 ; 2 2k2 i¼1;2 (11)
where U(J) is a penalty term designed to enforce incompressibility. Thus, the volumetric and isochoric terms are decoupled. In Gasser et al. (2006) it was required that the fibers contribute only when extended (and not in compression), strictly when either I41 > 1 and/or I42 > 1. Because of the formulation based on the volumetric/isochoric split these are typically replaced by I 4i^C : (a0i 5 a0i) > 1, i ¼ 1, 2, which have a slightly different interpretation and significantly different consequences (Vergori et al., 2013). However, in Abaqus, for example, the switch that turns the fibers on was instead taken as Ei > 0, which is equivalent + + > 1. It was stated in to I 4i > 1, or, for strict incompressibility, I4i Dassault Systemes Simulia Corp (2013) that this corresponds to the strain in the fiber direction being positive, which is not always the case (the strain in the mean fiber direction is actually I4i 1). Part of the purpose of the present paper is to demonstrate that this switch can give erroneous results. mes Simulia Corp, 2013) the In Abaqus 6.13 (Dassault Syste volumetric part of the strain-energy function has the form
UðJÞ ¼
1 J2 1 ln J ; D 2
(12)
where D is a material parameter (D=2 is the compressibility of the material in the reference configuration), and the notation J el is used instead of the J used here. They also considered N fiber families and the notation C10 for what we call m/2. 2. Problem motivation for simple tension
+ Ei ¼ Hi : C 1 ¼ I4i 1;
i ¼ 1; 2;
(8)
where Hi : C ¼ tr(HiC), with tr standing for the trace of a secondorder tensor. In the above we have defined the standard notation + I4i ¼ kI1 þ ð1 3kÞI4i ;
I4i ¼ C : ða0i 5a0i Þ;
i ¼ 1; 2:
(9)
We recall from the appendix of Gasser et al. (2006) that the strain energy (7) is a convex function of the stretches if Ei > 0. In fact, it is easy to see from (9) that Ei > 0 whenever I4i > 1 since for an isochoric deformation I1 3. However, as we shall elaborate in the following section, Ei > 0 may be associated with I4i < 1, which implies that the mean fiber direction is under compression! For purposes of finite element implementation of the model it is cast in compressible form by using the isochoric/volumetric split F ¼ J1/3F (Flory, 1961; Ogden, 1978), where J ¼ det F, det F ¼ 1 and an overbar on various kinematical quantities indicates their isochoric modification. Thus,
C ¼ J 2=3 C;
I 1 ¼ J 2=3 I1 ;
I 4i ¼ J 2=3 I4i ;
To illustrate the problem highlighted in the previous section we consider a single transversely isotropic distribution of fibers with the mean direction M and a dispersion of the fibers with parameter k. For purposes of illustration it suffices to restrict attention to strict incompressibility. We consider a homogeneous uniaxial deformation with stretch l1 ¼ l in the mean fiber direction (taken as the 1 direction) with lateral contraction. Then, by symmetry, the lateral stretches are l2 ¼ l3 ¼ l1/2, the deformation gradient has the component form
i h ½F ¼ diag l; l1=2 ; l1=2 ;
Note the model is designed strictly for use only in the incompressible case, although its numerical implementation is based on a slightly compressible version. Specific constraints have to be used in the numerical scheme in order to enforce incompressibility. One example is the Augmented Lagrangian method (see, for example, Simo and Taylor, 1991).
(13)
and M has components [M] ¼ [1, 0, 0]T so that the current fiber direction m ¼ FM gives [m] ¼ [l, 0, 0]T ¼ l[M]. The first invariant is I1 ¼ tr C, while the fourth invariant associated with the direction M is denoted by I4 ¼ M,(CM). For the considered deformation these are given by
I1 ¼ l2 þ 2l1 ;
Ei ¼ Hi : C 1; (10)
1
2.1. An invariant-based continuum model
I4 ¼ l2 :
(14)
We also adopt the notation
I4+ ¼ kI1 þ ð1 3kÞI4 ;
(15)
+ defined in (9) , analogously to that for the generalized invariant I4i 1 so that the switch for the single distribution of fibers is written I4+ > 1. Then, with (14),
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
I4+ ¼ k l2 þ 2l1 þ ð1 3kÞl2 ¼ l2 ð1 2kÞ þ 2kl1 :
(16)
The dependence of I4+ on l is illustrated in Fig. 1, for several values of k together with the line I4+ ¼ I4+ ð1Þ ¼ 1. As can be seen from this plot, I4+ is greater than 1 for l > 1; however, for I4+ > 1, l is not in general greater than 1. Hence, the condition I4+ > 1 does not necessarily imply that the mean fiber direction is under extension. In addition, from Fig. 1, one can conclude that I4+ has a minimum for each k2[0, 1/3], specifically
3k2=3 ð1 2kÞ1=3
for
l3 ¼
k 1; 1 2k
1 l2 þ l þ 1
¼ cos2 Qc ;
(20)
the right-hand equality in which defines the maximum value of Q, say Qc, for which the fibers are extended, the subscript c standing for ‘critical’. From this relation we see that the fibers should only contribute to the strain-energy function if cos2Q is larger than cos2Qc, i.e. only fibers for which 0 Q < Qc (equivalently pQc < Q p), with 0 < Qc < p/2, contribute. For fibers outside these ranges the anisotropic term in the strain-energy function should be dropped.
(17)
the right-hand equality holding only for the limiting case k ¼ 1/3. The fact that I4+ > 1 does not necessarily require l > 1 is the first important point to note in this section. The second important point is that, for a given distribution of fibers, there is a range of the angles of individual fibers for which the fibers are shortened when l, the stretch in the mean fiber direction, is larger that 1, as we now show. These fibers should not therefore contribute to the energy function because, as is commonly assumed, they do not support compression. Consider a general direction N in the reference configuration which makes an angle Q with M. In three dimensions its (Cartesian) components are given in terms of spherical polar angles Q2[0, p], F2[0, 2p] by
½N ¼ ½cosQ; sinQ cosF; sinQ sinFT ;
cos2 Q >
563
2.1.1. Transversely isotropic elasticity For a general incompressible transversely isotropic elastic material with preferred direction M, J is a function of four invariants, typically taken to be
I1 ¼ tr C; I2 ¼
(21) and the Cauchy stress tensor s is given by (see, for example, Holzapfel, 2000)
s ¼ pI þ 2j1 b þ 2j2 I1 b b2 þ 2j4 m5m þ 2j5 ðm 5 bm þ bm 5 mÞ;
(18)
but because of the symmetry of the deformation about the 1 direction the square of the stretch in the direction N, i.e. n,n ¼ N,(CN), where n ¼ FN, is independent of F. Thus, with (13), we obtain
N,ðCNÞ ¼ l2 cos2 Q þ l1 sin2 Q h i 1 ≡1 þ l ðl 1Þ l2 þ l þ 1 cos2 Q 1 ¼ I4 ðQÞ; (19) wherein I4(Q) is defined, this being the I4 invariant associated with the fiber direction N. It follows that for N,(CN) > 1 with l > 1, we require
1 2 I1 tr C2 ; I4 ¼ M,ðCMÞ; I5 ¼ M, C2 M ; 2
(22)
where b ¼ FFT and ji ¼ vJ/vIi, i ¼ 1,2,4,5. In the present situation we are restricting attention to dependence of J on I1 and I4, so that the above reduces to
s ¼ pI þ 2j1 b þ 2j4 m5m;
(23)
and we consider that the stress response coefficient j4 should be set to 0 for Q in the range [pQc, Qc]. With respect to the mean fiber direction and the transverse direction(s) the Cauchy stress components are calculated as
s11 ¼ 2ðj1 þ j4 Þl2 p;
s22 ¼ s33 ¼ 2j1 l1 p;
(24)
while the normal Cauchy stress in the deformed fiber direction n ¼ FN, denoted s(n), is given by
h sðnÞ ¼ 2 j1 l4 cos2 Q þ l2 sin2 Q i. þ j4 l4 cos2 Q l2 cos2 Q þ l1 sin2 Q p:
(25)
If, in particular, we take the lateral stress to vanish (s22 ¼ s33 ¼ 0), so we have a state of uniaxial stress, then we can eliminate p from (24) and (25) to obtain
s11 ¼ 2j1 l2 l1 þ 2j4 l2
(26)
and
. sðnÞ ¼ l2 s11 cos2 Q l2 cos2 Q þ l1 sin2 Q :
Fig. 1. Plot of the function I4+ ¼ l2 ð1 2kÞ þ 2kl1 for different values of the dispersion parameter k2[0,1/3].
(27)
Thus, for uniaxial tension s11 > 0 with l > 1, the stress in the deformed fiber direction is positive for all Q s p/2 even though I4(Q) < 1 for Qc < Q p/2, i.e. these fiber directions are shortened by the deformation. For the specialization of the model from Section 1 to one fiber family, but with the notation of the present section, we have
564
J¼
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
1 k mðI 3Þ þ 1 2 1 2k2
2 exp k2 I4+ 1 1 ;
(28)
and the derivatives of J with respect to I1 and I4 are given by
2 1 ; j1 ¼ m þ kk1 I4+ 1 exp k2 I4+ 1 2 2 ; j4 ¼ ð1 3kÞk1 I4+ 1 exp k2 I4+ 1
(29)
(30)
so that (26) leads to
s11 ¼ m l2 l1 þ 2k1 ð1 2kÞl2 kl1 2 ; I4+ 1 exp k2 I4+ 1
1 2
Z rðQÞWðI4 ðQÞÞsinQdQ;
where the definition of S is now changed so that it corresponds to an appropriate range of values of Q2[0, p], which is either [0, p] itself if all fibers are counted or [0, Qc) if fibers for which I4(Q) 1 are excluded from the energy, and hence from the stress. Note that this formulation is a special case of a formulation which will be given in Section 4 for a general deformation, and this has some similarities with the structural model of Lanir (1983) and the model of F-actin networks of Holzapfel et al. (2014). The energy of the system is assumed to be the uncoupled combination of an isotropic matrix energy depending only on I1, denoted Wiso(I1), and the fiber energy (32):
(31)
with I4+ given by (16)2. The curves in Fig. 2 show representative plots of the dimensionless s11, defined by s11 ¼ s11/m, versus l for k ¼ 0, 0.1, 1/3 with illustrative material constants k1 ¼ k1/m ¼ 2 and k2 ¼ 1.
(32)
S
Wiso ðI1 Þ þ
1 2
Z rðQÞWðI4 ðQÞÞsinQdQ:
(33)
S
The associated Cauchy stress is then given by
s ¼ pI þ 2W 0 iso ðI1 Þb þ
Z
rðQÞW 0 ðI4 ðQÞÞn5nsinQdQ;
(34)
S
where n has components lcosQ in the direction of M and l1/2 sinQ laterally, and a prime indicates differentiation with respect to the considered argument, I1 or I4. By using symmetry this can be expressed in the form
s ¼ pI þ ½2W 0 iso ðI1 Þ þ bb þ ða bÞm5m;
(35)
where
Z a¼ Fig. 2. Dimensionless plots of s11 ¼ s11/m versus l from (31) for the values 0, 0.1, 1/3 of k.
The angle Qc is not included explicitly in the term I4+ and therefore not in the strain-energy function. Thus, in order to account for the individual fibers within the distribution that are extended and to exclude those that are not, a modification of the model is required. At this point we emphasize that the correct tensionecompression switch is I4 > 1 and not I4+ > 1. With the current model it is not possible to exclude compressed fibers within a distribution for which the mean fiber direction is under extension. However, a possible modification of the model that allows for such an exclusion could be considered by including dependence on k and the critical angle Qc within the parameter k1, for example, bearing in mind that Qc depends on the state of deformation. We do not follow this approach, but instead we consider an alternative modification, based on the individual fibers within the distribution. This is discussed in detail in the following subsection for the same simple deformation, then for a general plane strain deformation in Section 3.2, and finally generalized for an arbitrary deformation in Section 4. 2.2. Distributed fiber model 2.2.1. An orientation dependent fiber distribution model We now use the notation I4(Q) defined in (19)1 and associate with each fiber in the distribution an energy function, say W(I4(Q)), per unit reference volume. The total fiber energy density with the orientation density r(Q) satisfying the normalization defined in (5) is then given by
S
1 b¼ 2
rðQÞW 0 ðI4 ðQÞÞcos2 QsinQdQ; Z
(36) rðQÞW 0 ðI4 ðQÞÞsin3 QdQ:
S
For uniaxial stress in the mean fiber direction we obtain
s11 ¼ 2W 0 iso l2 l1 þ al2 bl1
(37)
and
h i . sðnÞ ¼ ð2W 0 iso þ aÞl3 ð2W 0 iso þ bÞ lcos2 Q l2 cos2 Q þ l1 sin2 Q ; (38) which should be compared with (26) and (27), respectively. 2.2.2. The limiting case k ¼ 0 Here it is interesting to consider the special case in which all the fibers are in the same direction M, so there is no dispersion. Recall that (4) and (5) must hold in the limit k / 0, so that r(Q) effectively becomes a delta function. First we note that from the definition of I4(Q) in (19) we obtain
I4 ðQÞ I4 ð0Þ ¼ l2 l1 sin2 Q:
(39)
Let us now consider a ‘narrow’ distribution of fibers so that Q is ‘small’ and we can apply the Taylor expansion to W0 (I4(Q)) to obtain the approximation
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
h i 00 W 0 ðI4 ðQÞÞ ¼ W 0 ðI4 ð0ÞÞ þ I4 ðQÞ I4 ð0Þ W ðI4 ð0ÞÞ þ …
(40)
Substitution of this into the expressions for a and b in (36) and making use of equations (4), (5) and (39), we obtain the approximations
00 a ¼ 2ð1 2kÞW 0 ðI4 ð0ÞÞ 4 l2 l1 k A W ðI4 ð0ÞÞ þ …; (41) 00 b ¼ 2kW 0 ðI4 ð0ÞÞ 2 l2 l1 AW ðI4 ð0ÞÞ þ …;
565
function. Note that in the limit k / 1/3, b / 0 the limiting value of r(Q) is 1. Now we plot s11 in dimensionless form s11 ¼ s11/mf against l in Fig. 3 for the two cases (i) S ¼ (0, p) and (ii) S ¼ (0, Qc). We take the material parameters kf1 ¼ kf1/mf ¼ 2 and kf2 ¼ 1 for k ¼ 0, 0.1, 1/3, respectively. The integration over S ¼ [0, p] gives the solid curves, and the integration over S ¼ [0, Qc) gives the dashed curves. As can be seen in the figure, the model using Qc gives a weaker response. This is because the excluded fibers no longer contribute to the strain energy in the fibers. For the case k ¼ 0 the curves coincide because there is no dispersion (see the analysis in Section 2.2.2).
(42)
where
A¼
1 4
Zp
rðQÞsin5 QdQ:
(43)
0
Since A and k are non-negative then, with reference to (4), it is clear that A k and hence in the limit k / 0 we obtain A / 0 (and further terms in the Taylor expansion will similarly disappear). Thus, we obtain, in the limit, a ¼ 2W0 (I4(0)) and b ¼ 0, and the Cauchy stress tensor reduces to
s ¼ pI þ 2W 0 iso ðI1 Þb þ 2W 0 ðI4 Þm5m;
(44) 3. Plane strain problems
which is a particular example of a transversely isotropic model with just the two invariants I1 and I4 ¼ I4(0) given by (23). 2.2.3. A specific energy function Next we make the model more specific. We take
1 m ðI 3Þ; 2 f 1
Wiso ðI1 Þ ¼
(45)
where mf is a neo-Hookean constant, and
WðI4 Þ ¼
h i o kf1 n exp kf2 ðI4 1Þ2 1 ; 2kf2
a ¼ kf1
(46)
i h rðQÞðI4 ðQÞ 1Þexp kf2 ðI4 ðQÞ 1Þ2 cos2 QsinQdQ;
S
(47) b¼
1 k 2 f1
Z
i h rðQÞðI4 ðQÞ 1Þexp kf 2 ðI4 ðQÞ 1Þ2 sin3 QdQ:
S
(48) For specific calculations we take the orientation density function to be the standard p-periodic von Mises distribution of the form
rffiffiffiffiffiffiffi b 1 pffiffiffiffiffiffi exp 2bcos2 Q ; rðQÞ ¼ 4 2p erfi 2b
We recall that for a general incompressible transversely isotropic material the Cauchy stress is given by equation (22). When the deformation is restricted to plane strain we have I2 ¼ I1, as is well known, and also, with M lying in the considered plane, I5 ¼ (I1 1)I4 1, as shown in Merodio and Ogden (2002). Thus, in this case the strain energy is simply a function of I1 and I4 and the (planar) Cauchy stress tensor can be written as
s ¼ pI þ 2j1 b þ 2j4 m5m;
where the material constants kf1 and kf2 are taken to be the same for each fiber in the distribution. Note, in particular, that mf, kf1 and kf2 are not in general the same as the m, k1 and k2 in (28). Then, with I4(Q) as defined in (19), we obtain from (36)
Z
Fig. 3. Dimensionless plots of s11 ¼ s11/m versus l from (37) for the values 0, 0.1, 1/3 of k.
(49)
as we have used in Gasser et al. (2006), where b > 0 is the concentration parameter and erfi(x) ¼ i erf(ix) is the imaginary error
(50)
as in (23), except that all the tensors are now two dimensional (in the considered plane). We now consider a planar fiber dispersion with r(Q) ¼ r(Q) and Q2S, which is now the range [p/2, p/2], as distinct from [0, p] used in (5). The normalization condition is now
1 p
Z rðQÞdQ ¼ 1;
(51)
S
which is the counterpart for plane strain of the three-dimensional normalization (5). The corresponding planar generalized structure tensor, the counterpart of (1), is defined as (Ogden, 2009)
H¼
1 p
Z rðQÞN5NdQ;
(52)
S
so that, with planar C,
I4+ ¼ trðCHÞ ¼
1 p
Z rðQÞI4 ðQÞdQ;
(53)
S
where I4(Q) ¼ N,(CN). Relative to Cartesian axes E1 ¼ M and E2 we take N to have the form N ¼ cosQE1þsinQE2. The components of H in (52) are then
566
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
H11 ¼ 1 k;
H12 ¼ 0;
H22 ¼ k ¼
1 p
Z
rQsin2 QdQ;
(54)
S
which defines the dispersion parameter k for this situation. With the help of the planar identity tensor I ¼ E15E1 þ E25E2 and (54), we obtain.
H ¼ kI þ ð1 2kÞM5M;
(55)
c J I1 ; I4+ ¼ JðI1 ; I4 Þ:
b 4 ¼ n I + 1 ¼ nkc2 : j 4
b 1 ¼ 1 m; j 2
(64)
Substitution into (62) gives
(65)
s11 ¼ p þ m þ 2nk2 c2 1 þ c2 þ 2nkð1 2kÞc2 ;
(66)
s22 ¼ p þ m þ 2nk2 c2 :
(67)
(57) 3.2. Distributed fiber model
Then
b 1 þ kj b 4; j1 ¼ j
b 4; j4 ¼ ð1 2kÞ j
(58)
b 1 and j b 4 are the counterparts of j1 and j4 for c where j J. On substitution into (50), we obtain
b 1 þ kj b 4 b þ 2ð1 2kÞ j b 4 m5m: s ¼ pI þ 2 j
½n ¼ ½F½N ¼
Wiso ðI1 Þ þ
rðQÞWðI4 ðQÞÞdQ
(68)
S
I4 ðQÞ ¼ 1 þ csin2Q þ c2 sin2 Q;
(69)
where
n1 cosQ þ c sinQ ; ¼ n2 sinQ
S¼
where c is the amount of shear, [N] has components [cosQ, sinQ]T and [M] ¼ [1, 0]T, so that the mean fiber direction is in the direction of shear. Then, the invariants become
I4 ð0Þ ¼ 1;
I4+ ¼ 1 þ kc2 ; (61) the latter coming from (56). Note also that m ¼ M. In terms of the planar principal stretches l1 ¼ l, l2 ¼ l1, c ¼ ll1, with c > 0 corresponding to l > 1. Then, from (59),
b 1 þ kj b 4 c: s12 ¼ 2 j
Z
s ¼ pI þ 2W 0 iso ðI1 Þb þ 2FSFT ;
(60)
I1 ¼ 2 þ c2 ;
1 p
per unit area of the considered plane, where S represents the range of angles in [p/2, p/2] for which I4(Q) > 1. The Cauchy stress tensor is
Now let us consider a simple example of a deformation, that of simple shear with
c ; 1
Similarly to the discussion in Section 2.2, suppose now that the energy of the distribution embedded in an isotropic matrix characterized by Wiso is
(59)
3.1. Application to simple shear
1 ½F ¼ 0
with I4+ instead of I4, so that, with the help of (61)4,
The corresponding normal stresses are
(63)
(56)
where I4 ¼ I4(0), noting that I1 is now defined based on the twodimensional b (plane strain). Note, in particular, the factor 1 2k in the above compared with 1 3k in (15). If M is the mean fiber direction associated with a distribution then we may consider, instead of J(I1, I4),
2 1 1 mðI1 3Þ þ n I4+ 1 ; 2 2
s12 ¼ m þ 2nk2 c2 c:
the counterpart of (3), from which we obtain the planar I4+ as
I4+ ¼ kI1 þ ð1 2kÞI4 ;
J¼
(62)
In this example, the critical angle Qc for which I4(Q) ¼ 1 is given by tanQc ¼ 2/c. If we take c > 0 and sinQc < 0, for example, then Qc2(p/2, 0), and the fibers are extended if either p/2 < Q < Qc or 0 < Q < p/2 in which case the definition of S should be modified so that only the contributions of the extended fibers are included. Example. Let us now adopt a neo-Hookean material with a standard reinforcing model (see, for example, Qiu and Pence, 1997) of the form
1 p
Z
rðQÞW 0 ðI4 ðQÞÞN5NdQ:
(70)
S
The latter expression must have the form
S ¼ aM5M þ bðI M5MÞ þ g M5M⊥ þ M⊥ 5M ;
(71)
where M⊥ is a unit vector perpendicular to M in the considered plane, with the properties M,N ¼ cosQ, M⊥,N ¼ sinQ. The coefficients are given by
a¼
b¼
g¼
1 p 1 p 1 p
Z
rðQÞW 0 ðI4 ðQÞÞcos2 QdQ;
(72)
rðQÞW 0 ðI4 ðQÞÞsin2 QdQ;
(73)
rðQÞW 0 ðI4 ðQÞÞsinQcosQdQ;
(74)
S
Z S
Z S
which, in principle, can be calculated for any W and specified plane deformation. Note that, since r(Q) is an even function of Q but, in general, I4(Q) is not, then gs0 in general.
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
Thus, the Cauchy stress tensor (69) takes on the form
M ¼ coscE1 þ sincE2 ;
567
M0 ¼ coscE1 sincE2 ;
(84)
0
s ¼ pI þ 2ðW iso ðI1 Þ þ bÞb þ 2ða bÞm5m þ 2g m5m⊥ þ m⊥ 5m ;
(75)
where m⊥ ¼ FM⊥. Note that the final term here with the coefficient g does not have a counterpart in (59). We now introduce a specific constitutive model for Wiso, given by (45) with the standard reinforcing model for W(I4), so that
W 0 iso ¼
1 m; 2 f
(76)
(85)
(77)
where e1 ¼ FE1 and e2 ¼ FE2, which replace m and m⊥, respectively, in (75), while j1 ¼ W 0 iso ðI1 Þ. This has the same structure as the Cauchy stress (75) in the fiber distribution model with a, b and g replaced by ðj4 þ j6 Þcos2 c, ðj4 þ j6 Þsin2 c and ðj4 j6 Þsinccosc, respectively.
Z
nf p
rðQÞðI4 ðQÞ 1ÞdQ: S
3.2.1. Application to simple shear Now let us specialize to simple shear. Then, it is easy to show from (72)e(74), using (76)2 and (61)2, that 2 +
a ¼ nf c a ;
2 +
b ¼ nf c b ;
+
g ¼ 2nf ca ;
(78)
where
a+ ¼
1 p
Z
rðQÞsin2 Qcos2 QdQ;
b+ ¼
S
1 p
Z
rðQÞsin4 QdQ:
S
(79) +
s12 ¼ 2ðW
iso
4. Three-dimensional transversely isotropic distribution Without loss of generality we may choose Cartesian axes E1, E2, E3 so that the mean fiber direction M ¼ E3. Then the direction N of an arbitrary fiber within a distribution with mean direction M is given by
N ¼ sinQcosFE1 þ sinQsinFE2 þ cosQM;
N,M ¼ cosQ: (86)
+
The values of a and b depend on the orientation density distribution r(Q) and Qc only. Since we are only including extending R Qc R p=2 fibers the integration should be p=2 þ 0 , where we emphasize that Qc is contained in (p/2, 0). Then, from (75) with (76)1 and (78) 0
i
s ¼ pI þ 2 j1 þ ðj4 þ j6 Þsin2 c b þ 2ðj4 þ j6 Þcos2ce1 5e1 þ2ðj4 j6 Þsinccoscðe1 5e2 þ e2 5e1 Þ;
W 0 ðI4 Þ ¼ nf ðI4 1Þ;
and hence, by (51) and (53)2, we obtain, from (72) and (73),
aþb¼
where the (fixed) angle c is introduced to characterize the orientations of the mean fiber directions relative to E1, and let j1, j4 and j6 be the constitutive coefficients. The Cauchy stress tensor is then
þ bÞc þ 2g ¼ mf þ 4nf a
+
+ 3
c þ 2nf b c :
(80)
Let r(Q) be the orientation density, normalized according to (5). In general, I4 ¼ N,(CN) depends on both Q and F so that, using (86), we obtain
I4 ðQ; FÞ ¼ cos2 QM,ðCMÞ þ sin2 Q i h cos2 FE1 ,ðCE1 Þ þ sin2 FE2 ,ðCE2 Þ þ2sinQcosQðcosFE1 þ sinFE2 Þ,ðCMÞ
By comparing this formula with that in (65), we obtain the connections.
m ¼ mf þ 4nf a+ ;
nk2 ¼ nf b+
(81)
between the material parameters in the two models. Thus, in this case, the consequences of the two models are essentially identical, but only as far as the shear response is concerned. The corresponding normal stresses are
þ2sin2 QsinFcosFE1 ,ðCE2 Þ:
Referring to (1), the generalized structure tensor takes on the form
H¼
þðE1 cosF þ E2 sinFÞ5M o þsin2 QsinFcosF E1 5E2 þ E2 5E1 sinQdQdF:
(82) (83)
However, these do not match the normal stresses in (66) and (67). This is not surprising since the structures of the two models are slightly different by virtue of the term in g in (75). However, if one considers two families of distributed fibers in the plane instead of one, with mean fiber directions M and M′ symmetric with respect to the axes E1 and E2, then the structure of the distributed fiber model can be reproduced. In standard notation (see, for example, Holzapfel, 2000) let I4 ¼ M,(CM) and I6 ¼ M′,(CM′) be the squares of the stretches in the mean fiber directions with
Z n 1 rðQÞ cos2 QM5M þ sin2 Q 4p U cos2 FE1 5E1 þ sin2 FE2 5E2 þsinQcosQ½M5ðE1 cosF þ E2 sinFÞ
s11 ¼ p þ mf þ 2nf b+ c2 1 þ c2 þ 2nf c2 a+ b+ þ 8nf c2 a+ ;
s22 ¼ p þ mf þ 2nf b+ c2 :
(87)
(88)
By carrying out the integration over F2[0, 2p] we are left with
H¼
1 2
Z rðQÞcos2 QsinQdQM5M S
1 þ 4
Z
(89) rðQÞsin3 QdQðE1 5E1 þ E2 5E2 Þ:
S
With the properties (4) and (5) and since E15E1þE25E2 ¼ I M5M, we recover (3), the generalized structure tensor for a transversely isotropic distribution.
568
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
4.1. Distributed fiber model
The Cauchy stress tensor (92) can now be written as
Next, let W(I4(Q, F)) be the strain energy of the fiber in the direction N (per unit reference volume). Then, with weighting r(Q), the total strain energy of the fibers is
1 4p
Z rðQÞWðI4 ðQ; FÞÞsinQdQdF;
(90)
U
which extends the energy (32) to three dimensions. Strictly this should only include fibers for which I4 > 1. Let Wiso (I1) be the contribution of the matrix. Then, assuming that there is no coupling, the composite has strain energy
Wiso ðI1 Þ þ
1 4p
rðQÞWðI4 ðQ; FÞÞsinQdQdF;
(101)
where m ¼ FM, e1 ¼ FE1 and e2 ¼ FE2. Now suppose that the properties of the fibers are given by the standard reinforcing model (76)2 and note that
a þ b þ b0 ¼
1 4p
(91)
Z
rðQÞW 0 ðI4 ðQ; FÞÞsinQdQdF
U
Z
(102) rðQÞðI4 ðQ; FÞ 1ÞsinQdQdF:
U 0
U
For the coefficients a, b, b , g, g0 and d we obtain, after some detailed calculations,
a ¼ Anf ð3I4 ð0; 0Þ I1 Þ þ nf BðI4 ð0; 0Þ 1Þ þ knf ðI1 5I4 ð0; 0Þ þ 2Þ; (103)
(92)
where
1 4p
þdðe1 5e2 þ e2 5e1 Þ;
Z
s ¼ pI þ 2W 0 iso ðI1 Þb þ 2FSFT ;
S¼
þgðm5e1 þ e1 5mÞ þ g0 ðm5e2 þ e2 5mÞ
n ¼ f 4p
generalizing (33), and the corresponding Cauchy stress tensor is
Z
s ¼ pI þ 2W 0 iso ðI1 Þb þ 2½be1 5e1 þ b0 e2 5e2 þ am5m
b ¼ knf ðI4 ð0; 0Þ 1Þ þ rðQÞW 0 ðI4 ðQ; FÞÞN5NsinQdQdF;
Anf ð3e1 ,e1 þ e2 ,e2 4I4 ð0; 0ÞÞ; 4
(104)
(93)
U
generalizing the definition (70). The latter expression can be written in the form
S ¼ aM5M þ bE1 5E1 þ b0 E2 5E2 þ gðM5E1 þ E1 5MÞ þg0 ðM5E2 þ E2 5MÞ þ dðE1 5E2 þ E2 5E1 Þ;
b0 ¼ knf ðI4 ð0; 0Þ 1Þ þ
Anf ðe1 ,e1 þ 3e2 ,e2 4I4 ð0; 0ÞÞ; 4 (105)
g ¼ Cnf m,e1 ;
g0 ¼ Cnf m,e2 ;
d ¼ nf De1 ,e2 ;
(106)
(94) the choice of E1 and E2 perpendicular to M such that E1,E2 ¼ 0 being arbitrary. Then,
a¼
1 4p
1 b¼ 4p
b0 ¼
g¼
g0 ¼
1 4p 1 4p 1 4p
1 d¼ 4p
Z
where
A¼ rðQÞW 0 ðI4 ðQ; FÞÞcos2 QsinQdQdF;
(95)
1 4
U
Z
rðQÞW 0 ðI4 ðQ; FÞÞsin Qcos2 FdQdF; 3
(96)
C¼
1 2
U
Z
rðQÞW 0 ðI4 ðQ; FÞÞsin3 Qsin2 FdQdF;
(97) k¼
rðQÞW 0 ðI4 ðQ; FÞÞcosQsinQcosFdQdF;
(98)
B¼
S
Z
rðQÞsin2 Qcos2 QdQ;
1 2
Z rðQÞsinQdQ;
(107)
S
D¼
S
1 8
Z
rðQÞsin4 QdQ;
S
rðQÞW 0 ðI4 ðQ; FÞÞcosQsinQsinFdQdF;
(99)
1 4
Z
rðQÞsin3 QdQ:
(109)
S 0
0
The coefficients a, b, b , g, g and d can in principle be calculated for any W.
U
Z
rðQÞsin5 QdQ;
(108)
U
Z
Z
4.2. Identification of the critical angles
U
Z
The goal is now to find the range of angles for which 0
2
rðQÞW ðI4 ðQ; FÞÞsin QsinFcosFdQdF; U
(100)
and we note that cosFE1þsinFE2 is an arbitrary vector normal to M, so that interchanging cosF and sinF makes no essential difference. These are the generalizations of the a, b and g introduced in Section 3.2.
I4 ¼ n,n ¼ N,ðCNÞ > 1
(110)
for a fiber distribution with mean fiber direction M. We assume a transversely isotropic distribution. With N,M ¼ cosQ we want to find the critical angles Q and F such that
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
N,ðCNÞ ¼ 1:
(111)
Explicitly, with
N ¼ sinQðcosFE1 þ sinFE2 Þ þ cosQE3 ;
(112)
(111) gives
I4 ðNÞ ¼ N,ðCNÞ ¼ sin2 Q C11 cos2 F þ 2C12 sinFcosF þ C22 sin2 F þC33 cos2 Q þ 2sinQcosQðC13 cosF þ C23 sinFÞ ¼ 1; (113) where Cij are the components of C. For each Q2[0, p] this defines the critical values of F beyond which N,(CN) < 1. In the special case in which E1, E2, E3 are eigenvectors of C, Cij ¼ 0 for isj and Cii ¼ l2i , i ¼ 1, 2, 3, and the incompressibility l1l2l3 ¼ 1, and eq. (113) reduces to
h i tan2 Q l22 l21 sin2 F þ l21 1 ¼ 1 l23 ;
cos2F ¼
l21 l22
:
excluded. This model has been developed both for plane strain and for general three-dimensional deformations, and yields explicit expressions for the Cauchy stress deformation relation. These are compared with those for the models based on the continuum approach described in Gasser et al. (2006). In particular, for simple shear we have established connections between the material parameters in the two models for the shear stresses to coincide. For general three-dimensional deformations we have provided an algorithm for obtaining the critical angles that can be used within a general stress computation. The stress relation and the (explicit) expressions for the critical angles provide a basis for finite element implementation of the model.
Acknowledgment The authors are grateful to Dr. Kewei Li, Graz University of Technology, Austria, for valuable discussions.
(114) References
from which F is easily determined for any given Q in the form
2 1 l23 cot2 Q þ 2 l21 þ l22
569
(115)
5. Concluding remarks In this study we have analyzed the effect of using an inappropriate tensionecompression switch for excluding the contribution of fibers under compression within a soft fibrous composite solid. In particular, for a simple tension illustration we have shown that the stiffness of the composite is in general significantly overestimated by not excluding fibers under compression, and we have computed the critical angle of the individual fibers relative to the fiber direction that corresponds to the switch. The correct tensionecompression switch is I4 > 1 and not I4+ > 1. However, for the current model, shortened fibers within a distribution with mean fiber direction under extension cannot be excluded. We have proposed a possible modification of the model within which such an exclusion can be considered. In addition, we have introduced a distributed fiber model with a weighted energy function that also allows shortened fibers to be
mes Simulia Corp, 2013. Abaqus 6.13 Analysis User's Guide. In: SecDassault Syste tion 22.5.3: Anisotropic Hyperelastic Behavior. Flory, P.J., 1961. Thermodynamic relations for high elastic materials. Trans. Farad. Soc. 57, 829e838. Gasser, T.C., Ogden, R.W., Holzapfel, G.A., 2006. Hyperelastic modelling of arterial layers with distributed collagen fiber orientations. J. R. Soc. Interface 3, 15e35. Holzapfel, G.A., 2000. Nonlinear Solid Mechanics: a Continuum Approach for Engineering. John Wiley & Sons, Chichester. Holzapfel, G.A., Ogden, R.W., 2010. Constitutive modelling of arteries. Proc. R. Soc. Lond. A 466, 1551e1597. Holzapfel, G.A., Unterberger, M.J., Ogden, R.W., 2014. An affine continuum mechanical model for cross-linked F-actin networks with compliant linker proteins. J. Mech. Behav. Biomed. Mater. 38, 78e90. Lanir, Y., 1983. Constitutive equations for fibrous connective tissues. J. Biomech. 16, 1e12. Merodio, J., Ogden, R.W., 2002. Material instabilities in fiber-reinforced nonlinearly elastic solids under plane deformation. Arch. Mech. 54, 525e552. Ogden, R.W., 1978. Nearly isochoric elastic deformations: application to rubberlike solids. J. Mech. Phys. Solids 26, 37e57. Ogden, R.W., 2009. Anisotropy and nonlinear elasticity in arterial wall mechanics. In: Holzapfel, G.A., Ogden, R.W. (Eds.), Lecture Notes, CISM Course on the Biomechanical Modelling at the Molecular, Cellular and Tissue Levels, CISM Courses and Lectures Series, vol. 508. Springer, Wien, pp. 179e258. Qiu, G.Y., Pence, T.J., 1997. Remarks on the behavior of simple directionally reinforced incompressible nonlinearly elastic solids. J. Elast. 49, 1e30. Simo, J.C., Taylor, R.L., 1991. Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms. Comput. Meth. Appl. Mech. Eng. 85, 273e310. Vergori, L., Destrade, M., McGarry, P., Ogden, R.W., 2013. On anisotropic elasticity and questions concerning its Finite Element implementation. Comp. Mech. 52, 1185e1197.
Contents lists available at ScienceDirect
European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol
On the tensionecompression switch in soft fibrous solids Gerhard A. Holzapfel a, *, Ray W. Ogden b a b
Institute of Biomechanics, Graz University of Technology, Kronesgasse 5-I, 8010 Graz, Austria School of Mathematics and Statistics, University of Glasgow University Gardens, Glasgow G12 8QW, Scotland, UK
a r t i c l e i n f o
a b s t r a c t
Article history: Received 20 May 2014 Accepted 18 September 2014 Available online
Many soft fibrous solids consist of an isotropic matrix with strong embedded fibers. In continuum models of such materials it is often considered that the fibers do not support compression and that their contribution to the overall mechanical response should be neglected when they are under compression. A tensionecompression ‘switch’ is used to turn off the contribution of the fibers in the model, but it is apparent that an inappropriate switch is often adopted and the predictions of the mechanical response can then be significantly different for different switching criteria. On the basis of a simple illustration this paper draws attention to the consequences of using an inappropriate switch, and identifies possible solutions. One of the possible solutions for rectifying the problem is provided in some detail by introducing a modified fiber distribution model and applying it to both plane strain and three-dimensional transversely isotropic fiber distributions. © 2014 Elsevier Masson SAS. All rights reserved.
Keywords: Fiber distribution Fibrous solid Constitutive modeling
1. Introduction Many biological tissues and engineering materials are composites that consist of fibers embedded in a matrix, which in many cases can be treated as an isotropic material. In the application of continuum mechanics to such structures it is often assumed that the fibers do not support compression. The implementation of this condition requires the use of a tensionecompression ‘switch’ that turns off the (anisotropic) mechanical contribution of the fibers when they undergo a compressive strain. However, there is some confusion in the literature and in finite element software concerning the nature of this switch. For example, an incorrect switch is adopted for the frequently-used dispersed fiber model from Gasser et al. (2006) in the Abaqus finite element software (Dassault mes Simulia Corp, 2013) which does not exclude fibers under Syste compression. One aim of this paper, therefore, is to highlight and explain this problem and then to offer possible modifications that circumvent the problem. In the paper by Gasser et al. (2006) we developed a model for describing the elastic properties of soft fibrous materials with dispersed fiber orientations. This model was based on the notion of a generalized structure tensor H defined, in slightly different notation from that used in Gasser et al. (2006), by
* Corresponding author. E-mail address: [email protected] (G.A. Holzapfel). http://dx.doi.org/10.1016/j.euromechsol.2014.09.005 0997-7538/© 2014 Elsevier Masson SAS. All rights reserved.
H¼
1 4p
Z rðQ; FÞNðQ; FÞ5NðQ; FÞdU;
(1)
U
where Q and F are spherical polar angles, U is the unit sphere defined by U ¼ S[0, 2p], with Q2S ¼ [0, p], r(Q, F) is an orientation density function normalized according to
1 4p
Z rðQ; FÞsinQdQdF ¼ 1;
(2)
U
and N(Q, F) is a unit vector representing the direction of an arbitrary fiber within a distribution of fibers with diverse orientations and weight r(Q, F). When specialized to a transversely isotropic distribution with r independent of F, H can be reduced to the form
H ¼ kI þ ð1 3kÞM5M;
(3)
where I is the identity tensor, the unit vector M is the mean direction of the distribution, k2[0, 1/3] is a dispersion parameter defined by
k¼
1 4
Z S
rðQÞsin3 QdQ;
(4)
562
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
and Z the normalization of r reduces to
1 2
and
rðQÞsinQdQ ¼ 1:
(5)
S
Note, in particular, the limiting case k ¼ 1/3, for which r(Q)≡1. The other limiting case, k ¼ 0, will be discussed in detail in Section 2.2.2. In principle the upper limit of k could be 1/2, but as shown in Holzapfel and Ogden (2010) the range (1/3,1/2] of values of k can lead to unphysical results. In its incompressible form,1 the model in Gasser et al (2006) has a strain-energy (or free-energy) function J given by
JðC; Hi Þ ¼ Jg ðCÞ þ
X
Jf i ðC; Hi Þ;
(6)
i¼1;2
where C ¼ FTF is the right CauchyeGreen tensor, F is the deformation gradient with the incompressibility constraint detF ¼ 1, Jg is the strain energy of the isotropic groundmatrix, Jfi is the free energy associated with the i th family of fibers (i ¼ 1,2), Hi is a generalized structure tensor of the form (3), with k taken to be the same for each family of fibers and with M replaced by a0i, i ¼ 1,2, which define the mean fiber directions for families i ¼ 1,2 in the reference configuration. Note that k represents a weighting parameter between isotropy (in the limit k ¼ 1/3) and transverse isotropy (for k ¼ 0). More specifically, Jg and Jfi were taken to have the forms
Jg ¼
1 mðI 3Þ; 2 1
Jf i ¼
k1 2k2
exp k2 Ei2 1 ;
(7)
where the constant m denotes the neo-Hookean parameter, k1 is a parameter with dimension of stress, k2 is a dimensionless parameter, I1 ¼ tr C is the first principal invariant of C and
JðC; H1 ; H2 Þ ¼ UðJÞ þ
i 1 k Xh 2 m I1 3 þ 1 exp k2 Ei 1 ; 2 2k2 i¼1;2 (11)
where U(J) is a penalty term designed to enforce incompressibility. Thus, the volumetric and isochoric terms are decoupled. In Gasser et al. (2006) it was required that the fibers contribute only when extended (and not in compression), strictly when either I41 > 1 and/or I42 > 1. Because of the formulation based on the volumetric/isochoric split these are typically replaced by I 4i^C : (a0i 5 a0i) > 1, i ¼ 1, 2, which have a slightly different interpretation and significantly different consequences (Vergori et al., 2013). However, in Abaqus, for example, the switch that turns the fibers on was instead taken as Ei > 0, which is equivalent + + > 1. It was stated in to I 4i > 1, or, for strict incompressibility, I4i Dassault Systemes Simulia Corp (2013) that this corresponds to the strain in the fiber direction being positive, which is not always the case (the strain in the mean fiber direction is actually I4i 1). Part of the purpose of the present paper is to demonstrate that this switch can give erroneous results. mes Simulia Corp, 2013) the In Abaqus 6.13 (Dassault Syste volumetric part of the strain-energy function has the form
UðJÞ ¼
1 J2 1 ln J ; D 2
(12)
where D is a material parameter (D=2 is the compressibility of the material in the reference configuration), and the notation J el is used instead of the J used here. They also considered N fiber families and the notation C10 for what we call m/2. 2. Problem motivation for simple tension
+ Ei ¼ Hi : C 1 ¼ I4i 1;
i ¼ 1; 2;
(8)
where Hi : C ¼ tr(HiC), with tr standing for the trace of a secondorder tensor. In the above we have defined the standard notation + I4i ¼ kI1 þ ð1 3kÞI4i ;
I4i ¼ C : ða0i 5a0i Þ;
i ¼ 1; 2:
(9)
We recall from the appendix of Gasser et al. (2006) that the strain energy (7) is a convex function of the stretches if Ei > 0. In fact, it is easy to see from (9) that Ei > 0 whenever I4i > 1 since for an isochoric deformation I1 3. However, as we shall elaborate in the following section, Ei > 0 may be associated with I4i < 1, which implies that the mean fiber direction is under compression! For purposes of finite element implementation of the model it is cast in compressible form by using the isochoric/volumetric split F ¼ J1/3F (Flory, 1961; Ogden, 1978), where J ¼ det F, det F ¼ 1 and an overbar on various kinematical quantities indicates their isochoric modification. Thus,
C ¼ J 2=3 C;
I 1 ¼ J 2=3 I1 ;
I 4i ¼ J 2=3 I4i ;
To illustrate the problem highlighted in the previous section we consider a single transversely isotropic distribution of fibers with the mean direction M and a dispersion of the fibers with parameter k. For purposes of illustration it suffices to restrict attention to strict incompressibility. We consider a homogeneous uniaxial deformation with stretch l1 ¼ l in the mean fiber direction (taken as the 1 direction) with lateral contraction. Then, by symmetry, the lateral stretches are l2 ¼ l3 ¼ l1/2, the deformation gradient has the component form
i h ½F ¼ diag l; l1=2 ; l1=2 ;
Note the model is designed strictly for use only in the incompressible case, although its numerical implementation is based on a slightly compressible version. Specific constraints have to be used in the numerical scheme in order to enforce incompressibility. One example is the Augmented Lagrangian method (see, for example, Simo and Taylor, 1991).
(13)
and M has components [M] ¼ [1, 0, 0]T so that the current fiber direction m ¼ FM gives [m] ¼ [l, 0, 0]T ¼ l[M]. The first invariant is I1 ¼ tr C, while the fourth invariant associated with the direction M is denoted by I4 ¼ M,(CM). For the considered deformation these are given by
I1 ¼ l2 þ 2l1 ;
Ei ¼ Hi : C 1; (10)
1
2.1. An invariant-based continuum model
I4 ¼ l2 :
(14)
We also adopt the notation
I4+ ¼ kI1 þ ð1 3kÞI4 ;
(15)
+ defined in (9) , analogously to that for the generalized invariant I4i 1 so that the switch for the single distribution of fibers is written I4+ > 1. Then, with (14),
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
I4+ ¼ k l2 þ 2l1 þ ð1 3kÞl2 ¼ l2 ð1 2kÞ þ 2kl1 :
(16)
The dependence of I4+ on l is illustrated in Fig. 1, for several values of k together with the line I4+ ¼ I4+ ð1Þ ¼ 1. As can be seen from this plot, I4+ is greater than 1 for l > 1; however, for I4+ > 1, l is not in general greater than 1. Hence, the condition I4+ > 1 does not necessarily imply that the mean fiber direction is under extension. In addition, from Fig. 1, one can conclude that I4+ has a minimum for each k2[0, 1/3], specifically
3k2=3 ð1 2kÞ1=3
for
l3 ¼
k 1; 1 2k
1 l2 þ l þ 1
¼ cos2 Qc ;
(20)
the right-hand equality in which defines the maximum value of Q, say Qc, for which the fibers are extended, the subscript c standing for ‘critical’. From this relation we see that the fibers should only contribute to the strain-energy function if cos2Q is larger than cos2Qc, i.e. only fibers for which 0 Q < Qc (equivalently pQc < Q p), with 0 < Qc < p/2, contribute. For fibers outside these ranges the anisotropic term in the strain-energy function should be dropped.
(17)
the right-hand equality holding only for the limiting case k ¼ 1/3. The fact that I4+ > 1 does not necessarily require l > 1 is the first important point to note in this section. The second important point is that, for a given distribution of fibers, there is a range of the angles of individual fibers for which the fibers are shortened when l, the stretch in the mean fiber direction, is larger that 1, as we now show. These fibers should not therefore contribute to the energy function because, as is commonly assumed, they do not support compression. Consider a general direction N in the reference configuration which makes an angle Q with M. In three dimensions its (Cartesian) components are given in terms of spherical polar angles Q2[0, p], F2[0, 2p] by
½N ¼ ½cosQ; sinQ cosF; sinQ sinFT ;
cos2 Q >
563
2.1.1. Transversely isotropic elasticity For a general incompressible transversely isotropic elastic material with preferred direction M, J is a function of four invariants, typically taken to be
I1 ¼ tr C; I2 ¼
(21) and the Cauchy stress tensor s is given by (see, for example, Holzapfel, 2000)
s ¼ pI þ 2j1 b þ 2j2 I1 b b2 þ 2j4 m5m þ 2j5 ðm 5 bm þ bm 5 mÞ;
(18)
but because of the symmetry of the deformation about the 1 direction the square of the stretch in the direction N, i.e. n,n ¼ N,(CN), where n ¼ FN, is independent of F. Thus, with (13), we obtain
N,ðCNÞ ¼ l2 cos2 Q þ l1 sin2 Q h i 1 ≡1 þ l ðl 1Þ l2 þ l þ 1 cos2 Q 1 ¼ I4 ðQÞ; (19) wherein I4(Q) is defined, this being the I4 invariant associated with the fiber direction N. It follows that for N,(CN) > 1 with l > 1, we require
1 2 I1 tr C2 ; I4 ¼ M,ðCMÞ; I5 ¼ M, C2 M ; 2
(22)
where b ¼ FFT and ji ¼ vJ/vIi, i ¼ 1,2,4,5. In the present situation we are restricting attention to dependence of J on I1 and I4, so that the above reduces to
s ¼ pI þ 2j1 b þ 2j4 m5m;
(23)
and we consider that the stress response coefficient j4 should be set to 0 for Q in the range [pQc, Qc]. With respect to the mean fiber direction and the transverse direction(s) the Cauchy stress components are calculated as
s11 ¼ 2ðj1 þ j4 Þl2 p;
s22 ¼ s33 ¼ 2j1 l1 p;
(24)
while the normal Cauchy stress in the deformed fiber direction n ¼ FN, denoted s(n), is given by
h sðnÞ ¼ 2 j1 l4 cos2 Q þ l2 sin2 Q i. þ j4 l4 cos2 Q l2 cos2 Q þ l1 sin2 Q p:
(25)
If, in particular, we take the lateral stress to vanish (s22 ¼ s33 ¼ 0), so we have a state of uniaxial stress, then we can eliminate p from (24) and (25) to obtain
s11 ¼ 2j1 l2 l1 þ 2j4 l2
(26)
and
. sðnÞ ¼ l2 s11 cos2 Q l2 cos2 Q þ l1 sin2 Q :
Fig. 1. Plot of the function I4+ ¼ l2 ð1 2kÞ þ 2kl1 for different values of the dispersion parameter k2[0,1/3].
(27)
Thus, for uniaxial tension s11 > 0 with l > 1, the stress in the deformed fiber direction is positive for all Q s p/2 even though I4(Q) < 1 for Qc < Q p/2, i.e. these fiber directions are shortened by the deformation. For the specialization of the model from Section 1 to one fiber family, but with the notation of the present section, we have
564
J¼
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
1 k mðI 3Þ þ 1 2 1 2k2
2 exp k2 I4+ 1 1 ;
(28)
and the derivatives of J with respect to I1 and I4 are given by
2 1 ; j1 ¼ m þ kk1 I4+ 1 exp k2 I4+ 1 2 2 ; j4 ¼ ð1 3kÞk1 I4+ 1 exp k2 I4+ 1
(29)
(30)
so that (26) leads to
s11 ¼ m l2 l1 þ 2k1 ð1 2kÞl2 kl1 2 ; I4+ 1 exp k2 I4+ 1
1 2
Z rðQÞWðI4 ðQÞÞsinQdQ;
where the definition of S is now changed so that it corresponds to an appropriate range of values of Q2[0, p], which is either [0, p] itself if all fibers are counted or [0, Qc) if fibers for which I4(Q) 1 are excluded from the energy, and hence from the stress. Note that this formulation is a special case of a formulation which will be given in Section 4 for a general deformation, and this has some similarities with the structural model of Lanir (1983) and the model of F-actin networks of Holzapfel et al. (2014). The energy of the system is assumed to be the uncoupled combination of an isotropic matrix energy depending only on I1, denoted Wiso(I1), and the fiber energy (32):
(31)
with I4+ given by (16)2. The curves in Fig. 2 show representative plots of the dimensionless s11, defined by s11 ¼ s11/m, versus l for k ¼ 0, 0.1, 1/3 with illustrative material constants k1 ¼ k1/m ¼ 2 and k2 ¼ 1.
(32)
S
Wiso ðI1 Þ þ
1 2
Z rðQÞWðI4 ðQÞÞsinQdQ:
(33)
S
The associated Cauchy stress is then given by
s ¼ pI þ 2W 0 iso ðI1 Þb þ
Z
rðQÞW 0 ðI4 ðQÞÞn5nsinQdQ;
(34)
S
where n has components lcosQ in the direction of M and l1/2 sinQ laterally, and a prime indicates differentiation with respect to the considered argument, I1 or I4. By using symmetry this can be expressed in the form
s ¼ pI þ ½2W 0 iso ðI1 Þ þ bb þ ða bÞm5m;
(35)
where
Z a¼ Fig. 2. Dimensionless plots of s11 ¼ s11/m versus l from (31) for the values 0, 0.1, 1/3 of k.
The angle Qc is not included explicitly in the term I4+ and therefore not in the strain-energy function. Thus, in order to account for the individual fibers within the distribution that are extended and to exclude those that are not, a modification of the model is required. At this point we emphasize that the correct tensionecompression switch is I4 > 1 and not I4+ > 1. With the current model it is not possible to exclude compressed fibers within a distribution for which the mean fiber direction is under extension. However, a possible modification of the model that allows for such an exclusion could be considered by including dependence on k and the critical angle Qc within the parameter k1, for example, bearing in mind that Qc depends on the state of deformation. We do not follow this approach, but instead we consider an alternative modification, based on the individual fibers within the distribution. This is discussed in detail in the following subsection for the same simple deformation, then for a general plane strain deformation in Section 3.2, and finally generalized for an arbitrary deformation in Section 4. 2.2. Distributed fiber model 2.2.1. An orientation dependent fiber distribution model We now use the notation I4(Q) defined in (19)1 and associate with each fiber in the distribution an energy function, say W(I4(Q)), per unit reference volume. The total fiber energy density with the orientation density r(Q) satisfying the normalization defined in (5) is then given by
S
1 b¼ 2
rðQÞW 0 ðI4 ðQÞÞcos2 QsinQdQ; Z
(36) rðQÞW 0 ðI4 ðQÞÞsin3 QdQ:
S
For uniaxial stress in the mean fiber direction we obtain
s11 ¼ 2W 0 iso l2 l1 þ al2 bl1
(37)
and
h i . sðnÞ ¼ ð2W 0 iso þ aÞl3 ð2W 0 iso þ bÞ lcos2 Q l2 cos2 Q þ l1 sin2 Q ; (38) which should be compared with (26) and (27), respectively. 2.2.2. The limiting case k ¼ 0 Here it is interesting to consider the special case in which all the fibers are in the same direction M, so there is no dispersion. Recall that (4) and (5) must hold in the limit k / 0, so that r(Q) effectively becomes a delta function. First we note that from the definition of I4(Q) in (19) we obtain
I4 ðQÞ I4 ð0Þ ¼ l2 l1 sin2 Q:
(39)
Let us now consider a ‘narrow’ distribution of fibers so that Q is ‘small’ and we can apply the Taylor expansion to W0 (I4(Q)) to obtain the approximation
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
h i 00 W 0 ðI4 ðQÞÞ ¼ W 0 ðI4 ð0ÞÞ þ I4 ðQÞ I4 ð0Þ W ðI4 ð0ÞÞ þ …
(40)
Substitution of this into the expressions for a and b in (36) and making use of equations (4), (5) and (39), we obtain the approximations
00 a ¼ 2ð1 2kÞW 0 ðI4 ð0ÞÞ 4 l2 l1 k A W ðI4 ð0ÞÞ þ …; (41) 00 b ¼ 2kW 0 ðI4 ð0ÞÞ 2 l2 l1 AW ðI4 ð0ÞÞ þ …;
565
function. Note that in the limit k / 1/3, b / 0 the limiting value of r(Q) is 1. Now we plot s11 in dimensionless form s11 ¼ s11/mf against l in Fig. 3 for the two cases (i) S ¼ (0, p) and (ii) S ¼ (0, Qc). We take the material parameters kf1 ¼ kf1/mf ¼ 2 and kf2 ¼ 1 for k ¼ 0, 0.1, 1/3, respectively. The integration over S ¼ [0, p] gives the solid curves, and the integration over S ¼ [0, Qc) gives the dashed curves. As can be seen in the figure, the model using Qc gives a weaker response. This is because the excluded fibers no longer contribute to the strain energy in the fibers. For the case k ¼ 0 the curves coincide because there is no dispersion (see the analysis in Section 2.2.2).
(42)
where
A¼
1 4
Zp
rðQÞsin5 QdQ:
(43)
0
Since A and k are non-negative then, with reference to (4), it is clear that A k and hence in the limit k / 0 we obtain A / 0 (and further terms in the Taylor expansion will similarly disappear). Thus, we obtain, in the limit, a ¼ 2W0 (I4(0)) and b ¼ 0, and the Cauchy stress tensor reduces to
s ¼ pI þ 2W 0 iso ðI1 Þb þ 2W 0 ðI4 Þm5m;
(44) 3. Plane strain problems
which is a particular example of a transversely isotropic model with just the two invariants I1 and I4 ¼ I4(0) given by (23). 2.2.3. A specific energy function Next we make the model more specific. We take
1 m ðI 3Þ; 2 f 1
Wiso ðI1 Þ ¼
(45)
where mf is a neo-Hookean constant, and
WðI4 Þ ¼
h i o kf1 n exp kf2 ðI4 1Þ2 1 ; 2kf2
a ¼ kf1
(46)
i h rðQÞðI4 ðQÞ 1Þexp kf2 ðI4 ðQÞ 1Þ2 cos2 QsinQdQ;
S
(47) b¼
1 k 2 f1
Z
i h rðQÞðI4 ðQÞ 1Þexp kf 2 ðI4 ðQÞ 1Þ2 sin3 QdQ:
S
(48) For specific calculations we take the orientation density function to be the standard p-periodic von Mises distribution of the form
rffiffiffiffiffiffiffi b 1 pffiffiffiffiffiffi exp 2bcos2 Q ; rðQÞ ¼ 4 2p erfi 2b
We recall that for a general incompressible transversely isotropic material the Cauchy stress is given by equation (22). When the deformation is restricted to plane strain we have I2 ¼ I1, as is well known, and also, with M lying in the considered plane, I5 ¼ (I1 1)I4 1, as shown in Merodio and Ogden (2002). Thus, in this case the strain energy is simply a function of I1 and I4 and the (planar) Cauchy stress tensor can be written as
s ¼ pI þ 2j1 b þ 2j4 m5m;
where the material constants kf1 and kf2 are taken to be the same for each fiber in the distribution. Note, in particular, that mf, kf1 and kf2 are not in general the same as the m, k1 and k2 in (28). Then, with I4(Q) as defined in (19), we obtain from (36)
Z
Fig. 3. Dimensionless plots of s11 ¼ s11/m versus l from (37) for the values 0, 0.1, 1/3 of k.
(49)
as we have used in Gasser et al. (2006), where b > 0 is the concentration parameter and erfi(x) ¼ i erf(ix) is the imaginary error
(50)
as in (23), except that all the tensors are now two dimensional (in the considered plane). We now consider a planar fiber dispersion with r(Q) ¼ r(Q) and Q2S, which is now the range [p/2, p/2], as distinct from [0, p] used in (5). The normalization condition is now
1 p
Z rðQÞdQ ¼ 1;
(51)
S
which is the counterpart for plane strain of the three-dimensional normalization (5). The corresponding planar generalized structure tensor, the counterpart of (1), is defined as (Ogden, 2009)
H¼
1 p
Z rðQÞN5NdQ;
(52)
S
so that, with planar C,
I4+ ¼ trðCHÞ ¼
1 p
Z rðQÞI4 ðQÞdQ;
(53)
S
where I4(Q) ¼ N,(CN). Relative to Cartesian axes E1 ¼ M and E2 we take N to have the form N ¼ cosQE1þsinQE2. The components of H in (52) are then
566
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
H11 ¼ 1 k;
H12 ¼ 0;
H22 ¼ k ¼
1 p
Z
rQsin2 QdQ;
(54)
S
which defines the dispersion parameter k for this situation. With the help of the planar identity tensor I ¼ E15E1 þ E25E2 and (54), we obtain.
H ¼ kI þ ð1 2kÞM5M;
(55)
c J I1 ; I4+ ¼ JðI1 ; I4 Þ:
b 4 ¼ n I + 1 ¼ nkc2 : j 4
b 1 ¼ 1 m; j 2
(64)
Substitution into (62) gives
(65)
s11 ¼ p þ m þ 2nk2 c2 1 þ c2 þ 2nkð1 2kÞc2 ;
(66)
s22 ¼ p þ m þ 2nk2 c2 :
(67)
(57) 3.2. Distributed fiber model
Then
b 1 þ kj b 4; j1 ¼ j
b 4; j4 ¼ ð1 2kÞ j
(58)
b 1 and j b 4 are the counterparts of j1 and j4 for c where j J. On substitution into (50), we obtain
b 1 þ kj b 4 b þ 2ð1 2kÞ j b 4 m5m: s ¼ pI þ 2 j
½n ¼ ½F½N ¼
Wiso ðI1 Þ þ
rðQÞWðI4 ðQÞÞdQ
(68)
S
I4 ðQÞ ¼ 1 þ csin2Q þ c2 sin2 Q;
(69)
where
n1 cosQ þ c sinQ ; ¼ n2 sinQ
S¼
where c is the amount of shear, [N] has components [cosQ, sinQ]T and [M] ¼ [1, 0]T, so that the mean fiber direction is in the direction of shear. Then, the invariants become
I4 ð0Þ ¼ 1;
I4+ ¼ 1 þ kc2 ; (61) the latter coming from (56). Note also that m ¼ M. In terms of the planar principal stretches l1 ¼ l, l2 ¼ l1, c ¼ ll1, with c > 0 corresponding to l > 1. Then, from (59),
b 1 þ kj b 4 c: s12 ¼ 2 j
Z
s ¼ pI þ 2W 0 iso ðI1 Þb þ 2FSFT ;
(60)
I1 ¼ 2 þ c2 ;
1 p
per unit area of the considered plane, where S represents the range of angles in [p/2, p/2] for which I4(Q) > 1. The Cauchy stress tensor is
Now let us consider a simple example of a deformation, that of simple shear with
c ; 1
Similarly to the discussion in Section 2.2, suppose now that the energy of the distribution embedded in an isotropic matrix characterized by Wiso is
(59)
3.1. Application to simple shear
1 ½F ¼ 0
with I4+ instead of I4, so that, with the help of (61)4,
The corresponding normal stresses are
(63)
(56)
where I4 ¼ I4(0), noting that I1 is now defined based on the twodimensional b (plane strain). Note, in particular, the factor 1 2k in the above compared with 1 3k in (15). If M is the mean fiber direction associated with a distribution then we may consider, instead of J(I1, I4),
2 1 1 mðI1 3Þ þ n I4+ 1 ; 2 2
s12 ¼ m þ 2nk2 c2 c:
the counterpart of (3), from which we obtain the planar I4+ as
I4+ ¼ kI1 þ ð1 2kÞI4 ;
J¼
(62)
In this example, the critical angle Qc for which I4(Q) ¼ 1 is given by tanQc ¼ 2/c. If we take c > 0 and sinQc < 0, for example, then Qc2(p/2, 0), and the fibers are extended if either p/2 < Q < Qc or 0 < Q < p/2 in which case the definition of S should be modified so that only the contributions of the extended fibers are included. Example. Let us now adopt a neo-Hookean material with a standard reinforcing model (see, for example, Qiu and Pence, 1997) of the form
1 p
Z
rðQÞW 0 ðI4 ðQÞÞN5NdQ:
(70)
S
The latter expression must have the form
S ¼ aM5M þ bðI M5MÞ þ g M5M⊥ þ M⊥ 5M ;
(71)
where M⊥ is a unit vector perpendicular to M in the considered plane, with the properties M,N ¼ cosQ, M⊥,N ¼ sinQ. The coefficients are given by
a¼
b¼
g¼
1 p 1 p 1 p
Z
rðQÞW 0 ðI4 ðQÞÞcos2 QdQ;
(72)
rðQÞW 0 ðI4 ðQÞÞsin2 QdQ;
(73)
rðQÞW 0 ðI4 ðQÞÞsinQcosQdQ;
(74)
S
Z S
Z S
which, in principle, can be calculated for any W and specified plane deformation. Note that, since r(Q) is an even function of Q but, in general, I4(Q) is not, then gs0 in general.
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
Thus, the Cauchy stress tensor (69) takes on the form
M ¼ coscE1 þ sincE2 ;
567
M0 ¼ coscE1 sincE2 ;
(84)
0
s ¼ pI þ 2ðW iso ðI1 Þ þ bÞb þ 2ða bÞm5m þ 2g m5m⊥ þ m⊥ 5m ;
(75)
where m⊥ ¼ FM⊥. Note that the final term here with the coefficient g does not have a counterpart in (59). We now introduce a specific constitutive model for Wiso, given by (45) with the standard reinforcing model for W(I4), so that
W 0 iso ¼
1 m; 2 f
(76)
(85)
(77)
where e1 ¼ FE1 and e2 ¼ FE2, which replace m and m⊥, respectively, in (75), while j1 ¼ W 0 iso ðI1 Þ. This has the same structure as the Cauchy stress (75) in the fiber distribution model with a, b and g replaced by ðj4 þ j6 Þcos2 c, ðj4 þ j6 Þsin2 c and ðj4 j6 Þsinccosc, respectively.
Z
nf p
rðQÞðI4 ðQÞ 1ÞdQ: S
3.2.1. Application to simple shear Now let us specialize to simple shear. Then, it is easy to show from (72)e(74), using (76)2 and (61)2, that 2 +
a ¼ nf c a ;
2 +
b ¼ nf c b ;
+
g ¼ 2nf ca ;
(78)
where
a+ ¼
1 p
Z
rðQÞsin2 Qcos2 QdQ;
b+ ¼
S
1 p
Z
rðQÞsin4 QdQ:
S
(79) +
s12 ¼ 2ðW
iso
4. Three-dimensional transversely isotropic distribution Without loss of generality we may choose Cartesian axes E1, E2, E3 so that the mean fiber direction M ¼ E3. Then the direction N of an arbitrary fiber within a distribution with mean direction M is given by
N ¼ sinQcosFE1 þ sinQsinFE2 þ cosQM;
N,M ¼ cosQ: (86)
+
The values of a and b depend on the orientation density distribution r(Q) and Qc only. Since we are only including extending R Qc R p=2 fibers the integration should be p=2 þ 0 , where we emphasize that Qc is contained in (p/2, 0). Then, from (75) with (76)1 and (78) 0
i
s ¼ pI þ 2 j1 þ ðj4 þ j6 Þsin2 c b þ 2ðj4 þ j6 Þcos2ce1 5e1 þ2ðj4 j6 Þsinccoscðe1 5e2 þ e2 5e1 Þ;
W 0 ðI4 Þ ¼ nf ðI4 1Þ;
and hence, by (51) and (53)2, we obtain, from (72) and (73),
aþb¼
where the (fixed) angle c is introduced to characterize the orientations of the mean fiber directions relative to E1, and let j1, j4 and j6 be the constitutive coefficients. The Cauchy stress tensor is then
þ bÞc þ 2g ¼ mf þ 4nf a
+
+ 3
c þ 2nf b c :
(80)
Let r(Q) be the orientation density, normalized according to (5). In general, I4 ¼ N,(CN) depends on both Q and F so that, using (86), we obtain
I4 ðQ; FÞ ¼ cos2 QM,ðCMÞ þ sin2 Q i h cos2 FE1 ,ðCE1 Þ þ sin2 FE2 ,ðCE2 Þ þ2sinQcosQðcosFE1 þ sinFE2 Þ,ðCMÞ
By comparing this formula with that in (65), we obtain the connections.
m ¼ mf þ 4nf a+ ;
nk2 ¼ nf b+
(81)
between the material parameters in the two models. Thus, in this case, the consequences of the two models are essentially identical, but only as far as the shear response is concerned. The corresponding normal stresses are
þ2sin2 QsinFcosFE1 ,ðCE2 Þ:
Referring to (1), the generalized structure tensor takes on the form
H¼
þðE1 cosF þ E2 sinFÞ5M o þsin2 QsinFcosF E1 5E2 þ E2 5E1 sinQdQdF:
(82) (83)
However, these do not match the normal stresses in (66) and (67). This is not surprising since the structures of the two models are slightly different by virtue of the term in g in (75). However, if one considers two families of distributed fibers in the plane instead of one, with mean fiber directions M and M′ symmetric with respect to the axes E1 and E2, then the structure of the distributed fiber model can be reproduced. In standard notation (see, for example, Holzapfel, 2000) let I4 ¼ M,(CM) and I6 ¼ M′,(CM′) be the squares of the stretches in the mean fiber directions with
Z n 1 rðQÞ cos2 QM5M þ sin2 Q 4p U cos2 FE1 5E1 þ sin2 FE2 5E2 þsinQcosQ½M5ðE1 cosF þ E2 sinFÞ
s11 ¼ p þ mf þ 2nf b+ c2 1 þ c2 þ 2nf c2 a+ b+ þ 8nf c2 a+ ;
s22 ¼ p þ mf þ 2nf b+ c2 :
(87)
(88)
By carrying out the integration over F2[0, 2p] we are left with
H¼
1 2
Z rðQÞcos2 QsinQdQM5M S
1 þ 4
Z
(89) rðQÞsin3 QdQðE1 5E1 þ E2 5E2 Þ:
S
With the properties (4) and (5) and since E15E1þE25E2 ¼ I M5M, we recover (3), the generalized structure tensor for a transversely isotropic distribution.
568
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
4.1. Distributed fiber model
The Cauchy stress tensor (92) can now be written as
Next, let W(I4(Q, F)) be the strain energy of the fiber in the direction N (per unit reference volume). Then, with weighting r(Q), the total strain energy of the fibers is
1 4p
Z rðQÞWðI4 ðQ; FÞÞsinQdQdF;
(90)
U
which extends the energy (32) to three dimensions. Strictly this should only include fibers for which I4 > 1. Let Wiso (I1) be the contribution of the matrix. Then, assuming that there is no coupling, the composite has strain energy
Wiso ðI1 Þ þ
1 4p
rðQÞWðI4 ðQ; FÞÞsinQdQdF;
(101)
where m ¼ FM, e1 ¼ FE1 and e2 ¼ FE2. Now suppose that the properties of the fibers are given by the standard reinforcing model (76)2 and note that
a þ b þ b0 ¼
1 4p
(91)
Z
rðQÞW 0 ðI4 ðQ; FÞÞsinQdQdF
U
Z
(102) rðQÞðI4 ðQ; FÞ 1ÞsinQdQdF:
U 0
U
For the coefficients a, b, b , g, g0 and d we obtain, after some detailed calculations,
a ¼ Anf ð3I4 ð0; 0Þ I1 Þ þ nf BðI4 ð0; 0Þ 1Þ þ knf ðI1 5I4 ð0; 0Þ þ 2Þ; (103)
(92)
where
1 4p
þdðe1 5e2 þ e2 5e1 Þ;
Z
s ¼ pI þ 2W 0 iso ðI1 Þb þ 2FSFT ;
S¼
þgðm5e1 þ e1 5mÞ þ g0 ðm5e2 þ e2 5mÞ
n ¼ f 4p
generalizing (33), and the corresponding Cauchy stress tensor is
Z
s ¼ pI þ 2W 0 iso ðI1 Þb þ 2½be1 5e1 þ b0 e2 5e2 þ am5m
b ¼ knf ðI4 ð0; 0Þ 1Þ þ rðQÞW 0 ðI4 ðQ; FÞÞN5NsinQdQdF;
Anf ð3e1 ,e1 þ e2 ,e2 4I4 ð0; 0ÞÞ; 4
(104)
(93)
U
generalizing the definition (70). The latter expression can be written in the form
S ¼ aM5M þ bE1 5E1 þ b0 E2 5E2 þ gðM5E1 þ E1 5MÞ þg0 ðM5E2 þ E2 5MÞ þ dðE1 5E2 þ E2 5E1 Þ;
b0 ¼ knf ðI4 ð0; 0Þ 1Þ þ
Anf ðe1 ,e1 þ 3e2 ,e2 4I4 ð0; 0ÞÞ; 4 (105)
g ¼ Cnf m,e1 ;
g0 ¼ Cnf m,e2 ;
d ¼ nf De1 ,e2 ;
(106)
(94) the choice of E1 and E2 perpendicular to M such that E1,E2 ¼ 0 being arbitrary. Then,
a¼
1 4p
1 b¼ 4p
b0 ¼
g¼
g0 ¼
1 4p 1 4p 1 4p
1 d¼ 4p
Z
where
A¼ rðQÞW 0 ðI4 ðQ; FÞÞcos2 QsinQdQdF;
(95)
1 4
U
Z
rðQÞW 0 ðI4 ðQ; FÞÞsin Qcos2 FdQdF; 3
(96)
C¼
1 2
U
Z
rðQÞW 0 ðI4 ðQ; FÞÞsin3 Qsin2 FdQdF;
(97) k¼
rðQÞW 0 ðI4 ðQ; FÞÞcosQsinQcosFdQdF;
(98)
B¼
S
Z
rðQÞsin2 Qcos2 QdQ;
1 2
Z rðQÞsinQdQ;
(107)
S
D¼
S
1 8
Z
rðQÞsin4 QdQ;
S
rðQÞW 0 ðI4 ðQ; FÞÞcosQsinQsinFdQdF;
(99)
1 4
Z
rðQÞsin3 QdQ:
(109)
S 0
0
The coefficients a, b, b , g, g and d can in principle be calculated for any W.
U
Z
rðQÞsin5 QdQ;
(108)
U
Z
Z
4.2. Identification of the critical angles
U
Z
The goal is now to find the range of angles for which 0
2
rðQÞW ðI4 ðQ; FÞÞsin QsinFcosFdQdF; U
(100)
and we note that cosFE1þsinFE2 is an arbitrary vector normal to M, so that interchanging cosF and sinF makes no essential difference. These are the generalizations of the a, b and g introduced in Section 3.2.
I4 ¼ n,n ¼ N,ðCNÞ > 1
(110)
for a fiber distribution with mean fiber direction M. We assume a transversely isotropic distribution. With N,M ¼ cosQ we want to find the critical angles Q and F such that
G.A. Holzapfel, R.W. Ogden / European Journal of Mechanics A/Solids 49 (2015) 561e569
N,ðCNÞ ¼ 1:
(111)
Explicitly, with
N ¼ sinQðcosFE1 þ sinFE2 Þ þ cosQE3 ;
(112)
(111) gives
I4 ðNÞ ¼ N,ðCNÞ ¼ sin2 Q C11 cos2 F þ 2C12 sinFcosF þ C22 sin2 F þC33 cos2 Q þ 2sinQcosQðC13 cosF þ C23 sinFÞ ¼ 1; (113) where Cij are the components of C. For each Q2[0, p] this defines the critical values of F beyond which N,(CN) < 1. In the special case in which E1, E2, E3 are eigenvectors of C, Cij ¼ 0 for isj and Cii ¼ l2i , i ¼ 1, 2, 3, and the incompressibility l1l2l3 ¼ 1, and eq. (113) reduces to
h i tan2 Q l22 l21 sin2 F þ l21 1 ¼ 1 l23 ;
cos2F ¼
l21 l22
:
excluded. This model has been developed both for plane strain and for general three-dimensional deformations, and yields explicit expressions for the Cauchy stress deformation relation. These are compared with those for the models based on the continuum approach described in Gasser et al. (2006). In particular, for simple shear we have established connections between the material parameters in the two models for the shear stresses to coincide. For general three-dimensional deformations we have provided an algorithm for obtaining the critical angles that can be used within a general stress computation. The stress relation and the (explicit) expressions for the critical angles provide a basis for finite element implementation of the model.
Acknowledgment The authors are grateful to Dr. Kewei Li, Graz University of Technology, Austria, for valuable discussions.
(114) References
from which F is easily determined for any given Q in the form
2 1 l23 cot2 Q þ 2 l21 þ l22
569
(115)
5. Concluding remarks In this study we have analyzed the effect of using an inappropriate tensionecompression switch for excluding the contribution of fibers under compression within a soft fibrous composite solid. In particular, for a simple tension illustration we have shown that the stiffness of the composite is in general significantly overestimated by not excluding fibers under compression, and we have computed the critical angle of the individual fibers relative to the fiber direction that corresponds to the switch. The correct tensionecompression switch is I4 > 1 and not I4+ > 1. However, for the current model, shortened fibers within a distribution with mean fiber direction under extension cannot be excluded. We have proposed a possible modification of the model within which such an exclusion can be considered. In addition, we have introduced a distributed fiber model with a weighted energy function that also allows shortened fibers to be
mes Simulia Corp, 2013. Abaqus 6.13 Analysis User's Guide. In: SecDassault Syste tion 22.5.3: Anisotropic Hyperelastic Behavior. Flory, P.J., 1961. Thermodynamic relations for high elastic materials. Trans. Farad. Soc. 57, 829e838. Gasser, T.C., Ogden, R.W., Holzapfel, G.A., 2006. Hyperelastic modelling of arterial layers with distributed collagen fiber orientations. J. R. Soc. Interface 3, 15e35. Holzapfel, G.A., 2000. Nonlinear Solid Mechanics: a Continuum Approach for Engineering. John Wiley & Sons, Chichester. Holzapfel, G.A., Ogden, R.W., 2010. Constitutive modelling of arteries. Proc. R. Soc. Lond. A 466, 1551e1597. Holzapfel, G.A., Unterberger, M.J., Ogden, R.W., 2014. An affine continuum mechanical model for cross-linked F-actin networks with compliant linker proteins. J. Mech. Behav. Biomed. Mater. 38, 78e90. Lanir, Y., 1983. Constitutive equations for fibrous connective tissues. J. Biomech. 16, 1e12. Merodio, J., Ogden, R.W., 2002. Material instabilities in fiber-reinforced nonlinearly elastic solids under plane deformation. Arch. Mech. 54, 525e552. Ogden, R.W., 1978. Nearly isochoric elastic deformations: application to rubberlike solids. J. Mech. Phys. Solids 26, 37e57. Ogden, R.W., 2009. Anisotropy and nonlinear elasticity in arterial wall mechanics. In: Holzapfel, G.A., Ogden, R.W. (Eds.), Lecture Notes, CISM Course on the Biomechanical Modelling at the Molecular, Cellular and Tissue Levels, CISM Courses and Lectures Series, vol. 508. Springer, Wien, pp. 179e258. Qiu, G.Y., Pence, T.J., 1997. Remarks on the behavior of simple directionally reinforced incompressible nonlinearly elastic solids. J. Elast. 49, 1e30. Simo, J.C., Taylor, R.L., 1991. Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms. Comput. Meth. Appl. Mech. Eng. 85, 273e310. Vergori, L., Destrade, M., McGarry, P., Ogden, R.W., 2013. On anisotropic elasticity and questions concerning its Finite Element implementation. Comp. Mech. 52, 1185e1197.