## Abstract

The governing equations for the theory of poroelastic materials with hierarchical pore space architecture and compressible constituents undergoing small deformations are developed. These equations are applied to the problem of determining the exchange of pore fluid between the vascular porosity (PV) and the lacunar–canalicular porosity (PLC) in bone tissue due to cyclic mechanical loading and blood pressure oscillations. The result is basic to the understanding of interstitial flow in bone tissue that, in turn, is basic to understanding of nutrient transport from the vasculature to the bone cells buried in the bone tissue and to the process of mechanotransduction by these cells. A formula for the volume of fluid that moves between the PLC and PV in a cyclic loading is obtained as a function of the cyclic mechanical loading and blood pressure oscillations. Formulas for the oscillating fluid pore pressure in both the PLC and the PV are obtained as functions of the two driving forces, the cyclic mechanical straining and the blood pressure, both with specified amplitude and frequency. The results of this study also suggest a PV permeability greater than 10^{−9} m^{2} and perhaps a little lower than 10^{−8} m^{2}. Previous estimates of this permeability have been as small as 10^{−14} m^{2}.

## 1. Introduction

A poroelastic model for the interstitial fluid flow space in bone tissue, with a reasonably accurate anatomical model for the architecture of its pore space structure, is described. In order to characterize the special type of porous material pore structure considered, the phrase ‘hierarchical’ is used as an adjective to modify ‘poroelasticity’. Alternatively, it could be described as a set of nested porosities like a set of Russian nested dolls or babushka (matryoshki), a set of dolls of decreasing sizes placed one inside another; each doll but the smallest may be opened to reveal another doll of the same sort inside. The idea of a smaller structure within a larger, similarly shaped, structure is the idea that is to be transferred from a set of babushka to sets of different pore structures in a porous material. The body fluids in tissues reside in such nested, topologically similar, pore structures with different pore sizes in bone, tendon, meniscus and possibly other tissue types.

The animal vascular tree is an example of a pore structure with two such nested systems that are connected. In a microcirculatory bed, blood flows from arteries to arterioles, then to capillaries and on to venules and into the veins; in each of these pore structures the pore size is relatively uniform, but it varies monotonically between levels of porosity characterized by their pore size. The arterial system consists of the capillaries nested within the arterioles that are nested within the arteries. The venous system consists of the capillaries nested within the venules that are nested within the veins. The capillary plexuses of the two nested systems are connected.

Body fluids in tissues reside in such nested, topologically similar, pore structures with different pore sizes in bone, tendon, meniscus and possibly other tissue types. The nesting or ordering criterion is the porosity pore size. The nested porosities are connected so the pore fluid may flow easily through each and across the boundaries between the two nearest neighbour porosities, but any particular pore size porosity may only interchange its pore fluid with the next larger pore size porosity and the next smaller pore size porosity. The flow of interstitial fluid in tissues such as bone, tendon, meniscus and possibly other tissue types is similar to blood flow in the vascular system in the sense that the different pore size porosities are nested, but unsimilar in three important aspects: (i) there are only two levels of pore size porosity important for bone fluid flow; (ii) there is no flow out of, or in from, the smaller pore size porosity into any even smaller pore size porosity, only into, or in from, a larger pore size porosity; and (iii) the fluid flow direction reverses in the normal physiological function of these tissues.

With only the exception of the animal vascular tree, the applications of poroelasticity to fluid movement in biological tissues have simply transferred to tissues the models of the pore structure employed in geomechanics. Existing theories of the poroelasticity of materials with multiple connected porosities with different characteristic sizes and therefore different permeabilities do not address the case of nested porosities. These existing theories are appropriate for their intended use, fractured porous geological structures, but they are not appropriate for the biological tissues of interest; the nested porosities in biological tissue are hierarchical based on the average diameters of their fluid transport channels while the multiple porosity poroelasticity theories for fractured geological structures are democratic; their fluid transport channels of a particular size may exchange fluid with transport channels of any pore size. The primary objective of this work is to provide a model of a poroelastic pore structure appropriate for bone tissue; it is a model that is easily extended to other tissues such as the tendon and the meniscus. Concerning bone, a principal focus is the modelling of the mechanical and blood pressure load-driven movement of the interstitial bone fluid flow.

The absence of the assumption of incompressible constituents is a significant difference between the version of poroelasticity theory employed and the poroelasticity theory used for previous published solutions involving soft tissues. The assumption of incompressible constituents, while appropriate for soft tissues, is inaccurate for hard tissues. Gailani & Cowin (2008) have recently obtained the solution for the unconfined compression of an annular, transversely isotropic, poroelastic hollow cylinder with compressible constituents. On the basis of this solution a protocol has been devised for an experimental test procedure to determine tissue permeabilities for the smallest nested bone porosity, the osteonal lumen wall and the osteonal cement line. This protocol extends to bone tissue an experimental technique that has been very effective in determining soft tissue poroelastic properties (Gailani *et al*. in press).

Current theoretical and experimental evidence suggests that the bone cells in the lacunae (pores) of the lacunar–canalicular porosity (PLC) are the principal mechanosensory cells of bone, and that they are activated by the induced drag from fluid flowing through the PLC (Weinbaum *et al*. 1994; Cowin *et al*. 1995). The movement of bone fluid from the region of the bone vasculature through the canaliculi and the lacunae of the surrounding mineralized tissue accomplishes three important tasks. First, it transports nutrients to the cells in the lacunae buried in the mineralized matrix. Second, it carries away the cell waste. Third, the bone fluid exerts a force on the cell process, a force that is large enough for the cell to sense. This is thought to be the basic mechanotransduction mechanism in bone, the way in which bone senses the mechanical load to which it is subjected. Understanding bone mechanotransduction is fundamental to the understanding of how to treat osteoporosis, how to cope with microgravity in long-term manned space flight and how to design prostheses that are implanted in bone tissue to function for longer periods.

The porosities of bone tissue and their permeabilities are described in the following section. In §3 the innovation in poroelasticity theory to include the hierarchical structure to the pore space is described. In §4 the poroelastic boundary value problems employed to model the transport of interstitial pore fluid between porosity levels in bone tissue are formulated. In §§5 and 7 solutions to these boundary value problems are obtained and in §6 the fluid exchange between the porosity levels is determined. The implications for the determination of the permeabilities of these porosities are discussed in §8. These results are discussed in §9.

A slightly unconventional tensor–matrix notation is employed in this presentation. Its objective is to represent fourth-rank tensors as matrices that are composed of tensor components, something that the classical Voigt (1928) matrix notation for the anisotropic elasticity tensor does not achieve. In the notation employed here second- and fourth-rank tensors in three dimensions are represented as vectors and second-rank tensors, respectively, in six dimensions. Transformations in the six-dimensional space, corresponding to three-dimensional transformations, are six-by-six matrix multiplications that are easily entered and quickly computed with the symbolic algebra software (MAPLE, MATHEMATICA, MACSYMA and MATLAB). In particular, the three-dimensional fourth-rank elasticity tensor is represented as a second-rank tensor in a space of six dimensions. This notation has been developed and employed in a number of papers (Cowin 2003; Cowin & Doty 2007; Cowin & Mehrabadi 2007). In particular, Cowin & Mehrabadi (2007) employed this notation in the development of poroelasticity that is followed in the present paper.

## 2. Bone tissue porosity and permeability

