## Abstract

Recently, the mesoscale of cortical bone has been given particular attention in association with novel experimental techniques such as nanoindentation, micro-computed X-ray tomography and quantitative scanning acoustic microscopy (SAM). A need has emerged for reliable mathematical models to interpret the related microscopic and mesoscopic data in terms of effective elastic properties. In this work, a new model of cortical bone elasticity is developed and used to assess the influence of mesoscale porosity on the induced anisotropy of the material. Only the largest pores (Haversian canals and resorption cavities), characteristic of the mesoscale, are considered. The input parameters of the model are derived from typical mesoscale experimental data (e.g. SAM data). We use the method of asymptotic homogenization to determine the local effective elastic properties by modelling the propagation of low-frequency elastic waves through an idealized material that models the local mesostructure. We use a novel solution of the cell problem developed by Parnell & Abrahams. This solution is stable for the physiological range of variation of mesoscopic porosity and elasticity found in bone. Results are computed efficiently (in seconds) and the solutions can be implemented easily by other workers. Parametric studies are performed in order to assess the influence of mesoscopic porosity, the assumptions regarding the material *inside* the mesoscale pores (drained or undrained bone) and the shape of pores. Results are shown to be in good qualitative agreement with existing schemes and we describe the potential of the scheme for future use in modelling more complex microstructures for cortical bone. In particular, the scheme is shown to be a useful tool with which to predict the qualitative changes in anisotropy due to variations in the structure at the mesoscale.

## 1. Introduction

Cortical bone is a highly organized hard tissue that represents approximately 80% of the skeletal mass in the human adult. It can be conveniently described as a two-phase composite material: a soft phase, i.e. pores, containing fluid and soft tissues such as cells, blood vessels and nerves distributed inside a complex dense matrix phase. Porosity is distributed over several length scales: the diameter of pores ranges from a few micrometres to several hundreds of micrometres (Martin *et al*. 1998). Only the two largest pore types, resorption cavities (approx. 50–200 μm) and Haversian canals (approx. 50 μm), contribute to the so-called *mesoscale structure* (or mesostructure) and hence mesoscale porosity that is characteristic of the higher level of organization in bone. Mesoscopic pores are roughly aligned with the longitudinal direction of the bone (Currey 2002; Cooper *et al*. 2007). The cross section of Haversian pores is roughly circular, while that of the resorption cavities is often of a complex shape due to the association with bone remodelling activity.

The mesoscale porosity is embedded in a bony matrix, which is a mixture of the smaller pores such as Volkmann's canals, osteocyte lacunae, etc. (at least one order of magnitude smaller than the mesoscopic pores) and mineralized collagen fibrils. With standard techniques (micro-computed tomography, acoustic microscopy, microradiographs, etc.), given a typical resolution of approximately 10–20 μm, the smaller pores cannot be distinguished from the mineralized collagen matrix. Accordingly, the *local* mesoscale porosity in samples taken from cortical bone cross sections is that which is effectively assessed with standard techniques. It is defined as that due to the Haversian canals and resorption cavities and lies in the range of 5–15 per cent (Bousson *et al*. 2001; Dong & Guo 2004; Baron *et al*. 2007). However, note that the macroscopic porosity defined as an average of the mesoscopic porosity over local representative volume element (RVE) locations varies much less and is usually in the range of 7–9 per cent. In this paper, we consider the influence of *mesoscopic porosity* (i.e. that due to Haversian canals and resorption cavities) so that, in fact, the lower level porosity associated with smaller pores contributes to the overall (homogenized) elastic properties of the bony matrix phase.

Recently, the mesoscale of cortical bone has been given particular attention (Hoc *et al*. 2006; Verhulp *et al*. 2006; Grimal *et al*. 2008), in association with novel experimental techniques such as nanoindentation (Zysset *et al*. 1999) and quantitative scanning acoustic microscopy (SAM; Raum *et al*. 2003). These techniques give access to the elastic properties of the bony matrix mentioned above. Furthermore, SAM and various microscopic imaging techniques give access to the structural distribution of mesoscopic porosity.

There is a need for reliable methods to interpret SAM and nanoindentation data in terms of macroscopic elasticity. In particular, phenotyping, the investigation of structure–function relationships and remodelling, and numerical modelling of bone response to mechanical loads at various scales would all benefit from such methods. One of the main difficulties here is to account properly for the material anisotropy. This arises from two sources: the mesoscopic porosity and the *anisotropy of the matrix* itself, mainly due to the orientation of the collagen fibrils (Weiner *et al*. 1999). The respective contributions of the two sources of anisotropy have thus far not been made clear. Evidently, it is difficult to design experiments for this purpose, although there are particular bone tissues with specific organization, which can provide some information with regard to this issue (Currey & Zioupos 2001). Mathematical models using homogenization techniques and micromechanical methods can give excellent insight since they explicitly relate the effective properties to the meso/microstructural properties and in the final instance they give relatively accurate predictions of effective moduli.

