Nội dung

Vietnam Journal of Mechanics, VAST, Vol. 33, No. 4 (2011), pp. 225 – 243 SIMULATION OF TRANSIENT ULTRASONIC WAVE PROPAGATION IN FLUID-LOADED HETEROGENEOUS CORTICAL BONE Vu Hieu Nguyen, Salah Naili Université Paris-Est, France Abstract. This work deals with the ultrasonic wave propagation in the cortical layer of long bones which is known as being a functionally-graded anisotropic material coupled with fluids. The motivation arises from mechanical modeling of the ultrasound axial transmission technique in vivo for cortical long bone which is known as being a functionally-graded anisotropic material. The proposed method is based on a combined Laplace-Fourier transform which substitutes a problem defined by partial differential equations into a system of differential equations established in the frequency-wavenumber domain. In the spectral domain, as radiation conditions may be exactly introduced in the infinite fluid halfspaces, only the heterogeneous solid layer needs to be analyzed using finite element method. Several numerical tests are presented showing very good performance of the proposed approach. Keywords: Spectral finite element, transient wave, ultrasound, anisotropic, vibroacoustic, cortical bone, axial transmission. 1. INTRODUCTION In recent years, quantitative ultrasound (QUS) has demonstrated its promising potential in assessment of in vivo bone characteristics. An advantage of QUS over X-ray techniques is its ability to give some information about the elastic properties and defects of bones. Moreover, ultrasound is non-ionizing and the ultrasonic apparatus is relatively inexpensive and can be made portable. For measuring in vivo 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 on 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 can allow to the quantification of the geometrical information as well as mechanical characteristics of the cortical bone at the measured skeletal site [2, 3]. Mechanical modeling of this experiment deals with considering a model describing vibro-acoustic interactions of a solid waveguide (which represents the cortical bone) coupled with two fluid media (which represents soft tissues such as skin or marrow). The cortical bone may be described as plate-like or cylindrical-like structures. This technique of nondestructive testing used to evaluate the material properties requires careful analysis of 226 Vu Hieu Nguyen, Salah Naili the reflections, conversion modes and interferences of longitudinal and shear waves within the bone structure. Many studies have focused on the modeling of guided waves in long bones by using fluid-loaded (multilayer) plate models. The analysis of wave phenomena in the structure have been considered in the frequency-domain [4, 5, 6, 7] or in the time-domain [8, 9, 10, 11]. For the analysis of multilayer structures in the time-domain, there are mainly two approaches. The first one involves using (semi-)analytical methods such as the generalized ray/Cagniard-de-Hoop technique [9] or the direct stiffness matrix method [12]. The second one involves using numerical methods such as the finite difference method (FDM) [8, 13, 14] or the finite element method (FEM) [11, 13, 15]. 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 when considering unbounded domains [13, 16, 17]. For analyzing the waveguide with geometrical and mechanical properties which are constant along one or two directions, an efficient approach called the Hybrid Numerical Method (HNM, see e.g. [18, 19, 20]), alternatively called Spectral Finite Element Method (SFEM, see e.g. [21, 22]). The key point of this method consists of 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 problem into the frequency-wavenumber domain. Then, the wave equations in the spectral domain governed in a 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 finite element method [10, 18, 19, 22]. In comparison with classical FEM which requires extensive discretization due to small wavelengths for high frequency problems, the SFEM presents important advantages in time and memory saving. However, noting that when using the SFEM, finite element analysis should be performed at each value in the set of frequencies/wavenumbers, the computational cost increases rapidly with the size of the cross-sectional domain. In addition, when performing the inverse Fourier transform to obtain the time response of undamped structures, some numerical difficulties may arise because of singular-valued solutions at the poles of the frequency response solution. In order to improve the numerical stability, an appropriate artificial damping would be added into the model, removing the singularity or almost singular behavior near the real axis of the complex frequency response function [21]. In this work, we propose to use the spectral finite element method to study time response solve the transient wave propagation problem in a layered or functionally graded elastic plate loaded by two infinite fluid domains. This geometrical configuration is appropriate to study the radio-frequency signals obtained from the AT technique in order to recover the geometrical and mechanical characteristics of cortical bone. By assuming that the fluids are homogeneous, the general solutions of acoustic waves in fluids can be solved analytically in the spectral domain, so only the solid domain needs to be analyzed by finite element method. Simulation of transient ultrasonic wave propagation in fluid-loaded heterogeneous cortical bone 227 The paper is organized as follows. In Section 2, we present a description of the geometrical configuration of the coupled fluid-solid system. Then, the governing equations as well as the boundary and interface conditions of the considered problem are given. Section 3 presents the problem formulation in the Laplace-Fourier domain (i.e. frequencywavenumber domain). Section 4 provides the spectral finite element formulation and the procedure to obtain the solution in the spatio-temporal domain. Section 5 presents some numerical results. It begins with several validations which are performed in comparison with finite element results. Then, we show a case in which the influence of the variation of mechanical properties in the solid layer is considered. Finally, the conclusion and perspectives of this work are given in Section 6. 2. STATEMENT OF THE PROBLEM In this section, we describe an idealized model of the AT technique performed on cortical bone. The modeling consists of a plate of infinite extent (representing the cortical bone layer) and two semi-infinite ideal fluid media (representing soft tissues) (see Fig. 1). A source of pressure and a array of the receivers are placed in the upper fluid medium. The excitation represents an infinite line source perpendicular to the considered plane. As a consequence, a two-dimensional (2D) model is appropriate to modeling wave phenomena in the system. 2.1. Description of the problem Let R(O; e1 , e2 ) be the reference Cartesian frame where O is the origin of the space and (e1 , e2 ) is an orthonormal basis for this space. The coordinates of a point x in R are specified by (x1 , x2 ) and the time is denoted by t. Consider an infinite solid layer occupying the domain Ωb that represents the cortical bone plate with a constant thickness h (Ωb = {x(x1 , x2 ); 0 ≥ x2 ≤ −h}). This bone plate is loaded on its upper and lower surfaces by two fluid halfspaces. The upper fluid domain is denoted by Ωf1 (Ωf1 = {x(x1 , x2); x2 ≥ 0}) and the lower one is denoted by Ωf2 (Ωf2 = {x(x1 , x2 ); x2 ≤ −h}). The interfaces between the bone (Ωb) and the fluids (Ωf1 bf and Ωf2 ) are denoted by Γbf 1 and Γ2 , respectively (see Fig. 1). Γf∞ 1 Fluid 1 (Ωf1 ) e2 e1 source (xs1 , xs2 ) receivers Γbf 1 nbΓ1 O Solid (Ωb ) Fluid 2 (Ωf2 ) Γbf 2 nbΓ2 Γf∞ 2 Fig. 1. Geometrical configuration of the trilayer model for ultrasound axial transmission test 228 Vu Hieu Nguyen, Salah Naili Both fluid domains (Ωf1 and Ωf2 ) are considered as homogeneous fluid media and the solid plate is assumed to be an anisotropic elastic solid. Moreover, according to the cortical bone properties, we assume that mechanical properties of the solid may vary along the e2 -axis but they are homogeneous along the longitudinal direction described by the e1 -axis. The system is excited by an acoustic source located at a point xs = (xs1 , xs2 ) located in the upper fluid domain Ωf1 . 2.2. Governing equations The governing equations for the wave propagation in a coupled fluid-solid system can be found in many references, see e.g. [23]. Here, we only outline the main equations to solve the considered trilayer problem. In what follows, we denote respectively by ∇·, ∇ and ∇2 the divergence, gradient and Laplacian operators with respect to x and by the superposed dot, differentiation with respect to time t. Governing equations in the fluid domains (Ωf1 and Ωf2 ). Let us first consider the fluid domain Ωf1 . Since the considered compressible fluid is assumed to be inviscid, the equation of motion in the domain Ωf1 can be described by using a scalar field p1 , which satisfies the wave equation: 1 f p̈1 − ∇2 p1 = Q̇, ∀x ∈ Ω1 , (1) 2 c1 where p1 (x, t) is the acoustic pressure in the domain Ωf1 , Q(x, t) is the acoustic source density, the constant c1 is the wave velocity in Ωf1 at rest which depends on the bulk p modulus K1 of the fluid and its mass density ρ1 by the relation c1 = K1 /ρ1 . The impulsive infinite line source acting at xs = (0, xs2) is described as: Q̇ = ρ1 F (t)δ(x1 )δ(x2 − xs2 ), (2) where F (t) is a given real scalar function depending only on the time and δ(.) is Dirac’s delta function. Similarly to the equation of motion in the domain Ωf1 , the pressure field p2 in the fluid domain Ωf2 reads: 1 f p̈2 − ∇2 p2 = 0, ∀x ∈ Ω2 , (3) c22 p where c2 = K2 /ρ2 is the wave velocity in Ωf2 , K2 and ρ2 are respectively the bulk modulus of the fluid and its mass density at rest. Governing equations in the solid domain (Ωb). We denote by u(x, t) = {u1 , u2 }T the time-dependent vector of displacement at a point x located in the solid domain Ωb , where the superscript "T" designates the transpose operator. The associated linearized strain and stress tensors are denoted by (u) and σ(u), respectively. Neglecting the body force, the equation of motion in the solid domain Ωb is given by: ρü − ∇ · σ = 0, where ρ denotes the mass density of the solid. ∀x ∈ Ωb (4) Simulation of transient ultrasonic wave propagation in fluid-loaded heterogeneous cortical bone 229 For the finite element development afterwards, it would be more convenient to rewrite the stress and strain tensors in vectorial form: s = {σ11 , σ22 , σ12 }T and e = {11 , 22 , 212}T . Using this notation, the strain-displacement relation is given by:   ∂ 0   ∂x1  ∂    e = Lu, L =  0 (5) ,  ∂x2   ∂ ∂  ∂x2 ∂x1 and the constitutive equation describing the anisotropic elastic behavior of the solid reads:   c11 c12 c16 s = C e, with C = c12 c22 c26  , (6) c16 c26 c66 where C is the elasticity tensor written in the Voigt matrix notation. The equation of the motion (4) may be restated as follows: ρü − LT s = 0, (7) which becomes in substituting the constitutive equation (6): ρü − LT CLu = 0. (8) We recall that the material properties of the solid depend only on x2 , i.e. ρ = ρ(x2 ) and C = C(x2 ). Boundary and interface conditions. The interface conditions at the fluid-solid bf interfaces (Γbf 1 and Γ2 ) may be described by the same way. The condition required for the continuity of normal velocities at the interfaces reads: ∇pα · nfΓα = −ρα üΓα · nfΓα , ∀x ∈ Γbf α (α = 1, 2) (9) whereas the condition of the continuity of normal stresses at the interfaces requires: σnbΓα = −pα nbΓα , ∀x ∈ Γbf α (α = 1, 2) (10) f where nfΓα (α = 1, 2) is the unit normal vector at the interface Γbf α pointing out of Ωα and b nbΓα (α = 1, 2) is the unit normal vector of the interface Γbf α pointing out of Ω (see Fig. 1). Finally, the radiation conditions of two fluid halfspaces in the far field reads: pα → 0, ∀x → Γfα∞ (α = 1, 2). (11) Boundary value problem in terms of (p1 , u, p2 ). This paragraph rewrites all previous equations in a more compact form. Note that for the considered configuration (see Fig. 1) wherein nbΓ1 = −nfΓ1 = {0, 1}T and nbΓ2 = −nfΓ2 = {0, −1}T , the equations in terms of (p1 , u, p2 ) of the coupled fluid-solid problem may be represented as follows: 230 Vu Hieu Nguyen, Salah Naili ◦ Fluid domain Ωf1 : 1 p̈1 − ∇2 p1 = ρ1 F (t)δ(x1 )δ(x2 − xs2 ), 2 c1 1 ∂p1 = −ü2 , ∀x ∈ Γbf 1 , ρ1 ∂x2 p1 → 0, ∀x ∈ Ωf1 , ∀ x → Γf1 ∞ . (12) (13) (14) ◦ Solid domain Ωb : ρü − LT CLu = 0, T t = {0, −p1 } , t = {0, −p2 }T , ∀x ∈ Ωb , ∀x ∈ ∀x ∈ Γbf 1 , Γbf 2 , (15) (16) (17) where t := {σ12 , σ22 }T . ◦ Fluid domain Ωf2 : 1 p̈2 − ∇2 p2 = 0, ∀x ∈ Ωf2 , c22 1 ∂p2 = −ü2 , ∀x ∈ Γbf 2 , ρ2 ∂x2 p2 → 0, ∀x → Γf2 ∞ . (18) (19) (20) 3. EQUATIONS IN THE LAPLACE-FOURIER DOMAIN The problem presented in (12) - (20) deals with solving a system of linear partial differential equations in which the coefficients are homogeneous in the longitudinal direction given by x1 -axis. Here, we propose to solve the system as follows: (i) the system of equations is firstly transformed into frequency-wavenumber domain by using a Fourier transform with respect to x1 combined with a Laplace transform with respect to t; (ii) in the frequency-wavenumber domain, the equations in both fluid domains are solved analytically giving impedance boundary and interface 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. The general form of a Laplace-Fourier transform (LF) applied to a real-valued function y(x1 , x2 , t) denoted by ỹ(k1 , x2 , s) is defined as: ỹ(k1 , x2 , s) : = L{F[y(x1, x2 , t)]}  Z ∞ Z +∞ −ik1 x1 = y(x1 , x2 , t)e dx1 e−st dt, 0 √ (21) −∞ where i = −1, s ∈ C is the complex Laplace variable, k1 ∈ R is the real Fourier variable representing the wavenumber on the x1 -axis, R and C denote the set of all the real and complex numbers, respectively. Simulation of transient ultrasonic wave propagation in fluid-loaded heterogeneous cortical bone 231 In (s-k1 ) domain, the time derivative and the spatial derivative with respect to x1 ∂(.) ∂(.) can be replaced by → −s (.) and → ik1 (.), respectively. ∂t ∂x1 3.1. Transformed problem in (s - k1 ) domain for the fluids Ωf1 and Ωf2 By applying the Laplace-Fourier transform (21) to Eqs. (12) - (14), the problem so defined is reduced to be a boundary-value differential equation for p̃1 with respect only to x2 :  2  ∂ 2 p̃1 s 2 + k p̃ − = ρ1 F̃0 (s)δ(x2 − xs2 ), (22) 1 1 c21 ∂x22   1 ∂ p̃1 = −s2 ũ2 (k1 , 0, s), (23) ρ1 ∂x2 x2 =0 p̃1 → 0, ∀x2 → +∞, (24) of which the solution for p̃1 may be expressed in semi-explicit form as:  ρ 1 ρ1 ˜  −α1 (xs −x2 ) s 1 2 2 p̃1 = − F0 e + e−α1 (x2+x2 ) + s Ũ21 e−α1 x2 , for 0 ≤ x2 ≤ xs2 , (25) 2 α1 α1  ρ 1 ρ1 ˜  α1 (xs2 −x2 ) s 1 2 p̃1 = − F0 e + e−α1 (x2 +x2 ) + s Ũ21 e−α1 x2 , for x2 ≥ xs2 , (26) 2 α1 α1 q 2 where α1 := sc2 + k12 , Ũ21 is the LF transform of the vertical displacement of the solid 1 layer at the upper fluid-solid interface Γbf 1 which is defined by Ũ21 := ũ2 (k1 , 0, s). Similarly, applying the LF-transform to Eqs. (18) - (20), the general solution of p̃2 may be obtained: p̃2 = − where α2 = q s2 c22 ρ2 2 s Ũ22 eα2 (x2+h) , α2 ∀x2 ≤ −h, (27) + k12 and Ũ22 := ũ2 (k1 , −h, s). 3.2. Transformed problem in (s - k1 ) domain for the solid (Ωb ) By applying the LF-transform (21) to the system of equations associated with the solid given by Eqs. (15) - (17), we obtain: ρs2 ũ − L̃T CL̃ũ = 0, (28) where   ũ1 , ũ2 ∂ L̃ = ik1 L1 + L2 ∂x2     1 0 0 0 with L1 = 0 0 and L2 = 0 1 . 0 1 1 0 ũ = (29) (30) 232 Vu Hieu Nguyen, Salah Naili For each value of the couple (s, k1 ) ∈ C × R, the differential equations given by Eq. (28) lead to a system for ũ with respect only to x2 :