The pore sizes in cortical bone are of approximately three discrete magnitudes: the largest pore size (approx. 50 *μ*m diameter) is associated with the vascular porosity (PV), the second largest pore size (approx. 0.3 *μ*m diameter) with the canaliculi in the PLC and the smallest pore size (approx. 10 nm diameter) is in collagen-apatite porosity (PCA) (Cowin 1999). The typical pore size of the PV channels is not the blood vessel pore size; rather it is the size of the tubular tunnels (Haversian systems or osteons and Volkmann canals) containing the blood vessels, the arterioles and the venules, with the actual dimensions of these vessels subtracted from the volume of the tubular tunnels. The interstitial pore fluid in the PCA has been shown (Wehrli & Fernández-Seara 2005) to be bound to the solid structure, and the PCA is not of interest in the present considerations for that reason. Thus, the two pore size structures of interest in the flow of bone pore fluid are the PV and the PLC; the ratio of the pore sizes of these two porosities is approximately 167.

The measurement of the permeability of the PV has not been accomplished with sufficient accuracy to date, primarily because of the topological intertwining of the PV with the PLC. This difficulty and others were discussed in some detail by Johnson (1984), who mentioned experiments by Rouhana *et al*. (1981) that reported bovine cortical bone permeabilities of the order of 10^{−14} m^{2}. Johnson (1984) presented factors to suggest that the actual permeability could be much larger than that value. The PV permeabilities reported in the past are thought to be significant underestimates because these previous values represent a PLC and PV lumped measurement rather than a PV measurement alone. These lumped measurements also compromise the measurement of the permeability of the PLC as pointed out by Beno *et al.* (2006) in the discussion of the reported PLC permeabilities (Dillaman *et al.* 1991; Steck *et al.* 2003). There are four reported estimates of the permeability of the PLC. A theory-based estimate of the intrinsic permeability of the PLC was provided by Zhang *et al.* (1998) as *k*=1.47×10^{−20} m^{2}. An experimentally based estimate was provided later by Smit *et al.* (2002), *k*=2.2×10^{−22} m^{2}. Oyen (2008) employed a nanoindentation technique with a poroelastic analysis to provide an estimate of *k*=4.1×10^{−24} *m*^{2} (this corresponds to a hydraulic permeability of κ=*k*/*μ*=4.1×10^{−21} m^{4} Ns^{−1}). The assumption of an incompressibility in the analysis of experimental data renders the Oyen (2008) estimate smaller than it should be. A measurement of the PLC permeability by Gailani *et al.* (in press) based on the analytical model of Gailani & Cowin (2008) yielded values of the order of 10^{−24} m^{2} to 10^{−25} m^{2}. Each of these latter estimates is two orders of magnitude smaller than the previous one for the permeability of the PLC. Thus, the estimates for the ratio of the permeabilities of the PV to those of the PLC are as high as 10^{10}. These ratio estimates and the approximate ratio of the diameters as 167 are measures of the significant size difference, pore pressure difference and relaxation times in these two distinct pore size porosities.

The pore fluid pressures in these two pore size bone porosities are distinct and vary very differently with time under mechanical loading of the whole bone. Under physiologically possible rapid rise-time loadings of bone, the pore fluid pressure may rise considerably in the PLC (Zhang *et al.* 1998). The decay time for this pore pressure rise is much larger in the PLC than it is in the PV (Johnson *et al.* 1982; Johnson 1984). The PV is a low-pore-fluid-pressure domain because the PV permeability is sufficiently large to permit a rapid decay of a pressure pulse. This must be the case because the PV contains thin-walled blood vessels carrying blood with a pressure of 40–60 mm Hg (5.3–8 kPa); a pore fluid pressure significantly greater than 40–60 mm Hg will collapse these blood vessels and a prolonged increase in the pore fluid pressure significantly above 40–60 mm Hg would deprive the tissue of oxygen and nutrients. The thin capillary walls of the endothelium, the absence of a muscle layer and the sparse basement membrane lower the collapse pressure of blood vessel below that of a vessel in soft tissue.

These considerations suggest that the PV and PLC function almost independently, the prime mechanical influences for the two porosities being very different as are the time scales of their response. The mechanical loading of the whole bone moves the bone fluid in the PLC. When the bone is compressed, the bone fluid is driven from the PLC into the low-pressure PV and, when the bone is in tension, or the compression is reduced, the bone fluid is sucked from the PV into the PLC. These drainage and imbibing processes occur on a pressure time scale that is much larger than the short pressure adjustment relaxation time for the PV and therefore have minimum influence on the pressure in the PV. The change in interstitial pore fluid pressure in the PV owing to inflow or outflow of bone fluid from the PLC is insignificant because of the time period of pressure adjustment that is much shorter (estimates of these time periods are in Zhang *et al.* 1998) than the pressure adjustment time period for the PLC. While the bone fluid in the PLC is significantly affected for a longer time by the mechanical loading of the whole bone, the bone fluid pressure in the PV is hardly affected because the PV relaxes the pressure pulse by diffusion so rapidly.

## 3. One poroelastic porosity level in a hierarchy of poroelastic porosity levels

There are two objectives of this section. One is to present the model of a hierarchy of poroelastic porosity levels and the other is to summarize the governing equations of saturated, compressible and poroelasticity in the form that they will be employed in this paper. A porous material containing a hierarchical system of uniform porosity levels is considered. One typical porosity level is then focused upon. This porosity level is denoted by *P*, and *P*+1 represents the next porosity level with the larger uniform pore size and the porosity level with the next uniform smaller pore size is denoted by *P*−1. Because the pore size and, we assume, the total pore volume are different at the different porosity levels it follows that the permeability will be different at the different porosity levels as well as the elastic constants associated with the different porosity levels. The lowest uniform porosity level occurs for the matrix material level in which the porosity is zero; this is denoted by *P*_{0}=0 or *P*_{m}, where m indicates the matrix material. The elastic constants at the matrix material level, *P*_{m}, will be greater than those at the *P*_{1}=*P*+1 level because the elastic constants at the *P*_{1} level will be compromised by the porosity not existent at the lower level. It follows then that the elastic constants at each higher level will be further compromised because at each higher level the material contains an additional level of porosity. Thus, while the elastic constants at each higher level are reduced, the total permeability is increased at each level. Traditionally the material at this ‘no porosity’ level has been called the matrix material; however, in the hierarchical scheme of porosity levels, the porosity level below the porosity level under consideration plays the role of the ‘material’ for the particular porosity level under consideration. Thus, only in the case of the first level of porosity above the ‘no porosity’ level do the terminologies of matrix material level and ‘no porosity’ level coincide.

In this section the basic equations of poroelasticity are recorded for a continuum with a *P* level porosity admitting the possibility of fluid exchange (sources and/or sinks) from the larger uniform pore size porosity continuum at level *P*+1 and the next uniform smaller pore size porosity level continuum at *P*−1. Thus a separate continuum model represents each level of porosity *P*. The fluid transport connection between the different continuum model porosity levels occurs through sources (sinks) connecting two different porosity levels. In this section the sources (sinks) are treated in a general way because their exact nature will vary with the particular porous material being modelled. In subsequent sections addressing bone tissue specifically, examples of the porosity level connections will be made specific.

The traditional theory of poroelasticity (Biot 1941, 1955; Rice & Cleary 1976; Thompson & Willis 1991; Cowin & Mehrabadi 2007) is recovered from this hierarchical system of uniform porosity levels if the possibility of fluid exchange between porosity levels is excluded and the generic *P* level is taken to be *P*_{1}=*P*+1 level porosity with the *P*−1 level taken to be the matrix material level *P*_{0}=0 or *P*_{m}; that is to say, the structure proposed is a generalization of traditional poroelasticity to a hierarchical system of uniform porosity levels that allows fluid exchange between the different levels. Each porosity level is a separate continuum level model that interacts with the other continuum level porosity level models by sources and sinks that transfer fluid between the different levels.

In the remainder of this section the equations governing the generic *P* level porosity are recorded; they are equivalent to the traditional equations of poroelasticity except for the fluid supply terms from the two next nearest porosity levels. There are three sets of poroelastic constitutive equations, the stress–strain–pore pressure, the fluid content–stress–pore pressure and Darcy’s Law. The purpose of this section is now to record the form of these constitutive equations free from the incompressibility constraint on any constituent. In the following four subsections, these equations are described in the order indicated above for the constitutive equations and followed by a subsection on the governing equations.

### (a) Stress–strain–pore pressure relations

The displacement vector is denoted by **u** and second-order tensor representation in three dimensions of the strain tensor is **E**. The strain–displacement relations in three dimensions are written in the form
3.1
The symbol is used to represent the second-order strain tensor **E** as a vector in six dimensions,
The stress–strain–pore pressure constitutive relation for a saturated porous medium is that the strain in the saturated porous medium is linearly related, not only to the stress but also to the fluid pressure *p* in the fluid-filled pores; thus, one can write the stain–stress–pore pressure relation as
3.2
or the stress–pore pressure–strain relation as
3.3
where represents the drained anisotropic compliance elastic constants of the saturated porous medium at the porosity level *P* and is its reciprocal, the drained anisotropic elasticity tensor. They are called the drained elastic constants because they are measured in situations in which the fluid pressure *p*^{P} in the fluid-filled pores is negligible or zero. Draining all the pores before the test or doing the test so slowly causes the pores to drain from a negligible fluid pressure. In a porous medium the pores are assumed to be connected; there are no unconnected pores that prevent the flow of fluid through them. The six-dimensional vector (three-dimensional symmetric second rank tensor **A**^{P}), introduced by equation (3.2), is called the *Biot effective stress coefficient vector* (six-dimensional) or *tensor* (three-dimensional) at the porosity level *P*. The Biot effective stress coefficient vector is related to the difference between effective drained elastic constants and the elastic compliance tensor at the *P*−1 porosity level, , by the formula (Nur & Byerlee 1971; Carroll 1979; Cowin & Doty 2007; Cowin & Mehrabadi 2007)
3.4
where is the six-dimensional vector representation of the three-dimensional unit tensor **1**. The Biot effective stress coefficient vector is so named because it is employed in the definition of the *effective stress*,
3.5
This definition of effective stress reduces the stress–strain–pressure relation (3.2) to the same form as Hooke’s Law; thus,
3.6
This constitutive equation is a modification of Hooke’s Law to include the effect of the pore pressure. When *p*^{P}=0 the stress–strain–pore pressure relations (3.2) and (3.3) or (3.5) and (3.6) coincide with anisotropic Hooke’s Law. The advantage of the representation (3.6) is that the fluid-saturated porous material may be thought of as an ordinary elastic material with a compliance matrix , but one subjected to the ‘effective stress’ rather than (ordinary) stress . The drained elastic compliance tensor may be evaluated from knowledge of the pore structure of the medium and the elastic compliance tensor , or vice versa, using composite or effective medium theory. At the lowest level of porosity the elastic compliance tensor will coincide with the elastic compliance tensor of the (pore free) matrix material .

### (b) Fluid content–stress–pore pressure relations

The fluid content–stress–pore pressure constitutive relation involves all the basic field variables for poroelasticity, the total stress , the pore pressure *p*^{P}, the strain in the solid matrix and the variation in (dimensionless) fluid content *ζ*^{P}. The variation in fluid content *ζ*^{P} is the variation of the fluid volume per unit volume of the porous material owing to diffusive fluid mass transport; it is defined as the difference between the strain of the pore space and the strain of the fluid volume in the pore space. The variation in fluid content *ζ*^{P} is equal to the product of the porosity *ϕ*^{P} at the level *P* and the fluid density *ρ*_{f}, divided by a reference value of the fluid density *ρ*_{f0}; thus, *ζ*^{P}=*ϕ*^{P}*ρ*_{f}/*ρ*_{f0}. The variation in fluid content *ζ*^{P} depends not only upon the strain in the solid matrix but also upon the strain induced by the pore pressure *p*^{P}; thus,
3.7
where
3.8
and where *ϕ*^{P} is the porosity at the level *P*, *K*^{f} is the bulk modulus of the fluid and is the Reuss effective bulk modulus of an anisotropic elastic material (Cowin & Mehrabadi 2007), in this case the porous material at the *P*−1 level of porosity. In the case of isotropy this expression for *Λ*^{P} reduces to
3.9
The variation in fluid content *ζ*^{P} is also linearly related to both the stress and the pore pressure *p*^{P},
3.10
where
3.11
The relation between and *Λ*^{P} is
3.12

### (c) Darcy’s Law

The constitutive equations of poroelasticity developed thus far are the strain–stress–pore pressure relations (3.2) and the fluid content–stress–pore pressure relation (3.7). The other constitutive relation of poroelasticity is Darcy’s Law, which relates the fluid mass flow rate, *ρ*_{f}**v**, to the gradient (**∇***p*^{P}) of the pore pressure *p*^{P} (Darcy 1856),
3.13
The coefficient tensor **H**^{P} in Darcy’s Law may be represented by **H**^{P}=*ρ*_{f0}**K**^{P}/*μ*, where **K**^{P} is the intrinsic Darcy’s Law permeability tensor, *ρ*_{f0} is a reference value of the fluid density and *μ* is the fluid viscosity. The intrinsic permeability tensor **K**^{P} has units of length squared and is a function of the porous structure only, not the fluid in the pores; thus, Darcy’s Law (3.13) takes the form
3.14
where the mass flux **q** has the dimensions of velocity.

### (d) Mass and momentum conservation

The conservation of mass in the form expressed by the equation of continuity is
3.15
The form of the mass conservation equation (3.15) is altered to apply to the *P* level pore fluid volume by first replacing *ρ* by *ϕ*^{P}*ρ*_{f} in equation (3.15) and then dividing the equation through by *ρ*_{f0}; thus,
3.16
This relationship is easily modified to account for fluid influxes or outfluxes from the larger uniform pore size porosity level *P*+1 and the next uniform smaller pore size porosity level *P*−1; thus,
3.17
where the supply terms from the next nearest-neighbour porosities are denoted by *γ*^{P−1} and *γ*^{P+1}. These supply terms have the dimensions of density over time and represent the fluxes of pore fluids between porosity levels. The form of equation (3.17) is changed by noting that the variation in fluid content *ζ*^{P} is given by *ζ*^{P}=*ϕ*^{P}*ρ*_{f}/*ρ*_{f0} and the use of equation (3.14); thus,
3.18
The stress equations of motion in three dimensions at porosity level *P*,
3.19
have no simple representation in the six-dimensional vector notation, and the conventional notation is employed; represents the acceleration and **b** the body force.

### (e) Governing differential equations

There are many methods of approach to the solution of poroelastic problems for compressible media. The method selected depends upon the information that is provided and the fields that are to be calculated. One approach that has been effective is to solve for the variation in fluid content *ζ*^{P} if the stress or the strain field is known or may be calculated without reference to the variation in fluid content *ζ*^{P}. The diffusion equation for the variation in volume fraction is obtained by first substituting Darcy’s Law (3.14) into the expression (3.18) for the conservation of mass and subsequently eliminating the pore pressure by use of equation (3.7); thus,
3.20
where is the six-dimensional vector representation of the three-dimensional double gradient tensor defined by **O**=**∇**⊗**∇** (tr**O**=∇^{2}). The components of are
The differential equation (3.20) shows that the time rate of change of the fluid content *ζ*^{P} is due either to fluid flux or to volume changes that are caused by the strain field. It is possible to replace the strain on the right-hand side of equation (3.20) by stress, in which case equation (3.20) becomes
3.21
where and *Λ*^{P} are related by equation (3.21).

Diffusion equations for the pressure field are also employed in the solution of poroelastic problems. The first diffusion equation for the pore pressure field is obtained by first substituting Darcy’s Law (3.14) into the expression (3.18) for the conservation of mass and subsequently eliminating the variation in fluid content *ζ*^{P} by use of equation (3.10); thus,
3.22
The alternative diffusion equation for the pore pressure field is obtained by replacing the strain on the right-hand side of equation (3.22) by stress; thus,
3.23
The boundary conditions on the pore pressure field customarily employed in the solution of this differential equation are that: (i) the external pore pressure *p*^{P} is specified at the boundary (a lower pressure permits flow across the boundary); (ii) the pressure gradient **∇***p*^{P} at the boundary is specified (a zero pressure gradient permits no flow across the boundary); and (iii) some linear combination of (i) and (ii) is specified.

## 4. Outline of the coupled boundary value problems considered

Using the hierarchical scheme described in the previous section a model is formulated in this section for the transport of bone interstitial fluid between the PV and PLC porosity levels in osteonal cortical bone. A section of this bone is illustrated in figure 1. The osteon at the top of this figure is entirely PLC porosity except for its central lumen, called the osteonal canal or Haversian canal, which is part of the PV porosity. The PV porosity consists of the volume of all the tunnels in bones that contain blood vessels and includes all the osteonal canals and all the Volkmann canals with the volume of the tunnels actually occupied by the blood vessels removed.

In this paper the whole bone and the single osteon will both be modelled as poroelastic hollow circular cylinders. The typical osteon, which are numerous as suggested in figure 1, is represented by a poroelastic hollow circular cylinder. The hollow circular cylinder representing the typical osteon contains only the PLC. The numerous poroelastic hollow circular cylinder PLC models connect through their inner cylindrical wall to the PV; the hollow part of osteon is actually part of the PV. The inner surface of the typical cylinder representing the PLC is the surface across which the two porosities exchange pore fluids. The PLC is assumed to permit flow across its inner radial boundary, but not its outer radial boundary. While other assumptions with respect to the outer boundary are possible, an earlier study (Wang *et al.* 1999) showed that this is a reasonable assumption. The whole bone PV is also assumed to permit flow across only its internal radial boundary, into the medullary canal. The typical PLC hollow cylinder represents the typical osteon of figure 1. The PV hollow cylinder is the entire bone of figure 1 with the central lumen of the whole bone, the medullary canal, constituting the hollow part of the PV model. Both of these models are continuum models and the transport connection between the two continua is the outflow–influx that occurs across the osteonal or Volkmann inner wall between the PLC and the PV (Fritton *et al.* 2001). In the domain between the inner surface and outer surface of the PV cylinder there are areal or area sources–sinks that permit interchange of fluid between the two continuum models representing the PV and PLC.

In this model the fluid movement will be driven by three external force systems. The whole bone and its surrounding soft tissue structures are assumed to be cyclically strained in the long bone direction by the axial strain ϵ(*t*)=ϵ_{0}e^{iωt}. This straining occurs as a result of mechanical loading and muscle stimulation. At the endosteum, the wall of the central lumen of the whole bone, the medullary canal, the pore fluid pressure is assumed to be the same as the blood pressure, , where mm Hg. We assume that the periosteum is impermeable; it has been identified as a barrier to interstitial fluid flow (Li *et al.* 1987). We note that lymphatic circulation is unlikely to play a role in bone fluid transport in normal bone because lymphatic vessels are absent in bone. A chapter on the physiology of blood circulation in a book entitled *Blood vessels and lymphatics in organ systems* (Abramson & Dobrin 1984) contains no description of the bone lymphatic system. Vittas & Hainau (1989) and Edwards *et al.* (2008) determined by immunohistochemistry that lymphatics were not present in normal bone.

In the next section of this paper the generic poroelastic hollow circular cylinder model, considered in a general way in the present section, will be specialized to the PLC. In the third section following, the generic poroelastic hollow circular cylinder model will be specialized to the PV.

## 5. Poroelastic model of the PLC level in bone tissue

The general formulation of the mathematical problem for compressible poroelastic hollow circular cylinders has been described earlier (Cowin & Mehrabadi 2007; Gailani & Cowin 2008). As earlier, transversely isotropic material symmetry is assumed. The porosity level is designated by L as shorthand for PLC; later, V is shorthand for PV. The formulation is based on the following assumptions concerning the functional dependence of variables: cylindrical coordinates are employed and axial symmetry is assumed, and displacement field components , , the pressure *p*^{L} and the variation in fluid content *ζ*^{L} are assumed to have the functional dependencies indicated below (Cowin & Mehrabadi 2007, eqn. (47); Gailani & Cowin 2008, eqn. (1)):
5.1
where ϵ(*t*) is the applied axial loading that is ϵ(*t*)=ϵ_{0}e^{iωt} in this case. Thus, the key variables are taken to be *u*^{L}(*r*,*t*),*p*^{L}(*r*,*t*) and *ζ*^{L}(*r*,*t*) and the solutions are assumed to be of the form
5.2
where the frequency *Ω* is associated with the blood pressure, , and the specific solutions sought are for *u*^{Lω}(*r*),*u*^{LΩ}(*r*),*p*^{Lω}(*r*),*p*^{LΩ}(*r*), *ζ*^{Lω}(*r*) and *ζ*^{LΩ}(*r*). The non-zero stress components *T*_{rr},*T*_{θθ} and *T*_{zz} are determined to be (Cowin & Mehrabadi 2007, eqn. (60); Gailani & Cowin 2008, eqn. (20))
5.3
where the components of are the components of the drained elasticity tensor at the PLC porosity level and are those of the Biot effective stress coefficient vector at the PLC porosity level. The values of these coefficients and others used in plotting the curves in this paper are given in table 1. Table 2 contains formulas relating the various coefficients used in the calculations that follow.

The only equation of equilibrium that is not satisfied automatically is (Cowin & Mehrabadi 2007, eqn. (61); Gailani & Cowin 2008, eqn. (21))
5.4
Substitution of the non-zero stress components (5.3) into the equilibrium equation (5.4) yields a result that is subsequently integrated with respect to *r* to produce the following expression:
5.5
where is a function of time obtained in the integration with respect to *r* and is the *r* component of the Biot effective stress coefficient. As both the functions *u*^{L}(*r*,*t*) and *p*^{L}(*r*,*t*) in equation (5.2) are of the form of a product of a function of *r* times e^{iωt} plus another, but different, function of *r*, times *e*^{iΩt}, we assume that the function of time obtained in the integral (equation (5.5)) has the representation . The expression (5.5) for (*∂**u*^{L}/*∂**r*)+(*u*^{L}/*r*)=(1/*r*)(*∂*/*∂**r*)(*r**u*^{L}) may be used to calculate two expressions for the quantity ,
5.6
(Cowin & Mehrabadi 2007, eqn. (59); Gailani & Cowin 2008, eqn. (23)), one in terms of the variation in fluid content *ζ*^{L} and the other in terms of the pressure *p*^{L}. These expressions are needed in the diffusion equations (3.20) and (3.21) for *ζ*^{L} and *p*^{L}, respectively. To find in terms of *p*^{L}, substitute equation (5.5) in equation (5.6) to obtain
5.7
and to find in terms of *ζ*^{L}, eliminate *p*^{L} in equation (5.7) using equation (3.7) with *P* replaced by *L*; thus,
5.8