The purpose of this paper is to introduce a new model of cortical bone elasticity and investigate the relationships between the mesoscale information and the homogenized local elastic properties. The mesoscale information of interest here is the porosity, the pore shape and the distribution and elasticity of the mineralized matrix. Specific attention will be given to the sources of overall anisotropy. Pathologies and drugs produce effects that may in some cases *independently* alter the anisotropy of the matrix (e.g. due to modifications in the collagen and mineral properties) and the mesoscopic porosity due to perturbations of remodelling activity. For this reason and for the sake of clarity, it is proposed that anisotropy be investigated in two stages. In this paper, we shall consider the influence of the porosity when the matrix is considered to be *isotropic*. The influence of matrix anisotropy will be discussed in a forthcoming paper.

Strictly speaking, bone is a poroelastic material and therefore a relevant framework with which to study the effective behaviour of the medium is that presented by Cowin (1999) who applied the classical Biot poroelasticity theory to study deformation driven fluid motion. However, in many situations where the focus is on the effective elastic behaviour and not on the fluid motion or distribution of stresses inside the pores, the fluid can be considered to be static. Thus, we treat the fluid as a static phase within the elastic medium; alternatively, this is equivalent to considering *closed pores* within the material.

Mechanical models proposed in the past two decades have provided a great deal of insight into the micro–meso–macro relationships. Among others, Crolet *et al*. (1993) and Aoubiza *et al*. (1996) considered periodic microstructure in the framework of the asymptotic theory of homogenization; Sevostianov & Kachanov (2000) and Hellmich & Ulm (2004) proposed models in the spirit of micromechanics, based on the Eshelby tensor (Eshelby 1957), using Mori–Tanaka (Benveniste 1987) and self-consistent schemes (Sabina & Willis 1988; Kanaun & Levin 2003). All these models consider elementary shapes for the pores and simple mesostructure. Recently, Baron *et al*. (2007) and Grimal *et al*. (2008) have developed numerical models to compute the effective properties of cortical bone with realistic mesoscale structures derived from images. However, these latter models are not appropriate to use when performing systematic parametric studies due to their high computational cost. It appears that existing models do not readily lend themselves to a fast and reliable computation of macroscopic properties based on experimental mesoscale information. Furthermore, it appears that the following two fundamental points should be revisited and explored in greater depth.

Existing models assume very specific constitutive behaviour of the material

*inside*the pores. For example, Aoubiza*et al*. (1996) and Sevostianov & Kachanov (2000) considered empty pores (drained bone), which do not possess stiffness. On the other hand, for example, Hellmich & Ulm (2004) stated that because fluids saturate the pores, they must possess significant stiffness and should therefore be considered as*hard*pores (undrained bone). The influence of these hypotheses on the overall anisotropy is not clear. We note that*in vitro*measurements of bone properties have been performed on both undrained and also*dried*samples (e.g. Lang 1970; Yoon & Katz 1976*b*), where some of the fluids have been removed. Models of both undrained and drained bones are therefore relevant to interpret such data. Furthermore, pores are often considered to be of circular cross section but how does a non-circular shape affect predictions?The level of approximation of a modelling scheme is very often unclear and there is a great need for a comparison between the predictions given by alternative schemes. At present, it is unknown which of the schemes (Mori–Tanaka, self-consistent, asymptotic homogenization, numerical models, etc.) is most appropriate. Furthermore, the accuracy of the prediction may often depend on the level of porosity and the phase contrast. In particular, it is important that the model be stable when porosity and phase contrast vary within their physiological range.

In the present paper, we use the method of asymptotic homogenization (AS to denote asymptotic scheme), assuming that the pores in the local region of interest are periodically distributed within the matrix material, specifically on a hexagonal lattice, which leads to a transversely isotropic material for pores with circular cross section. This symmetry is thought to be a reasonable approximation of the effective behaviour of cortical bone (Yoon & Katz 1976*a*,*b*). Asymptotic homogenization is one of the most successful schemes for determining the effective moduli of periodic media (Bakhvalov & Panasenko 1989; Parton & Kudryavtsev 1993). The approach exploits the separation of scales within the composite material and uses the method of multiple scales in order to derive a homogenized equation governing the effective behaviour of the material in question. Some of the advantages of the scheme are as follows.

Given the meso/microstructure within the periodic cell, the resulting (induced) anisotropy arises

