INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN BIOMEDICAL ENGINEERING Int. J. Numer. Meth. Biomed. Engng. (2012) Published online in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/cnm.2462
Simulation of ultrasonic wave propagation in anisotropic poroelastic bone plate using hybrid spectral/finite element method Vu-Hieu Nguyen*,† and Salah Naili Laboratoire Modélisation et Simulation Multi Echelle, Université Paris-Est, MSME UMR 8208 CNRS, 61 Avenue du Général de Gaulle, 94010 Créteil, France
SUMMARY This paper deals with the modeling of guided waves propagation in in vivo cortical long bone, which is known to be anisotropic medium with functionally graded porosity. The bone is modeled as an anisotropic poroelastic material by using Biot’s theory formulated in high frequency domain. A hybrid spectral/finite element formulation has been developed to find the time-domain solution of ultrasonic waves propagating in a poroelastic plate immersed in two fluid halfspaces. The numerical technique is based on a combined Laplace–Fourier transform, which allows to obtain a reduced dimension problem in the frequency– wavenumber domain. In the spectral domain, as radiation conditions representing infinite fluid halfspaces may be exactly introduced, only the heterogeneous solid layer needs to be analyzed by using finite element method. Several numerical tests are presented showing very good performance of the proposed procedure. A preliminary study on the first arrived signal velocities computed by using equivalent elastic and poroelastic models will be presented. Copyright © 2012 John Wiley & Sons, Ltd. Received 6 July 2011; Revised 7 December 2011; Accepted 12 December 2011 KEY WORDS:
ultrasound; cortical bone; poroelastic; functionally graded plate; hybrid spectral/finite element
1. INTRODUCTION In recent years, quantitative ultrasound techniques for assessing the properties of bone have received much attention owing to its potential in estimating bone fragility and/or fracture risk. Ultrasound is non-ionizing, and moreover, ultrasonic devices are relatively inexpensive and can be made portable. From the point of view of mechanics, bone is a heterogeneous, anisotropic and porous material. There are two kinds of bone: cortical (compact) bone and cancellous (trabecular or spongy) bone. The porosity varies from 5% to 15% in cortical bone and from 50% to 90% in cancellous bone. For measuring in vivo material properties of cortical long bones, a so-called ‘axial transmission’ (AT) technique has been developed [1]. The axial transmission technique uses a set of ultrasonic transducers (transmitters and receivers) placed in the same side of the investigated skeletal site along a direction close to the long bone axis. Consequently, the transducers form a line in contact with the skin along the bone axial axis. The transmitter emits an ultrasound pulse wave (around 250 KHz 2 MHz) that propagates along the cortical layer of bones. The analysis of the signals received at the receivers serve for estimating the geometric parameters as well as mechanical characteristics of the cortical bone at the measured skeletal site [2, 3]. Mechanical modeling of this experiment deals with describing a vibro-acoustic problem of a solid waveguide (which represents cortical bone) coupled with two fluid media (which represents soft tissues such as skin or marrow). Although in vivo cortical bone may have slightly varied *Correspondence to: Vu-Hieu Nguyen, Université Paris-Est, Laboratoire Modélisation et Simulation Multi Echelle, MSME UMR 8208 CNRS, 61 Avenue du Général de Gaulle, 94010 Créteil, France. † E-mail: [email protected] Copyright © 2012 John Wiley & Sons, Ltd.
V.-H. NGUYEN AND S. NAILI
thickness along the bone axis, the bone’s thickness is usually assumed to be constant. Therefore, the cortical bone may be described as plate-like or cylindrical-like structures. Assessment of the bone’s material properties requires careful analysis of the reflections, conversion modes and interferences of longitudinal and shear waves within the bone structure. However, current development of ultrasound techniques is still limited because the information potentially contained in the transmitted ultrasonic waves is not fully exploited, and parameters such as bone material properties or microstructural parameters are still difficult to recover. The interaction between ultrasound and bone remains poorly understood from a physical point of view owing to the complex nature of bone, which has a viscoelastic porous microstructure spanning many length scales. At the macroscopic scale, porosity in the radial direction (which is defined in the bone’s crosssection plan) is heterogeneous. The observations at all ages and for both genders show that the mean porosity in the endosteal region (inner part of the bone) is significantly higher than the porosity in the periosteal region (outer part of the bone) [4–6]. This observation may be explained by the fact that cortical bone is affected by age-related bone resorption and osteoporosis, causing reduction of bone shell thickness as well as increase of porosity, namely in the endosteal region. Moreover, the macroscopic mechanical properties of bone have been shown to strongly depend on its porosity [7, 8]. As a consequence, cortical bone may naturally be considered as a functionally graded material. It has been shown that the longitudinal curvature of bone had minor effects in estimating the first arrived signal (FAS) (see e.g. [9]). Therefore, many of the past studies have focused on the modeling of guided waves in long bones by using fluid-loaded homogeneous/multilayer/functional graded plate models. The understanding of wave phenomena involved in the multilayer structures has been studied by many authors in the frequency domain [10–13] or in the time domain [14–17]. In these studies, the cortical bone material has been considered as an equivalent (visco-)elastic medium of which the effective macroscopic mass density and effective macroscopic elasticity tensor are estimated from its porosity. The presence of the interstitial fluid, which was considered in many works for the analysis of the behavior of cortical bone tissue under low frequency loading (e.g. [18–21]) or of ultrasonic wave propagation through cancellous bones (e.g. [22–25]), has usually been neglected when studying ultrasonic wave propagation in cortical bones. To solve the time-domain wave propagation problem in the multilayer structures, there are mainly two approaches. The first one involves using (semi-)analytical methods such as the generalized ray/Cagniard-de-Hoop technique [15] or the direct stiffness matrix method [26]. The second one involves using numerical methods such as the finite difference method [14, 27, 28] or the finite element method (FEM) [17, 29]. Although the analytical methods are attractive to obtain reliable transient responses of the structure, numerical methods are often more efficient to treat problems with inhomogeneous materials or complex geometries. However, most numerical methods require important computational costs, especially for problems in the high-frequency domain. Moreover, absorbing boundary conditions are required for considering unbounded domains [28, 30, 31]. Alternatively, for considering waveguides with geometrical and mechanical properties, which are constant only along one or two directions, the hybrid numerical method (see for example [32–34]), alternatively called spectral finite element method (see for example [35, 36]), have been employed. The key point of this method consists in using a hybrid algorithm, which begins by employing the Fourier transform (with respect to time and to the longitudinal direction of the waveguide) to transform the problem into the frequency-wavenumber domain. Then, the wave equations in the spectral domain governed in cross section (or even a 1D domain in the case of infinite plates or axisymmetric waveguides), which may actually have inhomogeneous material properties, can be easily handled using the FEM [16, 32, 33, 36]. Note that this technique uses the conventional finite elements in spectral domain and does not use high-order finite element families known as ‘spectral elements’. The objective of this paper is twofold. The first objective is to present a numerical procedure to simulate the ultrasonic wave propagation in bone plate immersed in fluid. The bone plate is not modeled as an equivalent elastic medium but as an anisotropic poroelastic medium by using the Biot’s theory in high frequency domain. The numerical method is based on a spectral finite element technique, which allows to consider the variation of bone properties along its thickness direction. As far as the authors know, no equivalent numerical procedure based on spectral finite element technique has been presented previously in the literature for an anisotropic poroelastic medium. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
The second objective is to check if using equivalent macroscopic elastic models [10, 11, 14–17] is sufficient to evaluate the bone properties for all porosities in physiological ranges. The paper is organized as follows. After this introduction on the context of this study, as well as on methods to solve the time-domain wave propagation problem in the multilayer structures, Section 2 provides the governing equations of an anisotropic poroelastic bone plate immersed in fluid. Next, Sections 3 and 4 present the problem governed in Laplace–Fourier domain and develop the finite element formulation in this domain. Then, some numerical tests are presented in Section 5 to validate the proposed procedure. We will also present some comparison of the FAS velocities computed by using equivalent elastic and poroelastic models. Finally, conclusion and discussion will be presented in Section 6. 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 two-dimensional (2D) space, respectively. The symbol ‘’ will denote the scalar product, and the symbol ‘W’ between tensors of any order denotes their double contraction. In this work, we will employ a 2D configuration, which has been shown to be an appropriate geometrical model for describing the ultrasound wave in long bone in the framework of the axial transmission technique [9]. Let R.OI e1 , e2 / be the reference Cartesian frame where O is the origin and .e1 , e2 / is an orthonormal basis for the bi-dimensional space. We denote the coordinates of a point x f f in R by .x1 , x2 / and the time by t . The fluid occupies both half-spaces 1 and 2 . The bone layer occupies the unbounded domain b . The two plane interfaces between the bone layer b and the f f fluid domains 1 and 2 are respectively denoted by 1bf and 2bf , as shown in Figure 1. The thickf ness of the bone plate is denoted by h. The acoustical source is located in the upper fluid domain 1 . f
f
2.1. Equations in the fluid domains 1 and 2
The fluid occupying the domain f1 is modeled by an acoustic fluid whose mass density and bulk modulus at equilibrium are denoted by 1 and K1 , respectively. By using the linear acoustic theory and by neglecting the body forces (other than inertia), the wave equation for the fluid in the domain f 1 may be expressed only in terms of the pressure field p1 .x, t / as follows: 1 P pR1 r 2 p1 D Q, c12
8x 2 f1 ,
(1) f
where Q.x, t / is p the acoustic source density; the constant c1 is the wave celerity in 1 which is defined by: c1 D K1 =1 . Here, we impose an impulsive point source acting at point x s D .0, x2s / f 1 f) 1
Fluid 1 ( e2 e1
source (xs1, xs2)
receivers bf 1
nb1
O Solid ( Fluid 2 (
b) f) 2
bf 2
nb2
f 2
Figure 1. The trilayer model used for representing ultrasound axial transmission test. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
which is defined by QP D 1 F .t / ı.x1 /ı.x2 x2s /,
(2)
where F .t / is a given real scalar function depending only on time t and ı../ is Dirac’s delta function. f f Similarly to the description in the domain 1 , the pressure field p2 .x, t / in the fluid domain 2 reads as follows: 1 pR2 r 2 p2 D 0, c22
f
8 x 2 2 ,
(3)
p f where c2 D K2 =2 is the wave celerity in 2 in which the bulk modulus and mass density of the fluid at equilibrium are denoted by K2 and 2 , respectively. In the domain fj (j D 1, 2), the velocity vector vj .x, t / is relied to the gradients of pressure fields pj .x, t / following the Euler equation: f
j vP j C rpj D 0
8 x 2 j .
(4)
2.2. Equations in the anisotropic poroelastic layer (b ) The anisotropic poroelasticity theory is used to describe the behavior of bone layer. The poroelastic equations employed here are based upon Biot’s original works [37–39] as well as the recent developments in anisotropic constitutive equations [40–42]. The bone material is assumed to consist of a solid skeleton (with mass density s ) and a connected pore network saturated by fluid (with mass density f ). Neglecting the body forces (other than inertia), the equations describing the wave propagation within the bone layer read: r D uR s C f w R, s
1
rp D f uR C k
(5) w P Cbw R,
(6)
where D f C .1 / s is the mixture density, is the porosity, .x, t / is the total stress tensor and p.x, t / is the interstitial fluid pressure in the pores; the vectors of displacement of the solid skeleton and of the fluid are denoted by us .x, t / and uf .x, t /, respectively; the vector of relative displacement between the fluid and the solid frame weighted by the porosity is denoted by w D .uf us /; the tensor k is the anisotropic permeability tensor and the tensor b is defined as b D .f =/ a with a is the tortuosity tensor. Note that both tensors k and b are symmetric second-order tensors. The constitutive equations for the anisotropic linear poroelastic material are given by D C W ˛p,
(7)
1 p D r wC˛W, M
(8)
where C.x/ is the elasticity fourth-order tensor of drained porous material; ˛, which is a symmetric second-order tensor, is the Biot effective tensor; M is the Biot scalar modulus; .x, t / is the infinitesimal strain tensor, which is defined as the symmetric part of rus . For the subsequent development, it is more convenient to represent the following equations in vector form. Let us then use the Voigt’s notation, which expresses the symmetric second-order tensors as vectors, so the stress is denoted s D f11 , 22 , 12 gT , the strain by e D f11 , 22 , 212 gT and the Biot’s effective tensor by ˛L D f˛11 , ˛22 , ˛12 gT , where the superscript ?T designates the transpose operator. We may also introduce an operator L defined as follows: 2 3 2 3 1 0 0 0 (9) L D L1 @1 C L2 @2 , where L1 D 4 0 0 5 , L2 D 4 0 1 5 , 0 1 1 0 where @1 and @2 denote the partial differentiation operators with respect to x1 and x2 , respectively. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
Using these notations, the balance equations of linear momentum (5) and (6) may be rewritten as: R LT s D 0, uR s C f w s
1
f uR C k
(10) T
w P C bw R C L mp D 0,
(11)
where m D f1, 1, 0gT . The constitutive equations for poroelastic material (7) and (8) read: s D C e ˛L p, p D M mT Lw C ˛L T Lus , where C.x/ is the drained elastic tensor using Voigt’s notation: 3 2 c11 c12 c16 C D 4 c21 c22 c26 5 . c61 c62 c66
(12) (13)
(14)
By noting that e D L u and by substituting Equation (13) into Equation (12), one has s D Cu L us C C˛ L w, mp D CM Lw C CT˛ us ,
(15) (16)
where the quantities Cu , C˛ and CM are defined as: Cu D C C M ˛L ˛L T ,
(17)
T
L , C˛ D M ˛m
(18)
T
(19)
CM D M mm .
The tensor Cu is called the undrained elasticity tensor, which may be considered as the rigidity of an equivalent elastic medium in which the relative movement between the interstitial fluid and solid skeleton is null (i.e. when w D 0). Moreover, one may notice that while Cu and CM are symmetric, C˛ is not. In fact, C˛ is symmetric only if the considered poroelastic medium is isotropic (i.e. when ˛11 D ˛22 and ˛12 D 0). 2.3. Boundary conditions At both interfaces 1bf and 2bf , the continuity of pressure and stress fields between the porous medium and the fluid domains requires the following: p D pj , bf nbf j D pj nj ,
8x 2 jbf .j D 1, 2/,
(20)
bf b where nbf j is the normal unit vector to j pointing from the porous domain toward the fluid f domain j (see Figure 1). In addition, open-pore condition at the interfaces jbf (j D 1, 2) is assumed, requiring
P s / nbf w P nbf j D .v j u j ,
8x 2 jbf .j D 1, 2/.
In view of the Euler Equation (4), the interface condition (21) may be rewritten as follows: 1 s nbf rpj C w R C uR 8x 2 jbf .j D 1, 2/. j D 0, j Copyright © 2012 John Wiley & Sons, Ltd.
(21)
(22)
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
Note that for the present configuration (see Figure 1), the normal unit vectors of two interfaces 1bf T bf and 2bf are defined by nbf 1 D n2 D f0, 1g . As a consequence, the boundary conditions may be expressed as follows: 9 1 > @2 pj D .uR 2 C wR 2 / > = j (23) , 8x 2 jbf , .j D 1, 2/, p D pj > > ; t D f 0, pj gT pj ! 0,
f1
8x ! j
,
.j D 1, 2/.
(24)
where us D fu1 , u2 gT and w D fw1 , w2 gT ; the traction vectors t at the interfaces are defined by t D LT2 s D f21 , 22 gT . 3. EQUATIONS IN THE LAPLACE–FOURIER (s k1 ) DOMAIN The boundary value problem given by Equations (1), (3), (10), (11), (23) and (24) deals with solving a system of linear partial differential equations of which the coefficients are constant in the longitudinal direction defined by e1 -axis. Here, we propose to solve the system as follows: (i) the system of equations is firstly transformed into wavenumber-frequency domain by using a Fourier transform with respect to x1 combining to a Laplace transform with respect to t ; (ii) in the wavenumberfrequency domain, the equations in both fluid domains are analytically solved giving impedance boundary conditions for the solid layer that may be solved by the finite element method; (iii) the space–time solution is finally obtained by performing two inverse transforms. We recall that the general form of a Laplace–Fourier transform applied to a real-valued function y.x1 , x2 , t / is defined as: Z 1 Z C1 y.x1 , x2 , t /e ik1 x1 dx1 e st dt , (25) y.k Q 1 , x2 , s/ WD 0
1
p where i is the imaginary unit defined by i D 1; s 2 is the complex Laplace variable; k1 2 is the real Fourier variable representing the wavenumber on the e1 -axis; and denote the set of all the real and complex numbers, respectively. In (s k1 ) domain, the time derivative and the spatial differential operation with respect to x1 can be replaced by @.?/ ! s.?/ and @1 .?/ ! i k1 .?/, respectively. @t f
f
3.1. Transformed problem in (s k1 ) domain for the fluids 1 and 2
By applying the Laplace–Fourier transform (25) to Equations (1), (23) and (24), we obtain a boundary value problem whose differential equation is defined for pQ1 with respect only to x2 : 2 s 2 C k1 pQ1 @2 pQ1 D 1 FQ0 .s/ı.x2 x2s /, 8x2 > 0, (26) c12 @2 pQ1 D 1 s 2 .uQ 2 C wQ 2 /, pQ2 ! 0,
at x2 D 0, for x2 ! C1.
(27) (28)
The solution of this system may be presented in a semi-explicit form as follows: s 2 GQ 1 FQ0 ˇ1 .x s x2 / s 1 21 ˇ1 x2 2 e C e ˇ1 .x2 Cx2 / C e , 2ˇ1 ˇ1 s 2 GQ 1 FQ0 ˇ1 .x s x2 / s 1 21 ˇ1 x2 2 pQ1 D C e ˇ1 .x2 Cx2 / C e , e 2ˇ1 ˇ1 r 2 where ˇ1 D cs 2 C k12 and GQ 21 D uQ 2 .k1 , 0, s/ C wQ 2 .k1 , 0, s/. pQ1 D
8 0 x2 x2s
(29)
8 x2 x2s
(30)
1
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
Similarly, the solution of pQ2 in (s k1 ) domain may be expressed as follows: pQ2 D r where ˇ2 D
s2 c22
2 s 2 GQ 22 ˇ2 .x2 Ch/ e , ˇ2
8x2 h,
(31)
C k12 and GQ 22 D uQ 2 .k1 , h, s/ C wQ 2 .k1 , h, s/.
At two interfaces x2 D 0 and x2 D h, the solutions of pQ1 and pQ2 are given by 1 FQ0 ˇ1 x s 1 s 2 GQ 21 2 C e , ˇ1 ˇ1 2 s 2 GQ 22 pQ2 .k1 , h, s/ D . ˇ2
pQ1 .k1 , 0, s/ D
(32) (33)
3.2. Transformed problem in (s k1 ) domain for the solid b Q D i k1 L1 C @2 L2 , the By noting that in the Fourier–Laplace domain, the operator L becomes L boundary valued problem defined by Equations (10) and (11) may be transformed into the (s k1 ) domain as follows: s 2 uQ s C f s 2 w Q i k1 LT1 sQ @2 LT2 sQ D 0 , 8x2 2 h, 0Œ, (34) Q C i k1 LT1 mpQ C LT2 @2 .mp/ Q D0 f s 2 uQ s C s k1 C s 2 b w or
s 2 1 f s 2 1 f s 2 1 s k1 C s 2 b
s LT1 sQ LT2 sQ uQ 0 @ D , i k1 2 0 w Q LT1 mpQ LT2 mpQ
(35)
where 1 is the 2 2 identity matrix and (see Equations (15) and (16)) sQ D i k1 Cu L1 uQ s C i k1 C˛ L1 w Q C Cu L2 @2 uQ s C C˛ L2 @2 w, Q mpQ D
i k1 CT˛ L1 uQ s
C i k1 CM L1 w Q C CT˛ L2 @2 uQ s
C CM L2 @2 w. Q
(36) (37)
The vectors concerning sQ and .mp/ Q in the system (35) may be evaluated in terms of uQ s and w Q as follows: 11
s 12
s C11 C12 Cu Cu uQ uQ LT1 sQ ˛ ˛ D i k1 C , (38) T 11 21 T 12 @2 w w Q Q LT1 mpQ C C C11 C ˛ ˛ M M 21
s 22
s C21 C22 Cu Cu uQ uQ LT2 sQ ˛ ˛ @2 D i k C , (39) T 1 21 22 T 22 w Q w Q LT2 mpQ C C C12 C ˛ ˛ M M ab ab where the 2 2 matrices Cab u , C˛ and CM with a, b D 1, 2 are defined by: T Cab u D La Cu Lb ,
T Cab ˛ D La C˛ Lb ,
T Cab M D La CM Lb
with a, b D 1, 2.
(40)
By replacing (38) and (39) into Equation (35) and then making some arrangements, we obtain a vectorial form of this problem: s 2 A1 Ú C sA2 Ú C k12 A3 Ú i k1 A4 @2 Ú @2 Ø D O,
(41)
where s Q Ú D uw , Q Copyright © 2012 John Wiley & Sons, Ltd.
LT2 sQ D i k1 AT4 Ú C A5 @2 Ú, ØD LT2 mpQ
(42)
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
and
0 0 1 f 1 A1 D , , A2 D f 1 b 0 k1
11 C11 Cu C12 ˛ u T , A D A3 D 11 4 21 T C11 C˛ C ˛ M
(43) 12
C22 A5 D 22u T C˛
C˛ , C12 M
22
C˛ . C22 M
(44)
The 4 4 matrices A1 , A2 , A3 , A4 and A5 depend only on the material properties of the elastic porous material. A physical interpretation of these matrices can be given: A1 relates to the mass density; A2 relates to the viscous property via the permeability; A3 , A4 and A5 relate to the poroelastic properties. It is worth noting that whereasA1 , A2 , A3 and A5 are symmetric, the matrix A4 is not because we have seen that the matrix C˛ defined in Equation (18) is nonsymmetric when the poroelastic medium is anisotropic. At two interfaces x2 D 0 and x2 D h, the continuity condition (23) between fluid pressures and stress traction require that LT2 sQ PQ 0 Q where P D (45) D Q pQ LT2 mpQ x D0,h P 2
where pQ D pQ1 at x2 D 0 and pQ D pQ2 at x2 D h. In view of Equations (32) and (33), the boundary conditions of the system (41) may be expressed by:
Ø.0/ D P1 Ú.0/ C F0 , Ø.h/ D P2 Ú.h/,
(46)
where 2
P1 D
1 s P, ˇ1
n F0 D 0,
1 FQ0 ˇ1
2
P2 D
s
2 s P, ˇ2
e ˇ1 x2 , 0,
1 FQ0 ˇ1
2 0 60 PD4 0 0 oT s
e ˇ1 x2
0 1 0 1
0 0 0 0
3 0 17 , 05 1
.
(47)
(48)
4. FINITE ELEMENT FORMULATION IN SPECTRAL DOMAIN AND TIME–SPACE SOLUTION 4.1. Weak formulation The weak formulation of the boundary value problem given by (41) and (46) may be now carried out using a usual procedure. Let C ad be the admissible function space of the problem constituted by the sufficiently differentiable complex-valued functions such as: x2 2 H b D h, 0Œ ! ı Ú.x2 / 2 . The conjugate transpose of ı Ú is denoted by ı Ú . Upon integrating (41) against a test vector function ı Ú and integrating by part, then applying the boundary condition (46), the weak formulation of the differential equation (41) reads: For all k1 fixed in and for all s fixed in , find Ú.k1 , x2 , s/ 2 C ad such that:
s 2 h ı Ú, A1 Ú iH b C s h ı Ú, A2 Ú iH b C k12 h ı Ú, A3 Ú iH b ˚˝ ˛ C i k1 @2 ı Ú, AT4 Ú H b h ı Ú, A4 @2 Ú iH b C h @2 ı Ú, A5 @2 Ú iH b
(49)
C ı Ú .0/ P1 Ú.0/ C ı Ú .h/ P2 Ú.h/ D ı Ú .0/ F0 ,
8 ı Ú 2 C ad ,
where h., .iH b denotes the inner product over H b , for instance, h ı Ú, A1 Ú iH b R0 h ı Ú.x2 / A1 Ú.x2 / dx2 . Copyright © 2012 John Wiley & Sons, Ltd.
D
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
4.2. Finite element formulation We proceed by introducing a finite element mesh of the domain Œh, 0, which contains nel eleS ments: Œh, 0 D e e with e D 1, : : : , nel . By using the Galerkin finite element method, both functions Ú and ı Ú in each element e are approximated using the same shape function:
Ú.x2 / D Ne Ve , ı Ú.x2 / D Ne ıVe , 8 x2 2 e ,
(50)
where Ne is the shape function, Ve and ıVe are the vectors of nodal solutions of Ú and ı Ú within the element e , respectively. Replacing Equation (50) into Equation (49) and assembling the elementary matrices lead to a static-like linear system of equations for each couple values (s, k1 ): (51) Kb C K V D F, where V is the global nodal solution vector; Kb is the global ‘stiffness matrix’ of the poroelastic layer; K represents the coupled operator between the fluid and poroelastic layers; the vector F is the external force vector due to the acoustical source. Noting the number of nodes by N , the size of vectors V and F is neq D 4N because each node has four degrees of freedom. The sizes of Kb and K are neq neq . For all couples .s, k1 / fixed in , these quantities may be expressed by:
Kjk D
8 ˆ ˆ < ˆ ˆ : (
Fj D
1 s 2 ˇ1
if .j , k/ D .2, 2/, .2, 4/, .4, 2/, .4, 4/
2 s 2 ˇ2
if .j , k/ D .neq 2, neq 2/, .neq 2, neq /, .neq , neq 2/, .neq , neq /
0
otherwise
1 FQ0 ˇ1
ˇ1 x2s
e
0
if j D 2, 4
(53)
otherwise
Kb D s 2 K1 C s K2 C k12 K3 C i k1 K4 C K5
(54)
where the matrices K1 , K2 , K3 , K4 and K5 are independent to (s, k1 ) and are defined by: [Z [Z K1 D NTe A1 Ne dx2 I K2 D NTe A2 Ne dx2 , e
K3 D
[Z e
K5 D
(52)
[Z e
e
e
e
NTe A3 Ne
dx2 I
K4 D
[Z e
(55)
e
n
.@2 Ne /T A4 Ne
e
.@2 Ne /T A5 .@2 Ne / dx2 ,
o a
dx2 ,
(56) (57)
e
in which the notation f?ga designates the anti-symmetric part of the f?g. 4.3. Computation of time–space solution For fixed values of .s, k1 / in the Laplace–Fourier transformed domain, the solution of Ú may be computed by solving the system of linear equations (51) in the complex domain. The solutions for pQ1 and pQ2 in two fluid domains may then be determined by using the Equations (29), (30) and (31), respectively. To obtain the spatio–temporal solution, we need to perform a numerical inverse Laplace–Fourier transform. In this paper, the inverse Fourier transform is computed by using the usual fast Fourier transform technique. The inverse Laplace transform is carried out using the Quadrature convolution method, which has proved to be a very efficient technique for computing the time response solution in many dynamic problems [43]. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
5. NUMERICAL TESTS This section presents some numerical tests describing an in vivo ultrasound test on human cortical long bones. The acoustic source (2) located at the position .x1s , x2s / D .0, 2/ .in millimeter/ in the upper fluid domain f1 has the time-history function given by: 2
F .t / D F0 e 4.fc t1/ sin.2 fc t /,
(58)
where F0 = 1 m.s2 and fc D 1 MHz which is the central frequency. Using a Fourier transform shows that this signal has a spectrum, which contains a frequency bandwidth about 0 2.5 MHz (the reader could refer to Figure 2 in [16]). For all tests presented later, the total duration of the simulation will be taken by T D 2 105 s. 5.1. Physical parameters f
f
Both fluid domains 1 and 2 are assumed to be identical and are considered as an idealized water. The mechanical properties of the fluids are given by 1 D 2 D 1000 kg.m3 and K1 D K2 D 2.25 GPa. Thus the wave velocities in two fluid domains are: c1 D c2 D 1500 m.s1 . To describe the behavior of the poroelastic bone plate, the drained elasticity tensor C as well as Biot’s effective coefficients ˛ and M used in (7) and (8) should be provided. For this study, these parameters are derived from the characteristics of the interstitial fluid and solid skeleton phases by using a continuum micromechanics model proposed by Hellmich et al. [44]. According to this model, the micropores at the microstructural scale are regarded as cylindrical pores with a circular cross section. In drained condition, the constitutive behavior of the material inside the pores does not possess stiffness. Hence, the estimated drained microstructural stiffness of the bone whose solid bone matrix’s elasticity tensor is C m reads: h i1 1 C D .1 /C m W .1 /I C I P cyl W C m
(59)
where I denotes the fourth-order identity tensor; the fourth-order tensor P cyl is the Hill’s tensor for materials with periodical cylindrical inclusions, which may be derived in closed form [44]. The tensor ˛ and the constant M can be then evaluated by [40, 42]: ˛ D I C W .Sm W I/ ,
1 D C ˛ W S W ˛, M
(60)
where S D .C/1 and Sm D .C m /1 are respectively the drained and solid material compliance tensors, I designates the second-order tensor identity and the scalar C denotes the effective compressibility of porous matrix material, which is given by 1 1 CD C K Km
1 1 m K Kf
,
(61)
where K D .I W S W I/1 , K f and K m D .I W Sm W I/1 are respectively the bulk moduli of the drained porous matrix, of the interstitial fluid and of the poroelastic material. For the numerical tests presented in this paper, the interstitial fluid is also assumed to be the same as the exterior fluid domain, that is f D 1000 kg.m3 and Kf D 2.25 GPa. The solid bone has the mass density s D 1722 kg.m3 and is assumed to be a transversely isotropic elastic material m [44]. The components of the elastic tensor are given by using the Voigt’s notations: c11 D 28.7 GPa, m m m m m c22 D 23.6 GPa, c12 D 9.9 GPa, c66 D 7.25 GPa, c16 D c26 D 0. The viscosity of interstitial fluid is D 1 103 Pa.s and the permeability tensor is roughly taken by 11 D 22 D 5 1013 m2 [45]. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
5.2. Validation This section will present some calculations to validate the proposed method. First, the hybrid spectral/FE solutions will be compared with the ones obtained by using a conventional finite element procedure. Then, the behavior of the solution vis-à-vis the parameters of numerical discretization will be studied. We consider a bone layer, which has a constant thickness h D 4 mm. The numerical parameters used for simulation have been chosen in a similar way to which required for a classical procedure of dynamic finite element analysis [24]. The space interval x is about the 1=10 shortest wavelength. The time step is chosen for satisfying the Courant–Friedrichs–Lewy condition. To obtain a boundedvalue solution in space on the e 1 -axis, the domain size should be chosen sufficiently large to avoid the fact that fastest waves pass over the boundaries after the duration T . Here, the space interval in e 1 -direction has been chosen to be x D 0.75 104 m. The time step used is t D 107 s. To carry out the finite element analysis on vertical direction e 2 , the bone plate’s thickness is discretized into 12 three-noded quadratic isoparametric Lagrangian elements. Figure 2 presents snapshots of the wave field in the coupled trilayer system at different instants t D 1 s, t D 3 s and t D 6 s. At each point x in the domain, the quantity presented for obtaining f f these snapshots is log.jp1 .x, t /j/ if x 2 1 , log.jp2 .x, t /j/ if x 2 2 and log.jp.x, t //j if x 2 b . In these graphs, we can clearly observe the transmission and reflection of ultrasonic waves through the interfaces. As the longitudinal waves propagating in the poroelastic layer are much faster than acoustic waves in the fluid, the wavefronts of two symmetrical head waves, which have truncated cone forms, are very well displayed in the last two graphs in Figure 2. Note that using the proposed method allows us to take into account the exact radiation conditions of two fluid domains, therefore the results involve only the waves diffracted by the bone layer. Figure 3 presents the numerical solution of p1 measured at two positions R1 .x1 , x2 / D .2, 2/ (mm) (upper) and R2 .x1 , x2 / D .20, 2/ (mm) (lower) by using the proposed spectral/finite element procedure. We also present numerical solutions obtained by using conventional time domain finite element analysis [24]. The comparison shows that S/FEM and FEM results perfectly match each other. While about 5 h are needed on a common desktop personal computer (using COMSOL
Figure 2. Snapshots of the fluid pressure (p1 , p and p2 ) in the trilayer. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
p1(x1,t) (Pa)
x1 = 0.002 m SFEM FEM
50 0 −50 0
0.5
1
1.5
x 10−5
time (s)
p1(x1,t) (Pa)
x1 = 0.02 m SFEM FEM
20 0 −20 0.6
0.8
1
1.2
1.4
1.6
1.8
2
x 10−5
time (s)
Figure 3. Comparison between spectral finite element and finite element solutions of p1 at x1 D 2 mm (upper graph) and at x1 D 20 mm (lower graph) from the source.
multiphysics software (COMSOL, Inc., Stockholm, Sweden) [46]) for simulation using the dynamic finite element method, the proposed spectral/finite element method needs only about 400 s on the same personal computer. Note that both methods use the same time step for the calculation. The computational advantage of analysis in spectral domain over the analysis using direct time FEM is mainly due to the fact that only very small matrices need to be factorized. Because the exact solution is not available for this example, we do not calculate the numerical p1 errors for convergence analysis, but investigate the behavior of a function L2,T .x/, which is defined as follows, with respect to the parameters of numerical dicretization (Nt , N1 , nel ): v u Nt u X p1 L2,T .x/ D t t .p1 .x, n t //2 , (62) nD0 1 Figure 4 presents a study with respect to different parameters by calculating Lp2,T at two points R1 and R2 . Figure 4(A) shows the effect of sampling number N1 , which is associated with the wavenumber k1 , by fixing the parameters Nt D 2048 and nel D 12. We recall that Nt is related to the total number of time step and nel is related to the element number. This last parameter combined with the parameter N1 is related to the discretization in space. Next, in Figure 4(B), the parameters N1 D 2048 and nel D 12 are fixed, and the effect on the convergence of the total number of
(A) Nt = 2048; nel = 12
(B) N1 = 2048; nel = 12 0.06
0.06 0.04 0.02
8
10
12
log2(N1)
0.04 0.02 0 8
x1 = 0.002m
x1 = 0.02m
x1 = 0.02m
0.04
Lp2,1T(x1)
x1 = 0.002m
x1 = 0.02m
Lp2,1T(x1)
p
L2,1T(x1)
x1 = 0.002m
0
(C) N1 = 2048; Nt = 2048 0.05
0.03 0.02
10
log2(Nt)
12
0.01
5
10
15
20
number of elements
1 Figure 4. Lp 2,T versus (A) number of space interval N1 , (B) number of time steps Nt and (C) the number of elements.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
time step Nt is shown. Note that the total number of time step Nt is related to time step size t . 1 Last, Figure 4(C) presents the values of Lp2,T versus the number of finite elements nel by fixing N1 D 2048 and Nt D 2048. These studies show that the proposed procedure using spectral finite elements may have fast convergence when refining parameters of numerical discretization. Moreover, it can be deduced that, for this numerical example and for the observations point located at R1 and R2 , convergence is reached for Nt D 2048, N1 D 2048 and nel D 12 which correspond well to the parameters used for obtaining results presented in Figure 3. 5.3. Influence of porosity on first arriving signal velocity When the axial transmission technique is used, the earliest event or wavelet usually called first arriving signal (FAS) of the multicomponent signal recorded at the receivers has been the most often investigated. In fact, the FAS velocity have been shown to be a relevant index of bone status [47–49]. In this section, we present some numerical results to study the influence of the bone’s porosity on the FAS velocity denoted by VF . For this purpose, the material properties of bone matrix (m , C m ) were fixed. The poroelastic properties were then determined by using the micromechanics analysis presented in Section 5.1. To determine VF , the p1 -signal is captured at an array of 14 sensors, which has a gap 0.8 mm from each other. The distance from the emitted source to the first sensor is 11 mm. At each sensor, the first zero-crossing location on time axis of measured p1 -signal may
sensor number #
15
10
5
0 0
0.5
1
1.5
2
x 10−5
time (s)
Figure 5. Case of a bone layer with thickness h D 4 mm and D 5%: seismograph of p1 .t / measured at 14 sensors. The gray line (red line if the color is available) is the linear regression of zero-crossing positions. bone thickness h = 0.6mm 3900
4000
3850
FAS (ms−1)
FAS (ms−1)
bone thickness h = 4mm 4050
3950 3900 3850
3800 3750 3700 poroelastic model undrained elastic model
poroelastic model undrained elastic model 3800
0
5
Porosity
10
(%)
15
3650 0
5
Porosity
10
15
(%)
Figure 6. Comparison of first arriving signal velocities for different bone porosities in two cases of bone’s thickness: h D 4 mm (left) and h D 0.8 mm (right). Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
be determined. Then VF is evaluated by the slope of the linear regression of these location over 14 sensors. An example of the results obtained at 14 sensors is shown in Figure 5. Figure 6 (left) shows the variation of the FAS velocity with respect to when considering a bone plate with thickness h D 4 mm. One may observe that the VF versus relation is practically linear and VF decreases when the porosity is higher. We also present in this graph the VF versus relation obtained by using equivalent elastic media. For a given porosity , the corresponding elastic models use the mass density taken as the mixture density given by (5), and the elasticity tensor is taken as the undrained elasticity tensor Cu defined by (17). For lower porosities, the VF obtained by using the poroelastic model are slightly different from the ones obtained using the elastic model. The difference becomes more significant for higher porosity because the fluid–solid movement is more important and would not be negligible. A similar discussion may be made when considering a thin bone layer (h D 0.8 mm) in Figure 6 (right). 6. CONCLUSION By using a hybrid semi-analytical/finite element technique, we have derived a finite element formulation in spectral domain for anisotropic poroelastic plate saturated by an interstitial fluid and immersed in a fluid. The procedure is based on a Laplace transform for time and a Fourier transform for the spatial variable in longitudinal direction defined by x1 . Performing these transforms leads to a one-dimensional differential equation with respect only to x2 , which may be solved by using conventional Galerkin finite element method. Note that in the Laplace–Fourier domain, the radiation conditions presenting two half-spaces of fluid may exactly be introduced. The validation has been carried out. Because only one-dimensional finite element problem need to be solved in spectral domain, the proposed procedure requires a very small memory size for computation. In addition, the computation time is very low, especially for the configuration used for ultrasonic test presented here, wherein we only need to compute the time-domain solution at some specific points. It has also been shown that the use of conjugated quadrature method for the inverse Laplace transform has very good stability when computing the time-domain solutions. Using the proposed method, numerical studies of ultrasonic wave propagation in poroelastic bone layer may be easily carried out, eventually in the very high frequency domain. This allows us to consider the anisotropy and the heterogeneity of the bone matrix material as well as bone porosity. Some preliminary results presented in this paper show that the FAS velocity is strongly influenced by the bone porosity. It seems that using of equivalent elastic model might lead to a poor estimation of FAS velocity in cortical bone with higher bone porosity. The detailed investigation of behavior of (heterogeneous) poroelastic cortical bone under broadband frequency ultrasound will be presented in a forthcoming paper. REFERENCES 1. Lowet G, Van der Perre G. Ultrasound velocity measurements in long bones: measurement method and simulation of ultrasound wave propagation. Journal of Biomechanics 1996; 29:1255–1262. 2. Nicholson P, Moilanen P, Karkkainen T, Timonen J, Cheng S. Guided ultrasonic waves in long bones: modelling, experiment and in vivo application. Physiological Measurement 2002; 23(4):755–768. 3. Moilanen P. Ultrasonic guided waves in bone. IEEE transactions of ultrasonics, ferroelectrics, and frequency control 2008; 55(6):1277–1286. 4. Bousson V, Meunier A, Bergot C, Vicaut E, Rocha M, Morais M, Laval-Jeantet A, Laredo J. Distribution of intracortical porosity in human midfemoral cortex by age and gender. Journal of Anatomy 2001; 16:1308–1317. 5. Thomas C, Feik SA, Clement J. Regional variation of intracortical porosity in the midshaft of the human femur: age and sex differences. Journal of Anatomy 2005; 206:115–125. 6. Sansalone V, Naili S, Bousson V, Bergot C, Peyrin F, Zarka J, Laredo JD, Haiat G. Determination of the heterogeneous anisotropic elastic properties of human femoral bone: from nanoscopic to organ scale. Journal of Biomechanics 2010; 43:1857–1863. 7. Dong NX, Guo EX. The dependence of transverse isotropic elasticity of human femoral cortical bone on porosity. Journal of Biomechanics 2004; 37(8):1281–1287. 8. Baron C, Talmant M, Laugier P. Effect of porosity on effective diagonal stiffness coefficients (ci i ) and anisotropy of cortical bone at 1 MHz: a finite difference time domain study. Journal of the Acoustical Society of America 2007; 122:1810–1817. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
9. Bossy E, Talmant M, Laugier P. Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models. Journal of the Acoustical Society of America 2004; 115:2314–2324. 10. Baron C, Naili S. Elastic wave propagation in a fluid-loaded anisotropic waveguide with laterally varying properties. Comptes Rendus Mecanique 2008; 336:772–730. 11. Baron C, Naili S. Propagation of elastic waves in a fluid-loaded anisotropic functionally graded waveguide: application to ultrasound characterization. Journal of the Acoustical Society of America 2010; 127(3):1307–1317. 12. Nayfeh HA. Wave Propagation in Layered Anisotropic Media. Elsevier: North-Holland, 1995. 13. Shuvalov A, Poncelet O, Deschamps M. Analysis of the dispersion spectrum of fluid-loaded anisotropic plates: leaky-wave branches. Journal of Sound and Vibration 2006; 296(3):494–517. 14. Bossy E, Talmant M, Laugier P. Effect of bone cortical thickness on velocity measurements using ultrasonic axial transmission: A 2D simulation study. Journal of the Acoustical Society of America 2002; 112:297–307. 15. Grimal Q, Naili S. A theoretical analysis in the time-domain of wave reflection on a bone plate. Journal of Sound and Vibration 2006; 298:12–29. 16. Desceliers C, Soize C, Grimal Q, Haiat G, Naili S. A time-domain method to solve transient elastic wave propagation in a multilayer medium with a hybrid spectral-finite element space approximation. Wave Motion 2008; 45(4):383–399. 17. Haiat G, Naili S, Grimal Q, Talmant M, Desceliers C, Soize C. Influence of a gradient of material properties on ultrasonic wave propagation in cortical bone: application to axial transmission. Journal of the Acoustical Society of America 2009; 125(6):4043–4052. 18. Cowin SC. Bone poroelasticity. Journal of Biomechanics 1999; 32:217–238. 19. Zhang D, Cowin SC. Oscillatory bending of a poroelastic beam. Journal of the Mechanics and Physics of Solids 1994; 42(10):1575–1599. 20. Nguyen VH, Lemaire T, Naili S. Poroelastic behaviour of cortical bone under harmonic axial loading: a finite element study at the osteonal scale. Medical Engineering & Physics 2010; 32(4):384–390. 21. Nguyen VH, Lemaire T, Naili S. Influence of interstitial bone microcracks on strain-induced fluid flow. Biomechanics and Modeling in Mechanobiology 2011; 10(6):963–972. 22. Williams JL. Ultrasonic wave propagation in cancellous and cortical bone: prediction of some experimental results by Biot’s theory. Journal of the Acoustical Society of America 1992; 91(2):1106–1112. 23. Hosokawa A. Simulation of ultrasound propagation through bovine cancellous bone using finite-difference timedomain methods. Journal of the Acoustical Society of America 2005; 118(3):1782–1789. 24. Nguyen VH, Naili S, Sansalone V. Simulation of ultrasonic wave propagation in anisotropic cancellous bones immersed in fluid. Wave Motion 2010; 47(2):117–129. 25. Nguyen VH, Naili S, Sansalone V. A closed-form solution for in vitro transient ultrasonic wave propagation in cancellous bone. Mechanics Research Communications 2010; 37(4):377–383. 26. Kausel E. Fundamental Solution in Elastodynamics. Cambridge University Press: New York, 2006. 27. Virieux J. P-SV-wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics 1986; 51(4):889–901. 28. Kampanis NA, Dougalis VA, Ekaterinaris JA. Effective Computational Methods for Wave Propagation. Chapman & Hall/CRC: New York, 2008. 29. Protopappas VC, Kourtis IC, Kourtis LC, Malizos KN, Massalas CV, Fotiadis DI. Three-dimensional finite element modeling of guided ultrasound wave propagation in intact and healing long bones. Journal of the Acoustical Society of America 2007; 121(6):3907–3921. 30. Thompson LL, Huan R. Implementation of exact non-reflecting boundary conditions in the finite element method for the time-dependent wave equation. Computer Methods in Applied Mechanics and Engineering 2000; 187(1-2): 137–159. 31. Givoli D. High-order local non-reflecting boundary conditions: a review. Wave Motion 2004; 39(4):319–326. 32. Han X, Liu GR, Lam KY. Transient waves in plates of functionally graded materials. International Journal for Numerical Methods in Engineering 2001; 52(8):851–865. 33. Han X, Liu GR, Xi ZC, Lam KY. Transient waves in a functionally graded cylinder. International Journal of Solids and Structures 2001; 38(17):3021–3037. 34. Liu G, Xi ZC. Elastic Waves in Anisotropic Laminates. CRC Press: Florida, 2002. 35. Gopalakrishnan S, Chakraborty A, Roy Mahapatra D. Spectral Finite Element Method. Springer-Verlag: London, 2008. 36. Marzani A. Time-transient response for ultrasonic guided waves propagating in damped cylinders. International Journal of Solids and Structures 2008; 45(25-26):6347–6368. 37. Biot MA. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics 1955; 26(2):182–185. 38. Biot MA, Willis DG. The elastic coefficients of the theory of consolidation. Journal of Applied Mechanics 1957; 79:594–601. 39. Coussy O. Poromechanics. John Wiley & Sons, Ltd: England, 2004. 40. Cheng AH. Material coefficients of anisotropic poroelasticity. International Journal of Rock Mechanics and Mining Sciences 1997; 34(2):199–205. 41. Cowin SC. A recasting of anisotropic poroelasticity in matrices of tensor components. Transport in Porous Media 2003; 50(1):35–56.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
42. Thompson M, Willis JR. A reformation of the equations of anisotropic poroelasticity. Journal of Applied Mechanics 1991; 58:612–616. 43. Schanz M, Antes H. Application of ’Operational Quadrature Methods’ in time domain boundary element methods. Meccanica 1997; 32:179–186. 44. Hellmich C, Ulm FJ. Microporodynamics of bones: prediction of the “Frenkel–Biot” slow compressional wave. Journal of Engineering Mechanics 2005; 131(9):918–927. 45. Cowin SC. Bone Mechanics Handbook (2nd edn). CRC Press: Boca Raton, FL, 2001. 46. COMSOL Multiphysics 3.5. User’s Guide, 2008, www.comsol.com 47. Barkmann R, Kantorovich E, Singal C, Hans D, Genant HK, Heller M, Gluer CC. A new method for quantitative ultrasound measurements at multiple skeletal sites. Journal Of Clinical Densitometry 2000; 3:1–7. 48. Hans D, Srivastav SK, Singal C, Barkmann R, Njeh CF, Kantorovich E, Gluer CC, Genant HK. Does combining the results from multiple bone sites measured by a new quantitative ultrasound device improve discrimination of hip fracture? Journal of Bone and Mineral Research 1999; 14(4):644–651. 49. Stegman MR, Heaney RP, Travers-Gustafson D, Leist J. Cortical ultrasound velocity as an indicator of bone status. Osteoporosis International 1995; 5(5):349–533.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
Simulation of ultrasonic wave propagation in anisotropic poroelastic bone plate using hybrid spectral/finite element method Vu-Hieu Nguyen*,† and Salah Naili Laboratoire Modélisation et Simulation Multi Echelle, Université Paris-Est, MSME UMR 8208 CNRS, 61 Avenue du Général de Gaulle, 94010 Créteil, France
SUMMARY This paper deals with the modeling of guided waves propagation in in vivo cortical long bone, which is known to be anisotropic medium with functionally graded porosity. The bone is modeled as an anisotropic poroelastic material by using Biot’s theory formulated in high frequency domain. A hybrid spectral/finite element formulation has been developed to find the time-domain solution of ultrasonic waves propagating in a poroelastic plate immersed in two fluid halfspaces. The numerical technique is based on a combined Laplace–Fourier transform, which allows to obtain a reduced dimension problem in the frequency– wavenumber domain. In the spectral domain, as radiation conditions representing infinite fluid halfspaces may be exactly introduced, only the heterogeneous solid layer needs to be analyzed by using finite element method. Several numerical tests are presented showing very good performance of the proposed procedure. A preliminary study on the first arrived signal velocities computed by using equivalent elastic and poroelastic models will be presented. Copyright © 2012 John Wiley & Sons, Ltd. Received 6 July 2011; Revised 7 December 2011; Accepted 12 December 2011 KEY WORDS:
ultrasound; cortical bone; poroelastic; functionally graded plate; hybrid spectral/finite element
1. INTRODUCTION In recent years, quantitative ultrasound techniques for assessing the properties of bone have received much attention owing to its potential in estimating bone fragility and/or fracture risk. Ultrasound is non-ionizing, and moreover, ultrasonic devices are relatively inexpensive and can be made portable. From the point of view of mechanics, bone is a heterogeneous, anisotropic and porous material. There are two kinds of bone: cortical (compact) bone and cancellous (trabecular or spongy) bone. The porosity varies from 5% to 15% in cortical bone and from 50% to 90% in cancellous bone. For measuring in vivo material properties of cortical long bones, a so-called ‘axial transmission’ (AT) technique has been developed [1]. The axial transmission technique uses a set of ultrasonic transducers (transmitters and receivers) placed in the same side of the investigated skeletal site along a direction close to the long bone axis. Consequently, the transducers form a line in contact with the skin along the bone axial axis. The transmitter emits an ultrasound pulse wave (around 250 KHz 2 MHz) that propagates along the cortical layer of bones. The analysis of the signals received at the receivers serve for estimating the geometric parameters as well as mechanical characteristics of the cortical bone at the measured skeletal site [2, 3]. Mechanical modeling of this experiment deals with describing a vibro-acoustic problem of a solid waveguide (which represents cortical bone) coupled with two fluid media (which represents soft tissues such as skin or marrow). Although in vivo cortical bone may have slightly varied *Correspondence to: Vu-Hieu Nguyen, Université Paris-Est, Laboratoire Modélisation et Simulation Multi Echelle, MSME UMR 8208 CNRS, 61 Avenue du Général de Gaulle, 94010 Créteil, France. † E-mail: [email protected] Copyright © 2012 John Wiley & Sons, Ltd.
V.-H. NGUYEN AND S. NAILI
thickness along the bone axis, the bone’s thickness is usually assumed to be constant. Therefore, the cortical bone may be described as plate-like or cylindrical-like structures. Assessment of the bone’s material properties requires careful analysis of the reflections, conversion modes and interferences of longitudinal and shear waves within the bone structure. However, current development of ultrasound techniques is still limited because the information potentially contained in the transmitted ultrasonic waves is not fully exploited, and parameters such as bone material properties or microstructural parameters are still difficult to recover. The interaction between ultrasound and bone remains poorly understood from a physical point of view owing to the complex nature of bone, which has a viscoelastic porous microstructure spanning many length scales. At the macroscopic scale, porosity in the radial direction (which is defined in the bone’s crosssection plan) is heterogeneous. The observations at all ages and for both genders show that the mean porosity in the endosteal region (inner part of the bone) is significantly higher than the porosity in the periosteal region (outer part of the bone) [4–6]. This observation may be explained by the fact that cortical bone is affected by age-related bone resorption and osteoporosis, causing reduction of bone shell thickness as well as increase of porosity, namely in the endosteal region. Moreover, the macroscopic mechanical properties of bone have been shown to strongly depend on its porosity [7, 8]. As a consequence, cortical bone may naturally be considered as a functionally graded material. It has been shown that the longitudinal curvature of bone had minor effects in estimating the first arrived signal (FAS) (see e.g. [9]). Therefore, many of the past studies have focused on the modeling of guided waves in long bones by using fluid-loaded homogeneous/multilayer/functional graded plate models. The understanding of wave phenomena involved in the multilayer structures has been studied by many authors in the frequency domain [10–13] or in the time domain [14–17]. In these studies, the cortical bone material has been considered as an equivalent (visco-)elastic medium of which the effective macroscopic mass density and effective macroscopic elasticity tensor are estimated from its porosity. The presence of the interstitial fluid, which was considered in many works for the analysis of the behavior of cortical bone tissue under low frequency loading (e.g. [18–21]) or of ultrasonic wave propagation through cancellous bones (e.g. [22–25]), has usually been neglected when studying ultrasonic wave propagation in cortical bones. To solve the time-domain wave propagation problem in the multilayer structures, there are mainly two approaches. The first one involves using (semi-)analytical methods such as the generalized ray/Cagniard-de-Hoop technique [15] or the direct stiffness matrix method [26]. The second one involves using numerical methods such as the finite difference method [14, 27, 28] or the finite element method (FEM) [17, 29]. Although the analytical methods are attractive to obtain reliable transient responses of the structure, numerical methods are often more efficient to treat problems with inhomogeneous materials or complex geometries. However, most numerical methods require important computational costs, especially for problems in the high-frequency domain. Moreover, absorbing boundary conditions are required for considering unbounded domains [28, 30, 31]. Alternatively, for considering waveguides with geometrical and mechanical properties, which are constant only along one or two directions, the hybrid numerical method (see for example [32–34]), alternatively called spectral finite element method (see for example [35, 36]), have been employed. The key point of this method consists in using a hybrid algorithm, which begins by employing the Fourier transform (with respect to time and to the longitudinal direction of the waveguide) to transform the problem into the frequency-wavenumber domain. Then, the wave equations in the spectral domain governed in cross section (or even a 1D domain in the case of infinite plates or axisymmetric waveguides), which may actually have inhomogeneous material properties, can be easily handled using the FEM [16, 32, 33, 36]. Note that this technique uses the conventional finite elements in spectral domain and does not use high-order finite element families known as ‘spectral elements’. The objective of this paper is twofold. The first objective is to present a numerical procedure to simulate the ultrasonic wave propagation in bone plate immersed in fluid. The bone plate is not modeled as an equivalent elastic medium but as an anisotropic poroelastic medium by using the Biot’s theory in high frequency domain. The numerical method is based on a spectral finite element technique, which allows to consider the variation of bone properties along its thickness direction. As far as the authors know, no equivalent numerical procedure based on spectral finite element technique has been presented previously in the literature for an anisotropic poroelastic medium. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
The second objective is to check if using equivalent macroscopic elastic models [10, 11, 14–17] is sufficient to evaluate the bone properties for all porosities in physiological ranges. The paper is organized as follows. After this introduction on the context of this study, as well as on methods to solve the time-domain wave propagation problem in the multilayer structures, Section 2 provides the governing equations of an anisotropic poroelastic bone plate immersed in fluid. Next, Sections 3 and 4 present the problem governed in Laplace–Fourier domain and develop the finite element formulation in this domain. Then, some numerical tests are presented in Section 5 to validate the proposed procedure. We will also present some comparison of the FAS velocities computed by using equivalent elastic and poroelastic models. Finally, conclusion and discussion will be presented in Section 6. 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 two-dimensional (2D) space, respectively. The symbol ‘’ will denote the scalar product, and the symbol ‘W’ between tensors of any order denotes their double contraction. In this work, we will employ a 2D configuration, which has been shown to be an appropriate geometrical model for describing the ultrasound wave in long bone in the framework of the axial transmission technique [9]. Let R.OI e1 , e2 / be the reference Cartesian frame where O is the origin and .e1 , e2 / is an orthonormal basis for the bi-dimensional space. We denote the coordinates of a point x f f in R by .x1 , x2 / and the time by t . The fluid occupies both half-spaces 1 and 2 . The bone layer occupies the unbounded domain b . The two plane interfaces between the bone layer b and the f f fluid domains 1 and 2 are respectively denoted by 1bf and 2bf , as shown in Figure 1. The thickf ness of the bone plate is denoted by h. The acoustical source is located in the upper fluid domain 1 . f
f
2.1. Equations in the fluid domains 1 and 2
The fluid occupying the domain f1 is modeled by an acoustic fluid whose mass density and bulk modulus at equilibrium are denoted by 1 and K1 , respectively. By using the linear acoustic theory and by neglecting the body forces (other than inertia), the wave equation for the fluid in the domain f 1 may be expressed only in terms of the pressure field p1 .x, t / as follows: 1 P pR1 r 2 p1 D Q, c12
8x 2 f1 ,
(1) f
where Q.x, t / is p the acoustic source density; the constant c1 is the wave celerity in 1 which is defined by: c1 D K1 =1 . Here, we impose an impulsive point source acting at point x s D .0, x2s / f 1 f) 1
Fluid 1 ( e2 e1
source (xs1, xs2)
receivers bf 1
nb1
O Solid ( Fluid 2 (
b) f) 2
bf 2
nb2
f 2
Figure 1. The trilayer model used for representing ultrasound axial transmission test. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
which is defined by QP D 1 F .t / ı.x1 /ı.x2 x2s /,
(2)
where F .t / is a given real scalar function depending only on time t and ı../ is Dirac’s delta function. f f Similarly to the description in the domain 1 , the pressure field p2 .x, t / in the fluid domain 2 reads as follows: 1 pR2 r 2 p2 D 0, c22
f
8 x 2 2 ,
(3)
p f where c2 D K2 =2 is the wave celerity in 2 in which the bulk modulus and mass density of the fluid at equilibrium are denoted by K2 and 2 , respectively. In the domain fj (j D 1, 2), the velocity vector vj .x, t / is relied to the gradients of pressure fields pj .x, t / following the Euler equation: f
j vP j C rpj D 0
8 x 2 j .
(4)
2.2. Equations in the anisotropic poroelastic layer (b ) The anisotropic poroelasticity theory is used to describe the behavior of bone layer. The poroelastic equations employed here are based upon Biot’s original works [37–39] as well as the recent developments in anisotropic constitutive equations [40–42]. The bone material is assumed to consist of a solid skeleton (with mass density s ) and a connected pore network saturated by fluid (with mass density f ). Neglecting the body forces (other than inertia), the equations describing the wave propagation within the bone layer read: r D uR s C f w R, s
1
rp D f uR C k
(5) w P Cbw R,
(6)
where D f C .1 / s is the mixture density, is the porosity, .x, t / is the total stress tensor and p.x, t / is the interstitial fluid pressure in the pores; the vectors of displacement of the solid skeleton and of the fluid are denoted by us .x, t / and uf .x, t /, respectively; the vector of relative displacement between the fluid and the solid frame weighted by the porosity is denoted by w D .uf us /; the tensor k is the anisotropic permeability tensor and the tensor b is defined as b D .f =/ a with a is the tortuosity tensor. Note that both tensors k and b are symmetric second-order tensors. The constitutive equations for the anisotropic linear poroelastic material are given by D C W ˛p,
(7)
1 p D r wC˛W, M
(8)
where C.x/ is the elasticity fourth-order tensor of drained porous material; ˛, which is a symmetric second-order tensor, is the Biot effective tensor; M is the Biot scalar modulus; .x, t / is the infinitesimal strain tensor, which is defined as the symmetric part of rus . For the subsequent development, it is more convenient to represent the following equations in vector form. Let us then use the Voigt’s notation, which expresses the symmetric second-order tensors as vectors, so the stress is denoted s D f11 , 22 , 12 gT , the strain by e D f11 , 22 , 212 gT and the Biot’s effective tensor by ˛L D f˛11 , ˛22 , ˛12 gT , where the superscript ?T designates the transpose operator. We may also introduce an operator L defined as follows: 2 3 2 3 1 0 0 0 (9) L D L1 @1 C L2 @2 , where L1 D 4 0 0 5 , L2 D 4 0 1 5 , 0 1 1 0 where @1 and @2 denote the partial differentiation operators with respect to x1 and x2 , respectively. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
Using these notations, the balance equations of linear momentum (5) and (6) may be rewritten as: R LT s D 0, uR s C f w s
1
f uR C k
(10) T
w P C bw R C L mp D 0,
(11)
where m D f1, 1, 0gT . The constitutive equations for poroelastic material (7) and (8) read: s D C e ˛L p, p D M mT Lw C ˛L T Lus , where C.x/ is the drained elastic tensor using Voigt’s notation: 3 2 c11 c12 c16 C D 4 c21 c22 c26 5 . c61 c62 c66
(12) (13)
(14)
By noting that e D L u and by substituting Equation (13) into Equation (12), one has s D Cu L us C C˛ L w, mp D CM Lw C CT˛ us ,
(15) (16)
where the quantities Cu , C˛ and CM are defined as: Cu D C C M ˛L ˛L T ,
(17)
T
L , C˛ D M ˛m
(18)
T
(19)
CM D M mm .
The tensor Cu is called the undrained elasticity tensor, which may be considered as the rigidity of an equivalent elastic medium in which the relative movement between the interstitial fluid and solid skeleton is null (i.e. when w D 0). Moreover, one may notice that while Cu and CM are symmetric, C˛ is not. In fact, C˛ is symmetric only if the considered poroelastic medium is isotropic (i.e. when ˛11 D ˛22 and ˛12 D 0). 2.3. Boundary conditions At both interfaces 1bf and 2bf , the continuity of pressure and stress fields between the porous medium and the fluid domains requires the following: p D pj , bf nbf j D pj nj ,
8x 2 jbf .j D 1, 2/,
(20)
bf b where nbf j is the normal unit vector to j pointing from the porous domain toward the fluid f domain j (see Figure 1). In addition, open-pore condition at the interfaces jbf (j D 1, 2) is assumed, requiring
P s / nbf w P nbf j D .v j u j ,
8x 2 jbf .j D 1, 2/.
In view of the Euler Equation (4), the interface condition (21) may be rewritten as follows: 1 s nbf rpj C w R C uR 8x 2 jbf .j D 1, 2/. j D 0, j Copyright © 2012 John Wiley & Sons, Ltd.
(21)
(22)
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
Note that for the present configuration (see Figure 1), the normal unit vectors of two interfaces 1bf T bf and 2bf are defined by nbf 1 D n2 D f0, 1g . As a consequence, the boundary conditions may be expressed as follows: 9 1 > @2 pj D .uR 2 C wR 2 / > = j (23) , 8x 2 jbf , .j D 1, 2/, p D pj > > ; t D f 0, pj gT pj ! 0,
f1
8x ! j
,
.j D 1, 2/.
(24)
where us D fu1 , u2 gT and w D fw1 , w2 gT ; the traction vectors t at the interfaces are defined by t D LT2 s D f21 , 22 gT . 3. EQUATIONS IN THE LAPLACE–FOURIER (s k1 ) DOMAIN The boundary value problem given by Equations (1), (3), (10), (11), (23) and (24) deals with solving a system of linear partial differential equations of which the coefficients are constant in the longitudinal direction defined by e1 -axis. Here, we propose to solve the system as follows: (i) the system of equations is firstly transformed into wavenumber-frequency domain by using a Fourier transform with respect to x1 combining to a Laplace transform with respect to t ; (ii) in the wavenumberfrequency domain, the equations in both fluid domains are analytically solved giving impedance boundary conditions for the solid layer that may be solved by the finite element method; (iii) the space–time solution is finally obtained by performing two inverse transforms. We recall that the general form of a Laplace–Fourier transform applied to a real-valued function y.x1 , x2 , t / is defined as: Z 1 Z C1 y.x1 , x2 , t /e ik1 x1 dx1 e st dt , (25) y.k Q 1 , x2 , s/ WD 0
1
p where i is the imaginary unit defined by i D 1; s 2 is the complex Laplace variable; k1 2 is the real Fourier variable representing the wavenumber on the e1 -axis; and denote the set of all the real and complex numbers, respectively. In (s k1 ) domain, the time derivative and the spatial differential operation with respect to x1 can be replaced by @.?/ ! s.?/ and @1 .?/ ! i k1 .?/, respectively. @t f
f
3.1. Transformed problem in (s k1 ) domain for the fluids 1 and 2
By applying the Laplace–Fourier transform (25) to Equations (1), (23) and (24), we obtain a boundary value problem whose differential equation is defined for pQ1 with respect only to x2 : 2 s 2 C k1 pQ1 @2 pQ1 D 1 FQ0 .s/ı.x2 x2s /, 8x2 > 0, (26) c12 @2 pQ1 D 1 s 2 .uQ 2 C wQ 2 /, pQ2 ! 0,
at x2 D 0, for x2 ! C1.
(27) (28)
The solution of this system may be presented in a semi-explicit form as follows: s 2 GQ 1 FQ0 ˇ1 .x s x2 / s 1 21 ˇ1 x2 2 e C e ˇ1 .x2 Cx2 / C e , 2ˇ1 ˇ1 s 2 GQ 1 FQ0 ˇ1 .x s x2 / s 1 21 ˇ1 x2 2 pQ1 D C e ˇ1 .x2 Cx2 / C e , e 2ˇ1 ˇ1 r 2 where ˇ1 D cs 2 C k12 and GQ 21 D uQ 2 .k1 , 0, s/ C wQ 2 .k1 , 0, s/. pQ1 D
8 0 x2 x2s
(29)
8 x2 x2s
(30)
1
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
Similarly, the solution of pQ2 in (s k1 ) domain may be expressed as follows: pQ2 D r where ˇ2 D
s2 c22
2 s 2 GQ 22 ˇ2 .x2 Ch/ e , ˇ2
8x2 h,
(31)
C k12 and GQ 22 D uQ 2 .k1 , h, s/ C wQ 2 .k1 , h, s/.
At two interfaces x2 D 0 and x2 D h, the solutions of pQ1 and pQ2 are given by 1 FQ0 ˇ1 x s 1 s 2 GQ 21 2 C e , ˇ1 ˇ1 2 s 2 GQ 22 pQ2 .k1 , h, s/ D . ˇ2
pQ1 .k1 , 0, s/ D
(32) (33)
3.2. Transformed problem in (s k1 ) domain for the solid b Q D i k1 L1 C @2 L2 , the By noting that in the Fourier–Laplace domain, the operator L becomes L boundary valued problem defined by Equations (10) and (11) may be transformed into the (s k1 ) domain as follows: s 2 uQ s C f s 2 w Q i k1 LT1 sQ @2 LT2 sQ D 0 , 8x2 2 h, 0Œ, (34) Q C i k1 LT1 mpQ C LT2 @2 .mp/ Q D0 f s 2 uQ s C s k1 C s 2 b w or
s 2 1 f s 2 1 f s 2 1 s k1 C s 2 b
s LT1 sQ LT2 sQ uQ 0 @ D , i k1 2 0 w Q LT1 mpQ LT2 mpQ
(35)
where 1 is the 2 2 identity matrix and (see Equations (15) and (16)) sQ D i k1 Cu L1 uQ s C i k1 C˛ L1 w Q C Cu L2 @2 uQ s C C˛ L2 @2 w, Q mpQ D
i k1 CT˛ L1 uQ s
C i k1 CM L1 w Q C CT˛ L2 @2 uQ s
C CM L2 @2 w. Q
(36) (37)
The vectors concerning sQ and .mp/ Q in the system (35) may be evaluated in terms of uQ s and w Q as follows: 11
s 12
s C11 C12 Cu Cu uQ uQ LT1 sQ ˛ ˛ D i k1 C , (38) T 11 21 T 12 @2 w w Q Q LT1 mpQ C C C11 C ˛ ˛ M M 21
s 22
s C21 C22 Cu Cu uQ uQ LT2 sQ ˛ ˛ @2 D i k C , (39) T 1 21 22 T 22 w Q w Q LT2 mpQ C C C12 C ˛ ˛ M M ab ab where the 2 2 matrices Cab u , C˛ and CM with a, b D 1, 2 are defined by: T Cab u D La Cu Lb ,
T Cab ˛ D La C˛ Lb ,
T Cab M D La CM Lb
with a, b D 1, 2.
(40)
By replacing (38) and (39) into Equation (35) and then making some arrangements, we obtain a vectorial form of this problem: s 2 A1 Ú C sA2 Ú C k12 A3 Ú i k1 A4 @2 Ú @2 Ø D O,
(41)
where s Q Ú D uw , Q Copyright © 2012 John Wiley & Sons, Ltd.
LT2 sQ D i k1 AT4 Ú C A5 @2 Ú, ØD LT2 mpQ
(42)
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
and
0 0 1 f 1 A1 D , , A2 D f 1 b 0 k1
11 C11 Cu C12 ˛ u T , A D A3 D 11 4 21 T C11 C˛ C ˛ M
(43) 12
C22 A5 D 22u T C˛
C˛ , C12 M
22
C˛ . C22 M
(44)
The 4 4 matrices A1 , A2 , A3 , A4 and A5 depend only on the material properties of the elastic porous material. A physical interpretation of these matrices can be given: A1 relates to the mass density; A2 relates to the viscous property via the permeability; A3 , A4 and A5 relate to the poroelastic properties. It is worth noting that whereasA1 , A2 , A3 and A5 are symmetric, the matrix A4 is not because we have seen that the matrix C˛ defined in Equation (18) is nonsymmetric when the poroelastic medium is anisotropic. At two interfaces x2 D 0 and x2 D h, the continuity condition (23) between fluid pressures and stress traction require that LT2 sQ PQ 0 Q where P D (45) D Q pQ LT2 mpQ x D0,h P 2
where pQ D pQ1 at x2 D 0 and pQ D pQ2 at x2 D h. In view of Equations (32) and (33), the boundary conditions of the system (41) may be expressed by:
Ø.0/ D P1 Ú.0/ C F0 , Ø.h/ D P2 Ú.h/,
(46)
where 2
P1 D
1 s P, ˇ1
n F0 D 0,
1 FQ0 ˇ1
2
P2 D
s
2 s P, ˇ2
e ˇ1 x2 , 0,
1 FQ0 ˇ1
2 0 60 PD4 0 0 oT s
e ˇ1 x2
0 1 0 1
0 0 0 0
3 0 17 , 05 1
.
(47)
(48)
4. FINITE ELEMENT FORMULATION IN SPECTRAL DOMAIN AND TIME–SPACE SOLUTION 4.1. Weak formulation The weak formulation of the boundary value problem given by (41) and (46) may be now carried out using a usual procedure. Let C ad be the admissible function space of the problem constituted by the sufficiently differentiable complex-valued functions such as: x2 2 H b D h, 0Œ ! ı Ú.x2 / 2 . The conjugate transpose of ı Ú is denoted by ı Ú . Upon integrating (41) against a test vector function ı Ú and integrating by part, then applying the boundary condition (46), the weak formulation of the differential equation (41) reads: For all k1 fixed in and for all s fixed in , find Ú.k1 , x2 , s/ 2 C ad such that:
s 2 h ı Ú, A1 Ú iH b C s h ı Ú, A2 Ú iH b C k12 h ı Ú, A3 Ú iH b ˚˝ ˛ C i k1 @2 ı Ú, AT4 Ú H b h ı Ú, A4 @2 Ú iH b C h @2 ı Ú, A5 @2 Ú iH b
(49)
C ı Ú .0/ P1 Ú.0/ C ı Ú .h/ P2 Ú.h/ D ı Ú .0/ F0 ,
8 ı Ú 2 C ad ,
where h., .iH b denotes the inner product over H b , for instance, h ı Ú, A1 Ú iH b R0 h ı Ú.x2 / A1 Ú.x2 / dx2 . Copyright © 2012 John Wiley & Sons, Ltd.
D
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
4.2. Finite element formulation We proceed by introducing a finite element mesh of the domain Œh, 0, which contains nel eleS ments: Œh, 0 D e e with e D 1, : : : , nel . By using the Galerkin finite element method, both functions Ú and ı Ú in each element e are approximated using the same shape function:
Ú.x2 / D Ne Ve , ı Ú.x2 / D Ne ıVe , 8 x2 2 e ,
(50)
where Ne is the shape function, Ve and ıVe are the vectors of nodal solutions of Ú and ı Ú within the element e , respectively. Replacing Equation (50) into Equation (49) and assembling the elementary matrices lead to a static-like linear system of equations for each couple values (s, k1 ): (51) Kb C K V D F, where V is the global nodal solution vector; Kb is the global ‘stiffness matrix’ of the poroelastic layer; K represents the coupled operator between the fluid and poroelastic layers; the vector F is the external force vector due to the acoustical source. Noting the number of nodes by N , the size of vectors V and F is neq D 4N because each node has four degrees of freedom. The sizes of Kb and K are neq neq . For all couples .s, k1 / fixed in , these quantities may be expressed by:
Kjk D
8 ˆ ˆ < ˆ ˆ : (
Fj D
1 s 2 ˇ1
if .j , k/ D .2, 2/, .2, 4/, .4, 2/, .4, 4/
2 s 2 ˇ2
if .j , k/ D .neq 2, neq 2/, .neq 2, neq /, .neq , neq 2/, .neq , neq /
0
otherwise
1 FQ0 ˇ1
ˇ1 x2s
e
0
if j D 2, 4
(53)
otherwise
Kb D s 2 K1 C s K2 C k12 K3 C i k1 K4 C K5
(54)
where the matrices K1 , K2 , K3 , K4 and K5 are independent to (s, k1 ) and are defined by: [Z [Z K1 D NTe A1 Ne dx2 I K2 D NTe A2 Ne dx2 , e
K3 D
[Z e
K5 D
(52)
[Z e
e
e
e
NTe A3 Ne
dx2 I
K4 D
[Z e
(55)
e
n
.@2 Ne /T A4 Ne
e
.@2 Ne /T A5 .@2 Ne / dx2 ,
o a
dx2 ,
(56) (57)
e
in which the notation f?ga designates the anti-symmetric part of the f?g. 4.3. Computation of time–space solution For fixed values of .s, k1 / in the Laplace–Fourier transformed domain, the solution of Ú may be computed by solving the system of linear equations (51) in the complex domain. The solutions for pQ1 and pQ2 in two fluid domains may then be determined by using the Equations (29), (30) and (31), respectively. To obtain the spatio–temporal solution, we need to perform a numerical inverse Laplace–Fourier transform. In this paper, the inverse Fourier transform is computed by using the usual fast Fourier transform technique. The inverse Laplace transform is carried out using the Quadrature convolution method, which has proved to be a very efficient technique for computing the time response solution in many dynamic problems [43]. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
5. NUMERICAL TESTS This section presents some numerical tests describing an in vivo ultrasound test on human cortical long bones. The acoustic source (2) located at the position .x1s , x2s / D .0, 2/ .in millimeter/ in the upper fluid domain f1 has the time-history function given by: 2
F .t / D F0 e 4.fc t1/ sin.2 fc t /,
(58)
where F0 = 1 m.s2 and fc D 1 MHz which is the central frequency. Using a Fourier transform shows that this signal has a spectrum, which contains a frequency bandwidth about 0 2.5 MHz (the reader could refer to Figure 2 in [16]). For all tests presented later, the total duration of the simulation will be taken by T D 2 105 s. 5.1. Physical parameters f
f
Both fluid domains 1 and 2 are assumed to be identical and are considered as an idealized water. The mechanical properties of the fluids are given by 1 D 2 D 1000 kg.m3 and K1 D K2 D 2.25 GPa. Thus the wave velocities in two fluid domains are: c1 D c2 D 1500 m.s1 . To describe the behavior of the poroelastic bone plate, the drained elasticity tensor C as well as Biot’s effective coefficients ˛ and M used in (7) and (8) should be provided. For this study, these parameters are derived from the characteristics of the interstitial fluid and solid skeleton phases by using a continuum micromechanics model proposed by Hellmich et al. [44]. According to this model, the micropores at the microstructural scale are regarded as cylindrical pores with a circular cross section. In drained condition, the constitutive behavior of the material inside the pores does not possess stiffness. Hence, the estimated drained microstructural stiffness of the bone whose solid bone matrix’s elasticity tensor is C m reads: h i1 1 C D .1 /C m W .1 /I C I P cyl W C m
(59)
where I denotes the fourth-order identity tensor; the fourth-order tensor P cyl is the Hill’s tensor for materials with periodical cylindrical inclusions, which may be derived in closed form [44]. The tensor ˛ and the constant M can be then evaluated by [40, 42]: ˛ D I C W .Sm W I/ ,
1 D C ˛ W S W ˛, M
(60)
where S D .C/1 and Sm D .C m /1 are respectively the drained and solid material compliance tensors, I designates the second-order tensor identity and the scalar C denotes the effective compressibility of porous matrix material, which is given by 1 1 CD C K Km
1 1 m K Kf
,
(61)
where K D .I W S W I/1 , K f and K m D .I W Sm W I/1 are respectively the bulk moduli of the drained porous matrix, of the interstitial fluid and of the poroelastic material. For the numerical tests presented in this paper, the interstitial fluid is also assumed to be the same as the exterior fluid domain, that is f D 1000 kg.m3 and Kf D 2.25 GPa. The solid bone has the mass density s D 1722 kg.m3 and is assumed to be a transversely isotropic elastic material m [44]. The components of the elastic tensor are given by using the Voigt’s notations: c11 D 28.7 GPa, m m m m m c22 D 23.6 GPa, c12 D 9.9 GPa, c66 D 7.25 GPa, c16 D c26 D 0. The viscosity of interstitial fluid is D 1 103 Pa.s and the permeability tensor is roughly taken by 11 D 22 D 5 1013 m2 [45]. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
5.2. Validation This section will present some calculations to validate the proposed method. First, the hybrid spectral/FE solutions will be compared with the ones obtained by using a conventional finite element procedure. Then, the behavior of the solution vis-à-vis the parameters of numerical discretization will be studied. We consider a bone layer, which has a constant thickness h D 4 mm. The numerical parameters used for simulation have been chosen in a similar way to which required for a classical procedure of dynamic finite element analysis [24]. The space interval x is about the 1=10 shortest wavelength. The time step is chosen for satisfying the Courant–Friedrichs–Lewy condition. To obtain a boundedvalue solution in space on the e 1 -axis, the domain size should be chosen sufficiently large to avoid the fact that fastest waves pass over the boundaries after the duration T . Here, the space interval in e 1 -direction has been chosen to be x D 0.75 104 m. The time step used is t D 107 s. To carry out the finite element analysis on vertical direction e 2 , the bone plate’s thickness is discretized into 12 three-noded quadratic isoparametric Lagrangian elements. Figure 2 presents snapshots of the wave field in the coupled trilayer system at different instants t D 1 s, t D 3 s and t D 6 s. At each point x in the domain, the quantity presented for obtaining f f these snapshots is log.jp1 .x, t /j/ if x 2 1 , log.jp2 .x, t /j/ if x 2 2 and log.jp.x, t //j if x 2 b . In these graphs, we can clearly observe the transmission and reflection of ultrasonic waves through the interfaces. As the longitudinal waves propagating in the poroelastic layer are much faster than acoustic waves in the fluid, the wavefronts of two symmetrical head waves, which have truncated cone forms, are very well displayed in the last two graphs in Figure 2. Note that using the proposed method allows us to take into account the exact radiation conditions of two fluid domains, therefore the results involve only the waves diffracted by the bone layer. Figure 3 presents the numerical solution of p1 measured at two positions R1 .x1 , x2 / D .2, 2/ (mm) (upper) and R2 .x1 , x2 / D .20, 2/ (mm) (lower) by using the proposed spectral/finite element procedure. We also present numerical solutions obtained by using conventional time domain finite element analysis [24]. The comparison shows that S/FEM and FEM results perfectly match each other. While about 5 h are needed on a common desktop personal computer (using COMSOL
Figure 2. Snapshots of the fluid pressure (p1 , p and p2 ) in the trilayer. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
p1(x1,t) (Pa)
x1 = 0.002 m SFEM FEM
50 0 −50 0
0.5
1
1.5
x 10−5
time (s)
p1(x1,t) (Pa)
x1 = 0.02 m SFEM FEM
20 0 −20 0.6
0.8
1
1.2
1.4
1.6
1.8
2
x 10−5
time (s)
Figure 3. Comparison between spectral finite element and finite element solutions of p1 at x1 D 2 mm (upper graph) and at x1 D 20 mm (lower graph) from the source.
multiphysics software (COMSOL, Inc., Stockholm, Sweden) [46]) for simulation using the dynamic finite element method, the proposed spectral/finite element method needs only about 400 s on the same personal computer. Note that both methods use the same time step for the calculation. The computational advantage of analysis in spectral domain over the analysis using direct time FEM is mainly due to the fact that only very small matrices need to be factorized. Because the exact solution is not available for this example, we do not calculate the numerical p1 errors for convergence analysis, but investigate the behavior of a function L2,T .x/, which is defined as follows, with respect to the parameters of numerical dicretization (Nt , N1 , nel ): v u Nt u X p1 L2,T .x/ D t t .p1 .x, n t //2 , (62) nD0 1 Figure 4 presents a study with respect to different parameters by calculating Lp2,T at two points R1 and R2 . Figure 4(A) shows the effect of sampling number N1 , which is associated with the wavenumber k1 , by fixing the parameters Nt D 2048 and nel D 12. We recall that Nt is related to the total number of time step and nel is related to the element number. This last parameter combined with the parameter N1 is related to the discretization in space. Next, in Figure 4(B), the parameters N1 D 2048 and nel D 12 are fixed, and the effect on the convergence of the total number of
(A) Nt = 2048; nel = 12
(B) N1 = 2048; nel = 12 0.06
0.06 0.04 0.02
8
10
12
log2(N1)
0.04 0.02 0 8
x1 = 0.002m
x1 = 0.02m
x1 = 0.02m
0.04
Lp2,1T(x1)
x1 = 0.002m
x1 = 0.02m
Lp2,1T(x1)
p
L2,1T(x1)
x1 = 0.002m
0
(C) N1 = 2048; Nt = 2048 0.05
0.03 0.02
10
log2(Nt)
12
0.01
5
10
15
20
number of elements
1 Figure 4. Lp 2,T versus (A) number of space interval N1 , (B) number of time steps Nt and (C) the number of elements.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
time step Nt is shown. Note that the total number of time step Nt is related to time step size t . 1 Last, Figure 4(C) presents the values of Lp2,T versus the number of finite elements nel by fixing N1 D 2048 and Nt D 2048. These studies show that the proposed procedure using spectral finite elements may have fast convergence when refining parameters of numerical discretization. Moreover, it can be deduced that, for this numerical example and for the observations point located at R1 and R2 , convergence is reached for Nt D 2048, N1 D 2048 and nel D 12 which correspond well to the parameters used for obtaining results presented in Figure 3. 5.3. Influence of porosity on first arriving signal velocity When the axial transmission technique is used, the earliest event or wavelet usually called first arriving signal (FAS) of the multicomponent signal recorded at the receivers has been the most often investigated. In fact, the FAS velocity have been shown to be a relevant index of bone status [47–49]. In this section, we present some numerical results to study the influence of the bone’s porosity on the FAS velocity denoted by VF . For this purpose, the material properties of bone matrix (m , C m ) were fixed. The poroelastic properties were then determined by using the micromechanics analysis presented in Section 5.1. To determine VF , the p1 -signal is captured at an array of 14 sensors, which has a gap 0.8 mm from each other. The distance from the emitted source to the first sensor is 11 mm. At each sensor, the first zero-crossing location on time axis of measured p1 -signal may
sensor number #
15
10
5
0 0
0.5
1
1.5
2
x 10−5
time (s)
Figure 5. Case of a bone layer with thickness h D 4 mm and D 5%: seismograph of p1 .t / measured at 14 sensors. The gray line (red line if the color is available) is the linear regression of zero-crossing positions. bone thickness h = 0.6mm 3900
4000
3850
FAS (ms−1)
FAS (ms−1)
bone thickness h = 4mm 4050
3950 3900 3850
3800 3750 3700 poroelastic model undrained elastic model
poroelastic model undrained elastic model 3800
0
5
Porosity
10
(%)
15
3650 0
5
Porosity
10
15
(%)
Figure 6. Comparison of first arriving signal velocities for different bone porosities in two cases of bone’s thickness: h D 4 mm (left) and h D 0.8 mm (right). Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
be determined. Then VF is evaluated by the slope of the linear regression of these location over 14 sensors. An example of the results obtained at 14 sensors is shown in Figure 5. Figure 6 (left) shows the variation of the FAS velocity with respect to when considering a bone plate with thickness h D 4 mm. One may observe that the VF versus relation is practically linear and VF decreases when the porosity is higher. We also present in this graph the VF versus relation obtained by using equivalent elastic media. For a given porosity , the corresponding elastic models use the mass density taken as the mixture density given by (5), and the elasticity tensor is taken as the undrained elasticity tensor Cu defined by (17). For lower porosities, the VF obtained by using the poroelastic model are slightly different from the ones obtained using the elastic model. The difference becomes more significant for higher porosity because the fluid–solid movement is more important and would not be negligible. A similar discussion may be made when considering a thin bone layer (h D 0.8 mm) in Figure 6 (right). 6. CONCLUSION By using a hybrid semi-analytical/finite element technique, we have derived a finite element formulation in spectral domain for anisotropic poroelastic plate saturated by an interstitial fluid and immersed in a fluid. The procedure is based on a Laplace transform for time and a Fourier transform for the spatial variable in longitudinal direction defined by x1 . Performing these transforms leads to a one-dimensional differential equation with respect only to x2 , which may be solved by using conventional Galerkin finite element method. Note that in the Laplace–Fourier domain, the radiation conditions presenting two half-spaces of fluid may exactly be introduced. The validation has been carried out. Because only one-dimensional finite element problem need to be solved in spectral domain, the proposed procedure requires a very small memory size for computation. In addition, the computation time is very low, especially for the configuration used for ultrasonic test presented here, wherein we only need to compute the time-domain solution at some specific points. It has also been shown that the use of conjugated quadrature method for the inverse Laplace transform has very good stability when computing the time-domain solutions. Using the proposed method, numerical studies of ultrasonic wave propagation in poroelastic bone layer may be easily carried out, eventually in the very high frequency domain. This allows us to consider the anisotropy and the heterogeneity of the bone matrix material as well as bone porosity. Some preliminary results presented in this paper show that the FAS velocity is strongly influenced by the bone porosity. It seems that using of equivalent elastic model might lead to a poor estimation of FAS velocity in cortical bone with higher bone porosity. The detailed investigation of behavior of (heterogeneous) poroelastic cortical bone under broadband frequency ultrasound will be presented in a forthcoming paper. REFERENCES 1. Lowet G, Van der Perre G. Ultrasound velocity measurements in long bones: measurement method and simulation of ultrasound wave propagation. Journal of Biomechanics 1996; 29:1255–1262. 2. Nicholson P, Moilanen P, Karkkainen T, Timonen J, Cheng S. Guided ultrasonic waves in long bones: modelling, experiment and in vivo application. Physiological Measurement 2002; 23(4):755–768. 3. Moilanen P. Ultrasonic guided waves in bone. IEEE transactions of ultrasonics, ferroelectrics, and frequency control 2008; 55(6):1277–1286. 4. Bousson V, Meunier A, Bergot C, Vicaut E, Rocha M, Morais M, Laval-Jeantet A, Laredo J. Distribution of intracortical porosity in human midfemoral cortex by age and gender. Journal of Anatomy 2001; 16:1308–1317. 5. Thomas C, Feik SA, Clement J. Regional variation of intracortical porosity in the midshaft of the human femur: age and sex differences. Journal of Anatomy 2005; 206:115–125. 6. Sansalone V, Naili S, Bousson V, Bergot C, Peyrin F, Zarka J, Laredo JD, Haiat G. Determination of the heterogeneous anisotropic elastic properties of human femoral bone: from nanoscopic to organ scale. Journal of Biomechanics 2010; 43:1857–1863. 7. Dong NX, Guo EX. The dependence of transverse isotropic elasticity of human femoral cortical bone on porosity. Journal of Biomechanics 2004; 37(8):1281–1287. 8. Baron C, Talmant M, Laugier P. Effect of porosity on effective diagonal stiffness coefficients (ci i ) and anisotropy of cortical bone at 1 MHz: a finite difference time domain study. Journal of the Acoustical Society of America 2007; 122:1810–1817. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
ULTRASONIC WAVE PROPAGATION IN ANISOTROPIC POROELASTIC BONE PLATE
9. Bossy E, Talmant M, Laugier P. Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models. Journal of the Acoustical Society of America 2004; 115:2314–2324. 10. Baron C, Naili S. Elastic wave propagation in a fluid-loaded anisotropic waveguide with laterally varying properties. Comptes Rendus Mecanique 2008; 336:772–730. 11. Baron C, Naili S. Propagation of elastic waves in a fluid-loaded anisotropic functionally graded waveguide: application to ultrasound characterization. Journal of the Acoustical Society of America 2010; 127(3):1307–1317. 12. Nayfeh HA. Wave Propagation in Layered Anisotropic Media. Elsevier: North-Holland, 1995. 13. Shuvalov A, Poncelet O, Deschamps M. Analysis of the dispersion spectrum of fluid-loaded anisotropic plates: leaky-wave branches. Journal of Sound and Vibration 2006; 296(3):494–517. 14. Bossy E, Talmant M, Laugier P. Effect of bone cortical thickness on velocity measurements using ultrasonic axial transmission: A 2D simulation study. Journal of the Acoustical Society of America 2002; 112:297–307. 15. Grimal Q, Naili S. A theoretical analysis in the time-domain of wave reflection on a bone plate. Journal of Sound and Vibration 2006; 298:12–29. 16. Desceliers C, Soize C, Grimal Q, Haiat G, Naili S. A time-domain method to solve transient elastic wave propagation in a multilayer medium with a hybrid spectral-finite element space approximation. Wave Motion 2008; 45(4):383–399. 17. Haiat G, Naili S, Grimal Q, Talmant M, Desceliers C, Soize C. Influence of a gradient of material properties on ultrasonic wave propagation in cortical bone: application to axial transmission. Journal of the Acoustical Society of America 2009; 125(6):4043–4052. 18. Cowin SC. Bone poroelasticity. Journal of Biomechanics 1999; 32:217–238. 19. Zhang D, Cowin SC. Oscillatory bending of a poroelastic beam. Journal of the Mechanics and Physics of Solids 1994; 42(10):1575–1599. 20. Nguyen VH, Lemaire T, Naili S. Poroelastic behaviour of cortical bone under harmonic axial loading: a finite element study at the osteonal scale. Medical Engineering & Physics 2010; 32(4):384–390. 21. Nguyen VH, Lemaire T, Naili S. Influence of interstitial bone microcracks on strain-induced fluid flow. Biomechanics and Modeling in Mechanobiology 2011; 10(6):963–972. 22. Williams JL. Ultrasonic wave propagation in cancellous and cortical bone: prediction of some experimental results by Biot’s theory. Journal of the Acoustical Society of America 1992; 91(2):1106–1112. 23. Hosokawa A. Simulation of ultrasound propagation through bovine cancellous bone using finite-difference timedomain methods. Journal of the Acoustical Society of America 2005; 118(3):1782–1789. 24. Nguyen VH, Naili S, Sansalone V. Simulation of ultrasonic wave propagation in anisotropic cancellous bones immersed in fluid. Wave Motion 2010; 47(2):117–129. 25. Nguyen VH, Naili S, Sansalone V. A closed-form solution for in vitro transient ultrasonic wave propagation in cancellous bone. Mechanics Research Communications 2010; 37(4):377–383. 26. Kausel E. Fundamental Solution in Elastodynamics. Cambridge University Press: New York, 2006. 27. Virieux J. P-SV-wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics 1986; 51(4):889–901. 28. Kampanis NA, Dougalis VA, Ekaterinaris JA. Effective Computational Methods for Wave Propagation. Chapman & Hall/CRC: New York, 2008. 29. Protopappas VC, Kourtis IC, Kourtis LC, Malizos KN, Massalas CV, Fotiadis DI. Three-dimensional finite element modeling of guided ultrasound wave propagation in intact and healing long bones. Journal of the Acoustical Society of America 2007; 121(6):3907–3921. 30. Thompson LL, Huan R. Implementation of exact non-reflecting boundary conditions in the finite element method for the time-dependent wave equation. Computer Methods in Applied Mechanics and Engineering 2000; 187(1-2): 137–159. 31. Givoli D. High-order local non-reflecting boundary conditions: a review. Wave Motion 2004; 39(4):319–326. 32. Han X, Liu GR, Lam KY. Transient waves in plates of functionally graded materials. International Journal for Numerical Methods in Engineering 2001; 52(8):851–865. 33. Han X, Liu GR, Xi ZC, Lam KY. Transient waves in a functionally graded cylinder. International Journal of Solids and Structures 2001; 38(17):3021–3037. 34. Liu G, Xi ZC. Elastic Waves in Anisotropic Laminates. CRC Press: Florida, 2002. 35. Gopalakrishnan S, Chakraborty A, Roy Mahapatra D. Spectral Finite Element Method. Springer-Verlag: London, 2008. 36. Marzani A. Time-transient response for ultrasonic guided waves propagating in damped cylinders. International Journal of Solids and Structures 2008; 45(25-26):6347–6368. 37. Biot MA. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics 1955; 26(2):182–185. 38. Biot MA, Willis DG. The elastic coefficients of the theory of consolidation. Journal of Applied Mechanics 1957; 79:594–601. 39. Coussy O. Poromechanics. John Wiley & Sons, Ltd: England, 2004. 40. Cheng AH. Material coefficients of anisotropic poroelasticity. International Journal of Rock Mechanics and Mining Sciences 1997; 34(2):199–205. 41. Cowin SC. A recasting of anisotropic poroelasticity in matrices of tensor components. Transport in Porous Media 2003; 50(1):35–56.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
V.-H. NGUYEN AND S. NAILI
42. Thompson M, Willis JR. A reformation of the equations of anisotropic poroelasticity. Journal of Applied Mechanics 1991; 58:612–616. 43. Schanz M, Antes H. Application of ’Operational Quadrature Methods’ in time domain boundary element methods. Meccanica 1997; 32:179–186. 44. Hellmich C, Ulm FJ. Microporodynamics of bones: prediction of the “Frenkel–Biot” slow compressional wave. Journal of Engineering Mechanics 2005; 131(9):918–927. 45. Cowin SC. Bone Mechanics Handbook (2nd edn). CRC Press: Boca Raton, FL, 2001. 46. COMSOL Multiphysics 3.5. User’s Guide, 2008, www.comsol.com 47. Barkmann R, Kantorovich E, Singal C, Hans D, Genant HK, Heller M, Gluer CC. A new method for quantitative ultrasound measurements at multiple skeletal sites. Journal Of Clinical Densitometry 2000; 3:1–7. 48. Hans D, Srivastav SK, Singal C, Barkmann R, Njeh CF, Kantorovich E, Gluer CC, Genant HK. Does combining the results from multiple bone sites measured by a new quantitative ultrasound device improve discrimination of hip fracture? Journal of Bone and Mineral Research 1999; 14(4):644–651. 49. Stegman MR, Heaney RP, Travers-Gustafson D, Leist J. Cortical ultrasound velocity as an indicator of bone status. Osteoporosis International 1995; 5(5):349–533.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