where the dimensionless quantity *J*^{L} is introduced,
5.9
and where
5.10
and
5.11
which are eqns (6.8) and (19C), respectively, in Cowin & Mehrabadi (2007). (There is an error in eqn. (6.8) in Cowin & Mehrabadi (2007), which is corrected in equation (5.10) above: specifically the last term in (5.10) is divided by 9 in Cowin & Mehrabadi (2007).) In the equations above, the kernel letter *K* stands for a bulk modulus, the superscripts d and f indicate drained material and the fluid, respectively, m indicates the matrix material and the subscript Reff indicates the effective Reuss lower bound on the indicated bulk moduli (Cowin & Mehrabadi 2007). In the compressible, transversely anisotropic case, the expression for *J*^{L} is given by Cowin & Mehrabadi (2007, eqn. (57)),
5.12
The diffusive differential equation for the pore pressure is obtained by substituting equation (5.7) and into equation (3.22) above, neglecting the flux terms from the other porosity levels; thus,
5.13
or, rearranging terms,
5.14
and introducing the notation
5.15
where equation (5.9) has been employed, it follows that
5.16
Introducing the dimensionless radius *λ*=*r*/*r*_{0}, the representation *p*^{L}(*r*,*t*)=*p*^{Lω}(*r*)e^{iωt}+*p*^{LΩ}(*r*)*e*^{iΩt} from equation (5.2) and the two dimensionless driving frequencies constant and , as well as a constant of dimension stress *Υ*^{L},
5.17
the differential equation takes the form
5.18
The solution to this equation is
5.19
This solution is subject to the following boundary conditions at *λ*=1 and *a*. The first boundary condition is that there is no flow across the cement line at *λ*=1. The second boundary condition is at the wall of the osteonal or Haversian canal, *λ*=*a*, where the PLC and the PV meet; thus,
5.20
where the superscript V indicates the PV. The second part of equation (5.20) is unusual in that it connects the two hierarchical porosity levels, and the appearance of *β* on the right-hand side of equation (5.20) is the connection to the vascular level. Recall that the symbols *r*, *r*_{0} and *r*_{ i} were used for radial variable, the external radius and internal radius of the PLC, respectively, and *R*,*R*_{0} and *R*_{i} will be used for the analogous quantities in the PV. The dimensionless radius was *λ*=*r*/*r*_{0} in the PLC and *β*=*R*/*R*_{0} is used for the dimensionless radius in the PV. Thus equation (5.20) states that the osteonal canal boundary PLC pore pressure depends upon where the osteon is located in the whole bone cross section. The first of the boundary conditions (5.20) shows that
5.21
which reduces the solution to
5.22
where
5.23
and where *ϖ* stands for either or , *λ* represents a dimensionless spatial variable and α represents a specific value of *λ*. The function *Π*(*ϖ*:*λ*,α) will appear in many of the subsequent equations and it is useful to note that
5.24
where
5.25
For future use we also introduce ,
5.26
In the special case when *λ*=α the functions (5.25) and (5.26) have the following values:
5.27
where is a Wronskian and the formula for its value employed in equation (5.27) is given by eqn. (9.6.15) of Abramowitz & Stegun (1964).

The second boundary condition in equation (5.20) may be used to eliminate and from equation (5.22). Thus,
5.28
and
5.29
where
5.30
In appendix A the boundary condition on the stress component *T*_{rr} is used to determine and and to eliminate the constants of integration, and . This is accomplished by requiring that the effective radial stress vanish at the outer surface of the osteon (*λ*=1); thus,
5.31
where equation (3.5) has been employed and the following dimensionless quantities have been defined:
5.32
The pore fluid pressure field (5.29) has the following form when the relations (5.31) are inserted into it:
5.33
where
5.34
and
5.35
where *Γ*(*ϖ*:*λ*,1) and are dimensionless. A plot of *p*^{L}(*λ*,*t*) in the case when *p*^{Vω}(*β*) and *p*^{VΩ}(*β*) vanish is shown in figure 2.

## 6. The fluid exchange between the PLC and the PV

The outflow–influx between the PLC level and the PV level will be calculated directly from the pore fluid pressure gradient at the osteonal wall of the PLC level. The pressure field is then given by equation (5.33) above where *p*^{V}=*p*^{Vω}(*β*)*e*^{iωt}+*p*^{VΩ}(*β*)e^{iΩt} is the fluid pore pressure in the PV at the osteonal wall, equation (5.20). The radial pressure gradient at the osteonal wall is then given by
6.1
where
6.2
Substituting equation (6.2) into equation (6.1), the radial pressure gradient at the osteonal wall is written
6.3
The flux *F*^{L} from the luminal surface or osteonal wall is calculated as follows:
6.4
where **q** is obtained from Darcy’s Law (3.14) and d**a**=(2*π**r*_{i}d*l*)**e**_{r} where d*l* is a differential length along the axis of the osteon and this integration is over a unit length *l*_{r}. Thus flux per unit length of an osteon is given by
6.5
where
and
6.6
The total volume of the pore fluid exchanged between the PLC and the PV, denoted by *V*_{total} in one loading cycle *π*/*ω*, is given by (the imaginary part, in this case)
6.7
If the pore pressure in the PV is zero (i.e. both *p*^{Vω}(*β*) and *p*^{VΩ}(*β*) are zero) there is only the pore pressure due to the mechanical loading driving the flow in and out; then from equation (6.7) and the first of equation (6.6) the total volume in the case of zero PV pore pressure is
6.8
and thus, using the first and last of equation (5.17), is given by
6.9
where the identity has been employed. This result shows that the total volume of the pore fluid exchanged between the PLC and the PV decreases as the frequency increases. The total volume of the pore fluid that may be exchanged must be less than or equal to the total volume of the pore fluid in a unit length of an osteon; thus,
6.10
Plots of *V*_{total} against *ω* are shown in figure 3 for values of *ϕ*^{L}, the PLC porosity, equal to 0.05, 0.1 and 0.15 for the case when the PV pressure is zero.