*naturally*in the governing homogenized equation.Asymptotic expansions of displacements allow the nature of meso/microstructural fields (stresses, etc.) to be understood.

It provides an algorithmic, rather than an

*ad hoc,*scheme with which one can determine effective properties.Effective properties are defined in terms of a solution to the so-called

*cell**problem*; solutions to this are, in general, found numerically but semi-analytical solutions can be derived in the case of two-dimensional fibre-reinforced composites (Sabina*et al*. 2002; Parnell & Abrahams 2006, 2008).

The current approach is novel and unique due to the solution scheme of the periodic cell problem. It is based on the work developed by Parnell & Abrahams (2006, 2008) and was applied successfully in the context of industrial fibre-reinforced composites. In particular, the method appears to be stable even for high porosities and high-contrast phases as in the limiting cases of rigid fibres and voids. We note that since solutions are found in terms of asymptotic expansions, we are aware of the order of error made by assuming a leading-order solution. Furthermore, the scheme does not introduce any *ad hoc* assumptions, which are required in many other methods.

In §2, we present the assumptions of our model. In particular, we introduce the idea of a periodic cell, which defines the local structure of the cross section of cortical bone. Next, in §3 we describe the method of asymptotic homogenization with reference to the works of Parnell & Abrahams (2006, 2008), in particular focusing on the cell problem and the form of effective elastic moduli. Sections 4 and 5 focus on analysing the two fundamental points mentioned above, i.e. comparing the effect of pore properties and shape on effective elastic properties and assessing the efficacy of current modelling schemes in predicting the effective moduli of cortical bone. Specifically, a comparison between the results obtained by the proposed AS and extant methods is presented in §4 in which we assess the influence of the pore properties, i.e. whether the bone is drained or undrained and also their shape, in particular by describing the situation when the pore cross section is elliptical rather than circular. We close with a discussion in §5 and conclusions and plans for future work in §6.

## 2. Assumptions of the model

We shall now discuss how we model cortical bone as a two-phase composite material. Let us employ Cartesian coordinates where the *x*_{3} axis runs parallel to the long axis of the bone. As a first model, the mesoscale pores (Haversian canals and resorption cavities) are modelled as identical cylinders of infinite extent in the *x*_{3} direction and we assume that they are distributed on a hexagonal lattice in the *x*_{1}*x*_{2} plane. Note that in this paper we consider pores of circular *and* elliptical cross section. These pores represent Haversian canals or resorption cavities that are roughly aligned in cortical bone (Cooper *et al*. 2007). Resorption cavities tend to be less circular and therefore the analysis of elliptical cross sections is intended to simulate this effect. Furthermore, resorption cavities tend to possess a conical shape, but since the important notion is the *local property* (small length in the longitudinal direction) the cylindrical assumption should be reasonable. This kind of local analysis is required in any case since samples change porosity, for example, in different regions of the bone shaft. Although here the method is presented for only one pore inside the periodic cell, it can, in fact, deal with many pores inside the cell (Parnell & Abrahams 2006, 2008), which would be ideal when dealing with complex microstructure. The motivation behind this work is to maintain simplicity, however, and thus the restriction to a single pore.

The mesoscale pores are surrounded by bony matrix tissue which we assume to be *isotropic*. In reality, the matrix material is an inhomogeneous, anisotropic material, but the focus here is on the anisotropy due to the mesostructure and not due to the tissue itself. The hexagonal periodic cell shown in figure 1 contains a single pore and defines the mesoscale structure in the *x*_{1}*x*_{2} plane and (for circular cross sections) ensures transverse isotropy, which is the material symmetry commonly observed in human bone (Yoon & Katz 1976*a*,*b*). This periodic model material of infinite extent is thus intended to give a snapshot of the (local) effective elastic properties of bone at a specific RVE location at some point in the cortical bone cross section, where the RVE is defined by its porosity. This is illustrated in figure 2 where we show how the mesostructure within the local region of the cortical bone cross section is first replaced by an idealized hexagonal structure and then homogenized by using the concept of an infinite periodic medium as described above.

The two phases are assumed to be purely elastic and isotropic and are defined by their respective Lamé constants *μ*_{0}, *λ*_{0} and *μ*_{1}, *λ*_{1}, where *μ* is the shear modulus and *λ*=*κ*−2*μ*/3, where *κ* is the bulk modulus. The subscript 0 refers to the matrix material while 1 refers to the pore. Note that fluid fills the pore in the undrained case. We discuss how we describe this case later on. Furthermore, we define the length scale of the mesoscale (for example, distance between pores) as *q*. Within the periodic cell, the domain of the fibre is denoted by *D*_{1} and the matrix by *D*_{0}, their respective areas being denoted by |*D*_{1}| and |*D*_{0}| and we define as the total cell area. The interface between matrix and pore is defined as *∂D*_{1}.

