Simulation of ultrasonic wave propagation in anisotropic cancellous bone immersed in fluid

Wave Motion 47 (2010) 117–129

Contents lists available at ScienceDirect

Wave Motion journal homepage: www.elsevier.com/locate/wavemoti

Simulation of ultrasonic wave propagation in anisotropic cancellous bone immersed in fluid Vu-Hieu Nguyen, Salah Naili *, Vittorio Sansalone Laboratoire Mécanique Physique, Université Paris 12 – Val de Marne, 61, avenue du Général de Gaulle, 94010 Créteil Cedex, France

a r t i c l e

i n f o

Article history: Received 22 January 2009 Received in revised form 31 August 2009 Accepted 14 September 2009 Available online 20 September 2009 Keywords: Anisotropic poroelasticity Biot’s model Ultrasonic Fluid-structure coupling Finite element Cancellous bone In vitro testing

a b s t r a c t Cancellous bone at the macroscopic level is known as being a two-phase anisotropic material composed of a solid rod-like or plate-like skeleton filled with viscous fluid. Quantifying the mechanical properties of cancellous bones by using ultrasound techniques should take into account these complexities. This paper proposes a model to investigate the transient ultrasonic wave propagation in a cancellous bone immersed in an acoustic fluid. A finite element time-domain model based on us  w formulation is developed to consider an anisotropic porous medium (Biot’s model) coupled with two acoustic fluids. Some numerical results are shown to illustrate the influence of the bone anisotropy to the reflection and transmission of plane waves in a human cancellous bone specimen. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction During the last decade, using of quantitative ultrasound (QUS) techniques for measuring mechanical properties of bone has received significant attention, particularly for its potential in estimating the bone fragility or fracture risk. The measurement is based on the estimation of the speed of sound and the attenuation of the wave transmitted through the bone. However, the interpretation of QUS measurements raises numerous problems because of the great complexity of the bone material which is a porous, anisotropic and heterogeneous medium. Two types of bone may be distinguished: cortical (compact) and cancellous (trabecular or spongy) bone. In cortical bone, porosity varies between 5% and 10%. Cancellous bone is much more porous with porosity ranging from 50% to 90%. Cancellous bone at the macroscopic level is known as to be a twophase anisotropic material composed of a solid skeleton (an assembly of rod-like or plate-like trabeculae) filled with viscous fluid (the marrow). The acoustic modeling of cancellous bone is relatively difficult in comparison with that of the cortical bone due to its greater complexity. For the in vitro studies of the ultrasonic wave propagation in cancellous bone, most ultrasound devices use the through-transmission measurement on the bone sample to determine broadband ultrasound attenuation (BUA) and speed of sound (SOS) [1–11]. Different theoretical models have been proposed: Biot’s model [1,2,4–6,10,12], multilayer fluid–solid model [7] or l-CT based model [11,13]. It has been shown that both mechanical properties and structural configuration of trabeculae strongly influence the wave characteristics. By examining experimentally the ultrasound wave propagation in bovine cancellous bone, two distinguished wave modes were observed propagating in the direction parallel to the main trabecular alignment [5]. The fast and slow

* Corresponding author. Tel.: +33 1 45 17 14 45; fax: +33 1 45 17 14 33. E-mail address: [email protected] (S. Naili). 0165-2125/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.wavemoti.2009.09.002

118

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

modes correspond to the in-phase and out-of-phase relative motions between the fluid and solid, respectively. On the other hand, when considering the wave propagation in the direction perpendicular to the main trabecular alignment, only a single wave mode was observed. It has also been shown that numerical models based on Biot’s model could predict the sound speed of the two waves in good agreement with experimental results [6]. By using Biot’s equations, other authors [1,4,12] investigated the influence of the porosity, permeability or tortuosity on the propagation mechanism of in vitro ultrasonic waves. In these studies, the beam axis was assumed to be parallel to one of the principal directions of the bone specimen material (either parallel or perpendicular to the main trabecular alignment) and the transmitted waves were considered as plane waves. As a consequence, it was sufficient to use different sets of isotropic material properties to study the wave propagations in the principal directions of an orthotropic medium. Moreover, a simplified geometry and a linear elastic homogeneous material model were considered, which allowed to use a semi-analytical [4] or finite-difference [6] method to simulate the transient ultrasound wave transmission. To the best of our knowledge, the use of finite element method (FEM) has not yet been employed for this problem. However, FEM could be very advantageous in considering complex geometries as well as nonlinearities. Some experimental studies have investigated the propagation of ultrasonic waves in an arbitrary direction through anisotropic cancellous bone [5,7]. The transmitted wave modes were shown to strongly depend on the angle between the wave direction and the main trabecular alignment. For numerical modeling of this problem, the Scheonberg’s multilayer model has been used [7,14]. To the best of our knowledge, there are no published studies on the in vitro transmission test based on anisotropic Biot’s model. The objective of this paper is twofold. The first objective is to present an efficient numerical method for the simulation of ultrasonic wave transmitted in cancellous bone immersed in fluid. A finite element time-domain model of porous medium is proposed for simulating the high-frequency wave propagation in anisotropic cancellous bone. The second objective of this paper is to investigate the influence of bone anisotropy on the wave propagation. In particular, cancellous bone is modeled as an orthotropic medium and an ultrasonic wave pulse is emitted in a direction tilted with respect to one of the principal axes of the bone material. The interest of this study relies on the fact that, when a bone specimen is taken out from bone volume, the principal directions of the cancellous bone materials are not known precisely but are only estimated. Thus, when testing the specimen in vitro, the wave direction is likely to be misaligned with the principal directions of the material and this should be taken into account when analyzing experimental data. The paper is organized as follows. After this introduction as background, in Section 2, the governing equations of a bone sample immersed in an ideal fluid are presented. The bone is modeled as an anisotropic saturated porous medium by using Biot’s model in high frequency domain. Section 3 presents the weak formulation and the finite element implementation for the transient problem of a saturated porous medium coupled with a fluid. Section 4 begins with a numerical validation of the finite element solution. Then, the influence of the misalignment described above on the reflected and transmitted wave signals is investigated. Finally, conclusions and perspectives of this work are drawn in Section 5.