∂ ũ ∂ 2 ũ s2 A1 + k12 A2 ũ − ik1 A3 + AT3 − A4 = 0, ∂x2 ∂x22 (31) where the matrices A1 , A2 , A3 and A4 , which depend only on the physical parameters, are defined by:     ρ 0 c11 c16 T A1 = ; A2 = L1 CL1 = ; (32) 0 ρ c16 c66     c c c c A3 = LT2 CL1 = 16 66 ; A4 = LT2 CL2 = 66 26 . (33) c12 c26 c26 c22 The interface conditions given by Eqs. (16) - (17) in (s - k1 ) domain read:     0 0 , , t̃(−h) = t̃(0) = −p̃1 (0) −p̃2 (−h) (34) bf in which the solution for p̃1 (0) at Γbf 1 and the solution for p̃2 (−h) at Γ2 can be determined by using Eqs. (26) and (27), respectively:  ρ1  s p̃1 (0) = − F̃0 e−α1 x2 − s2 Ũ21 , (35) α1 ρ2 p̃2 (−h) = − s2 Ũ22 . (36) α2 Thus, the interface conditions given by Eq. (34) can be rewritten as: t̃(0) = F0 − P1 ũ(0), t̃(−h) = P2 ũ(−h), (37) where ! 0 F0 = ρ1 F̃ e−α1 xs2 , 0 α " 1 # " 0 0 0 P1 = 0 ρ1 s2 , P2 = 0 α1 (38) # 0 ρ2 2 . s α2 (39) Moreover, it would be useful to note that the constitutive relation of the solid given by Eq. (6) in the (s k1 ) domain may be restated as: s̃ = CL̃ũ, (40) and the vector t̃ (which is given by t̃ = LT2 s̃) can be expressed by: t̃ = ik1 A3 ũ + A4 ∂ ũ . ∂x2 (41) Simulation of transient ultrasonic wave propagation in fluid-loaded heterogeneous cortical bone 233 4. SPECTRAL FINITE ELEMENT FORMULATION AND TIME-SPACE SOLUTION 4.1. Weak formulation in the solid domain The weak formulation of the boundary-value problem given by Eqs. (31) - (37) may be now introduced using the classic procedure (see e.g. [15]). Let Cad be the function space which comprises all sufficiently smooth complex-valued admissible functions: x2 ∈ ] − h, 0[→ δ ũ(x2 ) ∈ C × C. The conjugate transpose of δ ũ is denoted by δ ũ∗ . Upon integrating the Eq. (31) against a test function δ ũ∗ ∈ Cad and integrating by parts, we obtain:     Z 0 Z 0 ∂δ ũ∗ ∂ ∗ 2 2 T ∂ δ ũ s A1 + k1 A2 − ik1 A3 ik1 A3 + A4 ũ dx2 + ũ dx2 ∂x2 ∂x2 −h −h ∂x2    0 ∂ ũ = 0. (42) − δ ũ∗ ik1 A3 + A4 ∂x2 −h The last terms associated with boundary conditions given by Eq. (42) may be calculated by using the relations (37) and (41):  0    0 ∂ ∗ ũ = δ ũ∗ t̃ −h δ ũ ik1 A3 + A4 ∂x2 (43) −h ∗ ∗ ∗ = δ ũ (0)F0 − δ ũ (0)P1 ũ(0) − δ ũ (−h)P2 ũ(−h). The weak formulation of the differential equation (31) reads: For all k1 fixed in R and for all s fixed in C, find ũ(k1 , x2 , s) ∈ Cad such that:     Z 0 Z 0 ∂δ ũ∗ ∂ ∗ 2 2 T ∂ ũ dx2 + ik1 A3 + A4 ũ dx2 δ ũ s A1 + k1 A2 − ik1 A3 ∂x2 ∂x2 −h ∂x2 −h + δ ũ∗ (0) P1 ũ(0) + δ ũ∗ (0) P2 ũ(−h) = δ ũ∗ (0) F0 , ∀ δ ũ ∈ Cad . (44) 4.2. Finite element formulation We proceed by introducing aSfinite element mesh of the domain [−h, 0] which contains nel elements Ωe : [−h, 0] = e Ωe (e = 1, ..., nel). By the Galerkin finite element method, both functions ũ and δ ũ in each element Ωe are approximated using the same shape function: ũ(x2 ) = Ne Ue , δ ũ(x2 ) = Ne δUe , ∀ x2 ∈ Ωe , (45) where Ne is the shape function, Ue and δUe are the vectors of nodal solutions of ũ and δ ũ within the element Ωe , respectively. Substituting (45) into (44) and assembling the elementary matrices, we obtain a linear system of equations:
K + KΓ U = F, (46) where U is the global nodal displacement vector, K is the global "stiffness matrix" of the solid, KΓ represents the operator of coupling between the fluids and the solid and F is the external force vector. For all couples (s, k1 ) fixed in C × R, these quantities may be 234 Vu Hieu Nguyen, Salah Naili expressed by: K = s2 KA1 + k12 KA2 + ik1 KA3 + KA4 ,   ρ1 2 ρ2 2 Γ s , 0, ..., 0, s , K = Diag 0, α1 α2  T ρ1 −α1 xs2 F = 0, F̃0 e , 0, ..., 0 , α1 (47) (48) (49) where the matrices KA1 , KA2 , KA3 and KA4 are independent of s and k1 and are defined by: [Z [Z T KA1 = Ne A1 Ne dx2 , KA2 = NTe A2 Ne dx2 , (50) e KA3 = [Z e Ωe Ωe e n T N0 e A 3 Ne o a dx2 , Ωe KA4 = [Z e Ωe T N0 e A4 N0 e dx2 , (51) in which the notation {.}a is devoted for the skew part of the {.} and {.}0 means the differentiation with respect to x2 . Remark 1. We recall that the solid layer is assumed to be heterogeneous and the elementary matrices A1 , A2 , A3 and A4 (see Eqs. (32) - (33)) are not constant but vary with respect to x2 . Here, the elementary matrices (see Eqs. (50) - (51)) are numerically computed by using the Gauss quadrature technique. However, when the material properties of the solid layer are homogeneous, analytical expressions of elementary matrices may be determined. Remark 2. To reduce the computational cost, the matrices KA1 , KA2 , KA3 and KA4 have only to be computed once before performing a loop on s and k1 . For each value of s and k1 , the global matrix K is obtained by performing the summation (see Eq. (47)). Moreover, one may note that: K(−k1 , s) = KT (k1 , s). (52) 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 in the complex domain (46). The solutions for p̃1 and p̃2 in two fluid domains may be then determined by using equations (25) and (26), respectively. In order 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 FFT (Fast Fourier Transform) technique. The inverse Laplace transform is carried out using the Quadrature Convolution Method which has been shown to be a very efficient technique for computing the time response solution in many dynamic problems [24]. Let N1 and ∆k1 respectively the sampling number and sampling rate of the wavenumber k1 using for the FFT procedure. The space-solution is reconstructed upon on a broad