## 3. Asymptotic homogenization and the cell problem

We choose to use the theory developed by Parnell & Abrahams (2006, 2008) in order to determine the effective elastic properties by modelling the propagation of low-frequency elastic waves through the material. We shall summarize the findings of the theory here, but we refer the reader to these papers for more in-depth findings and details, particularly with regard to the solution of the cell problem. Neglecting body forces and on defining , Navier's equations for time harmonic waves of frequency *ω* and for * x*∉

*∂D*

_{1}are given by (Graff 1991)(3.1)where we adopt Einstein's convention of summation over repeated

*subscripts*. Dividing by the matrix shear modulus

*μ*

_{0}and non-dimensionalizing length scales on

*q*, (3.1) becomes(3.2)for

*∉*

**x***∂D*

_{1}, where(3.3)(3.4)and(3.5)

In (3.2), we have defined the non-dimensional parameter , where is the shear wavenumber in the matrix material. Since we assume that we are in the low-frequency regime, where the wavelength of the propagating waves is much larger than the mesoscale *q*, we can assume that *ϵ*≪1. Therefore, a so-called separation of scales exists and we can perform an analysis based on asymptotic homogenization. Note that we can always rescale on some other wavenumber, for example, *κ*_{0} defined by associated with longitudinal waves in the matrix, but the important notion is that the newly defined non-dimensional parameter still satisfies . We assume boundary conditions on the pore/matrix interface ∂*D*_{1} of continuity of displacement and traction, i.e.(3.6)where *σ*_{ij} is the Cauchy stress tensor; *n*_{j} is the *j*th component of the normal to ∂*D*_{1} (pointing *into* the matrix); and where we have used the notation to denote the jump in the function *f*(*z*) across ∂*D*_{1}.

Asymptotic homogenization proceeds by defining the multiple scale variables,(3.7)(3.8)where is some expansion in *ϵ* and *L*_{2}∈. On introducing and we note from (3.7) and (3.8) that these are short and long length scales, respectively. We see that no mesovariable is required in the long bone axis direction and thus (3.7) and (3.8) imply that(3.9)(3.10)Note that in the following, Greek indices can run from 1 to 2 whereas Latin indices run from 1 to 3.

An expansion in *ϵ* is performed for the displacements(3.11)and since the material is doubly periodic (with hexagonal symmetry) with respect to * ξ*, we insist that each

*u*

_{kj}is doubly periodic in

*on the hexagonal lattice.*

**ξ**Substituting (3.7)–(3.11) into the governing equation (3.2) and boundary condition (3.6) and equating orders in *ϵ* we obtain a hierarchy of problems, one associated with each order in *ϵ* (Bakhvalov & Panasenko 1989). As is usual in the theory of homogenization for elastodynamics, the *O*(1) problem shows that the leading-order displacement field is (explicitly) independent of * ξ*, so that(3.12)At

*O*(

*ϵ*), it turns out that it is advantageous to write the solution in the separable form(3.13)This form then gives rise to the resulting

*cell problem*for . As we shall show in (3.18)–(3.24) below, since the superscripts

*p*and

*m*merely alter the

*forcing*in the boundary conditions of the cell problem, it is useful to employ the bold font notation . For isotropic phases, on denoting the derivative of

**N**

_{k}with respect to

*ξ*

_{α}by

**N**

_{k,}

_{α}, it turns out that the governing equations for the cell problems (in the

*r*th phase) are of the form(3.14)(3.15)and(3.16)for

*r*=0,1. Boundary conditions are(3.17)and(3.18)(3.19)and(3.20)We have also defined(3.21)(3.22)(3.23)and(3.24)Note that the cell problem for

**N**

_{3}is decoupled from that of

**N**

_{1}and

**N**

_{2}. The former is associated with

*anti-plane*motion whereas the latter is associated with

*in-plane*motion.

### 3.1 The homogenized wave equation

The *O*(*ϵ*^{2}) problem is used to determine the effective wave equation governing the leading-order displacement *U*_{k}(* z*). This is achieved by integrating the governing equation at

*O*(

*ϵ*

^{2}) over the periodic cell, employing Green's theorem and imposing the necessary boundary conditions and (hexagonal) double periodicity in

*.*

**ξ***Restricting attention to circular pores*, this gives rise to the following effective wave equations (Parnell & Abrahams 2006, 2008):(3.25)(3.26)and(3.27)where(3.28)and where

*U*

_{k,j}denotes the derivative of

*U*