2. Statement of the problem In this section we describe a simplified mechanical model of an in vitro transmission test performed on a specimen of trabecular bone. For this test, a slice-shaped bone specimen is immersed in a fluid. An acoustical source and a receiver are located in the fluid on opposite sides of the slice. The source emits a wave pulse which is recovered by the receiver. Travelling from the source to the receiver, the wave encounters the bone specimen, where it is in part reflected and in part transmitted. The characteristics of the reflected and transmitted waves depend on the bone material. Thus, their analysis may be usefully used to characterize the bone material and infer information on its structure.

2.1. Description of the geometrical configuration A 2D schema of the test is represented in Fig. 1. We consider a bone layer embedded in two half-spaces occupied by fluid. Let RðO; x1 ; x2 Þ be the reference Cartesian frame where O is the origin and ðx1 ; x2 Þ is an orthonormal basis for the space. The coordinates of a point x in R are specified by ðx1 ; x2 Þ and the time by t. The fluid occupies the two half-spaces Xf1 and Xf2 . The bone layer occupies the unbounded domain Xb . The two plane interfaces between the bone ðXb Þ and the fluid domains bf (Xf1 and Xf2 ) are denoted by Cbf 1 and C2 , respectively. The thickness of the bone sample is denoted by L. The acoustical source f is located in the fluid domain X1 and the receiver is located in the fluid domain Xf2 . We assume that the beam axis coincides with the axis ðO; x1 Þ. The dashed rectangle ABCD represents a zone in which the wave field may be considered to be plane.

2.2. Governing equations In what follows, we will denote by a superposed dot the time derivative, by r and r the gradient and divergence operators in space, respectively. The symbol ‘‘” will denote the scalar product (between two vectors or tensors).

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

119

Fig. 1. Configuration of an in vitro ultrasonic transmission test.

2.2.1. Equations in the fluid domains (Xf1 and Xf2 ) The fluid occupying the domains Xf1 and Xf2 is modeled as an acoustic fluid whose mass density and bulk modulus at equilibrium are qf and K f , respectively. Let v 1 and q1 be the velocity and the mass density fields in the fluid domain Xf1 , respectively. In a linear acoustic theory, the mass conservation equation reads

q_ 1 þ qf r  v 1 ¼ Q in Xf1 ;

ð1Þ

where Q is the acoustic source density. Neglecting the body forces (other than inertia) and denoting by p1 the pressure field, the Euler equation reads

qf v_ 1 þ rp1 ¼ 0 in Xf1 :

ð2Þ

The constitutive equation for the acoustic fluid is given by

p1 ¼ c2f q1 ;

ð3Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffi where cf is the wave velocity in the fluid at equilibrium that is defined as cf ¼ K f =qf . By combining Eqs. (1)–(3), the wave equation for the fluid in the domain Xf1 may be expressed in terms of the pressure field only

1 1 2 1 €1  p r p1 ¼ Q_ Kf qf qf

in Xf1 :

ð4Þ

Following the same procedure, the wave equation for the fluid in the domain Xf2 (in which there is no acoustic source) reads

1 1 2 €2  r p2 ¼ 0 in Xf2 : p Kf qf

ð5Þ

2.2.2. Equations in the anisotropic poroelastic layer ðXb Þ We use the theory of anisotropic poroelasticity to describe the behavior of trabecular bone. The theory of isotropic poroelasticity was introduced by Biot [15] and was extended later to the anisotropic case by Biot [16] and Biot and Willis [17]. A detailed descriptions of the poroelastic theory can be found in [18]. The poroelastic equations employed here are based upon Biot’s original works as well as the recent developments in anisotropic constitutive equations [19–21]. The bone layer is assumed to consist of a solid skeleton (of density qs ) and of a connected pore network saturated by a fluid (of density qf ). Neglecting the body forces (other than inertia), the equations describing the wave propagation in the bone layer read

€ s þ qf w; € r  r ¼ qu €s

1

 rp ¼ q f u þ k

ð6Þ _ þ b w; € w

ð7Þ

where q ¼ / qf þ ð1  /Þ qs is the mixture density, / is the porosity, r is the stress tensor and p is the interstitial pore pressure. The displacement vectors of the solid skeleton and of the fluid are denoted by us and uf , respectively; the displacement vector of the fluid relative to the solid frame (weighted by the porosity) is defined by w ¼ /ðuf  us Þ. The two tensors k and b are symmetric second-order tensors. The first one designates the anisotropic permeability tensor. The second one is defined from the tortuosity tensor a as b ¼ ðqf =/Þ a . A more specific discussion on the determination of the tensors k and a for cancellous bone will be presented in Section 3.5. The constitutive equations for the anisotropic linear poroelastic material may be stated as

120

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

r ¼ C   a p; 1  p ¼ r  w þ a  ; M

ð8Þ ð9Þ

where C is the drained elasticity fourth-order tensor of the porous material, a is called Biot effective stress tensor, which is a symmetric second-order tensor and M is the Biot scalar modulus. The infinitesimal strain tensor  is the symmetric part of the gradient of us . 2.2.3. Boundary conditions – Exterior boundaries. Assuming the real box in which the in vitro test is realized to be sufficiently large with respect to the system depicted in Fig. 1, the waves reflected on the real boundary may be neglected. Thus, the boundary conditions for the system depicted in Fig. 1 may be expressed as

pi ¼ 0 on Cfi 1 p ¼ 0 on C

b1

for i ¼ 1; 2;

ð10Þ ð11Þ

;

us ¼ 0 on Cb1 :

ð12Þ

bf – Interfaces between the porous solid and the fluid domains. At the two interfaces Cbf 1 and C2 , the continuity of pressure and stress fields between the porous solid and the fluid domains requires

