- Email: [email protected]
Transverse isotropic poroelastic osteon model under cyclic loading
MECHANICS RESEARCH COMMUNICATIONS
Mechanics Research Communications 32 (2005) 645–651 www.elsevier.com/locate/mechrescom
Transverse isotropic poroelastic osteon model under cyclic loading Agne`s Re´mond, Salah Naili
*
Laboratoire de Me´canique Physique, CNRS (SPI) UMR 7052 B2OA, Universite´ Paris XII—Val-de-Marne, Faculte´ des Sciences et Technologie, 61 Avenue du Ge´ne´ral de Gaulle, 94010 Cre´teil, Ce´dex, France Available online 29 October 2004
Abstract The poroelastic problem associated with a hollow cylinder under cyclic loading is solved. This cylinder models an osteon, basic unit of cortical bone. Both fluid and solid phases are supposed compressible. Solid matrix is modeled as an elastic transverse isotropic material. An explicit close-form solution for the steady state is obtained. Fluid flow distribution as a function of poroelastic properties and cyclic loading is discussed as it could influence bone remodeling. Strain rate of loading is shown to play a significant role in mass flux in the porous material. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Poroelasticity; Transverse isotropy; Compact bone; Osteon
1. Introduction Cortical bone is a dense tissue found at the periphery of bones. Daily activities, such as walking, induce cyclic loading of bone tissue, which in turn adapts to its mechanical loading. Experiments from Lanyon (1984) showed that cyclic loading induces more bone adaptation than static loading. Turner (1988) deduced from experiments on cyclic loading of bone that the stimulus for bone remodeling is proportional to applied strain rate magnitude. Strain rate magnitude can directly be deduced from strain magnitude and frequency of loading. Cortical bone can be modeled as a fluid saturated porous material (Cowin, 2001). An osteon is the cylindrical basic unit of cortical bone and contains one of the level of cortical bone porosity, the lacuna canalicular porosity. At the cellular level, stimuli for bone adaptation has been hypothesized to be the shear stresses induced by interstitial fluid flow in the lacuna canalicular porosity (Burger and Klein-Nulend,
*
Corresponding author. Tel.: +33 1 45 17 14 45; fax: +33 1 45 17 14 33. E-mail addresses: [email protected] (A. Re´mond), [email protected] (S. Naili).
0093-6413/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2004.10.003
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
646
1999). Interstitial fluid flow is considered to occur at this porosity level and to result from mechanical loading of cortical bone. Poroelasticity theory, first developed by Biot, enables to take into account fluid–solid interaction and to model macroscopic averaged fluid flow through bone matrix (Biot, 1955). It allows to determine an average fluid flow in the porosity of a mechanically loaded material. Models of the osteon structure permits to gain insight in the fluid–solid interaction for this type of tissue. Fluid flow distribution and velocity are needed to understand fluid influence on cells lining within the porous matrix and part of bone remodeling mechanotransduction. Analytical solutions enable to study the relative influence of poroelastic characteristics and loading on the mechanical behavior of simple poroelastic models. Zeng et al. (1994) described an osteonÕs behavior modelling it as a hollow cylinder under cyclic loading. The corresponding elastic problem is solved to obtain the stress field. The stress field in the elastic case is used to calculate fluid pressure distribution. Partial fluid stresses have been neglected in this elastic model. Since only the porous matrix stress field was considered, fluid flow was not directly linked to mechanical loading of the osteon. Zhang et al. (1998) gave the close-form solution of a cyclic loaded poroelastic beam. This solution cannot be transposed to the cylindric geometry of the osteon. Furthermore, material isotropy of the tissue is supposed and bone tissue has different properties in longitudinal and radial directions (Turner et al., 1999). Abousleiman and Cui (1998) offered a solution method for poroelastic problems in cylindrical geometry, that could be used for the osteon. However, the time-dependent solution is not explicitly given and the geometry considered here is different. The explicit close-form solution of a poroelastic hollow cylinder is detailed in the following. This model is based on cylindrical geometry of the osteon and its poroelastic properties. Material is modeled as transverse isotropic. Both fluid and solid phases are supposed to be compressible, which is needed for bone tissue as compressibility of bone matrix and of fluid are significantly different. Cyclic loading is applied in the longitudinal direction, parallel to the cylinderÕs axis. This analytical steady state solution allows to determine mass flux through the porous matrix. Furthermore the solution can be used as a reference for numerical results. This introduction has given a brief review of some models used to study the osteonÕs mechanical behavior. Section 2 describes poroelasticity theory applied to an osteonÕs model. Equations of motion are given in the quasi-static case. Then, hypotheses and boundaries conditions corresponding to the osteon are described. The problem is solved and the mass flux is detailed. Results of the poroelastic model are then presented in Section 3. Their implication on the osteonÕs mechanical behavior are discussed in Section 4. Conclusion and prospects of this work are detailed in Section 5.
2. A poroelastic model for the osteon Poroelasticity theory is used to account for fluid–solid interactions in cortical bone. After stating the main equations for the theory, geometry and boundary conditions are specified to build a model for the osteon. OsteonÕs tissue is modeled as a porous material saturated with interstitial fluid. Porous matrix is supposed to be transverse isotropic and with an homogeneous porosity. Matrix material and fluid are assumed to be compressible. Constitutive laws for the solid matrix material and the saturating fluid are written: r
¼ Ce ap;
p
¼ Mðn trðaeÞÞ;
ð1Þ
where r is the total stress tensor, e is the apparent strain tensor of the solid matrix, p is the interstitial fluid pressure and n the fluid content variation. The fourth order drained stiffness tensor C is taken to be isotropic transverse. It depends only on Er and Ez, respectively radial and longitudinal Young modulus of drained
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
647
cortical bone, and on mr and mz, respectively radial and longitudinal drained PoissonÕs ratios. The Biot effective stress tensor a is taken to be isotropic transverse with the same principal directions as the compliance tensor. Biot modulus is noted M and links the fluid content variation to the pressure in the absence of solid matrix deformations. The trace operator is designated by tr. For a physiological range of activities excluding shocks, only low frequencies of loading occur. Thus, inertia terms are neglected in linear momentum conservation equation. No body forces are considered. In this case, linear momentum conservation equation becomes: divðCeÞ a gradðpÞ ¼ 0;
ð2Þ
where div and grad designate the usual gradient and divergence operator. DarcyÕs law is combined with the mass conservation law to obtain: o 1 p þ trðaeÞ ¼ KMp; ð3Þ ot M giving a relation between strains and interstitial fluid pressure. Permeability K ¼ jl is given by the intrinsic permeability j and the interstitial fluid viscosity l. The osteon is modeled as a transverse isotropic poroelastic hollow cylinder. Interior and exterior osteonÕs radii are designated by a and b, respectively. In cylindrical coordinates, the radial unit vector is written r and the longitudinal one z. The model is supposed to be axisymetric and the radial displacement Ur to depend only on the radius r and the longitudinal displacement Uz only on the longitudinal coordinate z. Fluid flow is hypothesized to be only radial since canals where fluid circulates are mostly oriented along the radial direction (Cowin, 2001). As fluid flow is considered to be only radial, the pressure term depends only on the radius r and on time t. Pressure inside the haversian canal is used as a reference pressure. Considering the low frequencies of loading, the haversian canal plays a reservoir role. Its dimensions allow the fluid to relax and its pressure is supposed to remain constant. The exterior surface of the osteon is supposed to be impermeable, meaning that there is no fluid flow through this surface. Full impermeability is an idealization of the existing boundary condition where few canaliculi can cross the cement line and thus allow leaking through the surface. This leaking is supposed to be sufficiently small to be neglected. Boundary conditions to take into account are then: pðr ¼ a; tÞ ¼ 0
op ðr ¼ b; tÞ ¼ 0: or
ð4Þ
Stresses on the haversian canal surface are supposed to be negligible compared to the stresses induced by the fluid flow. It leads to the boundary condition on stress rr = 0 for r = a. Cyclic loading in the longitudinal direction is applied. Applied loading results in constant longitudinal strains for the cylinder with an amplitude of ez0 and at the pulsation x. The pulsation is linked to the frequency f by the relation x = 2pf. Pressure distribution in the porous hollow cylinder is given by: pðr; tÞ ¼ Mðazz 2arr mz Þez0 gðr; tÞ; where g(r, t) is a function of radius and time given by: 1 1 00 qffiffiffiffi qffiffiffiffi qffiffiffiffi qffiffiffiffi ix ix ix ix b I r þ I b K r K 1 0 1 0 c c c c C C BB gðr; tÞ ¼ Re@@ qffiffiffiffi qffiffiffiffi qffiffiffiffi qffiffiffiffi 1Aeixt A: ix ix ix ix K1 b I0 a þ I1 b K0 a c c c c
ð5Þ
ð6Þ
The operator Re gives the real part of the complex number. Functions Ij et Kj are the modified Bessel function of jth order of the first and second kind respectively. The constant c is the consolidation constant for the considered poroelastic material and is given by:
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
648
c¼
j C 11 M ; l Ma2rr þ C 11
ð7Þ 2
z Er mz Þ where C 11 ¼ ð1þmr EÞðEr ðEz E 2 is a component of the stiffness tensor. z mr 2E r mz Þ As shear stresses on cell membranes are directly linked to fluid flow in the lacuno canalicular porosity, the mass flux qr is computed using DarcyÕs law:
qr ¼
j op ; l or
ð8Þ
developing into the following expression: 0 qffiffiffiffi qffiffiffiffi 0 qffiffiffiffi qffiffiffiffi 1 1 rffiffiffiffiffi ix ix ix ix K b I r I b K1 r 1 1 1 c c c c j C C B ixB q ffiffiffi ffi q ffiffiffi ffi q ffiffiffi ffi q qr ¼ Mðazz 2arr mz Þez0 Re@ ffiffiffiffi Aeixt A; @ l c K ix ix ix ix b I a þI b K a 1
c
0
c
1
0
c
ð9Þ
c
which has the dimensions of velocity.
3. Results Using elastic properties of cortical bone, mass flux qr is calculated with values from the literature and for osteonÕs radii a = 50 lm and b = 250 lm (Zhang et al., 1998; Cowin, 2001). Mass flux is proportional to cyclic loading amplitude ez0. It also depends on effective stress coefficients and PoissonÕs ratio. Dependency on frequency is less straightforward as mass flux has its amplitude being proportional to the square root of frequency. Its spatial distribution depends as well on the square root of frequency. Exterior surface is impermeable and qr is null at r = b and maximum at the interior radius for r = a. Fig. 1 shows qr versus time at different radii of the poroelastic cylinder for a constant applied strain e = 0.0002 and a 1 Hz loading frequency. At this frequency, qr varies mostly in amplitude for the different observed radii. Phase shift is not significant between location but is noticeable when compared to the
4 × 10 –11
qr
2 × 10 –11
0
– 2 × 10 –11
– 4 × 10 –11 0
0.2
0.4
0.6
0.8
1
1.2
t Fig. 1. Mass flux qr as a function of time t for a loading with a frequency f = 1 Hz and an applied strain ez0 = 0.0002 at different radii: inside radius r = a (solid line), r = 100 lm (- -), mid-radius r = 150 lm (– -), r = 200 lm (– –).
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
649
1
0.8
q /q r 0
0.6
0.4
0.2
0 0.00005
0.0001
0.00015
0.0002
0.00025
r Fig. 2. Normalized mass flux, mass flux qr divided by its value q0 at r = a, as a function of radius for an applied strain ez0 = 0.0002 for different frequency f = 1 Hz (solid line), f = 2 Hz (- -), f = 10 Hz (– -), f = 20 Hz (– –).
applied cyclic loading, which determinates the time scale. As expected, qr amplitude decreases when the radius increases and approaches the exterior radius where a no fluid flow boundary condition is imposed. As qr depends linearly on applied strain magnitude, its amplitude is directly scaled to the amplitude of cyclic loading. However, frequency dependency also influences spatial distribution of qr. To observe this influence, normalized mass flux, qr divided by its value q0 at the inside radius r = a, is plotted on Fig. 2 as a function of the cylinderÕs radius for different frequencies. At mid-radius, increasing the frequency of loading to f = 20 Hz yields to a qr half of what its normalized value is at the same radius for a frequency
1.2 × 10 –10 1 × 10 –10
qr
8 × 10 –11 6 × 10 –11 4 × 10 –11 2 × 10 –11
0 0
20
40
60
80
100
f Fig. 3. Mass flux qr as a function of frequency f for an applied strain ez0 = 0.0002 at different radii: inside radius r = a (solid line), r = 100 lm (- -), mid-radius r = 150 lm (– -), r = 200 lm (– –).
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
650
2 ×10 –10
qr
1.5 ×10 –10
1 ×10 –10
5 ×10 –11
0 0
5
10
15
20
f Fig. 4. Mass flux qr at the inside radius r = a as a function of frequency f for different applied strain rate e_ z0 ¼ 0:002 s1 (solid line), e_ z0 ¼ 0:001 s1 (- -), e_ z0 ¼ 0:0002 s1 (– -), e_ z0 ¼ 0:0001 s1 (– –).
f = 1 Hz. Mass flux qr decreases more rapidly along the radius at higher frequencies. For higher frequencies, exterior layers in the cylinder exhibit a very low mass flux. Fluid mean velocity is close to zero. Fig. 3 shows qr variation as a function of loading frequency for different radii. For low frequencies under 5 Hz, increasing the frequency of loading yields to an increase of qr for all considered radii. For higher loading frequencies, qr slightly decreases at the mid-radius and decreases more significantly at exterior radius. Loading is applied at a frequency that does not allow fluid flow to fully develop up to the exterior radius. Considering a loading resulting in a constant strain rate for the longitudinal strain, Fig. 4 shows the variation of qr at the inside radius r = a as a function of loading frequency for different strain rates. For example, a strain rate of e_ z0 ¼ 0:002 s1 (solid line) is obtained with an applied strain of ez0 = 0.002 at a frequency of 1 Hz and an applied strain of ez0 = 0.0002 at a frequency of 10 Hz. Mass flux qr peaks around f = 2.5 Hz frequency and the peak is proportional to the loading strain rate.
4. Discussion Considering a poroelastic hollow cylinder cyclically loaded in the longitudinal direction allowed to study the influence of poroelastic characteristics on pressure distribution and thus mass flux. Zeng et al. (1994) solution gave a similar pressure distribution but its magnitude is around 6% lower. Low porosity and pressure values explain this small difference and that fluid role in the poroelastic model is not very important. However since the pressure distribution is similar, these results allows to refine mechanotransduction signals in model Zeng et al. (1994) by using the exact magnitude. Furthermore this model permitted to link strains in the longitudinal direction to the mass flux qr, which could have an influence on mechanotransduction signals. Considering shear stresses induced by fluid flow are those signals, the distribution in the fluid flow velocity showed that shear stresses would be significantly higher near the Haversian canal than at the periphery of the osteon where the velocity is lower. This effect is even amplified at higher frequencies as the distributions of normalized mass flux as a function of radius showed. The permeability in the model is constant throughout the cylinder, so it does not take into account possible variations of the canaliculiÕs number at the exterior of the osteon. The boundary condition imposed at the exterior of the cylinder is weakly justified as the cement line could play a significant role in modulating fluid flow at the exterior
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
651
of the osteon structure. Furthermore, canaliculi can reach the cement line and enable fluid leaking through the corresponding surface. However, it was supposed that there are few canaliculi and that this would not influence greatly pressure building. Frequency influence has been observed on bone tissue as loading at small strains but high frequencies could lead to bone remodeling as well as lower frequencies higher strains (Turner, 1988). For this simplified model of an osteon, it is noticeable that frequency dependency gave an increase for qr at the inside radius but a decrease after the mid-radius. There might be a threshold value for frequency of loading after which the osteon exterior layers are not reached by fluid coming from the haversian canal. For higher frequencies, pressure in the haversian canal might not stay constant and modify fluid flow through the cylinder. Considering that bone remodeling stimuli could be directly proportional to strain rate, frequencies around f = 2.5 Hz at a constant strain rate of loading yield to a peak for qr. It could indicate that strain rate induce higher mechanotransduction signals at low frequencies. These results are also consistent with the observations of Warden and Turner on bone adaptation (Warden and Turner, 2004). They showed that cortical bone formation peaks around 5–10 Hz frequencies and seem to be less influenced at higher frequencies. These results show that fluid flow magnitude in bone increase with frequency till it reaches a plateau. They could explain partially the difference in bone adaptation when loading frequency goes over a threshold value. A more detailed model taking into account inhomogeneities could contribute to a better understanding of this phenomenon.
5. Conclusion Close-form solution of a poroelastic hollow cylinder modeling an osteon was obtained and used to study variations in the mass flux. Frequency of loading influence the distribution of mass flux as well at its amplitude during the loading cycle. Loading at a constant strain rate leads to a frequency range where mass flux peaks. This solution can also be used as a benchmark for numerical studies of models that could modeled more precisely the geometry and inhomogeneities of the osteon.
References Abousleiman, Y., Cui, L., 1998. Poroelastic solutions in transversely isotropic media for wellbore and cylinder. Int. J. Solids Struct. 35 (34–35), 4905–4929. Biot, M.A., 1955. Theory of elasticity and consolidation for a porous anisotropic solid. J. Appl. Phys. 26 (2), 182–185. Burger, E.H., Klein-Nulend, J., 1999. Mechanotransduction in bone: role of the lacuno-canalicular network. Faseb J. 13 (Suppl.), S101–S112. Cowin, S.C., 2001. Bone Mechanics Handbook, second ed. CRC Press, Boca Raton, FL. Lanyon, L.E., 1984. Functional strain as a determinant for bone remodeling. Calcif. Tissue Int. 36 (Suppl. 1), S56–S61. Turner, C.H., 1988. Three rules for bone adaptation to mechanical stimuli. Bone 23 (5), 399–407. Turner, C.H., Rho, J., Takano, Y., Tsui, T.Y., Pharr, G.M, 1999. The elastic properties of trabecular and cortical bone tissues are similar: results from two microscopic measurement techniques. J. Biomech. 32 (4), 437–441. Warden, S.J., Turner, C.H., 2004. Mechanotransduction in the cortical bone is most efficient at loading frequencies of 5–10 Hz. Bone 34 (2), 261–270. Zeng, Y., Cowin, S.C., Weinbaum, S., 1994. A fiber matrix model for fluid flow and streaming potentials in the canaliculi of an osteon. Ann. Biomed. Eng. 22 (3), 280–292. Zhang, D., Weinbaum, S., Cowin, S.C., 1998. On the calculation of bone pore water pressure due to mechanical loading. Int. J. Solids Struct. 35 (34–35), 4981–4997.
Mechanics Research Communications 32 (2005) 645–651 www.elsevier.com/locate/mechrescom
Transverse isotropic poroelastic osteon model under cyclic loading Agne`s Re´mond, Salah Naili
*
Laboratoire de Me´canique Physique, CNRS (SPI) UMR 7052 B2OA, Universite´ Paris XII—Val-de-Marne, Faculte´ des Sciences et Technologie, 61 Avenue du Ge´ne´ral de Gaulle, 94010 Cre´teil, Ce´dex, France Available online 29 October 2004
Abstract The poroelastic problem associated with a hollow cylinder under cyclic loading is solved. This cylinder models an osteon, basic unit of cortical bone. Both fluid and solid phases are supposed compressible. Solid matrix is modeled as an elastic transverse isotropic material. An explicit close-form solution for the steady state is obtained. Fluid flow distribution as a function of poroelastic properties and cyclic loading is discussed as it could influence bone remodeling. Strain rate of loading is shown to play a significant role in mass flux in the porous material. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Poroelasticity; Transverse isotropy; Compact bone; Osteon
1. Introduction Cortical bone is a dense tissue found at the periphery of bones. Daily activities, such as walking, induce cyclic loading of bone tissue, which in turn adapts to its mechanical loading. Experiments from Lanyon (1984) showed that cyclic loading induces more bone adaptation than static loading. Turner (1988) deduced from experiments on cyclic loading of bone that the stimulus for bone remodeling is proportional to applied strain rate magnitude. Strain rate magnitude can directly be deduced from strain magnitude and frequency of loading. Cortical bone can be modeled as a fluid saturated porous material (Cowin, 2001). An osteon is the cylindrical basic unit of cortical bone and contains one of the level of cortical bone porosity, the lacuna canalicular porosity. At the cellular level, stimuli for bone adaptation has been hypothesized to be the shear stresses induced by interstitial fluid flow in the lacuna canalicular porosity (Burger and Klein-Nulend,
*
Corresponding author. Tel.: +33 1 45 17 14 45; fax: +33 1 45 17 14 33. E-mail addresses: [email protected] (A. Re´mond), [email protected] (S. Naili).
0093-6413/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2004.10.003
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
646
1999). Interstitial fluid flow is considered to occur at this porosity level and to result from mechanical loading of cortical bone. Poroelasticity theory, first developed by Biot, enables to take into account fluid–solid interaction and to model macroscopic averaged fluid flow through bone matrix (Biot, 1955). It allows to determine an average fluid flow in the porosity of a mechanically loaded material. Models of the osteon structure permits to gain insight in the fluid–solid interaction for this type of tissue. Fluid flow distribution and velocity are needed to understand fluid influence on cells lining within the porous matrix and part of bone remodeling mechanotransduction. Analytical solutions enable to study the relative influence of poroelastic characteristics and loading on the mechanical behavior of simple poroelastic models. Zeng et al. (1994) described an osteonÕs behavior modelling it as a hollow cylinder under cyclic loading. The corresponding elastic problem is solved to obtain the stress field. The stress field in the elastic case is used to calculate fluid pressure distribution. Partial fluid stresses have been neglected in this elastic model. Since only the porous matrix stress field was considered, fluid flow was not directly linked to mechanical loading of the osteon. Zhang et al. (1998) gave the close-form solution of a cyclic loaded poroelastic beam. This solution cannot be transposed to the cylindric geometry of the osteon. Furthermore, material isotropy of the tissue is supposed and bone tissue has different properties in longitudinal and radial directions (Turner et al., 1999). Abousleiman and Cui (1998) offered a solution method for poroelastic problems in cylindrical geometry, that could be used for the osteon. However, the time-dependent solution is not explicitly given and the geometry considered here is different. The explicit close-form solution of a poroelastic hollow cylinder is detailed in the following. This model is based on cylindrical geometry of the osteon and its poroelastic properties. Material is modeled as transverse isotropic. Both fluid and solid phases are supposed to be compressible, which is needed for bone tissue as compressibility of bone matrix and of fluid are significantly different. Cyclic loading is applied in the longitudinal direction, parallel to the cylinderÕs axis. This analytical steady state solution allows to determine mass flux through the porous matrix. Furthermore the solution can be used as a reference for numerical results. This introduction has given a brief review of some models used to study the osteonÕs mechanical behavior. Section 2 describes poroelasticity theory applied to an osteonÕs model. Equations of motion are given in the quasi-static case. Then, hypotheses and boundaries conditions corresponding to the osteon are described. The problem is solved and the mass flux is detailed. Results of the poroelastic model are then presented in Section 3. Their implication on the osteonÕs mechanical behavior are discussed in Section 4. Conclusion and prospects of this work are detailed in Section 5.
2. A poroelastic model for the osteon Poroelasticity theory is used to account for fluid–solid interactions in cortical bone. After stating the main equations for the theory, geometry and boundary conditions are specified to build a model for the osteon. OsteonÕs tissue is modeled as a porous material saturated with interstitial fluid. Porous matrix is supposed to be transverse isotropic and with an homogeneous porosity. Matrix material and fluid are assumed to be compressible. Constitutive laws for the solid matrix material and the saturating fluid are written: r
¼ Ce ap;
p
¼ Mðn trðaeÞÞ;
ð1Þ
where r is the total stress tensor, e is the apparent strain tensor of the solid matrix, p is the interstitial fluid pressure and n the fluid content variation. The fourth order drained stiffness tensor C is taken to be isotropic transverse. It depends only on Er and Ez, respectively radial and longitudinal Young modulus of drained
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
647
cortical bone, and on mr and mz, respectively radial and longitudinal drained PoissonÕs ratios. The Biot effective stress tensor a is taken to be isotropic transverse with the same principal directions as the compliance tensor. Biot modulus is noted M and links the fluid content variation to the pressure in the absence of solid matrix deformations. The trace operator is designated by tr. For a physiological range of activities excluding shocks, only low frequencies of loading occur. Thus, inertia terms are neglected in linear momentum conservation equation. No body forces are considered. In this case, linear momentum conservation equation becomes: divðCeÞ a gradðpÞ ¼ 0;
ð2Þ
where div and grad designate the usual gradient and divergence operator. DarcyÕs law is combined with the mass conservation law to obtain: o 1 p þ trðaeÞ ¼ KMp; ð3Þ ot M giving a relation between strains and interstitial fluid pressure. Permeability K ¼ jl is given by the intrinsic permeability j and the interstitial fluid viscosity l. The osteon is modeled as a transverse isotropic poroelastic hollow cylinder. Interior and exterior osteonÕs radii are designated by a and b, respectively. In cylindrical coordinates, the radial unit vector is written r and the longitudinal one z. The model is supposed to be axisymetric and the radial displacement Ur to depend only on the radius r and the longitudinal displacement Uz only on the longitudinal coordinate z. Fluid flow is hypothesized to be only radial since canals where fluid circulates are mostly oriented along the radial direction (Cowin, 2001). As fluid flow is considered to be only radial, the pressure term depends only on the radius r and on time t. Pressure inside the haversian canal is used as a reference pressure. Considering the low frequencies of loading, the haversian canal plays a reservoir role. Its dimensions allow the fluid to relax and its pressure is supposed to remain constant. The exterior surface of the osteon is supposed to be impermeable, meaning that there is no fluid flow through this surface. Full impermeability is an idealization of the existing boundary condition where few canaliculi can cross the cement line and thus allow leaking through the surface. This leaking is supposed to be sufficiently small to be neglected. Boundary conditions to take into account are then: pðr ¼ a; tÞ ¼ 0
op ðr ¼ b; tÞ ¼ 0: or
ð4Þ
Stresses on the haversian canal surface are supposed to be negligible compared to the stresses induced by the fluid flow. It leads to the boundary condition on stress rr = 0 for r = a. Cyclic loading in the longitudinal direction is applied. Applied loading results in constant longitudinal strains for the cylinder with an amplitude of ez0 and at the pulsation x. The pulsation is linked to the frequency f by the relation x = 2pf. Pressure distribution in the porous hollow cylinder is given by: pðr; tÞ ¼ Mðazz 2arr mz Þez0 gðr; tÞ; where g(r, t) is a function of radius and time given by: 1 1 00 qffiffiffiffi qffiffiffiffi qffiffiffiffi qffiffiffiffi ix ix ix ix b I r þ I b K r K 1 0 1 0 c c c c C C BB gðr; tÞ ¼ Re@@ qffiffiffiffi qffiffiffiffi qffiffiffiffi qffiffiffiffi 1Aeixt A: ix ix ix ix K1 b I0 a þ I1 b K0 a c c c c
ð5Þ
ð6Þ
The operator Re gives the real part of the complex number. Functions Ij et Kj are the modified Bessel function of jth order of the first and second kind respectively. The constant c is the consolidation constant for the considered poroelastic material and is given by:
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
648
c¼
j C 11 M ; l Ma2rr þ C 11
ð7Þ 2
z Er mz Þ where C 11 ¼ ð1þmr EÞðEr ðEz E 2 is a component of the stiffness tensor. z mr 2E r mz Þ As shear stresses on cell membranes are directly linked to fluid flow in the lacuno canalicular porosity, the mass flux qr is computed using DarcyÕs law:
qr ¼
j op ; l or
ð8Þ
developing into the following expression: 0 qffiffiffiffi qffiffiffiffi 0 qffiffiffiffi qffiffiffiffi 1 1 rffiffiffiffiffi ix ix ix ix K b I r I b K1 r 1 1 1 c c c c j C C B ixB q ffiffiffi ffi q ffiffiffi ffi q ffiffiffi ffi q qr ¼ Mðazz 2arr mz Þez0 Re@ ffiffiffiffi Aeixt A; @ l c K ix ix ix ix b I a þI b K a 1
c
0
c
1
0
c
ð9Þ
c
which has the dimensions of velocity.
3. Results Using elastic properties of cortical bone, mass flux qr is calculated with values from the literature and for osteonÕs radii a = 50 lm and b = 250 lm (Zhang et al., 1998; Cowin, 2001). Mass flux is proportional to cyclic loading amplitude ez0. It also depends on effective stress coefficients and PoissonÕs ratio. Dependency on frequency is less straightforward as mass flux has its amplitude being proportional to the square root of frequency. Its spatial distribution depends as well on the square root of frequency. Exterior surface is impermeable and qr is null at r = b and maximum at the interior radius for r = a. Fig. 1 shows qr versus time at different radii of the poroelastic cylinder for a constant applied strain e = 0.0002 and a 1 Hz loading frequency. At this frequency, qr varies mostly in amplitude for the different observed radii. Phase shift is not significant between location but is noticeable when compared to the
4 × 10 –11
qr
2 × 10 –11
0
– 2 × 10 –11
– 4 × 10 –11 0
0.2
0.4
0.6
0.8
1
1.2
t Fig. 1. Mass flux qr as a function of time t for a loading with a frequency f = 1 Hz and an applied strain ez0 = 0.0002 at different radii: inside radius r = a (solid line), r = 100 lm (- -), mid-radius r = 150 lm (– -), r = 200 lm (– –).
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
649
1
0.8
q /q r 0
0.6
0.4
0.2
0 0.00005
0.0001
0.00015
0.0002
0.00025
r Fig. 2. Normalized mass flux, mass flux qr divided by its value q0 at r = a, as a function of radius for an applied strain ez0 = 0.0002 for different frequency f = 1 Hz (solid line), f = 2 Hz (- -), f = 10 Hz (– -), f = 20 Hz (– –).
applied cyclic loading, which determinates the time scale. As expected, qr amplitude decreases when the radius increases and approaches the exterior radius where a no fluid flow boundary condition is imposed. As qr depends linearly on applied strain magnitude, its amplitude is directly scaled to the amplitude of cyclic loading. However, frequency dependency also influences spatial distribution of qr. To observe this influence, normalized mass flux, qr divided by its value q0 at the inside radius r = a, is plotted on Fig. 2 as a function of the cylinderÕs radius for different frequencies. At mid-radius, increasing the frequency of loading to f = 20 Hz yields to a qr half of what its normalized value is at the same radius for a frequency
1.2 × 10 –10 1 × 10 –10
qr
8 × 10 –11 6 × 10 –11 4 × 10 –11 2 × 10 –11
0 0
20
40
60
80
100
f Fig. 3. Mass flux qr as a function of frequency f for an applied strain ez0 = 0.0002 at different radii: inside radius r = a (solid line), r = 100 lm (- -), mid-radius r = 150 lm (– -), r = 200 lm (– –).
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
650
2 ×10 –10
qr
1.5 ×10 –10
1 ×10 –10
5 ×10 –11
0 0
5
10
15
20
f Fig. 4. Mass flux qr at the inside radius r = a as a function of frequency f for different applied strain rate e_ z0 ¼ 0:002 s1 (solid line), e_ z0 ¼ 0:001 s1 (- -), e_ z0 ¼ 0:0002 s1 (– -), e_ z0 ¼ 0:0001 s1 (– –).
f = 1 Hz. Mass flux qr decreases more rapidly along the radius at higher frequencies. For higher frequencies, exterior layers in the cylinder exhibit a very low mass flux. Fluid mean velocity is close to zero. Fig. 3 shows qr variation as a function of loading frequency for different radii. For low frequencies under 5 Hz, increasing the frequency of loading yields to an increase of qr for all considered radii. For higher loading frequencies, qr slightly decreases at the mid-radius and decreases more significantly at exterior radius. Loading is applied at a frequency that does not allow fluid flow to fully develop up to the exterior radius. Considering a loading resulting in a constant strain rate for the longitudinal strain, Fig. 4 shows the variation of qr at the inside radius r = a as a function of loading frequency for different strain rates. For example, a strain rate of e_ z0 ¼ 0:002 s1 (solid line) is obtained with an applied strain of ez0 = 0.002 at a frequency of 1 Hz and an applied strain of ez0 = 0.0002 at a frequency of 10 Hz. Mass flux qr peaks around f = 2.5 Hz frequency and the peak is proportional to the loading strain rate.
4. Discussion Considering a poroelastic hollow cylinder cyclically loaded in the longitudinal direction allowed to study the influence of poroelastic characteristics on pressure distribution and thus mass flux. Zeng et al. (1994) solution gave a similar pressure distribution but its magnitude is around 6% lower. Low porosity and pressure values explain this small difference and that fluid role in the poroelastic model is not very important. However since the pressure distribution is similar, these results allows to refine mechanotransduction signals in model Zeng et al. (1994) by using the exact magnitude. Furthermore this model permitted to link strains in the longitudinal direction to the mass flux qr, which could have an influence on mechanotransduction signals. Considering shear stresses induced by fluid flow are those signals, the distribution in the fluid flow velocity showed that shear stresses would be significantly higher near the Haversian canal than at the periphery of the osteon where the velocity is lower. This effect is even amplified at higher frequencies as the distributions of normalized mass flux as a function of radius showed. The permeability in the model is constant throughout the cylinder, so it does not take into account possible variations of the canaliculiÕs number at the exterior of the osteon. The boundary condition imposed at the exterior of the cylinder is weakly justified as the cement line could play a significant role in modulating fluid flow at the exterior
A. Re´mond, S. Naili / Mechanics Research Communications 32 (2005) 645–651
651
of the osteon structure. Furthermore, canaliculi can reach the cement line and enable fluid leaking through the corresponding surface. However, it was supposed that there are few canaliculi and that this would not influence greatly pressure building. Frequency influence has been observed on bone tissue as loading at small strains but high frequencies could lead to bone remodeling as well as lower frequencies higher strains (Turner, 1988). For this simplified model of an osteon, it is noticeable that frequency dependency gave an increase for qr at the inside radius but a decrease after the mid-radius. There might be a threshold value for frequency of loading after which the osteon exterior layers are not reached by fluid coming from the haversian canal. For higher frequencies, pressure in the haversian canal might not stay constant and modify fluid flow through the cylinder. Considering that bone remodeling stimuli could be directly proportional to strain rate, frequencies around f = 2.5 Hz at a constant strain rate of loading yield to a peak for qr. It could indicate that strain rate induce higher mechanotransduction signals at low frequencies. These results are also consistent with the observations of Warden and Turner on bone adaptation (Warden and Turner, 2004). They showed that cortical bone formation peaks around 5–10 Hz frequencies and seem to be less influenced at higher frequencies. These results show that fluid flow magnitude in bone increase with frequency till it reaches a plateau. They could explain partially the difference in bone adaptation when loading frequency goes over a threshold value. A more detailed model taking into account inhomogeneities could contribute to a better understanding of this phenomenon.
5. Conclusion Close-form solution of a poroelastic hollow cylinder modeling an osteon was obtained and used to study variations in the mass flux. Frequency of loading influence the distribution of mass flux as well at its amplitude during the loading cycle. Loading at a constant strain rate leads to a frequency range where mass flux peaks. This solution can also be used as a benchmark for numerical studies of models that could modeled more precisely the geometry and inhomogeneities of the osteon.
References Abousleiman, Y., Cui, L., 1998. Poroelastic solutions in transversely isotropic media for wellbore and cylinder. Int. J. Solids Struct. 35 (34–35), 4905–4929. Biot, M.A., 1955. Theory of elasticity and consolidation for a porous anisotropic solid. J. Appl. Phys. 26 (2), 182–185. Burger, E.H., Klein-Nulend, J., 1999. Mechanotransduction in bone: role of the lacuno-canalicular network. Faseb J. 13 (Suppl.), S101–S112. Cowin, S.C., 2001. Bone Mechanics Handbook, second ed. CRC Press, Boca Raton, FL. Lanyon, L.E., 1984. Functional strain as a determinant for bone remodeling. Calcif. Tissue Int. 36 (Suppl. 1), S56–S61. Turner, C.H., 1988. Three rules for bone adaptation to mechanical stimuli. Bone 23 (5), 399–407. Turner, C.H., Rho, J., Takano, Y., Tsui, T.Y., Pharr, G.M, 1999. The elastic properties of trabecular and cortical bone tissues are similar: results from two microscopic measurement techniques. J. Biomech. 32 (4), 437–441. Warden, S.J., Turner, C.H., 2004. Mechanotransduction in the cortical bone is most efficient at loading frequencies of 5–10 Hz. Bone 34 (2), 261–270. Zeng, Y., Cowin, S.C., Weinbaum, S., 1994. A fiber matrix model for fluid flow and streaming potentials in the canaliculi of an osteon. Ann. Biomed. Eng. 22 (3), 280–292. Zhang, D., Weinbaum, S., Cowin, S.C., 1998. On the calculation of bone pore water pressure due to mechanical loading. Int. J. Solids Struct. 35 (34–35), 4981–4997.