_{k}with respect to

*z*

_{j}. The form of the coefficients will be given shortly.

The form of these effective wave equations, including the relationship (3.28) between and is identical to that which governs wave propagation in a transversely isotropic *homogeneous* elastic medium with elastic moduli (using the engineering notation for ; see Sokolnikoff 1956). The coefficients therefore define the corresponding *five* effective elastic moduli of the transversely isotropic material. In order to completely define all five coefficients, we require the solution to the cell problem associated with only *three* of the forcings in the boundary conditions above. These are the in-plane problems with *p*=*m*=1 and *p*=*m*=3 and also the anti-plane problem with *p*=3, *m*=1.

The five effective properties are defined as(3.29)(3.30)(3.31)(3.32)and(3.33)where is the mesoscopic porosity of the bone (i.e. the porosity of the local region of interest) or, alternatively, in the language of composite materials, the volume fraction of the fibre phase. The coefficients *M*_{j} are determined by line integrals of around the cylindrical fibre boundary ∂*D*_{1}(3.34)(3.35)(3.36)and(3.37)where firstly we recall the notation for the (*p*,*m*) cell problem defined in (3.13) and we note that the local polar coordinate system , centred on the pore within the periodic cell has been introduced. Finally, is the normal to the fibre boundary, and thus since here this boundary is circular, we have *n*_{1}=cos*θ*, *n*_{2}=sin*θ*.

### 3.2 The solution to the cell problem

As described above, the effective properties for a transversely isotropic material are defined completely if we know the solution to the appropriate anti-plane and in-plane cell problems. These were solved in the papers by Parnell & Abrahams (2006, 2008). The methods used are based on complex variable theory and multipole expansions of doubly periodic functions, specially constructed to enable complex meso/microstructure to be modelled. Local solution expansions are posed inside the pore and in the matrix and we match these on the pore boundary ∂*D*_{1} using the appropriate boundary conditions for the cell problem. The system of equations is closed using the extra conditions that arise from imposing the condition of double periodicity in * ξ*. This is achieved by constructing multipole expansions in the matrix material, in terms of doubly periodic basis functions. These are then matched to the local solution expansions in the matrix, thus providing the necessary extra equations required to close the system.

In the case of the in-plane problem, the above process results in a simple linear system, which must be solved for vectors (bold fonts indicating the dependence on the (*p*,*m*) cell problem) and the solution on the boundary of a *circular* fibre of radius *R* takes the form(3.38)and(3.39)where Re and Im define the real and imaginary parts, respectively, and where we have defined . The additional terms denoted by +… are higher order terms associated with cos *nθ* and sin *nθ*, but we see from the form of (3.34)–(3.36) that these do not contribute to the expressions for the effective moduli due to orthogonality. Note that this is only the case for *circular* cylindrical fibres.

Similarly, in the anti-plane case we have to solve a linear system for the vector and we obtain the simple form(3.40)for the anti-plane displacement on the boundary of the fibre ∂*D*_{1}.

Programs to solve these problems and thus find the necessary coefficients in **N**_{1}, **N**_{2} and **N**_{3} above were written by the authors in Mathematica. This package permits very simple and concise algorithms to be used.

### 3.3 Pores of elliptical cross section

If pores are no longer circular, the boundary ∂*D*_{1} being described by, say,(3.41)then the induced anisotropy will in general no longer be transversely isotropic. In particular, for the case of *elliptical* cross sections, for example, the material becomes *orthotropic* since now the expression (3.28) does *not* hold and , etc. (see Parnell & Abrahams (2006, 2008) for more details). For ease of exposition, let us focus on just a small number of results. In particular, we note that in this orthotropic case,(3.42)and(3.43)where(3.44)and(3.45)where is the normal to the fibre boundary and note that there is no addition over the *superscripts* *j* on the r.h.s. of (3.44) and (3.45). As such, for this more complicated material, additional cell problems (in this case, (*p*,*m*)=(1,1) and (2,2) problems) must be solved for the distinct elastic properties. Similar results follow for the other elastic properties (Parnell & Abrahams 2006, 2008). Local expansions for the solutions inside the pores may be posed in the same form as in §3.2 for each of these individual cell problems as in (3.38)–(3.40), but now with *R*=*R*(*θ*) as in (3.41) and also note that due to the more complicated normal in the line integrals (3.44) and (3.45), all of the terms in these local expansions now contribute to the effective properties (3.42) and (3.43) via the line integrals (3.44) and (3.45).

## 4. Results