p ¼ pi

)

on Cbf i

bf rnbf i ¼ pi ni

for i ¼ 1; 2;

ð13Þ

bf f b where nbf i is the normal unit vector to Ci pointing from the porous domain X toward the fluid domain Xi (see Fig 1). In ði ¼ 1; 2Þ is assumed, requiring addition, open-pore condition at the interfaces Cbf i

bf _s _  nbf w i ¼ ðv i  u Þ  ni

on Cbf i for i ¼ 1; 2:

ð14Þ

In view of Eq. (2), the interface condition (14) may be rewritten as

1

qf

!  nbf on Cbf i ¼ 0 i for i ¼ 1; 2:

€ þu €s rpi þ w

ð15Þ

3. Finite-element time-domain solution for anisotropic bone coupled with fluid We assume the bone material to be a fluid-saturated orthotropic porous medium. When the ultrasound wave is transmitted along an axis parallel to one of the principal directions of the bone material, and if only plane waves are considered, the wave propagation problem may be described by a one-dimensional model. For this case, we developed a semi-analytical solution of the wave problem following the work of Schanz and Cheng [22] and Buchanan et al [23]. The solution was obtained using the Laplace transform technique to transform the time-dependent equations in the Laplace domain, where a closed-form solution can be found. Then, the time-dependent solution was computed by numerical integration of the inverse Laplace transform. Details of this computation are omitted for sake of space and will be published elsewhere. This one-dimensional model is well suited when the incident wave is parallel to one of the principal directions of the (orthotropic) bone material. However, in practice, these directions are not known precisely but only estimated. Thus, it is necessary to develop a strategy to analyze wave propagation in bone in a more general setting. In this section, we present a finite element formulation of the wave propagation problem in the in vitro test. This formulation is well suited to deal with a large class of situations – namely, to account for a misalignment between the wave direction and the material principal axes. 3.1. Weak formulation in the fluid domains Xf1 and Xf2 The pressure field in the fluid occupying the domain Xf1 is described by Eq. (4). The weak form of this equation may be ~1 (with support in Xf1 ) and then integrating over the domain Xf1 . obtained by multiplying it by a scalar valued test function p By applying the Green–Gauss theorem and taking into account the boundary conditions (10)–(15), the weak formulation of the acoustic wave problem in Xf1 reads

Z

Xf1

1 ~1 p €1 dV þ p Kf

Z

Xf1

1

qf

rp~1  rp1 dV þ

Z

Cbf 1

€ þu € s Þ  nbf ~ 1 ðw p 1 dS ¼

Z

Xf1

1

qf

~1 Q_ dV; p

ð16Þ

where dV and dS are the differential surface and volume elements, respectively. Note that the nbf 1 is pointing toward the interior of the fluid domain Xf1 . A similar treatment may be reserved to Eq. (5) describing the fluid pressure field in the domain Xf2 , leading to the weak form of the acoustic wave problem in Xf2

121

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

Z Xf2

1 ~2 p €2 dV þ p Kf

Z

1

qf

Xf2

rp~2  rp2 dV þ

Z Cbf 2

€ þu € s Þ  nbf ~ 2 ðw p 2 dS ¼ 0 ;

ð17Þ

~2 is a scalar valued test function with support in Xf2 . Equalities (16) and (17) are required to hold for any test funcwhere p ~2 . ~1 and p tions p 3.2. Weak formulation in the bone domain Xb In order to resolve the time-dependent problem of a porous medium, different finite element formulations were developed in the past. The typical finite element models are based on us  p, us  w  p or us  w formulations (e.g. see [24,25]). The suitability of these models depends essentially on the frequency range of the problem in hand [25,26]. For the present study, which is a high frequency problem, the acceleration term cannot be neglected, then the finite element analysis developed herein will be based on the us  w formulation. The balance equations of linear momentum (6) and (7) are already written in terms of us and w. Accordingly, constitutive Eqs. (8) and (9) may be restated as

r ¼ Cu  þ M a r  w;

ð18Þ

p ¼ Mða   þ r  wÞ;

ð19Þ

where Cu ¼ C þ M a  a is the undrained elasticity tensor. (The symbol  designates the tensorial product between two tensors.) ~ be two vector valued test functions with support in Xb . The weak form of balance Eqs. (6) and (7) may be ~ s and w Let u ~ respectively, integrating over Xb and applying the Green–Gauss theorem ~ s and w, obtained by multiplying them by u

Z Xb

Z Xb

~ s  qu € s dV þ u

Z

~  ðbwÞ € dV þ w

~ s  qf w € dV þ u

Xb

Z Xb

Z

~  qf u € s dV þ w

Xb

~ s  r dV ru

Z Xb

Z Cb

~ s  ðr nCb Þ dS ¼ 0; u

~  ðk1 wÞ _ dV þ w

Z Xb

~  ðp IÞ dV þ rw

ð20Þ Z Cb

~  ðp nCb Þ dS ¼ 0; w

ð21Þ