In the PV level continuum model the exchange of fluid with the PLC is represented by a distributed source (sink) of the influx of fluid to the PV from the PLC or the outflow of fluid to the PV from the PLC. The flux *F*^{L} (equation (6.5)) across the osteonal wall from the PLC to the PV (or vice versa) is related to the distributed source (sink) *γ*^{L} by
6.11
where denotes the area of an osteon inside the cement line, including the luminal domain. Notice that *γ*^{L} has the dimensions of density over time and represents the flux of fluid in or out of the PV from or to the PLC per unit area of the whole bone cross section.

## 7. Poroelastic model of the PV level in bone tissue

The poroelastic model for the PV level is similar to that for the PLC level in several aspects. The governing poroelastic theory and the domain shape for the two boundary value problems are the same and the boundary conditions have the same form, but boundary values are different for one condition. However, the PV has a supply term while the PLC does not, and the poroelastic constants of the two continuum models are different. Because of the similarity between these boundary value problems we seek to minimize the repetitive effort by several strategies. First, when similar equations from the PLC derivation are repeated in the PV derivation the superscript L will be replaced by the superscript V if appropriate; it will always be appropriate for the poroelastic constants.

We begin with the development of the differential equations for the two components of the vascular pore pressure *p*^{V}(*β*,*t*)=*p*^{Vω}(*β*)e^{iωt}+*p*^{VΩ}(*β*)e^{iΩt}; recall from the text following equation (5.21) that the symbols *r*, *r*_{0} and *r*_{i} were used for the radial variable, the external radius and internal radius of the PLC, respectively, and *R*, *R*_{0} and *R*_{i} will be used for the analogous quantities in the PV. The dimensionless radius was *λ*=*r*/*r*_{0} in the PLC whereas *β*=*R*/*R*_{0} is used for the dimensionless radius in the PV. The ratio of the inner to the outer radii in the osteon was denoted by *a*=*r*_{i}/*r*_{0} and for whole bone the analogous ratio is *A*=*R*_{i}/*R*_{0}. The representation for the two dimensionless frequencies and is obtained by updating equation (5.2) from *L* to *V* , *r*_{0} to *R*_{0} and using the superposed tilde rather than the superposed bar employed to the osteon; thus,
7.1
Similarly the differential equation (5.18) is updated from *L* to *V* , *λ* to *β* and the supply term to the PV from the PLC is added, consistent with the formulation of equation (3.22),
7.2
where *γ*^{L} is given by equation (6.11). The supply term to the PV from the PLC, the term in which *γ*^{L} occurs, is of primary interest in the present calculation, as it did not appear in the equation for the pore pressure in the PLC. The form of this supply term changes from the form that appears in equation (3.22) to the form that appears in the equation above owing to the manipulation of the differential equation for the vascular pore pressure. The steps are not recorded here for the PV pore pressure, but they were recorded for the same manipulations of the PLC pore pressure in the steps outlined by equations (3.22)–(5.18). In the equivalent of the transition from equation (5.13) to equation (5.14) in the present situation, the differential equation for pore pressure is multiplied by ; in the non-dimensionalization of the coordinates that occurs in the equivalent to the step from equation (5.16) to equation (5.18) the differential equation for pore pressure is multiplied by . Thus the term *ϕ*^{V}*γ*^{L}/*Λ*^{V}*ρ*_{f0} in the PV pore pressure differential equation equivalent (5.13) takes the form , which appears in equation (7.2) after these multiplications. This supply term to the PV from the PLC may be rewritten as follows by using equations (6.5) and (6.7):
7.3
Substituting equation (7.3) into equation (7.2) the differential equation for the PV pore pressure field becomes
7.4
where using equation (5.15) for *c*^{L} and *c*^{V} we find that
7.5
The three terms introduced in equation (7.5) arise from the fact that the supply term from the PLC to the PV (or vice versa) modifies the two frequencies of pore pressure oscillation in the PLC and in the PV as well as modifying the factor multiplying the applied strain driving the motion. Note that *Φ*^{LVω} depends upon both the PV and PLC properties and is therefore called the coupling parameter. The coupling parameter *Φ*^{LVω} is small for bone tissue. A plot of *Φ*^{LVω} against *ω* for bone tissue is shown in figure 4. As one can see from figure 4, the value of *Φ*^{LVω} starts at 0.01 near zero frequency and quickly drops to 0.001 as frequency increases.

Thus new dimensionless frequency of oscillation in the PLC is *ω* and new dimensionless frequency of oscillation in the PV is . Using equations (5.15), (5.17), (5.27) and (7.1) these last two equations may be rewritten as
7.6
and where . The numerical value of and are very much less than 1 for the parameter values of the PLC and the PV as one can see from the plot in figure 5. It follows then from equation (7.6) that and . As and represent the altered frequencies of oscillation of the pore pressure in the PV from the driving frequencies *ω* and *Ω*, respectively, it follows that the alteration in the frequencies *ω* and *Ω* is quite small and near zero frequency and decreases with increasing frequency. Thus the principal frequencies of fluid pore pressure oscillation in the PV are little changed.

The solution to the differential equation for the pore pressure (7.2) is
7.7
This solution is subjected to the following boundary conditions at *β*=1 and *A*:
7.8
where the first condition requires that there be no flow across the periosteum and the second condition requires that the pore pressure at the endosteal surface equals the blood pressure in the medullary canal. The application of these boundary conditions to equation (7.7) yields the following formulas for the previously undetermined coefficients in equation (7.7):
7.9
where *Π*(*ϖ*:*λ*,α) is defined in equation (5.23) above and its properties are documented in equations (5.24), (5.25) and (5.26). Placing the results (7.9) back into equation (7.7) the formula for the pore pressure in the PV takes the form
7.10
where is defined in equation (5.30).

In appendix B the boundary conditions on the stress component *T*_{rr} are used to determine and to eliminate the constants of integration, *c*^{Vω} and *c*^{VΩ}. The values of the formerly unknown coefficients and are determined to be
7.11
where
and
7.12
where *ϖ* stands for either and the superscript LV indicates a dependence upon the properties of both porosities while the superscript V indicates a dependence upon the properties of only the PV
7.13
where equations (5.9) and (5.17) have been used to simplify *H*^{V}.

In order to evaluate *p*^{Vω}(1,*t*) and *p*^{VΩ}(1,*t*) the formulas (7.11) for and are substituted into equation (7.10) for the pore pressure in the PV and the pore pressure in the PV is then evaluated at *β*=1. The result leads to the following expressions for *p*^{Vω}(1,*t*) and *p*^{VΩ}(1,*t*):
7.14

The pore fluid pressure field in the PV (7.10) has the following form when the relations (7.11) and (7.14) are inserted into it:
where
and
7.15
These expressions relate the pore pressure field in the PV to the two driving forces, the cyclic mechanical straining with amplitude ϵ_{0} and a frequency *ω* and the blood pressure with amplitude and a frequency *Ω*. A plot of the pore pressure in the PV, *p*^{V}(*β*,*t*) given by equation (7.15) above, is shown in figure 6 as a function of the dimensionless radial coordinate of the model bone *β* at various temporal points in the cyclic oscillation. Figure 7 is again a plot of the same equation for the pore pressure in the PV, but with the PV permeability reduced by two orders of magnitude, 6.35×10^{−8} m^{2}, from that employed for figure 6, 6.35×10^{−10} m^{2}. Notice from figure 6 that the peak pressure is near 15 MPa for a PV permeability of 6.35×10^{−10} m^{2} and from figure 7 it is about 9.5 kPa for a PV permeability of 6.35×10^{−8} m^{2}. As 8 kPa is about 60 mm Hg, the peak blood pressure, a PV pore pressure in the range of 8 kPa is a more reasonable result for several reasons. First, as a comparison of figures 6 and 7 shows the gradient of the PV pore pressure across the bone cortex is smaller in the case of the higher PV permeability; the bone pressure is near uniform across the cortex. Second, since the gradient is small at the endosteal surface, there is less flow in and out of the medullary canal; there is no known physiological advantage for interstitial flow in and out of the medullary canal. Lastly the pore pressure in the PV cannot exceed a pore pressure that would collapse the blood vessels present for any significant length of time. Figure 8 is a plot of the pore pressure due to blood pressure in the PV, *p*^{VΩ}(*β*,*t*) given by the second part of equation (7.15) as a function of the dimensionless radial coordinate of the model bone *β*, at various temporal points in the cyclic oscillation. Figure 8*a* is at a scale that obscures the small spatial gradient in the blood pressure; figure 8*b* expands the scale for a single time point in the oscillation so that the small gradient is apparent. These considerations concerning the plots in figures 6 and 7 suggest a PV permeability greater than 10^{−9} m^{2} and perhaps a little lower than 10^{−8} m^{2}.

Our interpretation of results presented in a recent paper by Goulet *et al.* (2008) reinforces the point made at the end of the previous paragraph concerning a PV permeability greater than 10^{−9} m^{2} and perhaps a little lower than 10^{−8} m^{2}. In that work a finite element model whose geometry was generated from a quantitative computed tomography scan of a section of human tibia was employed. The fluid velocities and the pore fluid pressure in the PLC and in the PLC + PV were calculated for a specified applied loading consisting of combined axial loading and bending. Figure 4*b* in that paper shows PV pore pressures in the MPa range. We think that this pore pressure is too high because of likely constriction of the blood vessels in the PV and the significant volume of interstitial fluid that will be drained into the medullary canal. We think that their model is likely correct and that the high PV pore pressures in the MPa range are due to the authors’ use of a PV permeability of 10^{−14} m^{2}, which we now consider, on the basis of the analysis reported here, to be much too small for the PV. As we began work on this paper we thought that a PV permeability of 10^{−14} m^{2} was possible. Our original version of the figure that is now figure 6 employed a PV permeability of 10^{−14} m^{2} and our results were similar to those of Goulet *et al.* (2008) in predicting an unrealistically high pore pressure in the PV. It is the analysis of the last paragraph that changed our minds. A subsequent literature evaluation of the references reporting the low PV permeabilities revealed that they were very rough estimates. We return to this subject in the next section.

Please note that the functions *p*^{Vω}(*β*,*t*) and *p*^{VΩ}(*β*,*t*) specified in equation (7.15) are the two functions required to complete the solution for equation (5.33), for *p*^{L}(*λ*,*t*).

A plot of the cyclic pore pressure in the PLC, *p*^{L}(*λ*,*t*) given by equation (5.33), is shown in figure 9 as a function of time for the case when the pore pressure in the PV is neglected (set equal to zero) and the case when the pore pressure in the PV is given by equation (7.15). The difference in the peaks of the two cyclic curves depends upon the permeability differences. They are both plotted for the same frequency in figure 9 and the peaks would separate if the frequencies were different.

## 8. Implications for the determination of the permeabilities

The measurement of the permeability of the PV has not been accomplished with sufficient accuracy to date, primarily because of the topological intertwining of the PV with the PLC, as we noted in §2. The final formulas for the pore fluid pressure in the PLC, equation (5.33), and the PV, equation (7.15), show that these pore pressures both depend upon the permeability of the own porosity as well as upon the permeability of the other porosity. Note that the formulas (5.33) and (7.15) for the pore pressures depend upon different dimensionless frequencies, each of which is non-dimensionalized by use of a permeability (cf. equations (5.17), (7.1) and (7.5)). Thus these pore pressures depend upon the two permeabilities only through their dependence upon the various dimensionless frequencies.

The PV permeabilities reported in the past are thought to be significant underestimates because these previous values represent a PLC and PV lumped permeability measurement rather than a PV or PLC measurement alone. The fact that, in general, the pore pressure in each porosity depends upon the per-meability of the other porosity means that special strategies should be designed for the experimental determination of each porosity’s permeability. Traditional permeability measurement techniques based on Darcy’s technique of measuring the volume of flow per unit area per unit time across a porous layer and dividing by the pore pressure gradient across the layer will require modification for porous media structured hierarchically.

If the permeability is measured *in vitro*, as will likely be the case, there will be no blood pressure and and, by equation (7.15), *p*^{VΩ}(*β*,*t*) will be zero and all the frequencies denoted using an upper case *Ω* will be zero. However, as may be seen from (5.33) the pore pressure in the PLC still depends upon the pore pressure in the PV through *p*^{Vω}(*β*) given in equation (7.15) and the pore pressure in the PV still depends upon the permeability of the PLC through coupling parameter *Φ*^{LVω}.

The two porosities of bone tissue occupy the same three-dimensional volume of bone tissue. In the cross section of the bone shown in figure 1 the PV channels are the central lumens of the osteons and the region immediately surrounding the osteon is part of the PLC. Thus it is impossible to obtain a reasonable sized specimen of bone tissue that contains some PV porosity without also containing some PLC porosity. This topological fact makes it difficult to design an experiment to measure the permeability of the PV. The present work contains the analytical basis for a number of different approaches to this measurement technique for the measurement of the permeably of the PV, and we will be pursuing this topic in the near future. From this discussion it is clear that previous measurements of the permeability of the PV were not accurate as they were based on models that did not include the interchange of interstitial fluid with the PLC. That is one reason why we have not honoured the reported measurements of the permeability of the PV that are in the literature, permeabilities as low as 10^{−14} m^{2}. As noted in the previous section the results of this study suggest a PV permeability greater than 10^{−9} m^{2} and perhaps a little lower than 10^{−8} m^{2}. The second reason is that permeabilities as low as 10^{−14} m^{2} predict PV pore pressures that are too high as this work, and that of Goulet *et al.* (2008), show.

We turn now to the question of the measurement of the permeability of the PLC. The entire contents of a single osteon considered as a hollow circular cylinder is composed of the PLC porosity and no PV porosity, and thus may be used to measure the permeably of the PLC. This has been accomplished (Gailani *et al.* in press) by isolating single osteons, imposing a compressive stress on them and measuring the time dependent axial strain. The PLC permeability determined was of the order of 10^{−24} m^{2}. The analysis of these experimental data employed the poroelastic model of this situation analysed by Gailani & Cowin (2008).

## 9. Discussion

The long-term objective of the work, of which this paper is a part, is to develop a model of all fluid transport in bone, particularly that due to mechanical loading and blood pressure. The result presented here is the basic piece of this model. It addresses the question of pore fluid movement from the perivascular region in bone to the bone cells housed in the lacunae of the mineralized matrix. Knowledge of this fluid movement is necessary for the determination of the nutrient transport and the effect of the fluid drag on the cell processes of osteocytes in the canaliculi, which is thought to be a principal mechanotransduction mechanism. Quantification of the fluid drag on the cell processes is necessary for the evaluation of osteocytic mechanosensitivity (Weinbaum *et al.* 1994; Cowin *et al.* 1995).

In this contribution the governing equations for the theory of poroelastic materials with hierarchical pore space architecture and compressible constituents undergoing small deformations are developed. These equations are applied to the problem of determining the exchange of pore fluid between the PV and the PLC in bone tissue due to cyclic mechanical loading and blood pressure oscillations. A formula for the volume of fluid that moves between the PLC and PV in a cyclic loading is obtained as a function of the cyclic mechanical loading and blood pressure oscillations. Formulas for the fluid pore pressure in both the PLC and the PV are obtained as a function of the two driving forces, the cyclic mechanical straining and the blood pressure, both with specified amplitude and frequency. A general observation is that the PV and the PLC are greatly different. The PLC is a relatively high-pressure domain characterized by a relatively small pore size and no vascular vessels. On the other hand the PV is a relatively low-pressure domain characterized by a relatively large pore size that contains vascular vessels. The model shows that, although the two porosities exchange fluids, their influence on one another in terms of pressure change is small. The results of this study also suggest a PV permeability greater than 10^{−9} m^{2} and perhaps lower than 10^{−8} m^{2}.

A number of papers have considered the pore pressure oscillations in the PLC of an osteon without consideration of the activity in the PV. These contributions include those of Weinbaum *et al.* (1994), Cowin *et al.* (1995), Swan *et al.* (2003), Rémond & Naili (2005) and Rémond *et al.* (2008). Two recent papers have computationally addressed the question of the two interacting porosities. Fornells *et al.* (2007) considered a mixture theory-based dual porosity model of the PV and the PLC; the model was based on the continuum model of Khalili-Naghadesh & Valliappan (1996). In this model the two porosities are considered to be interacting on the same level as opposed to the model in this paper where the two porosities are arranged hierarchically. As the authors point out their averaging process is associated with a typical RVE length that averages over the details of local pore pressure fields. In the paper of Goulet *et al.* (2008), discussed in the previous section, a finite element model whose geometry was generated from a quantitative computed tomography scan of a section of human tibia was employed. The fluid velocities and the pore fluid pressure in the PLC and in the PLC+PV were calculated for a specified applied loading consisting of combined axial loading and bending. The values they obtained for the pore pressure in the PV were discussed in §7.

The limitations of the present work include the assumption of right circular cylindrical geometry for both the whole long bone cross section and the osteon and the inadequacy of the lymphatics modelling. The problem of the idealized modelling may be addressed with computational models and the improvement of the lymphatics modelling depends upon our achieving a greater understanding of bone lymphatics. It may be reasonable to try to model the readsorption of interstitial fluid into the venous system as those vessels leave the bone. It is possible that this work will lead to improvements in both of these topics. For the bone lymphatics problem the model provided could be used as the analytical basis for the determination of the location of all the interstitial bone fluid in a whole bone experiment at different times with unexpected imbalances leading to the identification of the bone lymphatics. Some computational approaches have modelled both the PV and PLC with some degree of anatomical accuracy, but computational complexity. The hierarchical approach employed here can easily be incorporated in a computational approach; one needs to model only the PV and develop a separate subprogram to determine the interstitial fluid supply interchange between the PLC to the PV depending upon the applied loading and the pore pressure differences between the two porosities.

Planned extensions of this work include the determination of the solute transport between the PV and the PLC and extension of the mechanical loading considered here, the cyclic purely axial loading of a long bone, to the consideration of the cyclic combined axial and bending loading of a long bone. The extension is not simple because differential pressure between PV locations in a bone cross section will cause blood to move within the vessels in the cross section, as well as the interstitial fluid. In the case of cyclic purely axial loading of a long bone considered here, the pore pressure in the PV was homogeneous at any time and did not induce blood movement within the vessels.

## Acknowledgements

The NSF (NSF/CUNY AGEP #0450360), the PSC-CUNY Research Award Program of the City University of New York (PSC-CUNY-60014-38-39) and CUNY-LSAMP supported this work.

## Appendix A. Determination of and

Substituting expression (5.30) for the pressure field into equation (5.5) permits the solution for the displacement field *u*^{L} that can then be used to solve for the stress components (5.3). Applying boundary conditions to the stress component *T*_{rr} will determine and the constants of integration, *c*^{Lω} and *c*^{LΩ}, that occur in the solution of equation (5.5). First equation (5.5) is rewritten in terms of the dimensionless radial variable, *λ*=*r*/*r*_{0}; thus,
A 1
which upon substitution from equation (5.29) for *p*^{L}(*λ*,*t*) yields
A 2
which takes the following form when integrated:
A 3
where is given by equation (5.26). Using the above results, the radial stress *T*_{rr} is now calculated. The formula (5.3) for *T*_{rr} is rewritten as
A 4
and we note that
A 5
and
A 6
The substitution of equations (1.5) and (1.6) into equation (1.4) yields the desired formula for the radial stress
A 7
Requiring that the radial stress vanish on the inner surface of the hollow cylinder then imposes the following restriction:
A 8