Following experimental and numerical work (Kabel *et al*. 1999; Zysset *et al*. 1999; Rho *et al*. 2002), we choose *ν*_{0}=0.3 for the Poisson ratio and *E*_{0}=22.5 GPa for the Young modulus of the bone matrix. The value for the Young modulus is consistent with experimental data obtained by nanoindentation (Zysset *et al*. 1999). Furthermore, this specific value was chosen here so that we can compare results with those obtained previously by Grimal *et al*. (2008). We shall consider the two alternative cases of hard and soft pores, corresponding to undrained and drained bone, respectively, as used in the literature by, for example, Sevostianov & Kachanov (2000) and Hellmich & Ulm (2004), respectively. In the case of undrained bone, we assume the pore is filled with static water. The properties of inviscid water are conveniently described by its bulk modulus *κ*=2.2 GPa and shear modulus *μ*=0 (Lide 1993). The values *E*_{1}=0.13 GPa and *ν*_{1}=0.49 correspond to an approximation of the above values for *κ* and *μ* for inviscid water and are the same values as those used in Grimal *et al*. (2008) with which we compare our results. For the case of drained bone (empty or soft pores), we take the limiting case of a cavity, i.e. .

We shall compare the results from the following three methods:

the AS presented above;

the Mori–Tanaka micromechanical estimates (Eshelby 1957; Benveniste 1987; Zaoui 1997); and

finite-element computations (Grimal

*et al*. 2008): effective properties were computed with prescribed displacement boundary conditions on local bone volumes (RVEs) with detailed real imaged mesostructures.

We also mention the results of Yoon & Katz (1976*a*,*b*), Ashman *et al*. (1984) and Rho (1996), which correspond to experiments on human cortical bone with contact transmission ultrasonic methods. Although we do not plot these directly (since their spread is large and porosity uncertain), it has been seen that they certainly lie within the range of the predictions by the above homogenization and micromechanical schemes.

We plot a selection of the local effective moduli of the transversely isotropic material as a function of mesoscale porosity, denoting those associated with undrained bone by and those for drained bone with a hat, i.e. . Note that each of the in the expressions (3.29)–(3.33) are scaled on *μ*_{0}. However, in order to plot results, we rescale the effective properties on their corresponding value in the *matrix* material so that all effective properties at *ϕ*=0 are unity. In §§4.1–4.4 we focus on pores with circular cross section and then in §4.5 we study the change in properties when the cross section becomes elliptical.

### 4.1 Comparison of the asymptotic scheme, finite-element computations and the Mori–Tanaka method

The predicted effective moduli are plotted in figures 3–8, focusing, in particular, on , and . Furthermore, is considered in §4.4 because it corresponds physically to the plane–strain shear modulus and therefore can be compared with the Hashin–Rosen bounds. Without exception, it is seen that the results of the AS coincide with the Mori–Tanaka method. It is evident that these schemes also agree relatively well with the finite-element computations, for both drained and undrained bone, especially given the assumptions made regarding the distribution of mesostructure.

We can provide a measure of the difference between the two schemes by calculating the following percentage difference:(4.1)where the superscripts ‘as’ and ‘fe’ denote AS and finite elements, respectively, and we evaluate this at three mesoscale porosity values, *ϕ*=0.05, 0.1 and 0.15. The resulting data are given in tables 1 and 2 for the two pore cases and we see that for the undrained case in particular, the percentage difference between the two schemes is small, except for . The reason for this latter point is known and is, in fact, related to the way that the finite-element computations are calculated.

### 4.2 Deviation from the arithmetic mean (the Voigt estimate)

The arithmetic mean approximation to the effective elastic moduli, for example, corresponds to the well-known Voigt (1889) estimation. In fact, this is also known as an *upper bound* on the elastic moduli as shown by Paul (1960). Note that the effective properties defined in (3.29)–(3.33) are written as an arithmetic mean plus an additional term. The magnitude of this additional term is of interest, since it allows us to determine the deviation of the effective property from its arithmetic mean. Thus, on referring to (3.32), (3.29), (3.30) and (3.28) and on defining(4.2)and(4.3)we assess results in the two cases of and relating to longitudinal extension and plain–strain shear, respectively, written as(4.4)and(4.5)and we plot _{3} and _{6} in figures 9 and 10, respectively. Note the similarity in discrepancy of the drained and undrained cases for the in-plane shear modulus. By contrast, they are further apart in the case of longitudinal extension, where the deviation from the arithmetic mean is much greater in the case of drained bone.

### 4.3 Influence of pore properties on anisotropy