bf b1 where Cb ¼ Cbf and I is the second-order identity tensor. By taking into account the boundary conditions (12) 1 [ C2 [ C and (13) and the constitutive Eqs. (18) and (19), the weak formulation of the wave propagation problem in the bone layer may be restated as

W ui þ W ue þ W ut ¼ 0;

ð22Þ

w w w Ww i þ W d þ W e þ W p ¼ 0;

ð23Þ

~ s and w, ~ and where which are required to hold for any test functions u

W ui ¼ W ue ¼

Z Z

~ s  qu € s dV þ u

Xb

Z Xb

~ s  qf w € dV; u

ð24Þ

~ s Þ  ½Cu  þ M a ðr  wÞ dV; ðru Z Z ~ s  ðp1 n bf Þ dS þ ~ s  ðp2 n bf Þ dS; u u W ut ¼ C C Cbf 1

Ww i ¼ Ww d ¼ Ww e ¼ Ww p ¼

ð25Þ

Xb

1

Z

ZX

~  bw € dV þ w

b

ZX

b

Xb

ð26Þ

2

~  qf u € s dV; w

ð27Þ

~  ðk1 wÞ _ dV; w

ð28Þ

~ ðr  wÞ½M ða   þ r  wÞ dV; Z ~  ðp1 n bf Þ dS þ ~  ðp2 n bf Þ dS: w w C C

ð29Þ

Xb

Z

Z

Cbf 2

Cbf 1

1

Cbf 2

ð30Þ

2

3.3. Finite element formulation For implementation purposes, it is more convenient to represent the problem using the Voigt notation (see e.g. [25]). For instance, constitutive Eqs. (18) and (19) read

 ¼ Cu  þ M a r  w; r

ð31Þ

T

p ¼ Mða  þ r  wÞ;

ð32Þ T

T

T

 ¼ ða11 ; a22 ; a33 ; a12 ; a23 ; a13 Þ ; moreover, Cu is the  ¼ ðr11 ; r22 ; r33 ; r12 ; r23 ; r13 Þ ;  where r  ¼ ð11 ; 22 ; 33 ; 212 ; 223 ; 213 Þ ; a 6  6 elasticity matrix whose components may be obtained from Eq. (18).

122

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

The finite element method may be applied to the spatial discretization of the weak equations stated in the previous sections. The domain Xd (either Xf1 , or Xf2 , or Xb ) is partitioned into N ed elements Xe . In each element, the relevant state variables (either p1 , or p2 , or us and w, respectively) are approximated in terms of their nodal values (either Pe1 , or Pe2 , or Ue and We , respectively)

p1 ðx; tÞ  Np ðxÞPe1 ðtÞ in Xe  Xf1 ; p2 ðx; tÞ 

Np ðxÞPe2 ðtÞ

e

in X  X

e

s

ð33Þ

f 2;

ð34Þ e

u ðx; tÞ  Nu ðxÞU ðtÞ;

e

b

wðx; tÞ  Nw ðxÞW ðtÞ in X  X ;

ð35Þ

where Np ; Nu and Nw are interpolation (or shape) functions. Collecting the local nodal-value vectors in their global counterparts (P1 , or P2 , or U and W, respectively) and using the Galerkin method, Eqs. (16), (17), (22) and (23) lead to the system of equations

€ þ CXðtÞ _ MXðtÞ þ KXðtÞ ¼ FðtÞ;

ð36Þ

T

T

where X ¼ ðP1 ; U; W; P2 Þ is the global vector of nodal values to be solved for; FðtÞ ¼ ðP0 ðtÞ; 0; 0; 0Þ is the global force vector; M, C and K are the global mass, damping and stiffness matrices, respectively. These quantities may be represented in matrix form as

2

Mp1

6 6 0 M¼6 6 0 4 0

Mpu 1

Mpw 1

Mu

Muw

Mwu

Mw

Mpu 2

Mpw 2

3

0

2

7 0 7 7; 0 7 5 Mp2

0 0

0

60 0 0 6 C¼6 4 0 0 Cw 0 0

3

2

07 7 7; 05

6 up 6 K1 K¼6 6 wp 4 K1 0

0

0

0

Kp1

0

0

Ku

Kuw

Kwu

Kw

0

0

0

3

7 7 Kup 2 7 ; wp 7 K2 5 Kp2

ð37Þ

in which the sub-matrices noted with an upper-bar are defined on the interfaces between two domains and represent coupling operators; the sub-matrices noted without the upper-bar are defined in one domain. These global sub-matrices are determined by assembling the elementary matrices as

[Z

Mp1 ¼ Mp2 ¼ Mu ¼

[Z

qNTu Nu dV;

Xe

e

[Z

Kp1 ¼ Kp2 ¼ Ku ¼

Xe

NTw bNw dV;

ð39Þ

qf NTu Nw dV;

ð40Þ

1

[Z

[Z e

Mw ¼

NTw k Nw dV;

Xe

e

Xe

e

ð38Þ [Z e

[Z

Muw ¼ ðMwu ÞT ¼ Cw ¼

1 T N Np dV; Kf p

Xe

e

X

e

Xe

1

BTu Du Bu dV; T

Kuw ¼ ðKwu Þ ¼

[Z Xe

e

pu Mpu 1 ¼ M2 ¼

[Z e

up Kup 1 ¼ K2 ¼

ðrNp ÞT rNp dV;

qf

e

[Z e

Ce

Ce

ð41Þ

Kw ¼

ð42Þ

[Z e

Xe

M ðr  Nw ÞT ðr  Nw Þ dV;

MBTu aðr  Nw Þ dV;

NTp Nu dS;

NTu Np dS;

ð44Þ

pw Mpw 1 ¼ M2 ¼

[Z Ce

e

wp Kwp 1 ¼ K2 ¼

ð43Þ

[Z e

Ce

NTp Nw dS;

NTw Np dS;

ð45Þ ð46Þ

where Bu is the interpolation matrix allowing to evaluate the strain in Xb from the nodal displacements:   ¼ Bu U. For sake of notation, in Eqs. (38)–(44) we did not specify the domain that the elements Xe belong to – this is made clear looking at the functions to be integrated. Moreover, in Eqs. (45) and (46), Ce denotes the interface between two adjacent elements belonging to two different domains – again, the domains are not specified. 3.4. Implementation The weak formulation presented above has been successfully implemented in the finite element software COMSOL Multiphysics [27] to study the 2D problem sketched in Fig. 2. This system roughly corresponds to the box ABCD in Fig. 1. Quadratic shape functions were used for all the fields. Solution in time domain was achieved using the COMSOL built-in direct time integration algorithm based on the generalized-a scheme. This integration method is based on an implicit second-order scheme.

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

123

  In this study, we consider that the bone sample is located in the plane containing the direction parallel x01 and one direc 0 tion orthogonal x3 to the main trabecular alignment. The angle h between x1 and x01 – to be called propagation angle – measures the misalignment between the main trabecular direction and the beam axis. The wave source is supposed to be far enough to consider only plane waves arriving on the system schematized by the box ABCD (see Fig. 2). Hence, the boundary conditions may be restated as

p1 ¼ p0 ðtÞ on Cfsrc ; @p1 ¼ 0 on Cf1 ; @x2 @us @w ¼ 0 and ¼ 0 on Cb ; @x2 @x2 @p2 ¼ 0 on Cf2 : @x2

ð47Þ ð48Þ ð49Þ ð50Þ

The corresponding weak terms may be easily recovered.

3.5. Orthotropic constitutive model for human cancellous bone This section is devoted to specify the constitutive model for human cancellous bone. Cancellous bone is assumed to be orthotropic. Poroelastic coefficients, permeability and tortuosity tensors are then modeled accordingly.

3.5.1. Biot’s poroelastic constants Elastic modulus of trabecular bone Eb is known to be related to the elastic modulus and volume fraction of the solid phase. Gibson [28] proposed an empirical formula: Eb ¼ Es ð1  /Þn , where Es is the Young’s modulus of the solid phase, / is the porosity and n is a parameter to be set on the basis of empirical data. Typical values of n range between 1 and 3, depending on the alignment and the type of the trabecular structure [28,1]. To consider the orthotropic property of bone materials, Gibson’s formula was generalized by other authors. We use in this paper the set of elastic constants proposed by Yang et al. [29] for human cancellous bone

E10 ¼ 1:240 Es /1:80 ; s

E20 ¼ 0:885 Es /1:89 ; s

E30 ¼ 0:5288 Es /1:92 ; s

2G20 30 ¼ 0:5333 Es /2:04 ; 2G10 30 ¼ 0:6333 Es /1:97 ; s s 0:09 0:19 0 0 ¼ 0:256/ 0 0 ¼ 0:306/ 0 0 ¼ ; ; 23 13 12 s s

m

m

m

2G10 20 ¼ 0:9726 Es /1:98 ; s 0:176/s0:25 ;

ð51Þ ð52Þ ð53Þ

where /s ¼ 1  /. The parameters entering these formulae were identified empirically from the results of finite element simulations carried out on a large number of human bone specimens. While obtaining these relations, the solid bone phase was assumed to be isotropic and the trabecular bone to be an orthotropic material. These elastic moduli are expressed in a local     frame R0 O; x01 ; x02 ; x03 oriented along the three principal directions x01 ; x02 ; x03 of the bone material, the first direction corresponding to the predominant trabecular direction. Note that the local frame R0 may not coincide with the frame RðO; x1 ; x2 ; x3 Þ where the wave equations are expressed. Eventually, knowing the elastic tensors of the solid material Cs and of the trabecular frame C, the Biot’s constants a and M (Eqs. (8) and (9)) could be easily constructed [19,21]. Note that the tensor a of an orthotropic porous material is represented, in the frame R0 , by the diagonal matrix a ¼ Diagða10 ; a20 ; a30 Þ, where a10 ; a20 and a30 are the Biot coefficients in the principal directions.

  Fig. 2. Bidimensional model for plane wave propagation in cancellous bone. The principal axes of the bone material x0i may be misaligned with respect to the wave direction ðx1 Þ.

124

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

3.5.2. Permeability and tortuosity The anisotropic permeability tensor k can be determined as k ¼ ðg F r ðxÞÞ1 j [30], where g is the pore fluid dynamic viscosity and j is the intrinsic permeability tensor. For orthotropic porous media, the intrinsic permeability tensor j is expressed, in the frame R0 , by the diagonal matrix j ¼ Diagðj10 ; j20 ; j30 Þ, where j10 ; j20 and j30 are the permeabilities in the principal directions. The correction factor F r is introduced to take into account the viscous resistance of the fluid flow at high frequencies. Namely, F r ðxÞ is defined as the real part of a function FðxÞ depending on the load frequency x

1 n TðnÞ FðnÞ ¼ ; 4 1  2TðnÞ

with TðnÞ ¼

in

0 0 ber ðnÞ þ i bei ðnÞ  berðnÞ þ i beiðnÞ

sffiffiffiffiffiffiffiffiffiffi

xqf ; g

and n ¼ a0

ð54Þ

where a0 is the average pore size, berð nÞ and beið nÞ are the real and the imaginary parts of the Kelvin function, and a prime denotes differentiation. The tortuosity tensor a of an orthotropic porous medium is expressed, in the frame R0 , by the diagonal matrix a ¼ Diagða10 ; a20 ; a30 Þ, where a10 ; a20 and a30 are the tortuosity coefficients in the principal directions. In this paper, the quantities ai ði ¼ 10 ; . . . ; 30 Þ are assumed to be equal and are evaluated by using Berryman’s relation [31]

ai ¼ 1 þ r



 1 1 ; /

ð55Þ

where r is a parameter related to the microscopic structure of the solid. 4. Numerical results This section presents some numerical simulations of the in vitro ultrasound test on cancellous human bone by using the methods presented in the previous sections. As cancellous human bone is usually assumed to be an orthotropic poroelastic material, its principal directions provide a favorite natural frame to study it. To this purpose, when preparing samples for an in vitro test, the bone specimen is cut out from the cancellous bone volume according to the principal directions of the bone material, i.e. cross-sections of the specimen are taken either parallel or perpendicular to the main trabecular alignment. However, it is not easy in practice to set correctly the specimen since the trabecular structure is not known precisely but only estimated. For example, in [7] the cross-section of the bone specimens were determined with the naked eye. Thus, the principal directions of the bone material are not really orthogonal to the cross-sections of the bone specimen. In the in vitro ultrasound test, the wave direction is set orthogonal to a cross-section of the specimen. Thus, the behavior of the waves transmitting through the system may be ill predicted when supposing the wave direction to be parallel to one material principal direction. As a consequence, it is interesting to investigate the influence of this misalignment on the reflected and transmitted waves. In this section, we present the results of some numerical simulations that we performed to address this point. First, the actual parameters of the problem in hand will be specified. Then, a preliminary analysis will be done assuming the wave direction to be parallel to a material principal direction. Finally, the influence of a misalignment between these directions will be studied. 4.1. Parameters for the computations In the in vitro test that we simulate, a human cancellous bone specimen is immersed in water. The physical characteristics 3 of the water are given by the bulk modulus K f ¼ 2:25 GPa and the mass density qf ¼ 1000 kg m . Thus, the sound velocity 1 in water is cf ¼ 1500 m s . Bone physical parameters – typical for cancellous human bone – are presented in Table 1. Only essential input data are provided, the poroelastic parameters being computed using the formulae given in Section 3.5. The left and right boundaries of the rectangular domain ABCD in Fig. 2 are placed at x1 ¼ 0:02 m and x1 ¼ 0:03 m, respectively. The thickness of the bone specimen is L ¼ 0:01 m. The distance between the source Cfsrc and the bone specimen 4 m, respectively. As only plane waves boundary Cbf 1 is 0.02 m. The top and bottom boundaries are placed at x2 ¼ 0:5  10 are supposed to propagate perpendicularly to the bone sample, all the fields (p1 ; p2 ; us ; w, and p) are homogeneous in the x2 direction, that is to say their derivative with respect to x2 does vanish. As a consequence, with the boundary conditions presented in Eqs. (50) and (49), the size of the domain in the x2 direction does no matter and may be safely reduced. The excitation is a planar source emitting a broadband pressure pulse in the x1 direction. The time-history function of the pressure pulse is assigned as

Table 1 Poroelastic material properties of cancellous human bone. / (–)

K f (GPa)

qf ðkg m3 )

Es (GPa)

ms (–)

0.8

2.25

1000

20

0.3

g (Pa s) 3

10

j10 ðm2 Þ 8

10

j30 ðm2 Þ 10

9

r (see Eq. (55)) (–)

a0 ðÞ 3

0:5  10

0.25

125

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129 2

p0 ðtÞ ¼ P 0 e4ðf0 t1Þ sinð2pf0 tÞ;

ð56Þ

where f0 is the central frequency and P0 is a constant. In the numerical tests, the values f0 ¼ 1 MHz (commonly used in practical in vitro ultrasonic tests) and P0 ¼ 1 Pa were used. We point out that, in practice, common values of the pulse amplitude are around 10 MPa [32], that is 107 higher than the value we used for P 0 . This should be taken into account when comparing our results with experimental measurements. Two receivers R1 and R2 are placed in the fluid domains Xf1 and Xf2 to record the radiofrequency signals of the reflected and transmitted waves, respectively. Receiver R1 is placed on x1 -axis at the position x1 ¼ 0:002 m (i.e. at a distance of 0.002 m from the left edge Cbf 1 of the bone specimen) to sample the pressure p1 . Receiver R2 is placed on x1 -axis at x1 ¼ 0:012 m (i.e. at a distance of 0.002 m from the right edge Cbf 2 of the bone specimen) to sample the pressure p2 . Both porous solid and fluid domains were meshed using  quadraticalrectangular elements. The typical size of the elements was taken as h ¼ 1  104 m over the whole domain Xf1 [ Xb [ Xf2 , corresponding to a total number of 4814 degrees of freedom. The duration of the phenomenon studied was T ¼ 4  105 s and the time step used for transient analysis was set to Dt ¼ 0:35  108 s. The total computational time for this problem by using a desktop computer (IntelÒ Core 2 Quad CPU 2.83 MHz) was around 200 s. The proposed element size h was chosen to be about one-tenth of the shortest wave length kmin ¼ cmin =fmax , where cmin is the slowest wave velocity in the domain, and fmax is the maximal value of frequency. The time step used in the simulations was calculated proportionally to the ratio between the element size h and the fastest velocity of wave cmax in the domain: Dt ¼ vh=cmax . Here, relative small time steps are required ðv ¼ 0:1Þ when using the generalized-a method. As the considered bone sample is an orthotropic material, the fastest and slowest wave velocities may be roughly estimated along two principal directions. TM

4.2. Plane wave propagation along the principal directions of the bone material As discussed before, when the plane wave direction is parallel to one of the principal directions of the orthotropic bone material, a semi-analytical solution was computed, which we used to validate the finite element computation. We present the analytical and finite element results for the transmitted waves in Fig. 3, when waves propagate in the direction parallel to the main trabecular alignment (x01 x1 , i.e. h ¼ 0 ). We can observe that the finite element solution perfectly matches the analytical one. Similar results (not shown here) are found when the emitted wave is orthogonal to the trabecular alignment (x1 x03 , i.e. h ¼ 90 ). Hence, the us  w finite element formulation turns out to be quite appropriate to this high-frequency transient problem. Fig. 3 presents the transmitted pressure wave measured by the receivers R2 ðx1 ¼ 0:012 mÞ when x01 x1 . Two separate transmitted waves can be clearly observed. These two waves are well known as the two Biot’s waves [33] corresponding to the in-phase (the fast wave) and out-of-phase (the slow wave) relative motion of the fluid and the solid in the bone. As observed experimentally [6], the amplitude of the fast wave is smaller than the slow wave’s one. The following group bf of four pulsations are the waves reflected backwards on Cbf 2 and then again forwards on C1 . The path of these waves is sketched in the figure. 4.3. Plane wave propagation in orthotropic cancellous human bone

Pressure at x1 = 0.012m (Pa)

In this section, we investigate the behavior of the transmitted and reflected waves when the wave direction is neither parallel nor perpendicular to the main trabecular alignment. To this purpose, the propagation angle h between x1 and either x01 or x03 was slightly varied simulating a misalignment of the wave direction with respect to the bone material principal directions. Figs. 4 and 5 refer to the case where the wave direction is approximately parallel to the main trabecular alignment. In these simulations, h is varied from 0° to 15° by steps of 5°. Fig. 4 shows the transmitted wave signals at receiver R2 . While the difference between the signals corresponding to h ¼ 0 and h ¼ 5 is negligible, the signal obtained for h ¼ 15 is quite

0.5

analytical FEM slow wave fast wave

0

reflected waves from 2 interfaces

transmitted waves

−0.5 0

0.5

1

1.5

2 time (s)

2.5

3

3.5

4 −5

x 10

Fig. 3. Transmitted waves when the wave direction is parallel to the trabecular alignment ðx01 x1 Þ.

126

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

different. In this case, the amplitudes of the fast and slow waves become nearly equivalent. Results concerning thereflected  waves areshown in Fig. 5. On the top and on the bottom, it is shown a zoom on the reflected wave from the left Cbf and 1 boundaries of the bone sample, respectively. One can observe that the reflected waves from the left boundary Cbf right Cbf 2 1 (on the top of Fig. 5) do not depend perceptibly on the propagation angle h. The velocities of the reflected waves from the  right boundary Cbf 2 (on the bottom of Fig. 5) change significantly for h ¼ 10 and 15°, while their amplitudes do not. Multi-dimensional effects appear in the system when the wave propagation direction is not coincident with one of the axes of material symmetry of the bone material. To illustrate these effects, we present in Fig. 6 the time history of the displacement vector of the solid phase at the center of the bone sample ðx1 ¼ 0:005 m;x2 ¼ 0Þ when the wave direction is nearly  parallel to the main trabecular alignment x01  x1 . On the top, we show the longitudinal (with respect to the wave axis) component, i.e. in direction x1 , while on the bottom we show the transversal component, i.e. in direction x2 . No transversal displacement can be observed when h ¼ 0 . As the angular misalignment h increases, increasing transversal displacement can be observed. Thus, when h–0, not only compressional waves but also shear waves do propagate in the bone sample. Appearance of shear waves (i.e. of transversal components of the displacement) is an effect related to the anisotropy of the material. However, the amplitudes of the shear waves are much smaller (two orders of magnitude) in comparison with those of the longitudinal waves. Figs. 7 and 8 refer to the case where the wave direction is approximately orthogonal to the main trabecular alignment. In these simulations, h is varied from 75° to 90° (i.e. the angle between x1 and x03 is varied from 0° to 15°) by steps of 5°. One can

Fig. 4. Transmitted signals when the wave direction is nearly parallel to the main trabecular alignment ðx01  x1 Þ.

bf Fig. 5. Reflected signals from the left bone surface Cbf 1 (on the top) and from the right bone surface C2 (on the bottom) when the wave direction is nearly parallel to the main trabecular alignment ðx01  x1 Þ.

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

127

Fig. 6. Displacement components of the solid phase ðx1 ¼ 0:005 m; x2 ¼ 0Þ when the wave direction is nearly parallel to the main trabecular alignment ðx01  x1 Þ.

observe that both the transmitted and reflected waves (Figs. 7 and 8, respectively) do not depend perceptibly on the propagation angle. In summary, when considering waves propagating through cancellous bone in a direction parallel to the trabecular alignment, it is important to be as accurate as possible in determining the material principal directions. Even small errors could be

Fig. 7. Transmitted signals when the wave direction is nearly orthogonal to the main trabecular alignment ðx03  x1 Þ.

Fig. 8. Reflected signals when the wave direction is nearly orthogonal to the main trabecular alignment ðx03  x1 Þ.

128

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

misleading in interpreting the signals of the reflected and transmitted waves. On the contrary, their effect seems not to be important when considering the wave propagation in the direction perpendicular to the trabecular alignment. 5. Conclusion The in vitro transmission test is a useful tool to estimate mechanical characteristics of bone. However, experimental results may be misleading when the bone anisotropy is not properly dealt with. Results of this work will hopefully provide guidelines in interpreting the results of this test. In particular, the dependency of the features of the transmitted and reflected wave on the misalignment with respect to the material principal directions was highlighted. To model the in vitro ultrasound test on cancellous bone, bone material was considered as a saturated orthotropic poroelastic medium immersed in an acoustic fluid. We presented a finite element formulation to simulate the propagation of ultrasonic plane waves through a cancellous bone specimen immersed in fluid. The finite element analysis based on the us  w formulation has been shown to be suitable to describe correctly the wave propagation in cancellous bones in high frequency domain. This strategy allows investigating the wave propagation when the wave direction is not parallel to one principal direction of the bone material. Moreover, it also allows to take into account inhomogeneous physical characteristics (e.g. porosity, rigidity, permeability, etc.) of the bone material or irregularities of the bone surface. Note that the later issue would hardly be considered by using analytical or finite difference methods. When emitting an ultrasound pulsation in-axis with one material principal direction, the numerical results show two separate transmitted and reflected waves. These observations agree well with the numerical and experimental results given in the literature. Moreover, the numerical results obtained by using the FE method show that the characteristics of the transmitted and reflected waves depend strongly on the anisotropy of cancellous bone. The out-of-axis effect has been shown to play a important role when the wave direction is nearly parallel to the main trabecular direction. Otherwise, when the wave propagates approximately perpendicularly to the main trabecular direction, this effect seems to be negligible. The first results of our model are promising. Qualitative features of transmitted and reflected waves are in good agreement with the experimental data available in the literature. However, in order to have a more reliable model, some improvements would be done. Namely: (i) the constitutive equations for the bone material and the fluid might be improved; (ii) a quantitative comparison with available experimental data should be done for some specific cases; (iii) the effect of the position of the source as well as of the specimen geometry should be investigated. We will address these points in a forthcoming work. References [1] J.L. Williams, Ultrasonic wave propagation in cancellous and cortical bone: prediction of some experimental results by Biot’s theory, J. Acoust. Soc. Am. 91 (2) (1992) 1106–1112. [2] P. Nicholson, R. Strelizki, R. Cleveland, M. Bouxesein, Scattering of ultrasound in cancellous bone: predictions from a theoretical model, J. Biomech. 33 (2000) 503–506. [3] L. Cardoso, F. Teboul, L. Sedel, C. Oddou, A. Meunier, In vitro acoustic waves propagation in human and bovine cancellous bone, J. Bone Miner. Res. 18 (10) (2003) 1803–1812. [4] Z.E.A. Fellah, J.Y. Chapelon, B.S. Lauriks, C. Depollier, Ultrasonic wave propagation in human cancellous bone: application of Biot theory, J. Acoust. Soc. Am. 116 (1) (2004) 61–73. [5] A. Hosokawa, T. Otani, Acoustic anisotropy in bovine cancellous bone, J. Acoust. Soc. Am. 103 (5) (1998) 2718–2722. [6] A. Hosokawa, Simulation of ultrasound propagation through bovine cancellous bone using finite-difference time-domain methods, J. Acoust. Soc. Am. 118 (3) (2005) 1782–1789. [7] E.R. Hughes, T.G. Leighton, G.W. Petley, P.R. White, Ultrasonic propagation in cancellous bone: a new stratified model, Ultrasound Med. Biol. 25 (5) (1999) 811–821. [8] E.R. Hughes, T.G. Leighton, P.R. White, G.W. Petley, Investigation of an anisotropic tortuosity in a Biot model of ultrasonic propagation in cancellous bone, J. Acoust. Soc. Am. 121 (1) (2007) 568–574. [9] M. Hakulinen, J. Day, J. Töyräs, M. Timonen, H. Kröger, H. Weinans, I. Kiviranta, J. Jurvelin, Prediction of density and mechanical properties of human trabecular bone in vitro by using ultrasound transmission and backscattering measurements at 0.2–6.7 MHz frequency range, Phys. Med. Biol. 50 (2005) 1629–1641. [10] M. Pakula, F. Padilla, P. Laugier, M. Kaczmazek, Application of Biot’s theory to ultrasonic characterization of human cancellous bones: determination of structural, material and mechanical properties, J. Acoust. Soc. Am. 123 (4) (2008) 2415–2423. [11] G. Haïat, F. Padilla, F. Peyrin, P. Laugier, Fast wave ultrasonic propagation in trabecular bone: numerical study of the influence of porosity and structural anisotropy, J. Acoust. Soc. Am. 123 (3) (2008) 1694–1705. [12] L. Cardoso, A. Meunier, C. Oddou, In vitro acoustic wave propagation in human and bovine cancellous bone as predicted by Biot’s theory, J. Mech. Med. Biol. 8 (2) (2008) 183–201. [13] F. Padilla, E. Bossy, G. Haïat, F. Jenson, P. Laugier, Numerical simulation of wave propagation in cancellous bone, Ultrasonics 44 (2006) e239–e243. [14] F. Padilla, P. Laugier, Phase and group velocities of fast and slow compressional waves in trabecular waves, J. Acoust. Soc. Am. 108 (4) (2000) 1949– 1952. [15] M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (2) (1941) 155–164. [16] M.A. Biot, Theory of elasticity and consolidation for a porous anisotropic solid, J. Appl. Phys. 26 (2) (1955) 182–185. [17] M.A. Biot, D.G. Willis, The elastic coefficients of the theory of consolidation, J. Appl. Mech. 79 (1957) 594–601. [18] O. Coussy, Poromechanics, John Wiley and Sons, Ltd., 2004. [19] A.H. Cheng, Material coefficients of anisotropic poroelasticity, Int. J. Rock Mech. Mining Sci. 34 (2) (1997) 199–205. [20] S.C. Cowin, A recasting of anisotropic poroelasticity in matrices of tensor components, Transp. Porous Media 50 (1) (2003) 35–56. [21] M. Thompson, J.R. Willis, A reformation of the equations of anisotropic poroelasticity, J. Appl. Mech. ASME 58 (1991) 612–616. [22] M. Schanz, A.H. Cheng, Transient wave propagation in a one-dimensional poroelastic column, Acta Mech. 145 (2) (2000) 1–18. [23] J.L. Buchanan, R.P. Gilbert, A. Wirgin, Y. Xu, Transient reflection and transmission of ultrasonic wave in cancellous bone, Appl. Math. Comput. 142 (2003) 561–573.

V.-H. Nguyen et al. / Wave Motion 47 (2010) 117–129

129

[24] R.W. Lewis, B.A. Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, Wiley, 1998. [25] O.C. Zienkiewicz, A.H.C. Chan, M. Pastor, B.A. Schrefler, T. Shiomi, Computational Geomechanics with Special Reference to Earthquake Engineering, Wiley, 2000. [26] B. Simon, J.S.-S. Wu, O.C. Zienkiewicz, D.K. Paul, Evaluation of u  w and u  p finite element methods for the dynamic response of saturated porous media using one-dimensional models, Int. J. Numer. Anal. Methods Geomech. 10 (1986) 461–482. [27] COMSOL Multiphysics 3.5 - User’s Guide. [28] L. Gibson, The mechanical behaviour of cancellous bone, J. Biomech. 18 (5) (1985) 317–328. [29] G. Yang, J. Kabel, B. VanRietbergen, A. Odgaard, R. Huiskes, S.C. Cowin, The anisotropic Hooke’s law for cancellous bone and wood, J. Elasticity 53 (1999) 125–146. [30] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range, J. Acoust. Soc. Am. 28 (2) (1956) 179–191. [31] J.G. Berryman, Confirmation of Biot’s theory, Appl. Phys. Lett. 37 (1980) 382–384. [32] W.D. O’Brien, Ultrasound-biophysics mechanisms, Prog. Biophys. Mol. Biol. 93 (2007) 212–255. [33] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range, J. Acoust. Soc. Am. 28 (2) (1956) 168–178.