IX
Finite Chapter Deformation
1. Introduction
One of the conclusions that comes from the study of uni-axial confined compression is that the flow of fluid, and its interaction with the solid matrix, can have a profound effect on the biorheological properties of soft tissue. For example, in Chapt. V we saw that in stress relaxation the local compressive strains that are predicted using a linearly elastic solid can exceed 100% even for strain-rates that are well within physiological limits. Large strains actually occur in vivo. To illustrate using articular cartilage, a 160 lb person standing upright produces a stress of about 1 MPa on the articular surface in each knee (assuming a contact area of 10 cm2 ). Since the elastic moduli are on the order of 1 MPa then the compressive strains, based on the linear theory, are around 100%. This situation is even more severe when the knee is in motion since the loads are greater (they can be up to 18 times greater) and the diffusive drag acts to increase the strains near the articular surface. The necessity of a large strain theory is even more apparent for tissues such as the skin, the vascular wall, and the lung. Therefore, to have a comprehensive description of the deformational behavior of soft tissue it is necessary to develop the model so it includes finite deformation. We will do
IX.1
7/9/95
this within the confines of simple mixture theory so the equations of motion are given in the summary to Chapt. II. One of the central problems in finite deformation of a biphasic mixture is to find the appropriate constitutive laws for the elastic stress and the permeability. As we saw in Chapt. V, the stress-strain relation for articular cartilage in uni-axial confined compression is nonlinear for compressive strain greater than about 20% (see Fig. V.3). One of the questions addressed in this chapter is what function should be used to describe this nonlinear behavior? There are many possible choices and examples are E = a + b3 , and E = e - 1 .
Both of these functions have been used in biomechanics. The first has the advantage that it is easy to manipulate since it is a polynomial, while the second can describe the large tensile strains observed in some tissues and has the feature that the modulus, , can be determined by measuring the slope of the stress-strain curve using logarithmic coordinates. However, along with their simplicity come certain drawbacks. For example, they both predict a finite stress when the tissue is compressed to the point that = -100%. In situations where the strains are controlled this may not be a problem because the strains may never reach levels where this is an issue. This was the case when the linear stress-strain law was studied in Chapt. V. However, if one wants a description that is capable of handling a wide range of motions such limitations should be minimized. There are other similar restrictions on the constitutive laws and these are discussed in more detail in this chapter. Our approach to develop the large strain theory will parallel the discussion in Chapt. III. So, rather than starting with the stress-strain relationship we will begin with the formulation of the free energy function. It is very simple to determine this function for the above two functions. For example, in the case of uni-axial motion as described in Section 9.6, one finds that E = s . One finds in this case that for the first function 1 1 s = 2 a 2 + 4 b4 and for the second s = (e - )/ . However, the limitations of these expressions have already been stated. Consequently, in the first part of the chapter we will develop the constitutive modeling of an elastic material in the case of a multidimensional strain field. After this we will study the uniaxial confined compression experiment. The constitutive law we will use, in the case of uniaxial motion, is the following
IX.2
7/9/95
s
e = a 2 ,
2
where = 1 + is the stretch. From this we will show that
2 2 - 1 1 E = 2 H M 2+1 e ( - 1) .
The reasons for choosing this particular function will be discussed later but note that it becomes unbounded if -1 or if . Also, the theory now contains two material parameters (H M, ) associated with the elastic stress in the solid matrix rather than one as in the small strain formulation (H A). The subject of finite deformation continuum mechanics is quite broad. This chapter is meant only as an introduction to the constitutive modeling and the simpler loading configurations that arise in tissue mechanics. There is also a brief review of kinematics and some of the other facets of finite deformation are given in the exercises. For those who would like a more extensive presentation of the subject, the books by Marsden and Hughes (1983), Truesdell and Noll (1992), and Ogden (1986) are worth considering.
2. Constitutive Relations for the Solid Phase
2.1 Introduction to Hyperelasticity We need to specify how the solid stress depends on the deformation. This is done, as before, by specifying how the free energy function of the solid phase depends on the deformation. This function generally depends on many variables including such quantities as the strain, the rate of strain, the strain gradient, the volume fractions, etc. Just which of these should be included in the constitutive law depends on the material properties of the tissue as well as the types of loading conditions to be considered. We will restrict our attention to materials and loadings where the free energy function for the solid depends only on the strain. To be precise, we will assume that the solid is hyperelastic (Truesdell and Noll, 1992), i.e., s = s (C) , (FD.1)
IX.3
7/9/95
where C is the right Cauchy-Green deformation tensor for the solid (the definition of C is given shortly). The choice of C as the measure of strain is common in hyperelasticity but it is by no means the only one that can be made and some others are discussed in the exercises. Whichever one is used, because s is dependent solely on a measure of the elastic strain it is sometimes referred to as the strain energy function or the elastic potential. It is assumed here that s is a smooth symmetric function of the components of C . By symmetric we mean that s (CT) = s (C) . (FD.2)
This does not restrict the choice of the energy function because C is symmetric. Also, the free energy function is required to satisfy the Principle of Material FrameIndifference. However, it is not hard to show that the constitutive law in (FD.1) satisfies this requirement (see Ex. 1). As a final comment, we have implicitly assumed that the solid is homogeneous since s is taken to be independent of the spatial coordinates. The ramifications of assuming the solid is hyperelastic will be discussed shortly but first it is best to briefly review kinematics and some of the properties of tensor functions (the motivation for this should be apparent from Eq. (FD.1)). 2.1.1 Review of Kinematics and Tensor Analysis To introduce the kinematics of mixtures suppose we have a fluid-solid mixture whose current configuration, at time t , is the region (see Fig. FD.0). The reference, f s or initial, configuration of each phase is designated as 0 and 0 , respectively. Using Lagrangian coordinates, let X be the position of a particle of the th phase in its reference configuration ( = s, f ). The current configuration can then be represented as
x = x(X, t) , for X 0 .
(FD.0a)
In what follows we will use the same reference axes for both x and X . In this case the current position of the particle can be written as x = X + U(X, t) , (FD.0b)
where U is the displacement vector of the th phase. With this the velocity V (X, t) of the particle is related to the displacement in the usual way, that is,
IX.4
7/9/95
f o
uf
Xf X3, x3 Xs X2, x2
s o
us
X1, x1
Fig. FD.0 Kinematics of a fluid-solid mixture. V = tU( X, t) . X It is assumed that the motion is smooth and invertiable. One consequence of this is that (FD.0a) can be inverted to yield X = X(x, t) , for x . (FD.0a)
In terms of a spatial, or Eulerian, description we then have u(x, t) = U(X(x, t), t) and v(x, t) = V(X(x, t), t) . Moreover, the material derivative relative to the th phase is Dt = t + v. . To simplify the presentation the superscript designating the constituent will be omitted in the remainder of this review since the discussion applies to either fluid or solid. So, letting X represent the initial position of a material point then the motion of that point is described as x = x(X, t) . The deformation gradient associated with this motion is the second order tensor F defined as F = Xx , or equivalently, in Cartesian component form xi Fij = X . j (FD.3b) (FD.3a)
IX.5
7/9/95
The right Cauchy-Green deformation tensor is defined as C = F TF . (FD.4)
This tensor has a lot of nice properties which will be useful later. For example, it is symmetric and it has only positive eigenvalues (in other words, C is a positive definite). The deformation gradient, on the other hand, is not necessarily symmetric and can have complex eigenvalues. It is possible to write C in terms of the derivative of the displacement. This is done by recalling that the displacement U(X, t) is defined through the following equation x(X, t) = X + U(X, t) , and so, from Eq. (FD.3), F = I + XU . (FD.5)
We now substitute this result into Eq. (FD.4) and multiply everything out to obtain the following complicated expression C = I + XU + ( XU)T + (XU)T XU , or equivalently, in Cartesian component form Cij = ij + X Ui + X Uj + X Ul X Ul . j i i j l (FD.6a)
(FD.6b)
As is apparent from the above expressions, C is a nonlinear function of the displacement gradient (it is a quadratic function). This is one source of nonlinearity in the problems to come but it is by no means the only one. We also need some information on the derivative of a function of a tensor. This is because we will soon substitute (FD.1) into the entropy inequality and this means we will need to be able to deal with the material derivative of s . If one takes the material derivative of a component of the deformation gradient (FD.5) and then writes the result in tensor form the result is DtF = LF , (FD.7)
where D t = t + v. is the material derivative and the second order tensor L is the velocity gradient. The latter is defined as IX.6 7/9/95
L = v , or equivalently, in Cartesian component form vi Lij = xj .
(FD.8)
Since the derivatives in the definition of L are with respect to the spatial coordinates the operator in Eq. (FD.8) is sometimes referred to as the spatial gradient to distinguish it from the material gradient that appears in Eq. (FD.3). From this result and Eq. (FD.4) it follows that DtC = 2F TDF , where D = 2 ( L + L T ) is the rate of deformation tensor. If = (C) then using the chain rule and the properties of the trace function (see Ex. II.10) D t = Cij D tCij i,j
1
= tr C D tC T = tr 2F C FTD . The second order tensor C is the called the gradient of with respect to C and the i,j th entry of this tensor in Cartesian coordinates is . Cij Note that the transpose in (FD.9) is not necessary because of our assumption that the strain energy function is a symmetric function of C . (FD.9)
IX.7
7/9/95
2.2 Anisotropic Solid It is now possible to relate the solid stress to the free energy function in Eq. (FD.1). The procedure is exactly the same as was used in Chapt. 3 to derive Eq. (III.6) from the reduced entropy inequality (III.5). Substituting Eq. (FD.9) into this inequality, f and using the fact that D tf = 0 (see Eq. (FD.28) below), one obtains s T tr s + s pI - 2 s F FT Ds + (f + fpI)Df C + (vs - vf) . (ps - s ) 0 . We assume s , f and s are independent of Ds . The Principle of Dissipation states that this inequality must hold for all values of Ds . The only way this can happen is that the following equation is satisfied s s = - spI + 2 sF C FT . (FD.11) (FD.10)
To use this expression we need a constitutive law that specifies the dependence of s on C and this is usually done based on experimental measurements of the stress-strain relationship. The generality contained in Eq. (FD.11) is not pursued here as we will restrict our attention to isotropic tissues. The assumption of isotropy serves as a good starting point for finite deformation and once it is understood it is then easier to study more complicated structural models. As a final comment on Eq. (FD.11), because we have assumed there is no moment of momentum supply, s must be symmetric (see Eq. (II.9)). That requires the gradient of s in Eq. (FD.11) to be symmetric. This condition is satisfied, however, because s is a symmetric function of C . As mentioned earlier other measures of strain can be used in hyperelasticity and some of the more commonly used strain tensors are listed in Table VIII.1. The corresponding expressions for the solid stress are given in the exercises below. An interesting question arises when considering these other tensors as to how it might be possible to generalize the constitutive law given in Eq. (FD.1b). One way is to assume that the energy function depends on F , rather than C . However, as shown in Ex. 1, the requirement that s is objective forces this function to depend on F through C and therefore it's just as well to start with C in the formulation. One might also try using the left Cauchy-Green (or, Finger) deformation tensor B , defined as
IX.8
7/9/95
B = FFT .
(FD.12)
Interesting enough, if one makes this assumption then the solid must be isotropic for s to be objective (Ex. 11). As a final comment, in nonlinear elasticity it is usually easier to solve a problem if one uses material rather than spatial coordinates and in Ex. 3 the material description for the solid is investigated. Exercises: 1. a) Letting x* = Qx + b , where Q is a rotation tensor (see Ex. III.1) then the deformation gradient in the x* system satisfies F* = QF . From this show that C is invariant (i.e., C * = C ) but B is not. However, show B is objective and C is not. b) Show the free energy function, when using the constitutive law in Eq. (FD.1), is objective. c) Suppose one assumes s = - s(F) , where F is the deformation gradient for the - solid phase. For objectivity it is required that s(F*) = - s(F) for all F and Q . - s depends on F through C , i.e., s = s (C) . (Hint: From this deduce that use the polar decomposition F = RU , where R is a rotation tensor and U is the right stretch tensor. Also, note C = U2 .) 2. The Cauchy deformation tensor is c = B-1 and the Eulerian (or, Almansi) strain tensor is e = 2 (I - c) .
1
(FD.13)
a) Show that D te = D - eL - LTe (Hint: first differentiate the equation cB = I and then use Eqs. (FD.7) and (FD.12)). b) Assume s = - s(e) , where - s is a symmetric function of the Eulerian strain tensor e for the solid phase. Show that - s s = - spI + s(I - 2e) e . - s Is it also necessary to assume that e e is symmetric? c) If u(x, t) is the displacement in spatial coordinates show that e = 2 [u + (u) T - (u)T u] .
1
IX.9
7/9/95
3. a) Show that in material coordinates the momentum equation for the solid phase, in the absence of external body forces, is 0 t Us = XTs + s , where Ts = det(F)s (FT) -1 , and s = det(F)s . (FD.14c) (FD.14b)
s 2
(FD.14a)
Here Ts is the first Piola-Kirchhoff (or, Lagrangian, or nominal, or Piola, or Piola-Kirchhoff) stress tensor and it gives the contact force per unit area in the reference configuration. It is, in general, not symmetric but S = F -1Ts is (the s latter is the second Piola-Kirchhoff stress tensor). Also, 0 is the initial value of the mass density of the solid. b) Show that the solidity equation (II.3s) in material coordinates is 0 s (X,t) = det(F) where 0 = s (X, 0) . c) Show that 1 -1 s tr[ ( s + s pI)Ds ] = 2det(F) tr[ (S + 0PC )D t C ] Use this to rewrite the reduced entropy inequality given in (II.17). From this conclude that for an incompressible hyperelastic solid (also see Exercise 9.3.15)
s s S = - 0PC-1 + 20 C .
,
(FD.14d)
d) Show, from (FD.11), that for an incompressible hyperelastic solid
s s Ts = - 0P(FT) -1 + 20 F C .
(FD.14e)
4. a) Show that if one assumes s = - s(E) , where E is the Lagrangian strain tensor (Table VIII.1),
IX.10
7/9/95
- s s = - spI + s F E FT . b) Show that if one assumes s = - s(c) , where c is the Cauchy deformation tensor (Ex. 2), then - s s = - spI + 2 sc c . What symmetry assumption must be used to derive this result? 5. If one assumes s = s (s , C) show that Eq. (FD.11) is still obtained but p is replaced with the scalar s p+ s
s
.
6. Show that the gradient of s with respect to C is not objective but the stress tensor given in Eq. (FD.11) is objective.
2.3 Isotropic Solid Because we will concentrate on isotropic models we need to introduce the principal invariants of a second order tensor. These are defined as I = tr(C) , II = 2 [ I2 - tr(C 2) ] , and III = det(C) . (FD.17)
1
(FD.15) (FD.16)
The reason for calling these quantities invariants is that their values are not affected by orthogonal changes of coordinates (e.g., rotations or translations of the coordinate system). For those unfamiliar with invariants, they come from the characteristic polynomial for C . This polynomial is used to calculate the eigenvalues of C and it is the left hand side of the characteristic equation det(C - I) = 0 . (FD.18)
IX.11
7/9/95
The coefficients of this cubic equation for are the invariants of C , that is, if one multiplies (FD.18) out they obtain 3 - I 2 + II - III = 0 . The number of principal invariants clearly depends on the number of spatial dimensions and in our case there are three. Using these same ideas one can also prove that the maximum number of independent invariants is three, in other words, any other invariant must be a function of the above three. Of particular interest for us is a "representation theorem" which states how the free energy can depend on C . It is not hard to prove that for an isotropic material, s is not an arbitrary function of the individual components of C but it only depends on its invariants. This result simply says that if the material is isotropic then the energy density must depend on C only through those quantities that are directionally independent (i.e., the invariants). Therefore, for an isotropic material we can rewrite Eq. (FD.1) as follows s = s (I, II, III) . (FD.19)
To complete the description for an isotropic solid we need to transform the derivative of s in Eq. (FD.11) into derivatives with respect to the invariants. To do this we need to determine the derivatives of the invariants. The first principal invariant (FD.15) is simply the sum of the diagonal entries of C and is independent of the other components. Thus, I = ij , Cij or equivalently, I = I . C Similarly, one can show that II = I I- C , C and, assuming C is invertible, III = II I - I C + C 2 . C Now, from (FD.19) and the chain rule IX.12 7/9/95 (FD.20c) (FD.20b) (FD.20a)
s I s II s III s = I C + II C + III C . Cij ij ij ij Substituting (FD.11) into this expression and collecting terms it follows that s s s s s s s = I + I II + II III I - II + I III C + III C 2 . C Introducing this into Eq. (FD.11) yields s = - spI + 2 s(a0B + a1B2 + a2B3) , (FD.22) (FD.21)
where the ai 's are functions of the invariants I, II, III and are given in Eq. (FD.42). It should be recalled that the invariants of C and B are the same (see Ex. 9). The fact that B is symmetric is useful because it allows us to find an orthogonal coordinate system, consisting of eigenvectors of B , in which B is diagonal. Not surprisingly this is called the Principal Axis Theorem and using this basis of eigenvectors we get that
B =
2 1 0 0
0 2 0
2
0 2 3
0
,
(FD.23)
where the i's are non-negative numbers known as the principal stretches. It is a simple matter, by introducing (FD.23) into the definition of the principal invariants, to show that I = 1 + 2 + 3 , II = 1 2 + 1 3 + 2 3 , and III = 1 2 3 .
2 2 2 2 2 2 2 2 2 2 2 2
(FD.24a) (FD.24b)
(FD.24c)
Moreover, with this coordinate system the stress tensor is diagonal and has the form s
=
1 - spI + 0 0
0 0 2 0 0 3
,
(FD.25)
IX.13
7/9/95
where the i 's are principal stresses of the elastic component of the solid phase. One finds from (FD.22) and (FD.23) that the latter are given as i = 2s (a0 + a1 i + a2 i ) i ,
2 4 2
(FD.26)
where the aj 's are the coefficients in Eq. (FD.21) and are given in (FD.42). Since the stretches are obtained rather easily from experiment it is not unusual to see constitutive laws which have the strain energy depending on them rather than on the principal invariants as in Eq. (FD.19). If this assumption is made then it must be remembered that because of the assumption of isotropy it is necessary to require that s does not change if there is a pairwise interchange of the stretches (note that the principal invariants satisfy this requirement). So, s = 1 + 2 satisfies this requirement but s = 1 + 22 does not. Exercises 7. This problem examines what are called homogeneous deformations. This means that x = F0X + x0 where F 0 is a constant tensor and x0 is a constant vector. a) For rigid body motion F 0 is a rotation. In this case show I = II = 3III = 3 and (FD.22) reduces to s = 0I where 0 = -s p + 2s ( Is + 2 II s + III s ) . b) For uni-axial stretching of the form x = (X1, X 2, X3) show that I = 2 + 2 , s II = 1 + 22 , and III = 2. Also, using (FD.22) show s is diagonal with 33 = s - s p + 2s ( Is + 2 II s + III s ) and ii = -s p + 2s (Is + (1 + 2) II s + 2 III s ) for i = 1, 2 . c) For uniform dilatation F0 = I . Show that I = 3 2 , II = 34 , and III = 6 . Also show (FD.22) reduces to s = I where = -s p + 2s 2 ( Is + 22 II s + 4 III s ) . d) For simple shear, where x = X + (X 2, 0, 0) , show that I = II = 3 + 2 , and s III = 1 . Also, using (FD.22), show that the shear stress is 12 = 2s ( Is + II s ) . 8. Show that the principal invariants of the left and right Cauchy-Green deformations tensors are the same but their eigenvectors may be different. 9. Show that I = 2F , F II F = 2I F - 2FC , and II T -1 F = III (F ) .
10. a) Suppose one assumes s = - s(B) . Show that the only way for s to be objective is for the solid to be isotropic (see Ex. 1). To do this recall that the IX.14 7/9/95
- scalar function - s(B) is isotropic if s(BQ) = - s(B) for all orthogonal transformations Q . b) If it is assumed that s = - s(F) , show that isotropy requires that - s depends on F through B . (Hint: Use the polar decomposition F = VR , where R is a rotation tensor and V is the left stretch tensor. Also, note B = V2 .) 11. Using a basis in which C is diagonal show that: i) II = det(C). tr(C-1) , and ii) C-1 = (II I - IC + C 2)/III . 12. From the result in Ex. 3b and Eq. (FD.20c) show that s 1 s -1 C = - 2 C . 13. a) Using the Cayley-Hamilton theorem show that (FD.22) can be rewritten as s = - spI + 2 s(b0I + b1B + b 2B2) , where b0 = III III s , b1 = Is + I II s , b2 = -II s . b) Suppose one assumes s = s (J 1, J 2, J 3) where J 1 = I , J2 = II/III , J 3 = this case, determine the a i's in (FD.22) and the bi 's in part (a). 14. Using (FD.21), and the results from Exercises 11 and 9.2.3, show that S = 0I + 1C + 2C 2 , where 0 = 20a0 - 0II P/III , 1 = 20a1 + 0I P/III , and 2 = 20a2 - 0P/III . A solid is said to be transversely isotropic if it has a reflective-axis of symmetry, h , in its undistorted state (Beatty, 1987). This condition translates into the requirement that s (C) = s(QTCQ) for all rotations that satisfy Qh = h . Show s depends only on I , II , III , J1 = hTCh , and J 2 = hTC 2h - J 12 . Suppose the reflective-axis is the X3-axis. In this case show s depends only on 2 2 I , II , III , C 33 and C31 + C 23 . Show that J 1 = (h )2 where h is the stretch in the h direction. Also, J 2 = ||P(Ch)|| 2 , where P(q) is the orthogonal projection of q = Ch onto the plane with normal h . Suppose the solid is transversely isotropicwith the reflective-axis is the X 3-axis (see Exercise 15). For this material, consider the constitutive law s = 0 e
(I, II, J1, J 2)
III .
In
15.
a) b) c)
16.
III
,
IX.15
7/9/95
where = 1(I - 3) + 2(II - 3) + 3(I - 3) 2 + 4(J 1 - 1) + 5(J 1 - 1) 2 + 6(I - 3)(J 1 - 1) + 7J 2 . So, consists of all linear combinations of the invariants involving powers of the stretches through order four. a) Show E = 2 s
III
[ b1B + b 2B2 + b 3B3 + b 4G + b 5(BG + GB) ]
where G = FM 3FT and M 3 = . b) Show that the requirement of no stress when the material is at rest results in = 1 + 22 - 7 and 7 = - 4 . c) Show that for the material to reduce to the small strain formulation it is necessary that 4 = 0 (II.25).
3. Constitutive Relations for the Fluid Phase and Diffusive Drag
As in Chapt. III, it is assumed that the fluid is inviscid, and so, f = f( f) . (FD.27)
This constitutive law was used earlier when deriving Eq. (FD.11). To explain how, recall f f that f = Tf where T is constant. Thus, by taking the material derivative of f , and assuming the fluid is incompressible (Eq. (II.3b)), one finds that
f Dtf
df f f = D d f t = df f f f D d f T t (FD.28)
= 0 .
With this result the reduced entropy inequality given in (FD.19) is obtained from (III.5). This inequality can be simplified even further. Assuming s is independent of D f then from (FD.19) and the Principle of Dissipation we have that f = - fpI . (FD.29)
IX.16
7/9/95
This is exactly the same expression we derived in Chapt. III. As stated before, there are clear similarities between this expression and the constitutive law for a single phase inviscid fluid. It is also important to note the connection between (FD. 29) and the constitutive law for the solid phase given in (FD. 11). Both contain a pressure that arises from the assumption that each phase is incompressible and this pressure is distributed through each constituent according to its volume fraction. The third constitutive relationship we need to specify is for the momentum supply s (recall f = - s ). As in Chapt. III, this vector can be written as s = p s - 0 , (FD.30)
where 0 depends on the relative velocity vs - vf but is independent of Ds and D f . From the reduced entropy inequality, 0 is restricted by the condition (vs - vf) . 0 0 . (FD.31)
In the present model 0 represents the diffusive drag that is caused when the fluid and solid move relative to one another in the tissue. We will assume that this drag is proportional to the relative velocity of the fluid phase with respect to the solid phase. In other words we will assume that 0 = K(vs - vf) , (FD.32)
where K is the diffusive drag coefficient and from inequality (FD.31) it is non-negative. Also, in this expression the assumption that the mixture is isotropic with respect to fluid flow has been used.
4. Restrictions on the Constitutive Relations
Two constitutive laws are needed to make the model complete, one for the strain energy function and another for the permeability. It is not a simple matter to determine these functions. Clearly it is desirable that the resulting model is capable of reproducing whatever experimental measurements are available, but these are only of limited help. For example, we want to make sure that the model predicts realistic stresses inside the tissue as well as in configurations that have not necessarily been studied experimentally. The same difficulty arises in finite elasticity and much of what has been derived there can be applied here.
IX.17
7/9/95
The first constraint we will consider concerns the equilibrium position, that is, the configuration of the tissue when it is at rest and no external forces are present. This position must be stable in the sense that the tissue will return to it if a load is applied and then removed. One way to help guarantee that this happens is to require that the potential energy has a minimum at the equilibrium position. This idea is easy to visualize if one considers a ball rolling over a surface. In this case the potential energy is determined from the gravitational force and the ball comes to rest at the bottom of a valley (which represents a local minimum in the potential energy). For a multi-dimensional continuum like we are dealing with, however, this idea is harder to express. One can try to approach this problem by requiring that the energy is a convex function of its arguments. Again, this is easy to understand from the rolling ball experiment since in the neighborhood of the minimum the surface should be convex. In linear elasticity the assumption of convexity, as a function of the elements of the infinitesimal strain tensor, leads to the inequalities on the Lam coefficients given in (III.12). As stated in Chapt. III this is equivalent to assuming that the strain energy function is a positive definite function of the components of the infinitesimal strain tensor. We can impose a similar condition on the free energy function for the solid phase. One possibility is the strong ellipticity, or the Legendre-Hadamard, condition which requires that 2s Fij Fkl ijk l > 0 , for all nonzero vectors , . This condition, which is in some sense a requirement that the second derivative is positive, is so cumbersome that it is essentially unusable. Instead we will use a simpler set of constraints. We will require that the greater principal Cauchy stress occurs in the direction of the greater principal stretch, and we will require that the Cauchy stress increases with the principal stretch (see Ex. 16). To satisfy these conditions the principal stresses must satisfy the Baker-Ericksen inequalities: i - j i - j > 0 , if i j , and > 0 for i = 1, 2, 3 , i i (FD.34) (FD.33)
where i, j are any two principal stretches and i, j are the corresponding principal elastic stresses defined in (FD.25). For an isotropic material these can be expressed in IX.18 7/9/95
terms of the Helmholtz free energy function s by using Eqs. (FD.26) and (FD.24). After some manipulation one finds that they can be rewritten as (see Ex. 17)
s s 2 I + i II > 0 , for each i ,
(FD.35)
and
s 2 i
2
> 0 , for each i .
(FD.36)
Note that these inequalities are special cases of the strong ellipiticity condition. A third condition we will impose is that it takes infinite energy to stretch the solid to infinite length or compress it to zero length (Kwan, 1985; Giaquinta, et al., 1989). In other words, s if i 0 or if i . (FD.37)
The restrictions on the free energy function that have been considered so far have been adapted from similar conditions for a single phase elastic material. There are also constraints that come from the fact that we are dealing with a mixture. For example, the volume fractions must be non-negative and not greater than one. Thus, we require 0 s 1 , (FD.38)
which we will refer to as the solidity condition. To see why this is so important we must first rewrite the continuity equation (II.3s) for the solidity using material coordinates for s s the solid. Recalling that dV s = dV 0 , where V s, V 0 are the current and initial III values of the solid volume, respectively, it follows that s (Xs, t) = 0
III
,
(FD.39)
where s is the material description of the solidity and 0 is its initial value. To satisfy (FD.38) we must have 0 0
II . I
(FD.40)
To illustrate the implications of this inequality consider the uni-axial confined compression experiments discussed in Chapt. IV. From Ex. 8b, III = 2 , and so, for relatively large compressive strains is close to zero. In this case inequality (FD.40) is IX.19 7/9/95
violated (assuming 0 > 0 ). Thus, it appears that we must impose more restrictions on the constitutive relations to guarantee that this does not happen. To determine which relationships are affected we start with the permeability. It is reasonable to assume that k is zero when f = 0 . If this holds then it would take infinite time for the tissue to reach the point where s = 1 (no matter what the strain-rate). The result of this is that even if the elastic stress-strain law allows compressions to the point where (FD.38) no longer holds the permeability would make it so they can never be reached. Therefore, it is assumed that k 0 whenever s 1 . (FD.41)
It is also possible to prevent unrealistic volume fractions by assuming the free energy function has a singular point at s = 1 . Although this fits naturally into the theory (Holmes, 1986) this additional assumption will not be imposed here since (FD.40) is sufficient to guarantee that the solidity condition is satisfied. Exercises: 14. This exercise deals exclusively with a linear isotropic solid. a) Using Tables 1,2 show that 1 1 1 3 s 0s = ( s + 2s )I 2 - sII - (3 s + 2s )I + (3 s + 2s ) . 8 2 4 8 b) Show that i = s (1 + 2 + 3) + 2si where i = 1 + i . c) Show that the Baker-Ericksen inequalities reduce to s > 0 and s + 2s > 0 . d) In the linear theory (see Sec. III.2.4) it is assumed that the Lam coefficients satisfy s > 0 and 3 s + 2s > 0 . Are these conditions stronger or weaker than those in (c)? 15. This problem concerns the first Piola-Kirchhoff stress tensor Ts , which is defined in Exercise 9.3. a) Assuming x = ( 1X 1, 2X 2, 3X 3) show Ts is diagonal with (check on p)
s s III s s Ti = - i0p + 20 i I + (I - i2) II + 2 III . i
b) For the deformation in part (a), show the Baker Ericksen inequalities imply T i increases with stretch but that the largest T i need not occur in the direction of largest principal stretch ( j ).
IX.20
7/9/95
16. One of the more common ways to derive restrictions on the constitutive laws is to consider simplified loading situations. This exercise examines the homogeneous deformations given in Exercise 9.7. a) If the elastic phase is stress free at equilibrium then show that one must require Is + 2 II s + III s = 0 when I = II = 3III = 3 . b) Using the results from uniform dilatation, explain why one should require that Is + 22 II s + 4 III s is negative for 0 < < 1 and is positive for 1 < . In this requirement, I = 3 2 , II = 34 , and III = 6 . c) Using the results from simple shear, explain why one should require Is + II s > 0 if I = II = 3 + 2 , and III = 1 . 17. Derive (FD.35) and (FD.36) from (FD.33) and (FD.34) using (FD.26).
IX.21
7/9/95
5. Summary of Constitutive Laws for an Isotropic Biphasic Mixture
Solid s = - spI + 2 s(a0B + a1B2 + a2B3) , where s s s a0 = I + I II + II III , s s a1 = - II - I III , and s a2 = III . Fluid f = - fpI Momentum Transfer = p s - K(vs - vf) , where K = (f) 2 k , (IV.9) (III.10) (FD.29) (FD.42c) (FD.42a) (FD.22)
(FD.42b)
and the permeability k is independent of D s and D f .
IX.22
7/9/95
6. Uni-Axial Motion
The question that is now addressed is how finite deformation affects the response of the tissue in the loading configuration we know the most about, that is, uni-axial motion. As in Chapters IV and V it is assumed that the tissue is homogeneous, isotropic and that external body forces and inertia can be ignored. Letting the motion be in the zdirection (see Figs. SSP.1 and C.1) then the solid stress in that direction is zz = - sp + E ,
s
(FD.CC.1)
where E is the elastic component and it is determined from the free energy function as given in Eq. (FD.26). For the fluid phase, from Eq. (FD.29), zz = - fp .
f
(FD.CC.1b)
With the above assumptions the momentum equation (II.5s) for the solid in the z direction reduces to f zE = K(vs - vf) , and for the fluid zE = zp , (FD.CC.2b) (FD.CC.2)
where vs (z, t) is the velocity of the solid phase and v f(z, t) is the velocity of the fluid phase in the z direction. Also, the continuity equation (II.6) can be integrated to obtain s vs + fvf = v 0 , (FD.CC.3)
where v0(t) is a constant of integration that is determined later once the boundary conditions have been specified. As before we will reduce the problem to only one equation, and it will be for the solid displacement. In preparation for this we reduce the above equations by solving (FD.CC.3) for the fluid velocity and then substitute the result into the momentum equation for the solid (FD.CC.2). One obtains zE = K (v s - v0) . (f) 2 (FD.CC.4)
IX.23
7/9/95
The momentum equation (FD.CC.2b) for the fluid is actually solvable and one finds that p(z, t) = E + p 0(t) , (FD.CC.4b)
where p0(t) is a constant of integration and, like v 0 , it is determined once the boundary conditions are specified. To continue it is convenient to introduce the material coordinate Z for the along solid with the associated solid displacement U(Z, t) and solidity s (Z, t) . Introducing the stretch = 1 + ZU (FD.CC.5)
then the principal invariants for C are I = 2 + 2 , II = 1 + 22 , and III = 2 (see Ex. 8b). The elastic component of the solid stress given in Eq. (FD.CC.1) is determined from the free energy function using Eq. (FD.26). For the uni-directional motion considered here this relation gives
s s s s E() = 2 0 + 2 II + III , I s s
(FD.CC.6)
where 0 = T0 and 0 is the initial value of the solidity. The momentum equation (FD.CC.4) in material coordinates takes the form ZE = (tU - v0) . k Also, the equation for the solidity (FD.39) is 0 s = . (FD.CC.8) (FD.CC.7)
The permeability k in Eq. (FD.CC.7) is defined as it was before (see Eq. (CC.9)), that is, k = (f) 2 K , (FD.CC.9)
where in the material description the porosity is 0 f = 1 - . (FD.CC.10)
The material coordinates have simplified the solidity equation (FD.CC.8) significantly but it still remains to solve the momentum equation for the solid (FD.CC.7).
IX.24
7/9/95
It is easier to solve this equation by first writing it in terms of the displacement and to do so recall that the invariants depend on the deformation through the single variable . Because of this we can use the chain rule as follows 2U ZE = E Z = E 2 . Z Therefore, from Eq. (FD.CC.5) and (FD.CC.7), the equation to be solved in uni-axial motion in the case of finite deformation is H(ZU)ZU = tU - v0 , where 1 H() = k E() , for = 1 + , (FD.CC.12)
2
0<Z<h,
(FD.CC.11)
and h is the thickness of the undeformed tissue. There are material parameters associated with the coefficients of Eq. (FD.CC.11) that have a profound effect on the response. We introduced two associated with the permeability in Chapt. IV. At that time the variable was the axial strain but now it is more convenient to use the stretch. The parameters in this case are the inherent permeability k0 = k . =1 (FD.CC.12b)
This is the value of the permeability in the absence of an imposed load, that is, when there is no deformation of the solid matrix. Another parameter associated with k is M0 = 1 d k0 d k
=1 ,
(FD.CC.12c)
which measures the sensitivity of the permeability on strain. For example, the greater M0 the faster the permeability decreases with compressive strain. There are several parameters that can be identified for the elastic stress and one of particular importance is HM = E , for = 1 , (FD.CC.12d)
which is the finite deformation analog of the aggregate elastic modulus H A used in the linear theory studied in Chapts. III and IV. In fact, all three parameters (k 0, M 0, H M) are IX.25 7/9/95
present in the small strain theory. To describe the finite deformational characteristics of tissue in uni-axial motion we will include two more parameters. These will be used to describe how the permeability and elastic stress behave for very large strains. Eq. (FD.CC.11) is the fundamental equation for uni-axial motion in the case of finite deformation of the solid matrix. It generalizes Eq. (CC.8) which was derived in Chapt. IV for the case of infinitesimal strains. To be able to solve the above nonlinear diffusion equation for U(Z, t) it is necessary to specify the boundary and initial conditions. These depend on the particular experiment of interest and we will consider two of those described in Chapt. IV. Two constitutive laws also have to be specified. These include the permeability as well as the free energy function for the solid (recall that the free energy function for the fluid was assumed to be constant). As they appear in Eq. (FD.CC.12) they are arbitrary functions of the stretch . We could specify them now but it is best to wait until the relevant experimental data is presented. Note that this is the same situation we found ourselves in the small strain case, although, there we only had to determine the permeability. One thing that will help in determining these constitutive relations is that the above nonlinear diffusion equation is similar, in terms of its mathematical structure, to the one studied in Chapt. IV. Namely, they both have coefficients that depend on the first spatial derivative of the displacement. Although we will eventually determine the exact functional dependence of H in Eq. (FD.CC.11), at the moment it is arbitrary, just as it was when we began the study for small strains. Consequently, most of the solutions derived in Chapt. IV can be adapted relatively easily to solve the present problem. Exercises: s In this Chapter the scalar E designates the elastic component of zz . In the following exercises, the tensor E is the elastic component of s . 18. Show that in uni-axial confined compression, the first Piola-Kirchhoff form of the elastic stress is 0 0 TE = 0 0 E , 0 0 1 and so, TE = E (see Ex. 3). 19. In the case of uniform dilatation, F = I (see Ex. 7c). In this case show that _ E = ()I ,
IX.26
7/9/95
where 2 0 s _ s s () = I + 22 II + 4 III .
s
IX.27
7/9/95
7. Creep
The creep test is considered first because it will give us information on the constitutive relationship we presently know the least about, that is, the free energy function for the solid. This was the same reason why steady state permeation was the first experiment studied for small strains in Chapt. IV since it lead to the determination of the dependence of the permeability on strain. We also need to extend that analysis to finite deformation and this will be considered in the next section. The creep test was described in Chapt. IV but, in brief, the tissue sample is supported on its lower surface so there is no motion at the boundary Z = h of either the fluid or solid phase (Fig. FD.C.1). At the same time the tissue is compressed on its upper surface by applying a static load across a rigid porous filter. The eventual steady state strain, and how it depends on the applied load, is of interest here because it is what will enable us to determine the free energy function. We will also investigate the transient nature of the nonuniform strain distribution in the tissue to see what effects finite deformation has on the creep response. Since the lower boundary is supported by a rigid impervious surface there is no
Axial Compression
Porous Filter
Z = U(0,t)
Tissue
Bone
Z=h
Fig. FD.C.1 Schematic depiction of the confined compression experiment. The rigid porous filter allows fluid to flow freely from the tissue and the confing chamber acts to make the motion uni-directional. IX.28 7/9/95
flow at this interface. So, the equation of motion (FD.CC.11) for uni-axial confined compression reduces to H(ZU)ZU = tU ,
2
0<Z<h,
(FD.C.1)
Due to the static load applied to the upper surface we require E
Z=0
= -F0 ,
(FD.C.2)
and at the lower surface U(h, t) = 0 . We will also assume the motion starts from rest so U(Z, 0) = 0 . 6.1 Equilibrium The first aspect of the creep response that we consider is the equilibrium state. With the displacement independent of t , Eq. (FD.C.1) is simply d2U dZ2 (FD.C.4) (FD.C.3)
= 0 .
If U 0 is the steady state surface displacement then the solution of this equation that also satisfies boundary condition (FD.C.3) is U = 0(h - Z) , where 0 is given as U0 0 = h . The value of U0 is determined from the applied load, that is, it must satisfy E( 0) = - F 0 , for 0 = 1 - 0 . (FD.C.7) (FD.C.5)
(FD.C.6)
This equation is important because we can use it and the experimental determination of the surface displacement to find how the elastic stress depends on . This is relatively
IX.29
7/9/95
easy to do because at equilibrium the stretch is constant, so, it equals the final length of the sample divided by its original length. To illustrate how to use Eq. (FD.C.7) to help determine the constitutive law, the surface stresses obtained from bovine and human articular cartilage are shown in Fig. FD.C.2 as a function of the measured stretch. This tells us how E depends on but, obviously, there are a lot of functions that can be used to fit this data. We will make a specific choice by stating how the free energy depends on the invariants. From this we can then use Eq. (FD.C.7) to fit the elastic stress to the data in Fig. FD.C.2. To make this choice we will require that inequalities (FD.18a,b) are satisfied and condition (FD.37) holds. Out of the infinite number of functions that satisfy these constraints we take the following 0 s = 0
s
e 1(I - 3) + 2(II - 3) , III
(FD.C.8)
where 0 , 1 , 2 are positive constants. We assume there is no pre-stress in the reference configuration from which the deformation is measured. To have E = 0 when 1 = 2 = 3 = 1 it must be that = 1 + 22 . In this case the Baker-Ericksen inequalities (FD.18a,b) are satisfied as long as the i 's are positive (Ex. 21). Also, in the case of small strains the above function produces the same stress-strain relationship used in the infinitesimal theory studied in Chapts. 3 and 4 (Ex. 20). The elastic stress in confined compression can now be written as
2 2 - 1 1 E = 2 H M 2+1 e ( - 1) ,
(FD.C.9a)
where HM = 4( 1+ 22) 0 . (FD.C.9b)
The elastic stress therefore contains two material constants (H M, ). Also, with the above information about the elastic stress the function H in the momentum equation (FD.CC.11) is now H = 2(2 - 1) 2 + 2 + 1 k E . 2( 2 - 1) (FD.C.9c)
IX.30
7/9/95
The coefficient H M in Eq. (FD.C.9a) is the finite strain analogue of the aggregate elastic modulus HA , defined in Eq. (IV.1). To show this let = 1 + , where << 1 . Thus, under infinitesimal conditions, (2 + ) (2 + ) 1 E = 2 H M e (1 + )2b+1 H M [ 1 - 2 + 2 (1 + 4) 2 ] H M . Therefore, H M is a modulus which represents the slope of the stress-strain curve at small strains. This is exactly what it was defined to be earlier in Eq. (FD.CC.12d). The parameter , on the other hand, measures how sharply the stress changes with increasing strain. It is a relatively simple matter to determine the two stress parameters from data like that given in Fig. 2. Using nonlinear regression to fit Eq. (FD.C.9) to the two sets of data points one finds that for bovine articular cartilage HM = 0.33 MPa and = 0.761 and for the human cartilage H M = 0.407 and = 1.105 . The resulting values for the surface stresses are plotted in these figures and it is seen that the assumed constitutive law for the free energy function gives very good agreement with the experimental results. Based on the success of Eq. (FD.C.9) to reproduce the experimental data it would seem that we are more or less through with the constitutive law for the free energy (at least for articular cartilage at equilibrium). This is, however, not the case. It should be pointed out that by requiring the free energy function to satisfy the Baker-Ericksen inequalities we have some assurance that whatever deformations come from these problems that the constitutive relation is capable of producing physically realistic strain fields. A second point that needs to be made concerns the assumption of zero pre-stress at the reference configuration. As discussed in Chapt. I, the equilibrium of unloaded tissues like articular cartilage and the stroma is achieved when there is a balance between the swelling pressure of the proteoglycans and the tensile stress in the collagen matrix. One manifestation of this is that if such a tissue is cut the regions near the cut will swell. In Chapt. 10 we will discuss a model which accounts for the swelling pressure.
1 1
IX.31
7/9/95
a) Stress (MPa)
0.0 -0.2 -0.4 -0.6 -0.8 -1.0 0.3 0.4 0.5 0.6 0.7 Stretch
Theory Experiment
0.8
0.9
1.0
b) Stress (MPa)
0.0 -0.2 -0.4 -0.6 -0.8 -1.0 0.3 0.4 0.5 0.6 0.7 Stretch
Theory Experiment
0.8
0.9
1.0
Fig. FD.C.2 Comparison between the constitutive law (FD.C.9a) for the elastic stress in the solid and the experimental values measured in confined compression for articular cartilage (Kwan, 1985). For both human, (a), and bovine, (b), the constitutive law gives a very good description of the observed stress corresponding to H M , equal to 0.407 MPa, 0.761 and 0.33 MPa, 1.105, respectively. Given the variety of functions that can be used for the free energy function it is worthwhile to consider some other possibilities and a few are given in Table VIII.2. Different invariants are used but it is possible to convert them using the results in Table VIII.1. It is not possible to give a detailed discussion of the merits and limitations of each of them but a few comments are in order. The first entry in Table VIII.2, which is usually credited to St. Venant and Kirchhoff (Truesdell and Nolls, 1992), is nothing more than the strain energy function used in the infinitesimal theory. The strain energy is a
IX.32
7/9/95
nonlinear function of the deformation because the strain tensor is a quadratic function of the deformation gradient (see Eq. (FD.4)). The second entry in Table VIII.2 is the basis of the Mooney-Rivlin theory of rubber and is distinguished from the others by the fact that it is a linear function of the invariants. Because of its simplicity, and its apparent accuracy for rubber, this function has been studied extensively both theoretically and experimentally. This function, like the first entry, has not been widely used to describe soft tissue. The primary reason for this is the simplicity of these functions since they are limited in their dependence on the stretches. For example, in uni-axial motion (either compression or tension) the Mooney-Rivlin model is limited to a quadratic function of the stretch. This is sufficient for perhaps small to moderate strains but at relatively large strains (in either compression or tension) biological tissues usually exhibit almost exponential behavior and in this case the Mooney-Rivlin model is not appropriate. A discussion of the application of this function to blood vessels, and some of its limitations, can be found in Canfield and Dobrin (1987). The third entry is a simplified version of some of the functions that have been proposed to describe the tensile properties of skin as well as aortic and lung tissue. The polynomial portion of this function is included to account for the relatively slow increase in the stress at small strains and the exponential is present to reproduce the very sharp increase in stress at large tensile strains. There are a number of other functions that have been used to describe the stretch of skin and a nice summary of them can be found in Lanir (1987). One of the problems with most of the functions that have been used for skin is that little, if any, attention has been paid to the question of convexity. As indicated in the discussion earlier, if the strain energy function does not have certain convexity properties then the stress can show unphysical characteristics, like having the greatest stress not occur in the direction of greatest stretch. Exercises: 20. a) Show that in the case of small strains the coefficients in (FD.C.8) are related to the Lam parameters in (III.7) as follows s = 42 0 and s = 2(1 + 2) 0 . b) Using the result in (a) show that the Poisson's ratio s s = 2( + ) , s s satisfies 1 0 < s < 3 .
IX.33
7/9/95
21. This problem concerns inequalities (FD.38) and the constitutive law (FD.C.8) for the free energy function for the solid phase. a) Show that
s s 2 2 + i II = (1 + 2 i )s , I
b) Show that
s 2 i
2
= 2 { 2[ ( i - 1)1 + (i j + i k - 2)2 ] 2
2 2 2 2 2
+ (i + 1)1 + (i j + i k + 2)2 } s / i .
2 2 2 2 2 2
22. a) Using constitutive law (FD.C.8) show that in the case of uniform dilatition E _ = ()I , where (see Ex. 19) ( 2 - 1)[1 + 2 2( 2 + 1)] 3( 2 - 1)[ + 2 ( 2 + 1)] _ 1 2 () = 2 0 e . 4+61+122 b) Using constitutive law (FD.C.8) show that in the case of simple shear (see Ex. 3c) ( 1 + 2) 0 0 ( 1 + 2)2 ( + ) 0 0 Ts = 20e 10 2 0 2 2
.
23. In the general development of constitutive laws for elastic materials Truesdell and Noll (1992) state several principles, based on everyday experience, which they expected to be satisfied. These include: i) tension is required to produce uniform dilatation, and pressure is required to produce uniform condensation (the P-C inequality), ii) if a cube of material is stretched along one axis then tension is necessary to keep the opposite lateral faces from moving (the T-E inequality). Show that (FD.8) satisfies these conditions. 24. The Poynting effect occurs in simple shear when there are unequal normal (Cauchy) stresses and the Kelvin effect refers to the situation when normal forces are necessary to maintain a state of simple shear. Using the result of Ex. 22b, show that both of these are present when constitutive law (FD.C.8) is used.
IX.34
7/9/95
6.2 Transient Response To determine the transient response of the tissue in finite deformation it is necessary, in general, to use numerical methods to solve the problem. However, there are special circumstances where asymptotic approximations can be found and one of these is the initial response of the tissue in a creep experiment. The motivation and explanation of this particular solution method is the same as given for the infinitesimal case (Chapt. IV). Using exactly the same arguments as in Chapt. IV one finds that for small values of time U(Z, t) ~
t
Z f . t
(FD.C.9d)
The function f(s), from Eq. (FD.C.1), satisfies the ordinary differential equation 2H(f )f + sf - f = 0 , for 0 < s < , (FD.C.10a)
where, from the applied load condition (FD.C.2), H Mf(0) = -F0 . Also, from the zero displacement conditions (FD.C.3) and (C.12) f() = 0 . (FD.C.10c) (FD.C.10b)
As before, a very important conclusion comes from this result and this states that the surface displacement in the initial moments is U(0, t) ~ f(0) t . (FD.C.11)
The constant f(0) is determined from the solution of Eq. (FD.C.10a). As in Chapt. IV, we can find an approximation of f(0) by making the assumption of a small strain. Doing so one finds that U(0,t) ~ t , where for small strains
s s F 0[ M0(1 - 0) + 20] 2F 0 ~ H 1 - s M HM(1 - 0)
for t << td ,
(FD.C.22a)
k0H M
.
(FD.C.22b)
IX.35
7/9/95
The parameters k0, M 0 appear in the expression for the permeability and are defined in Eqs. (FD.CC.12b,c). The constant t d is the diffusion time scale defined in Eq. (C.11). The assumption of a small strain that was used to derive (FD.C.22b) clearly limits our analysis of the finite deformation problem. However, what (FD.C.22b) does show is the contribution the other nonlinearities at small strains. By comparing the above result with s the small strain case given in (C.22b) it is seen that the difference is the 20 term in the numerator in Eq. (FD.C.22b). Because the initial value of the solidity for many soft hydrated tissues is from 0.2 to 0.3, and since M0 is in the neighborhood of 5 to 10, it is clear that the strain dependence of the permeability makes much more of a contribution at small strain levels than the other nonlinear mechanisms contained in the model.
IX.36
7/9/95
8. Ultrafiltration
The second uni-axial compression test we will consider in this chapter is ultrafiltration. The procedure for this test is the same as used for steady state permeation described in Chapt. IV, and a schematic representation of the experiment is illustrated in Fig. FD.U.1. A cylindrical plug of tissue is placed in a snugly fitting confining chamber filled with fluid. The tissue is supported on its downstream surface ( Z = h ) by a rigid porous filter and a second filter is used on the upstream side to impose a prescribed strain c = u c/h on the sample, where uc is the imposed surface displacement. At equilibrium, so there is no fluid flow, c is uniform through the plug of material. Subsequently, a hydrostatic pressure difference P is applied across the filters to attain a steady flow velocity v 0 of the filtrate. The measurements of the filtration velocity v0(c, P) are
vf = v
0
Porous Filter
Z = Uc
Tissue
p = P
Z=h
Porous Filter
v f = v0
Fig. FD.U.1. Schematic of the experimental configuration for a steady state permeation test in which the tissue sample is compressed between rigid porous filters and a hydrostatic pressure difference P is applied across the tissue. IX.37 7/9/95
used to separate the effects of direct compression (c) from the effects caused by the drag of fluid flow through the tissue (P), that is, these measurements are used to delineate the role of strain and pressure in the tissue. From this we intend to determine the constitutive law for the permeability. In the case of steady flow Eq. (FD.CC.11) reduces to H d2U = - v0 , for 0 < Z < h . dZ2 (FD.U.1)
Because of the clamping strain imposed on the tissue U(0) = uc , and due to the rigid filter at the downstream surface ( Z = h ) U(h) = 0 . Also, from the imposed pressure difference p(0) - p(h) = P . The last condition can be rewritten in terms of the elastic stress using Eq. (FD.CC.4b) and the result is E - E = P . = 0 =h (FD.U.4) (FD.U.3) (FD.U.2)
known.
The three boundary conditions, Eqs. (FD.U.2)-(FD.U.4), are needed because the differential equation is second order in Z and the velocity v 0 in Eq. (FD.U.1) is not
f
As a final comment, the filtration velocity is related to the velocity v 0 of the fluid in the tissue by the relation v0 = fv0 . Note that because of the nonuniform compression within the porous-permeable solid f matrix, v0 is not constant within the tissue.
f
IX.38
7/9/95
7.1 Solution of the Steady State Problem The finite deformation description for ultrafiltration has been reduced to a nonlinear two-point boundary value problem (FD.U.1) for the displacement for the solid phase subject to boundary conditions (FD.U.2)-(FD.U.4). It is relatively easy to solve this problem numerically, however, the analytical solution is also easily obtained since it is possible to integrate the differential equation. The function H in Eq. (FD.U.1) depends on the displacement through the single variable dU = dZ . This fact is going to be used to integrate the equation and with this in mind let the function F() be such that d H() F() = H . d M We can now write Eq. (FD.U.1) as d dU HM dZ F dZ = - v 0 . The procedure we are using to solve the problem is the same as that used for the small strain case studied in Chapt. IV. For example, the above differential equation is the finite deformation version of Eq. (SSP.8) The next steps in the solution are exactly the same as before and one finds that 1 F -1(c0 - km r)dr , U(Z) = - h Z/h
(FD.U.5)
where hv0 . km = H M The constant c 0 and the filtration velocity v0 are determined from the imposed strain condition (FD.U.2) and the imposed pressure condition (FD.U.4). An important special case of the above solution occurs when the imposed pressure is small. This limit was also considered in Chapt. IV, and it lead to the expression for the
IX.39
7/9/95
permeability given in Eq. (SSP.15). Carrying out the same procedure for the finite deformation problem (see Ex. 25) one finds that at low pressure levels k( c) v0 h P , for c = 1 - c , c (FD.U.6)
where c = u c/h is the clamping strain. This result shows that in the case of small imposed pressures the intrinsic and apparent permeabilities are the same. It should be emphasized that k a and k are generally different, particularly when there is a nontrival pressure gradient causing significant deformation in the tissue. Exercise: 25. Derive (FD.U.6) by adapting the argument given in Chapt. IV to the finite deformation problem using the following steps: - c a) Assuming k m is small let c0 ~ c0 + -1km + -2km 2 + ... . Substituting this, and c (FD.U.5), into (FD.U.2) show d -1 3 - c = F -1(c0) + k m(c1 - 0.5) c d F () = -0 + O(km ) .
From this conclude -0 = F(- c) and c1 = 1/2 . c b) Using the result in part (a) and Eq. (FD.U.5) show that
1 Z HM dU = - c - km ( 2 - h ) H(- ) + O(km 2) . dZ c c) Substitute the result of (b) into Eq. (FD.U.4) and expand for small km to derive Eq. (FD.U.6).
7.2 Response in Ultrafiltration The experiments described in Chapt. IV showed that the apparent permeability is a strongly nonlinear function of the strain and applied hydrostatic pressure. To see what the finite deformation biphasic theory predicts we need to determine the constitutive law for the permeability. In doing so it must be remembered that the solidity condition (FD.40) needs to be satisfied, or in other words, the permeability must satisfy (FD.41). There are many functions that will satisfy this condition but for small strains we found that the permeability of articular cartilage is an exponential function of the infinitesimal
IX.40
7/9/95
strain. Based on these observations, the constitutive law we will use in finite deformation is f 0 M(III - 1)/2 e k = k0 , s (1 - 0) (FD.U.12)
where k0 is the permeability when the tissue is unloaded (see Eq. (FD.CC.12b)), is a positive parameter that measures how fast the permeability approaches zero as s 1 , and M comes from the observation that the permeability is an exponential function of the compressive strain. The latter parameter is related to the slope M0 of the normalized permeability curve, defined in Eq. (FD.CC.12c), through the relation M0 = M + 1 - 0 . Also note that because of the denominator in (FD.U.12), k whenever s 0 . This comes from our assumption that the tissue becomes infinitely permeable in the case of a single phase incompressible fluid. From our assumption of homogeneity, k 0, M, and are constants. To show Eq. (FD.U.12) reduces to the constitutive law obtained in Chapt. IV, let = 1 + , so from Eqs. (FD.CC.8) and (FD.CC.10) 0 s = 1 + , and f = 1 + - 0 . 1+
Substituting these expressions into Eq. (FD.U.12), recalling that III = 2 , and assuming << 1 1 + - 0 k = k0 1 - 0
2 eM(2 + )/2
= k 0 1 + 1 - 0
2 eM( + /2)
IX.41
7/9/95
k0 1 + 1 - e M(1 + /2) 0 k0 eM . The exponential is retained in the last expression because M is on the order of 5 to 10. Comparing this with Eq. (SSP.25) it is seen that Eq. (FD.U.12) does indeed reduce to the constitutive law used in Chapt. IV in the case of small strains. In Chapt. IV several methods were discussed that can be used to determine the permeability from the measurement of the filtration velocity. We will only consider the first method, that is, the one that uses the slope of the velocity curve at P = 0 . This method will work in finite deformation because the formula on which it is based, namely Eq. (FD.U.6), is the direct extension of the one used for small strains, Eq. (SSP.15). The data obtained from these measurements for bovine articular cartilage are shown in Fig. SSP.5 and they are replotted in Fig. FD.U.2 except that allowance has now been made for the modified thickness hc in Eq. (FD.U.6). Fitting Eq. (FD.U.6) to these data points it is found that k 0 = 2.519x10 -15 m4/Ns , M = 4.638 , and = 0.0848 . It is assumed here that 0 = 0.2 (see Table VIII.1). The corresponding values from Eq. (FD.U.7) are also shown in Fig. FD.U.2 and clearly they are in very good agreement.
30
Experiment Theory
Permeability
20 10 0
0.0
0.1
0.2 Clamping Strain
0.3
0.4
Fig.FD.U.2. The apparent permeability ( x10-16 m4/Ns) obtained from the filtration velocity measurements on bovine articular cartilage (average of five samples). Also shown are the permeability values obtained from constitutive law (FD.U.12) in the case of k0 = 2.519x10 -15 m4/Ns , M = 4.638 , = 0.0848 and 0 = 0.2 . IX.42 7/9/95
To investigate the steady state permeation response in finite deformation we use parameter values typical for bovine articular cartilage and let HM = 0.33 MPa , = 0.76, k0 = 2.519x10 -15 m4/Ns , M = 4.638 , = 0.0848 , 0 = 0.2 , and h = 2 mm . With these values the apparent permeability as a function of the normalized applied hydrostatic pressure, defined as P p = HM , (FD.U.7)
and clamping strain is shown in Fig. FD.U.3. The definition of the apparent permeability that is used here is ka = h cv0 P . (FD.U.8)
The strongly nonlinear dependence of k a on c is clearly seen in this figure. For example, at lower pressures the permeability decreases by two orders of magnitude over the clamping strains shown. The nonlinear dependence on applied pressure is also evident and in this case the apparent permeability is seen to be more sensitive at smaller clamping strains. The reason for this comes from the relatively sharp increase in the compressive elastic stress at large compressive strains (Fig. FD.C.2). For example, if c = 50% it takes a 0.28 MPa increase in the compressive stress to achieve a change in strain of 5% that at c = 10% only takes 0.06 MPa. This is significant since the permeability depends on the strain and for this reason the drop in ka with applied pressure is not as great at c = 50% as it is at c = 10% (Fig. FD.U.3). The drag on the solid as the fluid moves through the tissue produces a nonuniform compaction of the solid and this is shown in Fig. FD.U.4a in the case of when c = 30% and p = 1 . The permeability and porosity in the tissue are given in Fig. FD.U.4b. As expected the displacement is greater than the linear distribution, i.e., the constant equilibrium strain distribution that is present when there is no fluid flow. This nonuniform compaction causes the permeability and porosity to decrease with depth, which is particularly significant in the vicinity of the downstream surface ( Z = h ). This behavior is due to the fact that the solid is pushed by the fluid, in the direction of the flow, but it is restricted from moving at Z = h . The result is that the solidity in the neighborhood of the downstream boundary increases and it decreases at the upper surface relative to the no flow case.
IX.43
7/9/95
a)
0.3 0.2 0.1
U(Z)/h Equilibrium Solution
U(Z)/h
0.0
0.0
0.2
0.4
Z/h
0.6
0.8
1.0 0.8
b) 100 Permeability
Porosity
Porosity
10
Permeability
0.7
1
0
0.2
0.4
Z/h
0.6
0.8
1
0.6
Fig. FD.U.4. (a) The displacement, solid curve, of the solid phase of the tissue due to the flow of fluid through the tissue. Also shown, by the dashed line, is the linear displacement that is present in the absence of fluid flow. (b) The permeability (x10-16 m4/Ns) and porosity in the tissue due to fluid flow. In both figures c = 30% and P = HM . The increase in the porosity at Z = 0 is accompanied by a decrease in the elastic component of the solid stress, as compared to the case of when there is no applied pressure. This happens because the diffusive drag force acts in a direction opposite to the restoring force of the elastic component of the solid. In fact, if P is large enough the tissue will separate from the filter and the assumed condition of a fixed clamping strain no longer applies. This situation was discussed in depth in Chapt. IV and when this occurs the boundary condition in Eq. (FD.U.2) must be replaced by
IX.44
7/9/95
E = 0 . = 0
(FD.U.9)
larger pressures they are separated. The pressure at which separation occurs increases with clamping strain and this is shown in Fig. FD.U.3 by the curve on the apparent permeability surface. It is worth pointing out that for the linear biphasic model separation occurs when p = 2c but the separation pressure is larger in the finite deformation theory. The fact that the solidity and permeability can change with pressure is potentially significant for the filtration properties of the tissue. To illustrate this, a measure of the effective pore radius, defined as (see Chapt. 5)
This condition is included, where necessary, in the calculation of the apparent permeability in Fig. FD.U.4. To illustrate this situation, in Fig. FD.U.5 the displacement of the upper surface is shown in the case of when a 30% clamping strain is imposed. For applied pressures up to about p = 1.3 the filter and tissue remain in contact but for any
0.4 0.3 U(0)/h 0.2 0.1 0.0 0.0 0.5 1.0 1.5 Normalized Applied Pressure 2.0
Fig. FD.U.5. Displacement of the upper surface, solid curve, of the solid as a function of the applied pressure in the case of when a 30% clamping strain is imposed. The surface stays at 0.3h until the applied pressure reaches about 1.3H M and then for any larger pressure the tissue and upper filter are no longer in contact. The dashed curve is the displacement in the case of when no clamping strain is imposed. IX.45 7/9/95
a=
_ 8k f f
(FD.U.10)
is plotted in Fig. FD.U.6 for three applied pressures using a clamping strain of 30%. In _ these calculations, = 2 and f = 0.01 gm/cm. s Between the lowest and highest pressure the effective radius of the pores at the upper surface ( Z = 0 ) increases by approximately a factor of two while it drops by about the same factor at the downstream surface. As a function of Z , at the larger pressure the average radius of the pore decreases by a factor of five through the specimen. This indicates that the filtration properties of the tissue can be quite different than expected. First, the size of the particles capable of passing through the tissue very clearly depends on pressure. It therefore appears possible to significantly reduce the size of those particles that get through by raising the pressure. Second, larger particles may be able to enter the tissue yet not pass all the way through, thus possibility providing a mechanism for trapping these particles within the deeper layers of the tissue. A condition that is implicitly assumed in the ultrafiltration experiment is that the fluid velocity is at a steady state. However, the viscoelastic nature of the dynamic response requires that a certain amount of time is allowed to pass before the measurement of the velocity is made. To illustrate the transient behavior of the velocity, in Fig. FD.U.7 the velocity of the permeating fluid is given as a function of time for two applied pressures. In both cases c = 30% and the fluid pressure is reached using a ramp loading function over a period of 10 sec. As might be expected, both curves show an increase in
10 8 Radius (nm) 6 4 2 0 0.0 0.2 0.4 Z/h 0.6
p = 0.310 p = 0.605 p = 1.20
0.8
1.0
Fig. FD.U.6. The effective pore radius, Eq. (FD.U.10), in the tissue at three applied pressure levels in the case of when c = 30% . IX.46 7/9/95
the filtration velocity as the pressure increases in the compression stage ( 0 t 10 ). During the relaxation stage ( 10 < t ) the applied pressure is constant but the solid matrix continues to move towards its eventual steady state position. In particular, the solid continues to compact near the downstream surface and this reduces the permeability (see Fig. FD.U.4). The result of this is a decrease in the filtration velocity and this is seen in both curves in Fig. FD.U.7. Thus care must be exercised in any filtration experiment to allow sufficient time for the filtration velocity v 0 to reach equilibrium even though the applied pressure has been reached. Without knowledge of this fact one can be lead to significant misunderstandings of the filtration and permeability properties of the tissue. The large velocities for the p = 1.2 curve in Fig. FD.U.8 enhance the diffusive drag effect. This acts to compress the matrix and this is magnified because the permeability depends on strain. One consequence of this is that the velocity reaches a steady state faster when a higher pressure is used. To illustrate this, the time it takes the velocity of the permeating fluid to reach 99% of its eventual steady state level is plotted as a function of the applied pressure in Fig. FD.U.8 (the clamping strain and the loading function are the same as for Fig. FD.U.7). It takes approximately 10 to 11 minutes to reach a steady state at this clamping strain and clearly the time decreases with increasing applied pressure. This result can be used as a guide as to when one may begin measurement of the steady state permeation velocities in filtration experiments.
3 Filtration Velocity 2 1 0 p = 1.20 p = 0.12
0.1
1
Time (sec)
10
100
Fig. FD.U.7. The transient response of the filtration velocity v0 ( x10-6 cm/sec) for two steady state pressure levels. In both cases, c = 30% and the steady state value of the imposed pressure is reached using a ramp loading function over a period of 10 sec. IX.47 7/9/95
680 660 Time (sec) 640 620 600 0.0 0.2 0.4 0.6 0.8 1.0 Normalized Applied Pressure 1.2
Fig. FD.U.8. The time it takes the filtration velocity to reach 99% of its steady state value as a function of applied pressure. As in Fig. FD.U.7, c = 30% and the fluid pressure is applied using a ramp loading function over a period of 10 sec.
IX.48
7/9/95