The usual parameter studied in order to assess the change in anisotropy for transversely isotropic materials is the anisotropy ratio (AR) . We plot this as a function of porosity in figure 11. This figure shows qualitative agreement between the finite-element results and the AS with particularly good agreement in the undrained case, as is also seen in table 3 where the percentage difference is evaluated. The average (global) porosity in human cortical bone has been determined as 8.26 and 9.72 per cent for males and females, respectively (Bousson *et al*. 2001). For these porosity values, the AS predicts ARs of 1.10 and 1.14 for undrained and drained male bones, respectively, and 1.12 and 1.16 for female specimens. Note that these predictions are in line with those of Crolet *et al*. (2005). However, they are at odds with experimental results that indicate that the AR is usually in the range 1.3–1.5. This discrepancy may be fully explained by the fact that the anisotropy of the matrix material has been neglected (Crolet *et al*. 2005; Grimal *et al*. 2008).

### 4.4 Comparison with the Hashin–Rosen bounds

The physical relevancy of the homogenization scheme can be assessed upon comparison with the well-known variational bounds on fibre-reinforced composites found by Hashin & Rosen (1964). The closest bounds are obtained on the engineering moduli such as Young's moduli, Poisson's ratio, etc., rather than individual themselves, due to the way that these bounds are constructed via strain energy. We shall illustrate with two examples relating to the in-plane deformation, the plane–strain shear modulus and bulk modulus , defined by(4.6)We plot the results for in the cases of undrained and drained bones in figures 12 and 13, respectively, and for again in these cases in figures 14 and 15. It is clear that the finite-element computations will not always lie inside the bounds since these are calculations on local samples that do not have the perfect hexagonal periodicity. We see that the AS remains *inside* the bounds for all porosities. This is important since this means it gives physically valid predictions for a material of this symmetry. Note that for , since the bounds are so close, the bounds themselves give reasonable predictions of the effective properties. In the case of , however, the distance between the upper and lower bounds is much greater and thus the AS becomes particularly useful as a predictive scheme.

### 4.5 Influence of pore shape; pores of elliptical cross section

Real pores at the mesoscale level do not always have perfectly circular cross sections, particularly in the case of resorption cavities (as we indicated in figure 2). Thus, let us consider a single pore of elliptical cross section within a hexagonal periodic cell. We wish to assess how much this change in pore shape affects the local elastic properties. However, we must also mention that this specific case considered here will lead to a greater change in anisotropy than would be expected in cortical bone since, in a *true* sample, several of these non-circular cross-sectional pores would occur within a RVE, each having its own orientation and thus the induced orthotropy would then be fairly weak. Hence we can consider the cases discussed here to be fairly extreme situations and therefore we can usually assume that the true local effect will be somewhere between the case of circular pores and the elliptical case considered in this section. Let us first consider the results for before going on to assess the change in the ARs and . We analyse two specific cases of ellipse at three different porosities, *ϕ*=0.05, 0.1 and 0.15. Ellipse A has aspect ratio 4/3 and ellipse B has aspect ratio 2, where the semi-major axis is in the *x*_{1} direction. Note that ellipse B leads to a higher degree of anisotropy than ellipse A, as should be expected by the higher eccentricity of the elliptical pore shape. These results are tabulated in tables 4–7.

In the case of ellipse A, for the average porosity case of 8.26 per cent for males, we obtain ARs of 1.07 and 1.13 for undrained bone ( and , respectively) and 1.09 and 1.18 for drained bone. For female samples, where the average porosity is 9.72 per cent, we find 1.08 and 1.15 for undrained bone and 1.10 and 1.21 for drained bone.

## 5. Discussion

Note firstly that the asymptotic and Mori–Tanaka schemes coincide and therefore we conclude that these are complementary methods with which we can describe the effective properties of cortical bone. It is certainly the case, however, that asymptotic homogenization permits the incorporation of complex microstructure which is not as easy to include via the Mori–Tanaka scheme. For example, the local region indicated in figure 2 could be used as a periodic cell in the AS (Parnell & Abrahams 2006, 2008). Note also that these schemes lie within the Hashin–Rosen bounds.

As has been seen for the case of circular pores, the results obtained with the AS are in good quantitative and qualitative agreement with the finite-element results performed on real two-dimensional microstructures obtained from images. The relative difference between finite-element computations and the AS results (tables 1 and 2) is larger in the case of drained bone. This may be due to the fact that the AS is stable for a wide range of phase contrasts while the finite-element computational homogenization method may not be as stable since it depends on boundary conditions.

As should be the case (due to an increase in phase contrast), the AR is larger in the case of empty pores (figure 11). What is most interesting to note from this figure is that the increase of AR with porosity shows the same trend for the AS and finite-element results in both the undrained and drained bone cases.