Multiplying equation (1.8) by *a*^{2}/*λ*^{2} and subtracting it from equation (1.7), the unknown constants *c*^{Lω} and *c*^{LΩ} are eliminated from equation (1.7) resulting in the following representation for :
A 9
If we now require that the effective radial stress vanishes as the outer surface of the osteon (*λ*=1), the values of the formerly unknown coefficients and are determined and recorded as equation (5.31).

## Appendix B. Determination of and

The solution for the PV pressure field (7.10) now permits the solution for the displacement field *u*^{V} that can then be used to solve for the PV stress components (5.3). Applying boundary conditions to the stress component *T*_{rr} will determine and the constants of integration, *c*^{Vω} and *c*^{VΩ}, that occur in the solution. First equation (5.5) is rewritten in terms of the dimensionless radial variable *β*; thus,
B 1
which upon substitution from equation (7.10) for *p*^{V}(*β*,*t*) yields the following:
B 2
Integration of the above equation using equations (5.24) and (5.26) yields
B 3
where is given by equation (5.26) with the small ‘*a*’ replaced by capital ‘*A*’ and where *c*^{Vω} and *c*^{VΩ} are the constants of integration. Using the above results, the radial stress *T*_{rr} is now calculated. The formula (5.3) for *T*_{rr} is rewritten as
B 4
and we note that
B 5
and
B 6
The substitution of equations (2.5) and (2.6) into equation (2.4) yields the desired formula for the radial stress
B 7
Requiring that the radial stress vanishes on the inner surface of the hollow cylinder, *β*=*A*, then imposes the following restriction:
B 8
where the result , following from equation (5.30), has been employed. Multiplying the result (2.8) by *A*^{2}/*β*^{2} and subtracting it from equation (2.7) the unknown constants *c*^{Vω} and *c*^{VΩ} are eliminated from equation (2.7) resulting in the following representation for :
B 9
If we now set *p*^{V}(*β*,*t*)=*p*^{Vω}(*β*)e^{iωt}+*p*^{VΩ}(*β*)e^{iΩt} and require the effective radial stress *T*_{rr} to vanish as the outer surface of the periosteum (*β*=1), and then employ the results and that follow from equation (5.27), the following condition for the values of the formerly unknown coefficients and is obtained:
B 10
The values of the coefficients and are recorded as equation (7.11).

## Footnotes

One contribution of 12 to a Theme Issue ‘Mechanics in biology: cells and tissues’.

- © 2009 The Royal Society