We therefore conclude that the idealized pore model considered here could be used as a reasonable estimation of the local effective properties of cortical bone. However, this conclusion should be balanced by the fact that the finite-element computations (Grimal *et al*. 2008) have been performed with displacement boundary conditions alone and it could not be formally demonstrated that the size of the RVE used was large enough to obtain well-defined effective properties. In addition, Grimal *et al*. (2008) have questioned the validity of their computations for the case of *c*_{44}, which may explain the lack of agreement between the various methods in the case of and .

The study of the deviation from the arithmetic mean exhibits the value of homogenization formulae for the effective properties rather than using gross assumptions such as the Voigt average. Furthermore, it can be shown that this deviation is significantly larger (at least one order of magnitude) for certain elastic moduli and also depends strongly on the material inside the pores.

The AS lies inside the Hashin–Rosen bounds for the effective properties and (figures 12–15) and it is seen from these figures that, in fact, the bounds themselves can give reasonable predictions in the case of here but not for .

For elliptical cross sections with high eccentricity, corresponding to an aspect ratio of two (ellipse B above), the ARs are 1.06 and 1.25 compared to 1.16 with circular pores (9.72% porosity, drained bone, corresponding to the average porosity in female cortical bone). This result and the others in tables 4–7 are, as far as we know, the first direct quantification of the effect of pore shape on the anisotropy of cortical bone. The overall anisotropy in real bone (where pores have random orientation) is expected to deviate from the circular pore case less than our results for the elliptical cross section indicate here, since this latter situation is an extreme case where all pore cross sections are aligned.

The above analysis allows us to deduce the magnitude of the contribution to the anisotropy of the mesoscale porosity alone; the AR lies in the range 1–1.3, which is consistent with the results in Crolet *et al*. (2005) and Grimal *et al*. (2008). This approximate range of AR (1–1.3) was found for both drained and undrained bone and also for pores of elliptical and circular cross section.

Note that these ARs are lower than actual anisotropy in cortical bone since in the present work the bony matrix anisotropy has been neglected; the latter has a significant influence on cortical bone anisotropy, as was pointed out by Crolet *et al*. (2005) and Grimal *et al*. (2008).

## 6. Conclusions

We have applied the method of asymptotic homogenization in order to predict the local effective elastic properties of cortical bone. In particular, we have used the method developed by Parnell & Abrahams (2006, 2008), which is ideally suited to materials with high phase contrast and complex microstructure. Furthermore, as far as the authors are aware, this is the first time that comparisons have been made between effective moduli obtained on idealized mesostructures (using asymptotic homogenization) and on real images of cortical bone (using finite element methods; Grimal *et al*. 2008).

The benefits of the AS are numerous but most importantly it provides fast, consistent predictions of the effective properties of cortical bone that can be used for the interpretation of measurements on drained or undrained bone samples at all mesoscale porosities within the physiological range.

The proposed method is semi-analytical and given the correct expansions of the required doubly periodic potential functions, the only numerical computation involved is the solution of a linear system of equations. Hence, its implementation is straightforward and it is hoped that the scheme will be easily implemented by others. This is possible by downloading the appropriate Mathematica codes available online (Parnell 2008; www.maths.manchester.ac.uk/∼wparnell/research). The computation time for the entire elastic coefficient matrix is of the order of a second. As a consequence, the method is extremely well suited to parametric studies.

Note that the matrix has been modelled as a homogeneous material, whereas in reality the matrix of cortical bone possesses properties that vary from point to point, in particular, due to different degrees of mineralization. However, our assumption of homogeneity is based on the hypothesis that the inhomogeneity of matrix properties is expected to have a small effect on the overall properties because the elastic contrast between different regions in the matrix is small when compared with the contrast between mean matrix elasticity and the elasticity of the mesoscopic pores. Results that tend to confirm this hypothesis have been presented in Grimal *et al*. (2007). Note also that since the *local* effective properties have been determined, this means that the resulting *homogenized* material properties may vary within the cross section of the cortical bone (as depicted in figure 2), depending upon the mesoscale porosity within that local region.

The AS used in this work for a simple hexagonal lattice where pores are described as cylinders can be used straightforwardly to model matrix anisotropy and complex microstructures that more closely resemble that of cortical bone by considering patterns with several cells including different pore shapes. This work is underway and will be reported on in future articles. This will lend itself to validation when compared with experimental data at the meso- and macroscales.

## Acknowledgments

W.J.P. would like to acknowledge the support of the EPSRC via his Postdoctoral Research Fellowship, grant GR/S92151/01. Q.G. is grateful to Université Pierre et Marie Curie-Paris VI for the financial support BQR ‘Modélisation mécanique du tissu osseux’.

## Footnotes

- Received June 12, 2008.
- Accepted June 18, 2008.

- © 2008 The Royal Society