Determination of the parameters of cancellous ... - Semantic Scholar
ARTICLE IN PRESS
Mathematical and Computer Modelling (
)
– www.elsevier.com/locate/mcm
Determination of the parameters of cancellous bone using high frequency acoustic measurements James L. Buchanan a , Robert P. Gilbert b,∗ a Mathematics Department, United States Naval Academy, Annapolis, MD, United States b Department of Mathematical Sciences, University of Delaware, Newark, DE, United States
Received 3 March 2006; accepted 9 May 2006
Abstract The Biot model is widely used to model poroelastic media. Several authors have tested its applicability to cancellous bone, but to do so requires a priori estimation of the parameters of the Biot model, which is an uncertain and expensive endeavor. A method of computing acoustic pressure in the low 100 kHz range is developed. c 2006 Elsevier Ltd. All rights reserved.
Keywords: Osteoporosis; Cancellous bone; Poroelastic media; Biot model; Inverse problem; Parallel processing; Simplex method
1. Introduction Cancellous bone is a two component material consisting of a calcified bone matrix with interstial fatty marrow. Faced with the clinical need to understand better the effects of aging and disease on the elastic and strength properties of trabecular bone, researchers have begun to investigate whether factors such as a architecture and tissue quality can account for 20%–40% of unexplained variance in these mechanical properties [10]. Quantitative ultrasound QUS attempts to interogate the bone structure and determine whether the subject has osteoporosis; however, QUS is still in its infancy, Chaffai et al. [6]. It relies on assessment of two fundamental parameters: broadband ultrasound attenuation BUA and speed of sound SOS. The magnitude and frequency dependence of ultrasound attenuation is a complicated function of the composition and the micro-architecture of the propagating medium. However, one needs to note a lack of theoretical model that describes how the the attenuation is influenced by properties of bone. Hodgskinson et al. [7] studied the ability of ultrasound velocity to predict the Young’s modulus of elasticity of cancellous bone. QUS measurements have been successfully used recently to predict fracture risk in osteoporosis. This technique is now incorporated into commercial devices enabling attenuation to be measured. Hence mathematical models of poroelastic media have been applied to the study of bone architecture. Indeed, McKelvie and Palmer [9], Williams [11], and Hosokawa and Otani [8] discuss the application of the Biot model for a poroelastic medium
∗ Corresponding author. Tel.: +1 302 831 2315; fax: +1 302 831 4511.
E-mail address: [email protected] (R.P. Gilbert). c 2006 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2006.05.007
ARTICLE IN PRESS 2
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
to cancellous bone. In the experiments of McKelvie and Palmer and Hosokawa and Otani a small specimen of cancellous bone was placed in a tank of water between an acoustic source and receiver. An attempt was made to align the trabecular direction of the specimen along the source–receiver axis so that the isotropic Biot model would be applicable. In order to compare the predictions of the Biot model the experimenters needed estimates of the model’s parameters. These values were taken from the literature, obtained by experimental measurement, or in the case of permeability characterized as estimates without elaboration. In view of the acquisition expense and uncertainty of such estimates, it is worthwhile to explore whether these parameters can be ascertained from the measurements of the acoustic field generated by such experiments. In [5] and [4] Buchanan, Gilbert and Khashanah investigated the extent to which the most important parameters of the Biot model could be recovered by acoustic interrogation in a numerical experiment which simulated in two dimensions the physical experiment of Hosokawa and Otani. The finite element method was used to simulate both the target pressure data. This had the virtue that boundary conditions could be imposed on all four sides of the water tank used in the experiment, but the drawback that the long computation times necessitated the use of frequencies around 10 kHz and unrealistically small tank sizes. In this article we use a different mathematical technique, contour integration of the Green’s function, for the simulation. This has the virtue of permitting use of frequencies in the low 100 kHz range, which is closer to the frequencies used in the experiments and of more realistic tank dimensions. Its drawback is that the tank and bone specimen must be taken to be infinite in one direction. 2. The Biot model applied to cancellous bone The Biot model treats a poroelastic medium as an elastic frame with interstial pore fluid. Cancellous bone is anisotropic, however, as pointed out by Williams, if the acoustic waves passing through it travel in the trabecular direction an isotropic model may be acceptable. We will simulate a two dimensional version of the experiments described in McKelvie and Palmer and Hosokawa and Otani. The motion of the frame and fluid within the bone are tracked by position vectors u = [u, v] and U = [U, V ]. The constitutive equations used by Biot are those of a linear elastic material with terms added to account for the interaction of the frame and interstial fluid σx x = 2µex x + λe + Q, σ yy = 2µe yy + λe + Q, σx y = µex y ,
(1)
σ yx = µe yx ,
σ = Qe + R where the solid and fluid dilatations are given by e =∇·u=
∂v ∂u + , ∂x ∂y
=∇·U =
∂U ∂V + . ∂x ∂y
∂u ∂v + , ∂y ∂x
e yy =
(2)
The stress–strain relations are ex x =
∂u , ∂x
ex y = e yx =
∂v . ∂y
(3)
The parameter µ, the complex frame shear modulus can be measured. The other parameters λ, R and Q occurring in the constitutive equations are calculated from the measured or estimated values of the parameters given in Table 1 using the formulas 2 (K r − K b )2 − 2β K r (K r − K b ) + β 2 K r2 λ = Kb − µ + 3 D − Kb 2 2 β Kr R= D − Kb β K r ((1 − β)K r − K b ) Q= D − Kb
(4)
where D = K r (1 + β(K r /K f − 1)).
(5)
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
3
–
Table 1 Parameters in the Biot model Symbol
Parameter
ρf ρr Kb µ Kf Kr β η k α a
Density of the pore fluid Density of frame material Complex frame bulk modulus Complex frame shear modulus Fluid bulk modulus Frame material bulk modulus Porosity Viscosity of pore fluid Permeability Structure constant Pore size parameter
The bulk and shear moduli K b and µ are often given imaginary parts to account for frame inelasticity. Eqs. (1)–(3) and an argument based upon Lagrangian dynamics are shown in [2] to lead to the following equations of motion for the displacements u, U and dilatations e, ∂ ∂2 (ρ11 u + ρ12 U) + b (u − U) ∂t ∂t 2 (6) ∂2 ∂ ∇[Qe + R] = 2 (ρ12 u + ρ22 U) − b (u − U). ∂t ∂t Here ρ11 and ρ22 are density parameters for the solid and fluid, ρ12 is a density coupling parameter, and b is a dissipation parameter. These are calculated from the inputs of Table 1 using the formulas µ∇ 2 u + ∇[(λ + µ)e + Q] =
ρ11 = (1 − β)ρr − β(ρ f − mβ) ρ12 = β(ρ f − mβ) ρ22 = mβ 2 p F a ωρ f /η β 2 η b= k where αρ f m= β and the multiplicative factor F(ζ ), which was introduced in [3] to correct for the invalidity of the assumption of Poiseuille flow at high frequencies, is given by F(ζ ) =
ζ T (ζ ) 1 4 1 − 2T (ζ )/iζ
(7)
where T is defined in terms of Kelvin functions T (ζ ) =
ber0 (ζ ) + ibei0 (ζ ) . ber(ζ ) + ibei(ζ )
For time-harmonic vibrations u(x, y, t) = u(x, y)e−iωt , U(x, y, t) = U(x, y)e−iωt the equations become µ∇ 2 u + ∇[(λ + µ)e + Q] = −ω2 (ρ11 u + ρ12 U) − ib(u − U)
(8)
∇[Qe + R] = −ω (ρ12 u + ρ22 U) + ib(u − U).
(9)
2
Taking the divergence of these equations yields the system ∇ 2 ((λ + 2µ)e + Q) + p11 e + p12 = 0 ∇ 2 (Qe + R) + p12 e + p22 = 0
(10)
ARTICLE IN PRESS 4
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Table 2 Estimated values of some Biot parameters at different porosities taken from McKelvie and Palmer or Hosokawa and Otani β
k
a
0.72 0.75 0.81 0.83 0.95
5 × 10−9
4.71 × 10−4
7 × 10−9 2 × 10−8 3 × 10−8 5 × 10−7
8.00 × 10−4 1.20 × 10−3 1.35 × 10−3 2.20 × 10−3
α
Re K b
Re µ
1.10 1.08 1.06 1.05 1.01
3.18 × 109
1.30 × 109 1.10 × 109 7.38 × 108 6.27 × 108 1.05 × 108
2.69 × 109 1.80 × 109 1.55 × 109 2.57 × 108
where p11 := ω2 ρ11 + iωb,
p12 := ω2 ρ12 − iωb,
p22 := ω2 ρ22 + iωb.
(11)
The article of McKelvie and Palmer contains estimates of the Biot parameters of cancellous bone in the human os calcis (heel bone) for the normal (β = 0.72) and severely osteoporotic (β = 0.95) cases while the article of Hosokawa and Otani has estimates for bovine femoral bone for porosities of β = 0.75, 0.81 and 0.83. Table 2 contains estimates of these six Biot parameters for five bone specimens. In obtaining them we followed the estimation procedures used by Williams, McKelvie and Palmer, and Hosokawa and Otani. In generating test problems for a parameter recovery algorithm an estimate of the range of values a parameter might take is needed. Here is a discussion of how the values of the Biot parameters in Table 2 were calculated and our estimate of the range of values for each parameter: • The real parts of K b and µ were calculated using the formulas of Williams E Vn 3(1 − 2ν) f E Vn Re µ = 2(1 + ν) f
Re K b =
•
•
•
•
(12)
used by Hosokawa and Otani. Here V f = 1 − β is the bone volume fraction. Theoretically n = 1 for waves travelling in the trabecular direction and is between 2 and 3 for transverse waves, however there is enough randomness in the trabecular direction in cancellous bone that authors have empirically adjusted the exponent to agree better with experiment. Williams arrived at a value of n = 1.23 based on comparing the Biot predictions for Type I compressional and shear wave velocity assuming the form (12) to the measured speeds obtained from experiments conducted on samples taken from bovine tibia. Hosokawa and Otani found that n = 1.46 agreed well with their data from experiments on bone specimens from bovine femora. We followed Hosokawa and Otani who used the values E = 2.2 × 1010 , ν = 0.32 for the Young’s modulus and Poisson ratio of solid bone. For our parameter range we assumed that an attempt was made to align the source–receiver axis with the trabecular direction and took 1.1 ≤ n ≤ 1.6, an interval which symmetrically encompasses Hosokawa and Otani’s and William’s values of n. The imaginary parts of K b and µ were calculated using a log decrement `: Im K b∗ = `Re K b∗ /π, Im µ = `Re µ∗ /π with a value ` = 0.1 which is typical of that used in underwater acoustics. There appears to be little sensitivity to these parameters, however. The structure factor was calculated using the formula of Berryman [1] α = 1 − s(1 − 1/β). Williams and Hosokawa and Otani used s = 0.25. A value s = 0.227 was suggested by the finite element analysis of Yavari and Bedford [12]. We will calculate α assuming 0.2 ≤ s ≤ 0.3. The pore size parameter was estimated by McKelvie and Palmer using electron microscopy and by Hosokawa and Otani using x-ray examination. Fig. 1 shows that the measurements suggest a linear relationship between pore size and porosity. We took the uncertainty in this parameter to be the regression line value ±0.0002. Permeability is a difficult parameter to estimate. Fig. 2 shows the values for the five specimens of McKelvie and Palmer and Hosokawa and Otani. The estimates indicate that permeability is approximately a log-linear function of porosity. However McKelvie and Palmer characterize their values without elaboration as “estimates” Hosokawa and Otani state that their estimates are based upon those of McKelvie and Palmer. Hence the apparent log-linear
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
5
Fig. 1. Values of the pore size parameter measured by McKelvie and Palmer and Hosokawa and Otani. Solid line: Regression line for the measured values. Dashed lines: uncertainty interval assigned to this parameter.
Fig. 2. Values of the pore size parameter estimated by McKelvie and Palmer and Hosokawa and Otani. Solid line: Regression line for the estimated values. Solid curve: Estimate of permeabilty given by the Kozeny–Carmen equation. Dashed lines: uncertainty interval assigned to this parameter.
relation should be regarded circumspectly. Indeed according the Kozeny–Carmen formula βa 2 , (13) 4K where K ≈ 5 is an empirical constant, the relation is not log-linear. Fig. 2 shows that if pore size is indeed a linear function of porosity as indicated by Fig. 1, then permeability, as predicted by (13), will deviate significantly from log-linearity. We took the uncertainty in this parameter to be ±100.75 so that the Kozeny–Carmen values would be encompassed as well as the regression values. k=
Table 3 gives the values we shall use for the other Biot parameters. The value for ρ f is from McKelvie and Palmer, but the value used by Hosokawa and Otani, 930, is similar. Likewise both sets of authors used about the same value for viscosity η. The fluid bulk modulus K f is from Hosokawa and Otani. The frame material densities used by McKelvie and Palmer and Hosokawa and Otani were somewhat different. Williams reports a range of estimates for bovine
ARTICLE IN PRESS 6
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Table 3 Parameters for cancellous bone to be used for all specimens Parameter
Symbol
Value
Pore fluid density Fluid bulk modulus Pore fluid viscosity Frame material density Frame material bulk modulus
ρf Kf η ρr Kr
950 2.00 × 109 1.5 1960 2.00 × 1010
cortical bone of ρr = 1930 → 2000. We use Hosokawa and Otani’s value ρr = 1960. The frame material bulk modulus was calculated from (12) with V f = 1. 3. Mathematical simulation of the experiment 3.1. A Green’s function for acoustic pressure in water We formulate a mathematical model for calculating acoustic pressure in an experiment reminiscent of Hosokawa and Otani’s. A two-dimensional bone slab of width w is infinite in the y-direction and is in a tank of water with rigid walls at x = wl and x = wr . A sound source is located at the origin, the left edge of the bone slab is at x = bl and thus the right edge is at br = bl + w. In the water pressure p(x, y)e−iωt and horizontal displacement Uxw (x, y)e−iωt satisfy ∇ 2 p + k02 p = −δ(x)δ(y)
(14)
∂p + ρ0 ω2 Uxw = 0 ∂x
(15)
where k0 = ω/c0 , and c0 and ρ0 are the speed of sound in water and the density of water. To construct a Green’s function we start by seeking solutions of the form p(x, y, κ) = G x (x, κ)G y (y, κ)
(16)
where G 00y + k02 κG y = −δ(y)
(17) 2 G 00x + aw G x = −δ(x) √ with aw := k0 1 − κ. √ If the branch cut of the square root function is chosen so that Im κ > 0, then in order for G y → 0 as y → ±∞ we must have ( √ κy , y≥0 Aeik0 √ G y (y, κ) = −ik0 κ y Be , y < 0. To satisfy (17) we require that G y (0+, κ) − G y (0−, κ) = 0 G 0y (0+, κ) − G 0y (0−, κ) = −1. This gives √ 1 √ eik0 κ y , y ≥ 0 2ik0 κ G y (y, κ) = √ 1 − √ e−ik0 κ y , y < 0. 2ik0 κ
−
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
7
Since pressure is symmetric about y = 0, it suffices to construct the solution for y ≥ 0. We henceforth take p(x, y, κ) = −
G x (x, κ) ik0 √κ y √ e 2ik0 κ
Assuming that the walls of the tank are rigid, Uxw = 0 at x = wl , wr . From (15) we then have G 0x = 0 at x = wl , wr . This gives C cos aw (x − wl ), wl ≤ x < 0 1 sin aw x G x (x, κ) = A cos aw x + B , 0 < x < bl aw C8 cos aw (x − wr ), bl + w < x ≤ wr . The conditions G x (0+, κ) − G x (0−, κ) = 0 G 0x (0+, κ) − G 0x (0−, κ) = −1 give A = C1 cos aw wl ,
B = C1 aw sin aw wl − 1
and hence G x (x, κ) = C1 cos aw wl cos aw x + (C1 aw sin aw wl − 1) = C1 cos aw (x − wl ) −
sin aw x , aw
sin aw x aw
0 < x < bl .
3.2. Solution of the Biot equations In (10) we make the change of dependent variables τ := (λ + 2µ)e + Q,
σ := Qe + R,
(18)
the inverse transformation for which is e = a11 τ − a12 σ,
= −a12 τ + a22 σ
(19)
where a11 := R/d,
a22 := (λ + 2µ)/d
a12 := Q/d,
with d := R(λ + 2µ) − Q 2 . In terms of the new variables (10) becomes ∇ 2 τ + B11 τ + B12 σ = 0
(20)
∇ 2 σ + B21 τ + B22 σ = 0 where B11 := a11 p11 − a12 p12 ,
B12 := −a12 p11 + a22 p12 ,
B21 := a11 p12 − a12 p22 ,
B22 := −a12 p12 + a22 p22 .
From (9) the fluid displacement vector is U=−
1 (∇σ + p12 u). p22
(21)
ARTICLE IN PRESS 8
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Substituting this into (8) gives a partial differential equation for the frame displacement vector ∇ 2 u + A31 ∇τ + A32 ∇σ + B33 u = 0
(22)
where A31 :=
1 − µa11 , µ
A32 := a12 −
p12 , µp22
B33 :=
2 p11 p22 − p12 . µp22
In order to get the solutions to match across interfaces, we assume that all unknowns have the same dependence on y as p(x, y, κ) Uxw (x, y) = Uxw (x)eik0 τ (x, y) = τ (x)eik0
√ κy
√
κy , √ ik0 κ y
u x (x, y) = u x (x)e
σ (x, y) = σ (x)eik0 ,
√
κy
u y (x, y) = u y (x)eik0
√ κy
and so forth. Substituting this into (20) gives τ 00 (x) + aτ2 τ (x) + B12 σ (x) = 0 σ 00 (x) + B21 τ (x) + aσ2 σ (x) = 0 q q where aτ := B11 − k02 κ and aσ := B22 − k02 κ. Elimination of σ in (23) gives
(23)
τ 0000 + (aτ2 + aσ2 )τ 00 + (aτ2 aσ2 − B12 B21 )τ = 0 which has the general solution τ (x) = C2 cos m + (x − bl ) + C3
sin m − (x − bl ) sin m + (x − bl ) + C4 cos m − (x − bl ) + C5 m+ m−
(24)
where s m ± := s
aτ2 + aσ2 ±
p
(aτ2 + aσ2 )2 − 4(aτ2 aσ2 − B12 B21 ) 2
(aτ2 − aσ2 )2 + 4B12 B21 2 s p B11 + B22 − 2k02 κ ± (B11 − B22 )2 + 4B12 B21 = 2 s p B11 + B22 ± (B11 − B22 )2 + 4B12 B21 = − k02 κ. 2 =
aτ2 + aσ2 ±
p
From (24) and (23) we have σ (x) = −
1 (τ 00 (x) + aτ2 τ (x)). B12
The horizontal displacements can be calculated from (21) and (22) u 00x (x) + A31 τ 0 (x) + A32 σ 0 (x) + au2 u x (x) = 0 q where au := B33 − k02 κ and Ux (x) = −
1 (σ 0 (x) + p12 u x (x)). p22
(25)
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
9
The solution of (25) contains two new arbitrary constants which we denote by C6 and C7 . It is convenient for what ensues to put m ± and au in the same form as aw . Hence we write p √ au = k0 au0 − κ m ± = k0 m 0± − κ, where m 0± =
B11 + B22 ±
p
(B11 − B22 )2 + 4B12 B21 , 2k02
√
Since e(x)eik0 κ y = ∇ · hu x , u y i = computation for the fluid dilatation,
∂u x ∂x
+
∂u y ∂y
au0 =
= u 0x (x)eik0
√
B33 . k02
κy
√ √ + u y (x)ik0 κeik0 κ y we have, after a similar
1 √ (e(x) − u 0x (x)) ik0 κ 1 U y (x) = √ ((x) − Ux0 (x)). ik0 κ u y (x) =
The frame normal and tangential stress in the bone specimen are given by σx x (z) = λe(x) + 2µex x (x) + Q(x) = λe(x) + 2µu 0x (x) + Q(x) ∂u y ∂u x σx y (x) = µ + ∂y ∂x √ = µ(ik0 κu x (x) + u 0y (x)). Finally from (14) 2 we have Uxw (x) = −
1 G 0 (x, κ). ρ0 ω 2 x
3.3. Interface conditions At an interface x = xd between the water and the bone slab the following conditions are required Uxw (xd −) = βUx (xd +) + (1 − β)u x (xd +) G x (xd −) = σx x (xd +) + σ (xd +) G x (xd −) = σ (xd +)/β
(26)
0 = σx y (xd +). These conditions derive from continuity of flux, aggregate normal pressure, pore fluid pressure and tangential frame stress respectively. The conditions (26) applied at x = bl and x = bl + w lead to a system of eight linear equations in the unknowns C1 , . . . , C8 C1 .. M . = v. (27) C8 The general solution has the form p(x, y) =
k02 2πi
I Cκ
p(x, y, κ)dκ
(28)
ARTICLE IN PRESS 10
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
where Cκ is a suitably chosen contour in the κ-plane. Substituting this into (14) gives ∇ 2 p + k02 p =
k02 2πi
I Cκ
k4 + 0 2π i
I
k2 =− 0 2π i
I
k02 2π i
I
−
= −δ(x)
G 00x (x, κ)G y (y, κ)dκ +
Cκ
Cκ
k02 2πi
I Cκ
G x (x, κ)G 00y (y, κ)dκ
G x (x, κ)G y (y, κ)dκ
[k02 κG x (x, κ) + δ(x)]G y (y, κ)dκ G x (x, κ)[k02 (1 − κ)G y (y, κ) + δ(y)]dκ +
Cκ I 2 k0
2π i
k2 G y (y, κ)dκ − δ(y) 0 2π i Cκ
I Cκ
k04 2πi
I Cκ
G x (x, κ)G y (y, κ)dκ
G x (x, κ)dκ.
Since k02 2π i
I Cκ
G x (x, κ)dκ = δ(x),
for a positively oriented contour Cκ which encloses the singularities of G x , p is a solution of (14) providing Cκ is chosen to exclude the singularities of G y and enclose those of G x . Substituting the results of the previous section into (28) gives I k02 G x (x, κ) ik0 √κ y p(x, y) = − e dκ (29) √ 2π i Cκ 2ik0 κ I k0 G x (x, κ) ik0 √κ y = e dκ √ 4π Cκ κ k0 I C1 (κ) cos aw (x − wl ) ik0 √κ y e dκ, wl < x < 0 √ 4π Cκ κ k I C1 (κ) cos aw (x − wl ) − sinaaww x ik √κ y 0 = (30) e 0 dκ, 0 < x < bl √ 4π Cκ κ I C8 (κ) cos aw (x − wr ) ik0 √κ y k0 e dκ, br < x < wr . √ 4π Cκ κ √ The only singularity of G y is the κ, thus it can be excluded by taking Cκ to be a cut enclosing the positive κ-axis. Thus the first two entries in (30) are the same since sinaaww x is an entire function of κ. Observing that all of the functions used in the matrix M in (27) are analytic except possibly on the positive κ-axis, it follows that the only singularities of G x are the solutions {κn } of ∆ := det M = 0, which are then the eigenvalues of the problem. According to Cramer’s ∆ (κ) rule C1,8 = ∆1,8(κ) and thus from the residue theorem I k0 ∆1 (κ) cos aw (x − wl ) ik0 √κ y e dκ, wl < x < bl √ 4π ∆(κ) κ Cκ I p(x, y) = k0 ∆8 (κ) cos aw (x − wr ) ik0 √κ y e dκ, br < x < wr . √ 4π Cκ ∆(κ) κ √ ik0 X ∆1 (κn ) cos aw (κn )(x − wl )eik0 κn y , wl < x < bl √ d∆ 2 κn dκ (κn ) n = X √ ik0 ∆8 (κn ) cos aw (κn )(x − wr )eik0 κn y , br < x < wr . √ d∆ 2 κn (κn ) n dκ
(31)
(32)
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
11
Fig. 3. Level curves for ∆(κ) at 100 kHz.
3.4. Numerical implementation of the pressure equation The first of the four interface conditions involves displacement, whereas the remaining three involve pressure. Typically the two quantities differ by more than 10 orders of magnitude. This can lead to inaccuracies in calculating the characteristic determinant ∆(κ), but the problem is easily resolved by scaling the pressure equations by some characteristic pressure. We used ρ0 ω2 , the compressional Lam´e coefficient for water. We are interested in determining the pressure field in the ultrasound regime, say for frequencies of 100 kHz and higher. For our initial explorations we will use wl = −1, bl = 0.01, br = 0.02, wr = 1.02 m. For the Biot parameters for bone we will use those for cancellous bone of porosity 0.72. Fig. 3, which shows the level curves Cn : ∆(κ) = 10n , n = 0, 50, 100, . . . , 300 for 100 kHz with the stated parameters, indicates a more formidable difficulty in implementing the Eqs. (31) and (32). The magnitude of ∆(κ) ascends rapidly as κ moves away from the segment Re κ ≤ 1. The curve C300 is nearly the largest that can be computed without rescaling since the largest representable real number in IEEE floating point is Ω ≈ 10308 . The problem is the presence of terms of the form sin ad and cos ad where a is one of the coefficients aw , m ± or au and d is a distance. For instance the outer two columns of the characteristic matrix are of the form M sin a d · · · 0 11
w l
M21 cos aw dl M31 cos aw dl 0 M= 0 0 0 0
···
0 0 0 M58 sin aw dr M68 cos aw dr M78 cos aw dr 0
where dl = bl − wl , dr = wr − br and the coefficients M jk depend upon κ but do not grow exponentially as κ moves away from Re κ ≤ 1. Likewise the interior columns of M contain exponentially growing terms 0 M ··· M 0 13
M2...7
M22 M32 0 = M52 sin m + w M62 cos m + w M cos m w + 72 M82 sin m + w
0 0 M43 M53 cos m + w M63 sin m + w M73 sin m + w M83 cos m + w
16
···
0 0 M46 M56 cos au w M66 sin au w 0 M86 cos au w
M27 0 0 . M57 sin au w M67 cos au w 0 M87 sin au w
ARTICLE IN PRESS 12
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 4. Level curves for Γ (κ) at 100 kHz.
Fig. 5. Level curves for Γ (κ) at 500 kHz.
√ To deal with such terms we let a := k0 a0 − κ = ξ + ηi, η ≥ 0 and rewrite them as eiξ d e−2ηd − e−iξ d 2i iξ d −2ηd + e−iξ d ηd e e cos ad = e 2 sin ad = eηd
(33)
and factor the leading exponential out of the determinant. This gives ∆ = eηw dl eηw dr e2η+ w e2η− w e2ηu w Γ where aw = ξw + ηw i, m ± = ξ± + η± i, au = ξu + ηu i. Figs. 4–6 show the level curves of Γ at frequencies 100, 500 and 1000 kHz. The value of |Γ | decreases as κ moves away from Re κ ≤ 1. The stability of the computation also decreases in some regions of the κ-plane, more noticeably at the higher frequencies.
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
13
Fig. 6. Level curves for Γ (κ) at 1000 kHz.
The right hand side of (27) is cos aw bl − ρ0 ω2 sin a b w l ρ0 ω2 aw sin aw bl . v= ρ0 ω2 aw 0 .. . 0 iξ d −2η d w rx e w r x +e−iξw dr x Thus we have ∆8 = eηbl eηw dl e2η+ w e2η− w e2ηu w Γ8 . Finally cos aw (x − wr ) = eηw dr x e with 2 dr x = wr − x. Thus the second integrand in (31) is iξw dr x e−2ηw dr x + e−iξw dr x Γ ∆8 (κ) cos aw (x − wr ) ik0 √κ y 8 ik0 √κ y ηw (bl +dr x −dr ) e =e . (34) e e √ √ ∆(κ) κ 2Γ κ Since bl + dr x − dr = bl + br − x, the leading exponential factor will grow for br < x < bl + br and decay for bl + br < x < wr , but what the integrand does of course depends upon the factor Γ8 /Γ also. A similar analysis for wl < x < bl gives iξw dlx e−2ηw dlx + e−iξw dlx Γ ∆1 (κ) cos aw (x − wr ) ik0 √κ y 1 ik0 √κ y ηw (bl +dlx −dl ) e =e e . (35) e √ √ 2Γ κ ∆(κ) κ We now consider the residue expansion (32). Introducing the operators ∂ 1 ∂ ∂ ∂ 1 ∂ ∂ := −i , := +i , κ = κ1 + iκ2 , ∂κ 2 ∂κ1 ∂κ2 ∂κ 2 ∂κ1 ∂κ2 √ we have for a term of the form a = k0 a0 − κ = ξ + iη ∂ √ k0 ∂ξ ∂η k 0 a0 − κ = − √ = +i . ∂κ ∂κ ∂κ 2 a0 − κ
ARTICLE IN PRESS 14
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 7. Eigenvalues at 50 kHz for a receiver positioned at x = 0.03, y = 0 m.
Since a is an analytic function of κ, 0=
∂ξ ∂η dξ ∂η ∂a = +i ⇒ =i ∂κ ∂κ ∂κ dκ ∂κ
which gives k2i ∂η = 0 . ∂κ 4a This gives d∆ = dκ
"
k02 i 4
# dl + dr 1 1 1 ∂Γ ηw dl ηw dr 2η+ w 2η− w 2ηu w Γ+ e e e e e . + 2w + + aw m+ m− au ∂κ
The term ddκΓ is complicated and will be computed numerically. Thus for br < x < wr we have for the terms in (32) √ √ eηw (bl +dr x −dr ) Γ8 eiξw dr x e−2ηw dr x + e−iξw dr x eik0 κn y ∆8 (κn ) ik0 κn y i = √ h cos aw (κn )(x − wr )e √ d∆ 1 1 1 ∂Γ r κn dκ (κn ) 2 κn k40 i dla+d + 2w + + Γ + m+ m− au ∂κ w whereas for wl < x < bl √
∆1 (κn ) ∆ κn ddκ (κn )
cos aw (κn )(x − wl )e
√ ik0 κn y
√ eηw (bl +dlx −dl ) Γ1 eiξw dlx e−2ηw dlx + e−iξw dlx eik0 κn y i. = √ h 1 1 1 ∂Γ r 2 κn k40 i dla+d + 2w + + Γ + m m a ∂κ w + − u
3.5. Computation of acoustic pressure We can compute the acoustic pressure from either (31) or (32). In order to check the accuracy of our computations we shall use both equations. Figs. 7–11 show the eigenvalues at the frequencies 50 and 100 kHz. The parameters were wl = −1, bl = 0.01, br = 0.02, wr = 1. The porosity of the bone was 0.72 and the receiver was at x = 0.03, y = 0. As frequency increases the eigenvalues migrate toward κ = 1. The approximate contribution to pressure of an individual eigenvalue is indicated by its color. The decibel scale used was pdB = 20 log(10| p|). As can be seen, as frequency increases, the eigenvalues lying away from Re κ ≤ 1 make less of a contribution to the total pressure as do those lying on the axis out toward Re κ = −∞. Clearly regions to the left of Re κ = −10 would have to be included to accurately compute pressure at the frequencies shown. As Fig. 12 shows, the eigenvalues near the right
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
15
Fig. 8. Eigenvalues near Re κ ≤ 1 at 50 kHz for a receiver positioned at x = 0.03, y = 0 m.
Fig. 9. Eigenvalues at 100 kHz for a receiver positioned at x = 0.03, y = 0 m.
most part of Re κ ≤ 1 increasingly dominate the computation of pressure as the receiver position is moved away from the x-axis. The eigenvalues shown in Figs. 7–12 were computed using a recursive scheme which started with a subdivision of an inner parabolic slit-contour into quadrilaterals and a subdivision of the region between the inner parabola and an outer one into quadrilaterals (Fig. 13). The parabolas used approximated the level curves of kHz1 (see Fig. 3). For each quadrilateral the pressure pq was calculated using (31). If the pressure was less than some minimum tolerance τm , then the quadrilateral was dropped from consideration. If pq ≥ τm an attempt was made to find an eigenvalue in the quadrilateral using the secant method. If this attempt was successful then the pressure pe due to that eigenvalue was computed using (32). If pc and pe agreed to within a tolerance τa , then the eigenvalue was accepted as the unique eigenvalue in the quadrilateral. If no eigenvalue was found or pq and pe did not agree, then the quadrilateral was subdivided into four sub-quadrilaterals and the considerations just described were repeated for each sub-quadrilateral.
1 ∆(κ) could be computed up to about 10500 by rescaling the terms in the characteristic matrix M.
ARTICLE IN PRESS 16
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 10. Eigenvalues near Re κ ≤ 1 at 100 kHz for a receiver positioned at x = 0.03, y = 0 m.
Fig. 11. Eigenvalues at 250 kHz for a receiver positioned at x = 0.03, y = 0 m.
The procedure is obviously time-consuming, requiring around 40 min at 100 kHz and more at higher frequencies.2 It was effective up to around 750 kHz. By 1 MHz it was ineffective because the integrands around even a small quadrilateral were sometimes too oscillatory to compute accurately. See Fig. 14. The tolerances used for computing the eigenvalues in Figs. 7–12 were τm = 10−9 , τa = 5 × 10−9 . As mentioned above finding the eigenvalues is computationally intensive. Clearly integration along a contour that encloses all of the eigenvalues that make a significant contribution to the pressure will be more efficient, however knowing the location of such eigenvalues is necessary to do this. Figs. 15–17 show the paths of individual eigenvalues as frequency increases. The tolerance used for detecting eigenvalues was τm = 10−11 . In the frequency range 1–91 kHz (Fig. 15) there are two singular paths, one which starts near κ = 1 and circles counterclockwise, and one which moves toward the axis Re κ ≤ 1 from the first quadrant. Figs. 16 and 17 however give no indication that subsequent eigenvalues follow these trajectories. This is important since Figs. 4–6 indicate that Γ (κ) is difficult to compute far to the right of κ = 1. Also as frequency increases the trajectories tend to be closer to the axis 2 Unless otherwise stated the timed computations were done with a mobile Intel Pentium 4 processor running at 2.4 GHz.
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
17
Fig. 12. Eigenvalues at 100 kHz for a receiver positioned at x = 0.03, y = 0.01 m.
Fig. 13. Subdivision of inner (Top) and region between inner and outer parabolas (Bottom).
Re κ ≤ 1. Thus it seems that in the 100 kHz and higher range the parabolic contours shown in the figures should include all of the strongly propagating eigenvalues if they are extended far enough to the left. Consequently we will use slit-contours of the form shown in Fig. 18 which are based on the parabolic approximation of the level curve ∆(κ) = 10500 at 100 kHz. The slit runs from κ = −δ ±δi out to the parabolic portion of the contour. Unless otherwise stated δ = 10−10 . A major difficulty in the contour integration is the integration along the slit between κ = −δ ± δi and κ = 1 ± δi. Fig. 19 shows that the magnitude of the integrand fluctuates considerably along these lines with the degree of oscillation increasing with increasing frequency. To the right of κ = 1 the magnitude of the integrand decays. The numerical integrator used was MATLAB 6.5’s quad, which uses an adaptive Simpson quadrature scheme. The spike in the integrands near κ = 0 (Fig. 20) seemed to cause quad no difficulty. In the regions where the integrand was oscillatory (Figs. 21 and 22) quad would issue a warning that the result might be unreliable for very oscillatory integrands, but the absence of such a warning did not necessarily mean that the results were accurate to the specified tolerance. The sequence of computations shown here illustrates this:
ARTICLE IN PRESS 18
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 14. Integrands (real red, imag green) for the quadralateral 7.9436e−001−2.0027e−002i → 8.6311e−001−1.6589e−002i → 8.6274e−001− 9.7362e−003i → 7.9399e−001−1.1754e−002i at frequency 1 MHz. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 15. Paths of eigenvalues for frequencies 1, 16, 31, . . . , 91 kHz. The color sequence is black (1 kHz), red, magenta, yellow, green, cyan, blue (91 kHz). The parabolas (inner to outer) shown are parabolic approximations to the level curves ∆(κ) = 10−50 , 1, 1050 , . . . , 10500 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Frequency: 100 kHz, Rel. Tol. = 5.000000e−005, Int. Tol. = 1.466741e−006 |plast| = 2.933482e−002, time = 4.13
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
19
Fig. 16. Paths of eigenvalues for frequencies 200, 210, . . . , 260 kHz. The color sequence is black (200 kHz), red, magenta, yellow, green, cyan, blue (260 kHz). The parabolas (inner to outer) shown are parabolic approximations to the level curves ∆(κ) = 10−50 , 1, 1050 , . . . , 10500 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 17. Paths of eigenvalues for frequencies 470, 480, . . . , 530 kHz. The color sequence is black (470 kHz), red, magenta, yellow, green, cyan, blue (530 kHz). The parabolas (inner to outer) shown are parabolic approximations to the level curves ∆(κ) = 10−50 , 1, 1050 , . . . , 10500 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
N = 1, | p| = 2.929160e−002, rel. err. = 3.028707e−002, time = 4.90 N = 2, | p| = 3.213897e−002, rel. err. = 9.410061e−002, time = 5.36 N = 4, | p| = 2.955957e−002, rel. err. = 1.010306e−001, time = 6.34 N = 8, | p| = 2.913781e−002, rel. err. = 1.531679e−002, time = 7.50 N = 16, | p| = 2.911204e−002, rel. err. = 9.478357e−004, time = 8.65 N = 32, | p| = 2.914829e−002, rel. err. = 1.432803e−003, time = 10.20 N = 64, | p| = 2.915810e−002, rel. err. = 3.482070e−004, time = 12.03 N = 128, | p| = 2.915833e−002, rel. err. = 2.928985e−005, time = 13.91 N = 256, | p| = 2.915802e−002, rel. err. = 2.201273e−005, time = 17.49.
ARTICLE IN PRESS 20
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 18. Typical contour used for computing pressure.
Fig. 19. Magnitude of integrands along the line κ = −δ + δi to κ = 1 + δi.
In these computations rel. err. is relative to the previous value of N . The interval (−δ, 1) was subdivided into N + 2 intervals (−δ, pl ), (κn , κn+1 ) where κn = pl + (1 − pl ) log n/ log(N + 1), n = 1, . . . , N + 1. The logarithmic subdivision was used because the degree of oscillation tended to be greater near κ = 1; pl = 0.05 was used. Since it was found that quad’s value plast for the integration using the entire interval (−δ, 1) and a tolerance of 10−6 gave about the right order of magnitude, it was used to set the tolerance for the integration over each interval: τi = |plast|τr /(N + 1), where τr was the target relative error. Despite the fact that quad issued no warnings, four digit accuracy was not attained until N = 64 for a bone specimen of porosity β = 0.72 and a receiver position (x, y) = (0.03, 0). Table 4 shows the number of subdivisions required to attain four digit precision for the integration along the branch cut from (−δ, 1) for bone specimens of porosities β = 0.72, 0.83 and 0.95 when the receiver position was (x, y) = (0.03, 0). The formulas (31) and (32) can easily be vectorized in either of the receiver position variables x and y. MATLAB’s integrator quad was modified to accept vector arguments in the y-variable and require the
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
21
Fig. 20. Spike near κ = 0. Real (red) and imaginary (green) parts of the integrand. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 21. Real (red) and imaginary (green) parts of the integrand. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
acceptance criterion to be met for all y-values. The result of doing the branch-cut integration for x = 0.03 and y = 0, 0.001, 0.002, . . . , 0.02 is shown below: Frequency: 100 kHz, Rel. Tol. = 5.000000e−005, Int. Tol. = 1.444573e−006 time = 4.37 N = 1, rel. err. = 3.018602e−002, time = 5.15 N = 2, rel. err. = 1.204851e−001, time = 5.67
ARTICLE IN PRESS 22
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 22. Real (red) and imaginary (green) parts of the integrand. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Table 4 Number of intervals required to achieve a tolerance of 5e–5 Frequency
N (β = 0.72)
Time
N (β = 0.83)
Time
N (β = 0.95)
Time
100 250 500 750
64 128 512 2048
12.03 18.02 47.44 121.13
32 256 128 512
14.27 34.42 59.10 76.12
32 32 64 256
21.56 41.27 80.34 110.82
N = 4, rel. err. = 1.034279e−001, time = 6.65 N = 8, rel. err. = 1.531365e−002, time = 7.89 N = 16, rel. err. = 9.465125e−004, time = 9.12 N = 32, rel. err. = 1.619185e−003, time = 10.88 N = 64, rel. err. = 1.932502e−004, time = 12.78 N = 128, rel. err. = 9.186522e−005, time = 14.54 N = 256, rel. err. = 2.126272e−005, time = 18.32 N = 512, rel. err. = 8.372496e−006, time = 25.49. As can be seen by comparison to the scalar computation given above, the times required to do the branch-cut integration increased only slightly. Table 5 shows the number of subdivisions required to attain four digit precision for the branch-cut computation for the vectorized integrator with x = 0.03 and y = 0, 0.001, 0.002, . . . , 0.02. Comparison with Table 4 show the results of using the vectorized integrator can be unpredictable. The number of subdivisions required at 500 and 750 kHz went down, presumably because requiring the termination criterion to be met at multiple y-values inhibited premature termination. Tables 6 and 7 compare the results of computing the contribution to pressure of all eigenvalues detected to the right of Re κ = −10 using eigenvalues and contour integration respectively. The agreement between the two is close to the target precision of four digits.
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
23
–
Table 5 x = 0.03; y = 0: .001: .02 Frequency
N (β = 0.72)
Time
N (β = 0.83)
Time
N (β = 0.95)
Time
100 250 500 750
128 64 128 512
14.54 25.41 60.44 211.38
32 256 128 128
15.35 41.22 55.00 83.59
32 64 64 64
23.05 50.09 83.32 131.79
Table 6 x = 0.03; y = 0: .001: .02 Frequency
β
# Eigs
p
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
894 896 896 2214 2232 2230 4299 4439
2.9696e−002 + 2.5700e−004i 2.8956e−002 + 1.7255e−002i 2.4116e−002 + 3.8046e−002i 7.4316e−004 − 7.6031e−003i 1.1638e−002 + 1.7442e−004i 9.1995e−003 + 2.4580e−002i 1.0134e−003 − 9.0914e−005i 4.4549e−003 + 3.1661e−004i 3.6933e−003 + 1.5300e−002i
Frequency
β
p
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
2.9697e−002 + 2.5531e−004i 2.8955e−002 + 1.7256e−002i 2.4115e−002 + 3.8046e−002i 7.4262e−004 − 7.6030e−003i 1.1638e−002 + 1.7478e−004i 9.1994e−003 + 2.4581e−002i 1.0134e−003 − 9.0860e−005i 4.4550e−003 + 3.1657e−004i 3.6936e−003 + 1.5300e−002i
24.4 34.1 53.6 39.8 60.8 118.9 81.8 132.3 223.5
250
500
Table 7 x = 0.03; y = 0: .001: .02
250
500
Table 8 f = 100 kHz, y = 0 Re κl
p
Time
−20 −40 −80 −160 −320 −640 −1280 −2560
3.0638e−002 + 3.8728e−003i 2.5556e−002 + 3.8674e−003i 2.8747e−002 + 3.4166e−003i 2.9791e−002 + 3.0984e−003i 2.8229e−002 + 3.2654e−003i 2.8450e−002 + 3.2619e−003i 2.8666e−002 + 3.2362e−003i 2.9291e−002 + 3.1504e−003i
24.8 26.0 28.5 33.0 42.3 62.4 99.9 172.7
As indicated in Figs. 9–11 eigenvalues lying along the negative real κ-axis may make a significant contribution to the total pressure when the receiver position is at y = 0. Table 8 shows that contour integration along the outermost parabolic contour in Fig. 3 back to Re κ = −2560 results in only about two digit accuracy. The situation improves
ARTICLE IN PRESS 24
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Table 9 f = 100 kHz, y = 0.001 Re κl
p
Time
−20 −40 −80 −160 −320 −640 −1280
2.9157e−002 + 3.5317e−003i 2.8671e−002 + 3.4693e−003i 2.8904e−002 + 3.4328e−003i 2.8918e−002 + 3.4249e−003i 2.8910e−002 + 3.4257e−003i 2.8911e−002 + 3.4256e−003i 2.8911e−002 + 3.4256e−003i
24.7 25.8 28.3 32.7 43.3 60.6 98.1
Table 10 x = 0.03, y = 0.001 Frequency
β
Re κ
| p|
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
−320 −320 −320 −80 −80 −80 −40 −40 −20
2.9113e−002 3.5905e−002 4.6836e−002 6.5416e−003 1.2102e−002 2.6703e−002 9.4016e−004 3.7096e−003 1.5470e−002
43.3 51.6 69.7 40.8 59.2 109.1 65.4 93.0 156.4
250
500
Table 11 x = 0.03, y = 0.001: 0.001: 0.02 Frequency
β
κ
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
−320 −320 −320 −80 −80 −80 −40 −40 −20
48.2 59.7 81.2 57.3 73.1 123.4 122.6 125.0 224.2
Frequency
β
κ
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
−320 −320 −320 −80 −80 −80 −20 −20 −20
59.0 68.3 86.1 95.4 93.8 129.4 164.7 149.7 276.1
250
500
Table 12 x = 0.0, y = 0.001: 0.001: 0.02
250
500
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
25
Table 13 x = −0.01, y = 0.001: 0.001: 0.02 Frequency
β
κ
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
−320 −320 −320 −80 −80 −40 −20 −20 −20
53.7 60.0 77.4 73.3 75.7 112.5 139.6 130.9 211.4
250
500
Fig. 23. Magnitude of pressure at x = 0.03 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
if y > 0 since the exponential terms in (31) and (32) diminish the contribution of eigenvalues lying near Re κ < 0. Table 9 shows that four digit accuracy is attainable by integrating out to Re κ = −320 when x = 0.03, y = 0.001. Tables 10 and 11 show how far to the left the parabolic contour had to be extended to obtain four digit accuracy when x = 0.03 and y = 0.001 and y = 0.001, 0.002, . . . , 0.02 for several different bone specimens and frequencies. As can be seen the times required to compute the pressure at twenty y-values increased only moderately over the time required for one y-value. Tables 12 and 13 give the corresponding results when the horizontal coordinate of the receiver was x = 0.0 and −0.01 respectively. Figs. 23–28 show the magnitude of pressure for three different bone specimens at different frequencies and receiver positions. As can be seen the magnitude varies more with the bone specimen for receiver positions to the right of the specimen. For receiver positions left of the specimen the pressures were more differentiated for larger y-values.
ARTICLE IN PRESS 26
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 24. Magnitude of pressure at x = 0.0 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 25. Magnitude of pressure at x = −0.01 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
27
Fig. 26. Magnitude of pressure at x = −0.01 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 27. Magnitude of pressure at y = 0.001 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
ARTICLE IN PRESS 28
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 28. Magnitude of pressure at y = 0.001 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
References [1] J.G. Berryman, Confirmation of biot’s theory, Appl. Phys. Lett. 37 (1980) 382–384. [2] M.A. Biot, Theory of propogation of elastic waves in a fluid-saturated porous solid. i. lower frequency range, J. Acoust. Soc. Am. 28 (2) (1956) 168–178. [3] M.A. Biot, Theory of propogation of elastic waves in a fluid-saturated porous solid. ii. higher frequency range, J. Acoust. Soc. Am. 28 (2) (1956) 179–191. [4] J.L. Buchanan, R.P. Gilbert, K. Khashanah, Determination of the parameters of cancellous bone using low frequency acoustic measurements, J. Comput. Acoust. 12 (2) (2004). [5] J.L. Buchanan, R.P. Gilbert, K. Khashanah, Recovery of the poroelastic parameters of cancellous bone using low frequency acoustic interrogation, in: A. Wirgin (Ed.), Acoustics, Mechanics, and the Related Topics of Mathematical Analysis, World Scientific, 2002, pp. 41–47. [6] S. Chaffai, F. Padilla, G. Berger, P. Languier, In vitro measurement of the frequency dependent attenuation in cancellous bone between 0.2 and 2 MHz, J. Acoust, Soc. Amer. 108 (2000) 1281–1289. [7] R. Hodgskinson, C.F. Njeh, J.D. Currey, C.M. Langton, The ability of ultrasound velocity to predict the stiffness of cancellous bone in vitro, Bone 21 (1997) 183–190. [8] A. Hosokawa, T. Otani, Ultrasonic wave propagation in bovine cancellous bone, J. Acouist. Soc. Am. 101 (1997) 558–562. [9] M.L. McKelvie, S.B. Palmer, The interaction of ultrasound with cancellous bone, Phys. Med. Biol. 10 (1991) 1331–1340. [10] E.F. Morgon, O.C. Yeh, W.C. Chang, T.M. Keaveny, Nonlinear behavior of trabecular bone at small strains, J. Biomechanical Engineering (2001) 1–9. [11] J.L. Williams, Prediction of some experimental results by biot’s theory, J. Acoust. Soc. Am. 91 (1992) 1106–1112. [12] B. Yavari, A. Bedford, Comparison of numerical calculations of two biot coefficients with analytical solutions, J. Acoust. Soc. Am. 90 (1991) 985–990.
Mathematical and Computer Modelling (
)
– www.elsevier.com/locate/mcm
Determination of the parameters of cancellous bone using high frequency acoustic measurements James L. Buchanan a , Robert P. Gilbert b,∗ a Mathematics Department, United States Naval Academy, Annapolis, MD, United States b Department of Mathematical Sciences, University of Delaware, Newark, DE, United States
Received 3 March 2006; accepted 9 May 2006
Abstract The Biot model is widely used to model poroelastic media. Several authors have tested its applicability to cancellous bone, but to do so requires a priori estimation of the parameters of the Biot model, which is an uncertain and expensive endeavor. A method of computing acoustic pressure in the low 100 kHz range is developed. c 2006 Elsevier Ltd. All rights reserved.
Keywords: Osteoporosis; Cancellous bone; Poroelastic media; Biot model; Inverse problem; Parallel processing; Simplex method
1. Introduction Cancellous bone is a two component material consisting of a calcified bone matrix with interstial fatty marrow. Faced with the clinical need to understand better the effects of aging and disease on the elastic and strength properties of trabecular bone, researchers have begun to investigate whether factors such as a architecture and tissue quality can account for 20%–40% of unexplained variance in these mechanical properties [10]. Quantitative ultrasound QUS attempts to interogate the bone structure and determine whether the subject has osteoporosis; however, QUS is still in its infancy, Chaffai et al. [6]. It relies on assessment of two fundamental parameters: broadband ultrasound attenuation BUA and speed of sound SOS. The magnitude and frequency dependence of ultrasound attenuation is a complicated function of the composition and the micro-architecture of the propagating medium. However, one needs to note a lack of theoretical model that describes how the the attenuation is influenced by properties of bone. Hodgskinson et al. [7] studied the ability of ultrasound velocity to predict the Young’s modulus of elasticity of cancellous bone. QUS measurements have been successfully used recently to predict fracture risk in osteoporosis. This technique is now incorporated into commercial devices enabling attenuation to be measured. Hence mathematical models of poroelastic media have been applied to the study of bone architecture. Indeed, McKelvie and Palmer [9], Williams [11], and Hosokawa and Otani [8] discuss the application of the Biot model for a poroelastic medium
∗ Corresponding author. Tel.: +1 302 831 2315; fax: +1 302 831 4511.
E-mail address: [email protected] (R.P. Gilbert). c 2006 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2006.05.007
ARTICLE IN PRESS 2
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
to cancellous bone. In the experiments of McKelvie and Palmer and Hosokawa and Otani a small specimen of cancellous bone was placed in a tank of water between an acoustic source and receiver. An attempt was made to align the trabecular direction of the specimen along the source–receiver axis so that the isotropic Biot model would be applicable. In order to compare the predictions of the Biot model the experimenters needed estimates of the model’s parameters. These values were taken from the literature, obtained by experimental measurement, or in the case of permeability characterized as estimates without elaboration. In view of the acquisition expense and uncertainty of such estimates, it is worthwhile to explore whether these parameters can be ascertained from the measurements of the acoustic field generated by such experiments. In [5] and [4] Buchanan, Gilbert and Khashanah investigated the extent to which the most important parameters of the Biot model could be recovered by acoustic interrogation in a numerical experiment which simulated in two dimensions the physical experiment of Hosokawa and Otani. The finite element method was used to simulate both the target pressure data. This had the virtue that boundary conditions could be imposed on all four sides of the water tank used in the experiment, but the drawback that the long computation times necessitated the use of frequencies around 10 kHz and unrealistically small tank sizes. In this article we use a different mathematical technique, contour integration of the Green’s function, for the simulation. This has the virtue of permitting use of frequencies in the low 100 kHz range, which is closer to the frequencies used in the experiments and of more realistic tank dimensions. Its drawback is that the tank and bone specimen must be taken to be infinite in one direction. 2. The Biot model applied to cancellous bone The Biot model treats a poroelastic medium as an elastic frame with interstial pore fluid. Cancellous bone is anisotropic, however, as pointed out by Williams, if the acoustic waves passing through it travel in the trabecular direction an isotropic model may be acceptable. We will simulate a two dimensional version of the experiments described in McKelvie and Palmer and Hosokawa and Otani. The motion of the frame and fluid within the bone are tracked by position vectors u = [u, v] and U = [U, V ]. The constitutive equations used by Biot are those of a linear elastic material with terms added to account for the interaction of the frame and interstial fluid σx x = 2µex x + λe + Q, σ yy = 2µe yy + λe + Q, σx y = µex y ,
(1)
σ yx = µe yx ,
σ = Qe + R where the solid and fluid dilatations are given by e =∇·u=
∂v ∂u + , ∂x ∂y
=∇·U =
∂U ∂V + . ∂x ∂y
∂u ∂v + , ∂y ∂x
e yy =
(2)
The stress–strain relations are ex x =
∂u , ∂x
ex y = e yx =
∂v . ∂y
(3)
The parameter µ, the complex frame shear modulus can be measured. The other parameters λ, R and Q occurring in the constitutive equations are calculated from the measured or estimated values of the parameters given in Table 1 using the formulas 2 (K r − K b )2 − 2β K r (K r − K b ) + β 2 K r2 λ = Kb − µ + 3 D − Kb 2 2 β Kr R= D − Kb β K r ((1 − β)K r − K b ) Q= D − Kb
(4)
where D = K r (1 + β(K r /K f − 1)).
(5)
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
3
–
Table 1 Parameters in the Biot model Symbol
Parameter
ρf ρr Kb µ Kf Kr β η k α a
Density of the pore fluid Density of frame material Complex frame bulk modulus Complex frame shear modulus Fluid bulk modulus Frame material bulk modulus Porosity Viscosity of pore fluid Permeability Structure constant Pore size parameter
The bulk and shear moduli K b and µ are often given imaginary parts to account for frame inelasticity. Eqs. (1)–(3) and an argument based upon Lagrangian dynamics are shown in [2] to lead to the following equations of motion for the displacements u, U and dilatations e, ∂ ∂2 (ρ11 u + ρ12 U) + b (u − U) ∂t ∂t 2 (6) ∂2 ∂ ∇[Qe + R] = 2 (ρ12 u + ρ22 U) − b (u − U). ∂t ∂t Here ρ11 and ρ22 are density parameters for the solid and fluid, ρ12 is a density coupling parameter, and b is a dissipation parameter. These are calculated from the inputs of Table 1 using the formulas µ∇ 2 u + ∇[(λ + µ)e + Q] =
ρ11 = (1 − β)ρr − β(ρ f − mβ) ρ12 = β(ρ f − mβ) ρ22 = mβ 2 p F a ωρ f /η β 2 η b= k where αρ f m= β and the multiplicative factor F(ζ ), which was introduced in [3] to correct for the invalidity of the assumption of Poiseuille flow at high frequencies, is given by F(ζ ) =
ζ T (ζ ) 1 4 1 − 2T (ζ )/iζ
(7)
where T is defined in terms of Kelvin functions T (ζ ) =
ber0 (ζ ) + ibei0 (ζ ) . ber(ζ ) + ibei(ζ )
For time-harmonic vibrations u(x, y, t) = u(x, y)e−iωt , U(x, y, t) = U(x, y)e−iωt the equations become µ∇ 2 u + ∇[(λ + µ)e + Q] = −ω2 (ρ11 u + ρ12 U) − ib(u − U)
(8)
∇[Qe + R] = −ω (ρ12 u + ρ22 U) + ib(u − U).
(9)
2
Taking the divergence of these equations yields the system ∇ 2 ((λ + 2µ)e + Q) + p11 e + p12 = 0 ∇ 2 (Qe + R) + p12 e + p22 = 0
(10)
ARTICLE IN PRESS 4
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Table 2 Estimated values of some Biot parameters at different porosities taken from McKelvie and Palmer or Hosokawa and Otani β
k
a
0.72 0.75 0.81 0.83 0.95
5 × 10−9
4.71 × 10−4
7 × 10−9 2 × 10−8 3 × 10−8 5 × 10−7
8.00 × 10−4 1.20 × 10−3 1.35 × 10−3 2.20 × 10−3
α
Re K b
Re µ
1.10 1.08 1.06 1.05 1.01
3.18 × 109
1.30 × 109 1.10 × 109 7.38 × 108 6.27 × 108 1.05 × 108
2.69 × 109 1.80 × 109 1.55 × 109 2.57 × 108
where p11 := ω2 ρ11 + iωb,
p12 := ω2 ρ12 − iωb,
p22 := ω2 ρ22 + iωb.
(11)
The article of McKelvie and Palmer contains estimates of the Biot parameters of cancellous bone in the human os calcis (heel bone) for the normal (β = 0.72) and severely osteoporotic (β = 0.95) cases while the article of Hosokawa and Otani has estimates for bovine femoral bone for porosities of β = 0.75, 0.81 and 0.83. Table 2 contains estimates of these six Biot parameters for five bone specimens. In obtaining them we followed the estimation procedures used by Williams, McKelvie and Palmer, and Hosokawa and Otani. In generating test problems for a parameter recovery algorithm an estimate of the range of values a parameter might take is needed. Here is a discussion of how the values of the Biot parameters in Table 2 were calculated and our estimate of the range of values for each parameter: • The real parts of K b and µ were calculated using the formulas of Williams E Vn 3(1 − 2ν) f E Vn Re µ = 2(1 + ν) f
Re K b =
•
•
•
•
(12)
used by Hosokawa and Otani. Here V f = 1 − β is the bone volume fraction. Theoretically n = 1 for waves travelling in the trabecular direction and is between 2 and 3 for transverse waves, however there is enough randomness in the trabecular direction in cancellous bone that authors have empirically adjusted the exponent to agree better with experiment. Williams arrived at a value of n = 1.23 based on comparing the Biot predictions for Type I compressional and shear wave velocity assuming the form (12) to the measured speeds obtained from experiments conducted on samples taken from bovine tibia. Hosokawa and Otani found that n = 1.46 agreed well with their data from experiments on bone specimens from bovine femora. We followed Hosokawa and Otani who used the values E = 2.2 × 1010 , ν = 0.32 for the Young’s modulus and Poisson ratio of solid bone. For our parameter range we assumed that an attempt was made to align the source–receiver axis with the trabecular direction and took 1.1 ≤ n ≤ 1.6, an interval which symmetrically encompasses Hosokawa and Otani’s and William’s values of n. The imaginary parts of K b and µ were calculated using a log decrement `: Im K b∗ = `Re K b∗ /π, Im µ = `Re µ∗ /π with a value ` = 0.1 which is typical of that used in underwater acoustics. There appears to be little sensitivity to these parameters, however. The structure factor was calculated using the formula of Berryman [1] α = 1 − s(1 − 1/β). Williams and Hosokawa and Otani used s = 0.25. A value s = 0.227 was suggested by the finite element analysis of Yavari and Bedford [12]. We will calculate α assuming 0.2 ≤ s ≤ 0.3. The pore size parameter was estimated by McKelvie and Palmer using electron microscopy and by Hosokawa and Otani using x-ray examination. Fig. 1 shows that the measurements suggest a linear relationship between pore size and porosity. We took the uncertainty in this parameter to be the regression line value ±0.0002. Permeability is a difficult parameter to estimate. Fig. 2 shows the values for the five specimens of McKelvie and Palmer and Hosokawa and Otani. The estimates indicate that permeability is approximately a log-linear function of porosity. However McKelvie and Palmer characterize their values without elaboration as “estimates” Hosokawa and Otani state that their estimates are based upon those of McKelvie and Palmer. Hence the apparent log-linear
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
5
Fig. 1. Values of the pore size parameter measured by McKelvie and Palmer and Hosokawa and Otani. Solid line: Regression line for the measured values. Dashed lines: uncertainty interval assigned to this parameter.
Fig. 2. Values of the pore size parameter estimated by McKelvie and Palmer and Hosokawa and Otani. Solid line: Regression line for the estimated values. Solid curve: Estimate of permeabilty given by the Kozeny–Carmen equation. Dashed lines: uncertainty interval assigned to this parameter.
relation should be regarded circumspectly. Indeed according the Kozeny–Carmen formula βa 2 , (13) 4K where K ≈ 5 is an empirical constant, the relation is not log-linear. Fig. 2 shows that if pore size is indeed a linear function of porosity as indicated by Fig. 1, then permeability, as predicted by (13), will deviate significantly from log-linearity. We took the uncertainty in this parameter to be ±100.75 so that the Kozeny–Carmen values would be encompassed as well as the regression values. k=
Table 3 gives the values we shall use for the other Biot parameters. The value for ρ f is from McKelvie and Palmer, but the value used by Hosokawa and Otani, 930, is similar. Likewise both sets of authors used about the same value for viscosity η. The fluid bulk modulus K f is from Hosokawa and Otani. The frame material densities used by McKelvie and Palmer and Hosokawa and Otani were somewhat different. Williams reports a range of estimates for bovine
ARTICLE IN PRESS 6
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Table 3 Parameters for cancellous bone to be used for all specimens Parameter
Symbol
Value
Pore fluid density Fluid bulk modulus Pore fluid viscosity Frame material density Frame material bulk modulus
ρf Kf η ρr Kr
950 2.00 × 109 1.5 1960 2.00 × 1010
cortical bone of ρr = 1930 → 2000. We use Hosokawa and Otani’s value ρr = 1960. The frame material bulk modulus was calculated from (12) with V f = 1. 3. Mathematical simulation of the experiment 3.1. A Green’s function for acoustic pressure in water We formulate a mathematical model for calculating acoustic pressure in an experiment reminiscent of Hosokawa and Otani’s. A two-dimensional bone slab of width w is infinite in the y-direction and is in a tank of water with rigid walls at x = wl and x = wr . A sound source is located at the origin, the left edge of the bone slab is at x = bl and thus the right edge is at br = bl + w. In the water pressure p(x, y)e−iωt and horizontal displacement Uxw (x, y)e−iωt satisfy ∇ 2 p + k02 p = −δ(x)δ(y)
(14)
∂p + ρ0 ω2 Uxw = 0 ∂x
(15)
where k0 = ω/c0 , and c0 and ρ0 are the speed of sound in water and the density of water. To construct a Green’s function we start by seeking solutions of the form p(x, y, κ) = G x (x, κ)G y (y, κ)
(16)
where G 00y + k02 κG y = −δ(y)
(17) 2 G 00x + aw G x = −δ(x) √ with aw := k0 1 − κ. √ If the branch cut of the square root function is chosen so that Im κ > 0, then in order for G y → 0 as y → ±∞ we must have ( √ κy , y≥0 Aeik0 √ G y (y, κ) = −ik0 κ y Be , y < 0. To satisfy (17) we require that G y (0+, κ) − G y (0−, κ) = 0 G 0y (0+, κ) − G 0y (0−, κ) = −1. This gives √ 1 √ eik0 κ y , y ≥ 0 2ik0 κ G y (y, κ) = √ 1 − √ e−ik0 κ y , y < 0. 2ik0 κ
−
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
7
Since pressure is symmetric about y = 0, it suffices to construct the solution for y ≥ 0. We henceforth take p(x, y, κ) = −
G x (x, κ) ik0 √κ y √ e 2ik0 κ
Assuming that the walls of the tank are rigid, Uxw = 0 at x = wl , wr . From (15) we then have G 0x = 0 at x = wl , wr . This gives C cos aw (x − wl ), wl ≤ x < 0 1 sin aw x G x (x, κ) = A cos aw x + B , 0 < x < bl aw C8 cos aw (x − wr ), bl + w < x ≤ wr . The conditions G x (0+, κ) − G x (0−, κ) = 0 G 0x (0+, κ) − G 0x (0−, κ) = −1 give A = C1 cos aw wl ,
B = C1 aw sin aw wl − 1
and hence G x (x, κ) = C1 cos aw wl cos aw x + (C1 aw sin aw wl − 1) = C1 cos aw (x − wl ) −
sin aw x , aw
sin aw x aw
0 < x < bl .
3.2. Solution of the Biot equations In (10) we make the change of dependent variables τ := (λ + 2µ)e + Q,
σ := Qe + R,
(18)
the inverse transformation for which is e = a11 τ − a12 σ,
= −a12 τ + a22 σ
(19)
where a11 := R/d,
a22 := (λ + 2µ)/d
a12 := Q/d,
with d := R(λ + 2µ) − Q 2 . In terms of the new variables (10) becomes ∇ 2 τ + B11 τ + B12 σ = 0
(20)
∇ 2 σ + B21 τ + B22 σ = 0 where B11 := a11 p11 − a12 p12 ,
B12 := −a12 p11 + a22 p12 ,
B21 := a11 p12 − a12 p22 ,
B22 := −a12 p12 + a22 p22 .
From (9) the fluid displacement vector is U=−
1 (∇σ + p12 u). p22
(21)
ARTICLE IN PRESS 8
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Substituting this into (8) gives a partial differential equation for the frame displacement vector ∇ 2 u + A31 ∇τ + A32 ∇σ + B33 u = 0
(22)
where A31 :=
1 − µa11 , µ
A32 := a12 −
p12 , µp22
B33 :=
2 p11 p22 − p12 . µp22
In order to get the solutions to match across interfaces, we assume that all unknowns have the same dependence on y as p(x, y, κ) Uxw (x, y) = Uxw (x)eik0 τ (x, y) = τ (x)eik0
√ κy
√
κy , √ ik0 κ y
u x (x, y) = u x (x)e
σ (x, y) = σ (x)eik0 ,
√
κy
u y (x, y) = u y (x)eik0
√ κy
and so forth. Substituting this into (20) gives τ 00 (x) + aτ2 τ (x) + B12 σ (x) = 0 σ 00 (x) + B21 τ (x) + aσ2 σ (x) = 0 q q where aτ := B11 − k02 κ and aσ := B22 − k02 κ. Elimination of σ in (23) gives
(23)
τ 0000 + (aτ2 + aσ2 )τ 00 + (aτ2 aσ2 − B12 B21 )τ = 0 which has the general solution τ (x) = C2 cos m + (x − bl ) + C3
sin m − (x − bl ) sin m + (x − bl ) + C4 cos m − (x − bl ) + C5 m+ m−
(24)
where s m ± := s
aτ2 + aσ2 ±
p
(aτ2 + aσ2 )2 − 4(aτ2 aσ2 − B12 B21 ) 2
(aτ2 − aσ2 )2 + 4B12 B21 2 s p B11 + B22 − 2k02 κ ± (B11 − B22 )2 + 4B12 B21 = 2 s p B11 + B22 ± (B11 − B22 )2 + 4B12 B21 = − k02 κ. 2 =
aτ2 + aσ2 ±
p
From (24) and (23) we have σ (x) = −
1 (τ 00 (x) + aτ2 τ (x)). B12
The horizontal displacements can be calculated from (21) and (22) u 00x (x) + A31 τ 0 (x) + A32 σ 0 (x) + au2 u x (x) = 0 q where au := B33 − k02 κ and Ux (x) = −
1 (σ 0 (x) + p12 u x (x)). p22
(25)
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
9
The solution of (25) contains two new arbitrary constants which we denote by C6 and C7 . It is convenient for what ensues to put m ± and au in the same form as aw . Hence we write p √ au = k0 au0 − κ m ± = k0 m 0± − κ, where m 0± =
B11 + B22 ±
p
(B11 − B22 )2 + 4B12 B21 , 2k02
√
Since e(x)eik0 κ y = ∇ · hu x , u y i = computation for the fluid dilatation,
∂u x ∂x
+
∂u y ∂y
au0 =
= u 0x (x)eik0
√
B33 . k02
κy
√ √ + u y (x)ik0 κeik0 κ y we have, after a similar
1 √ (e(x) − u 0x (x)) ik0 κ 1 U y (x) = √ ((x) − Ux0 (x)). ik0 κ u y (x) =
The frame normal and tangential stress in the bone specimen are given by σx x (z) = λe(x) + 2µex x (x) + Q(x) = λe(x) + 2µu 0x (x) + Q(x) ∂u y ∂u x σx y (x) = µ + ∂y ∂x √ = µ(ik0 κu x (x) + u 0y (x)). Finally from (14) 2 we have Uxw (x) = −
1 G 0 (x, κ). ρ0 ω 2 x
3.3. Interface conditions At an interface x = xd between the water and the bone slab the following conditions are required Uxw (xd −) = βUx (xd +) + (1 − β)u x (xd +) G x (xd −) = σx x (xd +) + σ (xd +) G x (xd −) = σ (xd +)/β
(26)
0 = σx y (xd +). These conditions derive from continuity of flux, aggregate normal pressure, pore fluid pressure and tangential frame stress respectively. The conditions (26) applied at x = bl and x = bl + w lead to a system of eight linear equations in the unknowns C1 , . . . , C8 C1 .. M . = v. (27) C8 The general solution has the form p(x, y) =
k02 2πi
I Cκ
p(x, y, κ)dκ
(28)
ARTICLE IN PRESS 10
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
where Cκ is a suitably chosen contour in the κ-plane. Substituting this into (14) gives ∇ 2 p + k02 p =
k02 2πi
I Cκ
k4 + 0 2π i
I
k2 =− 0 2π i
I
k02 2π i
I
−
= −δ(x)
G 00x (x, κ)G y (y, κ)dκ +
Cκ
Cκ
k02 2πi
I Cκ
G x (x, κ)G 00y (y, κ)dκ
G x (x, κ)G y (y, κ)dκ
[k02 κG x (x, κ) + δ(x)]G y (y, κ)dκ G x (x, κ)[k02 (1 − κ)G y (y, κ) + δ(y)]dκ +
Cκ I 2 k0
2π i
k2 G y (y, κ)dκ − δ(y) 0 2π i Cκ
I Cκ
k04 2πi
I Cκ
G x (x, κ)G y (y, κ)dκ
G x (x, κ)dκ.
Since k02 2π i
I Cκ
G x (x, κ)dκ = δ(x),
for a positively oriented contour Cκ which encloses the singularities of G x , p is a solution of (14) providing Cκ is chosen to exclude the singularities of G y and enclose those of G x . Substituting the results of the previous section into (28) gives I k02 G x (x, κ) ik0 √κ y p(x, y) = − e dκ (29) √ 2π i Cκ 2ik0 κ I k0 G x (x, κ) ik0 √κ y = e dκ √ 4π Cκ κ k0 I C1 (κ) cos aw (x − wl ) ik0 √κ y e dκ, wl < x < 0 √ 4π Cκ κ k I C1 (κ) cos aw (x − wl ) − sinaaww x ik √κ y 0 = (30) e 0 dκ, 0 < x < bl √ 4π Cκ κ I C8 (κ) cos aw (x − wr ) ik0 √κ y k0 e dκ, br < x < wr . √ 4π Cκ κ √ The only singularity of G y is the κ, thus it can be excluded by taking Cκ to be a cut enclosing the positive κ-axis. Thus the first two entries in (30) are the same since sinaaww x is an entire function of κ. Observing that all of the functions used in the matrix M in (27) are analytic except possibly on the positive κ-axis, it follows that the only singularities of G x are the solutions {κn } of ∆ := det M = 0, which are then the eigenvalues of the problem. According to Cramer’s ∆ (κ) rule C1,8 = ∆1,8(κ) and thus from the residue theorem I k0 ∆1 (κ) cos aw (x − wl ) ik0 √κ y e dκ, wl < x < bl √ 4π ∆(κ) κ Cκ I p(x, y) = k0 ∆8 (κ) cos aw (x − wr ) ik0 √κ y e dκ, br < x < wr . √ 4π Cκ ∆(κ) κ √ ik0 X ∆1 (κn ) cos aw (κn )(x − wl )eik0 κn y , wl < x < bl √ d∆ 2 κn dκ (κn ) n = X √ ik0 ∆8 (κn ) cos aw (κn )(x − wr )eik0 κn y , br < x < wr . √ d∆ 2 κn (κn ) n dκ
(31)
(32)
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
11
Fig. 3. Level curves for ∆(κ) at 100 kHz.
3.4. Numerical implementation of the pressure equation The first of the four interface conditions involves displacement, whereas the remaining three involve pressure. Typically the two quantities differ by more than 10 orders of magnitude. This can lead to inaccuracies in calculating the characteristic determinant ∆(κ), but the problem is easily resolved by scaling the pressure equations by some characteristic pressure. We used ρ0 ω2 , the compressional Lam´e coefficient for water. We are interested in determining the pressure field in the ultrasound regime, say for frequencies of 100 kHz and higher. For our initial explorations we will use wl = −1, bl = 0.01, br = 0.02, wr = 1.02 m. For the Biot parameters for bone we will use those for cancellous bone of porosity 0.72. Fig. 3, which shows the level curves Cn : ∆(κ) = 10n , n = 0, 50, 100, . . . , 300 for 100 kHz with the stated parameters, indicates a more formidable difficulty in implementing the Eqs. (31) and (32). The magnitude of ∆(κ) ascends rapidly as κ moves away from the segment Re κ ≤ 1. The curve C300 is nearly the largest that can be computed without rescaling since the largest representable real number in IEEE floating point is Ω ≈ 10308 . The problem is the presence of terms of the form sin ad and cos ad where a is one of the coefficients aw , m ± or au and d is a distance. For instance the outer two columns of the characteristic matrix are of the form M sin a d · · · 0 11
w l
M21 cos aw dl M31 cos aw dl 0 M= 0 0 0 0
···
0 0 0 M58 sin aw dr M68 cos aw dr M78 cos aw dr 0
where dl = bl − wl , dr = wr − br and the coefficients M jk depend upon κ but do not grow exponentially as κ moves away from Re κ ≤ 1. Likewise the interior columns of M contain exponentially growing terms 0 M ··· M 0 13
M2...7
M22 M32 0 = M52 sin m + w M62 cos m + w M cos m w + 72 M82 sin m + w
0 0 M43 M53 cos m + w M63 sin m + w M73 sin m + w M83 cos m + w
16
···
0 0 M46 M56 cos au w M66 sin au w 0 M86 cos au w
M27 0 0 . M57 sin au w M67 cos au w 0 M87 sin au w
ARTICLE IN PRESS 12
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 4. Level curves for Γ (κ) at 100 kHz.
Fig. 5. Level curves for Γ (κ) at 500 kHz.
√ To deal with such terms we let a := k0 a0 − κ = ξ + ηi, η ≥ 0 and rewrite them as eiξ d e−2ηd − e−iξ d 2i iξ d −2ηd + e−iξ d ηd e e cos ad = e 2 sin ad = eηd
(33)
and factor the leading exponential out of the determinant. This gives ∆ = eηw dl eηw dr e2η+ w e2η− w e2ηu w Γ where aw = ξw + ηw i, m ± = ξ± + η± i, au = ξu + ηu i. Figs. 4–6 show the level curves of Γ at frequencies 100, 500 and 1000 kHz. The value of |Γ | decreases as κ moves away from Re κ ≤ 1. The stability of the computation also decreases in some regions of the κ-plane, more noticeably at the higher frequencies.
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
13
Fig. 6. Level curves for Γ (κ) at 1000 kHz.
The right hand side of (27) is cos aw bl − ρ0 ω2 sin a b w l ρ0 ω2 aw sin aw bl . v= ρ0 ω2 aw 0 .. . 0 iξ d −2η d w rx e w r x +e−iξw dr x Thus we have ∆8 = eηbl eηw dl e2η+ w e2η− w e2ηu w Γ8 . Finally cos aw (x − wr ) = eηw dr x e with 2 dr x = wr − x. Thus the second integrand in (31) is iξw dr x e−2ηw dr x + e−iξw dr x Γ ∆8 (κ) cos aw (x − wr ) ik0 √κ y 8 ik0 √κ y ηw (bl +dr x −dr ) e =e . (34) e e √ √ ∆(κ) κ 2Γ κ Since bl + dr x − dr = bl + br − x, the leading exponential factor will grow for br < x < bl + br and decay for bl + br < x < wr , but what the integrand does of course depends upon the factor Γ8 /Γ also. A similar analysis for wl < x < bl gives iξw dlx e−2ηw dlx + e−iξw dlx Γ ∆1 (κ) cos aw (x − wr ) ik0 √κ y 1 ik0 √κ y ηw (bl +dlx −dl ) e =e e . (35) e √ √ 2Γ κ ∆(κ) κ We now consider the residue expansion (32). Introducing the operators ∂ 1 ∂ ∂ ∂ 1 ∂ ∂ := −i , := +i , κ = κ1 + iκ2 , ∂κ 2 ∂κ1 ∂κ2 ∂κ 2 ∂κ1 ∂κ2 √ we have for a term of the form a = k0 a0 − κ = ξ + iη ∂ √ k0 ∂ξ ∂η k 0 a0 − κ = − √ = +i . ∂κ ∂κ ∂κ 2 a0 − κ
ARTICLE IN PRESS 14
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 7. Eigenvalues at 50 kHz for a receiver positioned at x = 0.03, y = 0 m.
Since a is an analytic function of κ, 0=
∂ξ ∂η dξ ∂η ∂a = +i ⇒ =i ∂κ ∂κ ∂κ dκ ∂κ
which gives k2i ∂η = 0 . ∂κ 4a This gives d∆ = dκ
"
k02 i 4
# dl + dr 1 1 1 ∂Γ ηw dl ηw dr 2η+ w 2η− w 2ηu w Γ+ e e e e e . + 2w + + aw m+ m− au ∂κ
The term ddκΓ is complicated and will be computed numerically. Thus for br < x < wr we have for the terms in (32) √ √ eηw (bl +dr x −dr ) Γ8 eiξw dr x e−2ηw dr x + e−iξw dr x eik0 κn y ∆8 (κn ) ik0 κn y i = √ h cos aw (κn )(x − wr )e √ d∆ 1 1 1 ∂Γ r κn dκ (κn ) 2 κn k40 i dla+d + 2w + + Γ + m+ m− au ∂κ w whereas for wl < x < bl √
∆1 (κn ) ∆ κn ddκ (κn )
cos aw (κn )(x − wl )e
√ ik0 κn y
√ eηw (bl +dlx −dl ) Γ1 eiξw dlx e−2ηw dlx + e−iξw dlx eik0 κn y i. = √ h 1 1 1 ∂Γ r 2 κn k40 i dla+d + 2w + + Γ + m m a ∂κ w + − u
3.5. Computation of acoustic pressure We can compute the acoustic pressure from either (31) or (32). In order to check the accuracy of our computations we shall use both equations. Figs. 7–11 show the eigenvalues at the frequencies 50 and 100 kHz. The parameters were wl = −1, bl = 0.01, br = 0.02, wr = 1. The porosity of the bone was 0.72 and the receiver was at x = 0.03, y = 0. As frequency increases the eigenvalues migrate toward κ = 1. The approximate contribution to pressure of an individual eigenvalue is indicated by its color. The decibel scale used was pdB = 20 log(10| p|). As can be seen, as frequency increases, the eigenvalues lying away from Re κ ≤ 1 make less of a contribution to the total pressure as do those lying on the axis out toward Re κ = −∞. Clearly regions to the left of Re κ = −10 would have to be included to accurately compute pressure at the frequencies shown. As Fig. 12 shows, the eigenvalues near the right
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
15
Fig. 8. Eigenvalues near Re κ ≤ 1 at 50 kHz for a receiver positioned at x = 0.03, y = 0 m.
Fig. 9. Eigenvalues at 100 kHz for a receiver positioned at x = 0.03, y = 0 m.
most part of Re κ ≤ 1 increasingly dominate the computation of pressure as the receiver position is moved away from the x-axis. The eigenvalues shown in Figs. 7–12 were computed using a recursive scheme which started with a subdivision of an inner parabolic slit-contour into quadrilaterals and a subdivision of the region between the inner parabola and an outer one into quadrilaterals (Fig. 13). The parabolas used approximated the level curves of kHz1 (see Fig. 3). For each quadrilateral the pressure pq was calculated using (31). If the pressure was less than some minimum tolerance τm , then the quadrilateral was dropped from consideration. If pq ≥ τm an attempt was made to find an eigenvalue in the quadrilateral using the secant method. If this attempt was successful then the pressure pe due to that eigenvalue was computed using (32). If pc and pe agreed to within a tolerance τa , then the eigenvalue was accepted as the unique eigenvalue in the quadrilateral. If no eigenvalue was found or pq and pe did not agree, then the quadrilateral was subdivided into four sub-quadrilaterals and the considerations just described were repeated for each sub-quadrilateral.
1 ∆(κ) could be computed up to about 10500 by rescaling the terms in the characteristic matrix M.
ARTICLE IN PRESS 16
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 10. Eigenvalues near Re κ ≤ 1 at 100 kHz for a receiver positioned at x = 0.03, y = 0 m.
Fig. 11. Eigenvalues at 250 kHz for a receiver positioned at x = 0.03, y = 0 m.
The procedure is obviously time-consuming, requiring around 40 min at 100 kHz and more at higher frequencies.2 It was effective up to around 750 kHz. By 1 MHz it was ineffective because the integrands around even a small quadrilateral were sometimes too oscillatory to compute accurately. See Fig. 14. The tolerances used for computing the eigenvalues in Figs. 7–12 were τm = 10−9 , τa = 5 × 10−9 . As mentioned above finding the eigenvalues is computationally intensive. Clearly integration along a contour that encloses all of the eigenvalues that make a significant contribution to the pressure will be more efficient, however knowing the location of such eigenvalues is necessary to do this. Figs. 15–17 show the paths of individual eigenvalues as frequency increases. The tolerance used for detecting eigenvalues was τm = 10−11 . In the frequency range 1–91 kHz (Fig. 15) there are two singular paths, one which starts near κ = 1 and circles counterclockwise, and one which moves toward the axis Re κ ≤ 1 from the first quadrant. Figs. 16 and 17 however give no indication that subsequent eigenvalues follow these trajectories. This is important since Figs. 4–6 indicate that Γ (κ) is difficult to compute far to the right of κ = 1. Also as frequency increases the trajectories tend to be closer to the axis 2 Unless otherwise stated the timed computations were done with a mobile Intel Pentium 4 processor running at 2.4 GHz.
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
17
Fig. 12. Eigenvalues at 100 kHz for a receiver positioned at x = 0.03, y = 0.01 m.
Fig. 13. Subdivision of inner (Top) and region between inner and outer parabolas (Bottom).
Re κ ≤ 1. Thus it seems that in the 100 kHz and higher range the parabolic contours shown in the figures should include all of the strongly propagating eigenvalues if they are extended far enough to the left. Consequently we will use slit-contours of the form shown in Fig. 18 which are based on the parabolic approximation of the level curve ∆(κ) = 10500 at 100 kHz. The slit runs from κ = −δ ±δi out to the parabolic portion of the contour. Unless otherwise stated δ = 10−10 . A major difficulty in the contour integration is the integration along the slit between κ = −δ ± δi and κ = 1 ± δi. Fig. 19 shows that the magnitude of the integrand fluctuates considerably along these lines with the degree of oscillation increasing with increasing frequency. To the right of κ = 1 the magnitude of the integrand decays. The numerical integrator used was MATLAB 6.5’s quad, which uses an adaptive Simpson quadrature scheme. The spike in the integrands near κ = 0 (Fig. 20) seemed to cause quad no difficulty. In the regions where the integrand was oscillatory (Figs. 21 and 22) quad would issue a warning that the result might be unreliable for very oscillatory integrands, but the absence of such a warning did not necessarily mean that the results were accurate to the specified tolerance. The sequence of computations shown here illustrates this:
ARTICLE IN PRESS 18
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 14. Integrands (real red, imag green) for the quadralateral 7.9436e−001−2.0027e−002i → 8.6311e−001−1.6589e−002i → 8.6274e−001− 9.7362e−003i → 7.9399e−001−1.1754e−002i at frequency 1 MHz. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 15. Paths of eigenvalues for frequencies 1, 16, 31, . . . , 91 kHz. The color sequence is black (1 kHz), red, magenta, yellow, green, cyan, blue (91 kHz). The parabolas (inner to outer) shown are parabolic approximations to the level curves ∆(κ) = 10−50 , 1, 1050 , . . . , 10500 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Frequency: 100 kHz, Rel. Tol. = 5.000000e−005, Int. Tol. = 1.466741e−006 |plast| = 2.933482e−002, time = 4.13
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
19
Fig. 16. Paths of eigenvalues for frequencies 200, 210, . . . , 260 kHz. The color sequence is black (200 kHz), red, magenta, yellow, green, cyan, blue (260 kHz). The parabolas (inner to outer) shown are parabolic approximations to the level curves ∆(κ) = 10−50 , 1, 1050 , . . . , 10500 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 17. Paths of eigenvalues for frequencies 470, 480, . . . , 530 kHz. The color sequence is black (470 kHz), red, magenta, yellow, green, cyan, blue (530 kHz). The parabolas (inner to outer) shown are parabolic approximations to the level curves ∆(κ) = 10−50 , 1, 1050 , . . . , 10500 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
N = 1, | p| = 2.929160e−002, rel. err. = 3.028707e−002, time = 4.90 N = 2, | p| = 3.213897e−002, rel. err. = 9.410061e−002, time = 5.36 N = 4, | p| = 2.955957e−002, rel. err. = 1.010306e−001, time = 6.34 N = 8, | p| = 2.913781e−002, rel. err. = 1.531679e−002, time = 7.50 N = 16, | p| = 2.911204e−002, rel. err. = 9.478357e−004, time = 8.65 N = 32, | p| = 2.914829e−002, rel. err. = 1.432803e−003, time = 10.20 N = 64, | p| = 2.915810e−002, rel. err. = 3.482070e−004, time = 12.03 N = 128, | p| = 2.915833e−002, rel. err. = 2.928985e−005, time = 13.91 N = 256, | p| = 2.915802e−002, rel. err. = 2.201273e−005, time = 17.49.
ARTICLE IN PRESS 20
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 18. Typical contour used for computing pressure.
Fig. 19. Magnitude of integrands along the line κ = −δ + δi to κ = 1 + δi.
In these computations rel. err. is relative to the previous value of N . The interval (−δ, 1) was subdivided into N + 2 intervals (−δ, pl ), (κn , κn+1 ) where κn = pl + (1 − pl ) log n/ log(N + 1), n = 1, . . . , N + 1. The logarithmic subdivision was used because the degree of oscillation tended to be greater near κ = 1; pl = 0.05 was used. Since it was found that quad’s value plast for the integration using the entire interval (−δ, 1) and a tolerance of 10−6 gave about the right order of magnitude, it was used to set the tolerance for the integration over each interval: τi = |plast|τr /(N + 1), where τr was the target relative error. Despite the fact that quad issued no warnings, four digit accuracy was not attained until N = 64 for a bone specimen of porosity β = 0.72 and a receiver position (x, y) = (0.03, 0). Table 4 shows the number of subdivisions required to attain four digit precision for the integration along the branch cut from (−δ, 1) for bone specimens of porosities β = 0.72, 0.83 and 0.95 when the receiver position was (x, y) = (0.03, 0). The formulas (31) and (32) can easily be vectorized in either of the receiver position variables x and y. MATLAB’s integrator quad was modified to accept vector arguments in the y-variable and require the
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
21
Fig. 20. Spike near κ = 0. Real (red) and imaginary (green) parts of the integrand. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 21. Real (red) and imaginary (green) parts of the integrand. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
acceptance criterion to be met for all y-values. The result of doing the branch-cut integration for x = 0.03 and y = 0, 0.001, 0.002, . . . , 0.02 is shown below: Frequency: 100 kHz, Rel. Tol. = 5.000000e−005, Int. Tol. = 1.444573e−006 time = 4.37 N = 1, rel. err. = 3.018602e−002, time = 5.15 N = 2, rel. err. = 1.204851e−001, time = 5.67
ARTICLE IN PRESS 22
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 22. Real (red) and imaginary (green) parts of the integrand. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Table 4 Number of intervals required to achieve a tolerance of 5e–5 Frequency
N (β = 0.72)
Time
N (β = 0.83)
Time
N (β = 0.95)
Time
100 250 500 750
64 128 512 2048
12.03 18.02 47.44 121.13
32 256 128 512
14.27 34.42 59.10 76.12
32 32 64 256
21.56 41.27 80.34 110.82
N = 4, rel. err. = 1.034279e−001, time = 6.65 N = 8, rel. err. = 1.531365e−002, time = 7.89 N = 16, rel. err. = 9.465125e−004, time = 9.12 N = 32, rel. err. = 1.619185e−003, time = 10.88 N = 64, rel. err. = 1.932502e−004, time = 12.78 N = 128, rel. err. = 9.186522e−005, time = 14.54 N = 256, rel. err. = 2.126272e−005, time = 18.32 N = 512, rel. err. = 8.372496e−006, time = 25.49. As can be seen by comparison to the scalar computation given above, the times required to do the branch-cut integration increased only slightly. Table 5 shows the number of subdivisions required to attain four digit precision for the branch-cut computation for the vectorized integrator with x = 0.03 and y = 0, 0.001, 0.002, . . . , 0.02. Comparison with Table 4 show the results of using the vectorized integrator can be unpredictable. The number of subdivisions required at 500 and 750 kHz went down, presumably because requiring the termination criterion to be met at multiple y-values inhibited premature termination. Tables 6 and 7 compare the results of computing the contribution to pressure of all eigenvalues detected to the right of Re κ = −10 using eigenvalues and contour integration respectively. The agreement between the two is close to the target precision of four digits.
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
23
–
Table 5 x = 0.03; y = 0: .001: .02 Frequency
N (β = 0.72)
Time
N (β = 0.83)
Time
N (β = 0.95)
Time
100 250 500 750
128 64 128 512
14.54 25.41 60.44 211.38
32 256 128 128
15.35 41.22 55.00 83.59
32 64 64 64
23.05 50.09 83.32 131.79
Table 6 x = 0.03; y = 0: .001: .02 Frequency
β
# Eigs
p
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
894 896 896 2214 2232 2230 4299 4439
2.9696e−002 + 2.5700e−004i 2.8956e−002 + 1.7255e−002i 2.4116e−002 + 3.8046e−002i 7.4316e−004 − 7.6031e−003i 1.1638e−002 + 1.7442e−004i 9.1995e−003 + 2.4580e−002i 1.0134e−003 − 9.0914e−005i 4.4549e−003 + 3.1661e−004i 3.6933e−003 + 1.5300e−002i
Frequency
β
p
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
2.9697e−002 + 2.5531e−004i 2.8955e−002 + 1.7256e−002i 2.4115e−002 + 3.8046e−002i 7.4262e−004 − 7.6030e−003i 1.1638e−002 + 1.7478e−004i 9.1994e−003 + 2.4581e−002i 1.0134e−003 − 9.0860e−005i 4.4550e−003 + 3.1657e−004i 3.6936e−003 + 1.5300e−002i
24.4 34.1 53.6 39.8 60.8 118.9 81.8 132.3 223.5
250
500
Table 7 x = 0.03; y = 0: .001: .02
250
500
Table 8 f = 100 kHz, y = 0 Re κl
p
Time
−20 −40 −80 −160 −320 −640 −1280 −2560
3.0638e−002 + 3.8728e−003i 2.5556e−002 + 3.8674e−003i 2.8747e−002 + 3.4166e−003i 2.9791e−002 + 3.0984e−003i 2.8229e−002 + 3.2654e−003i 2.8450e−002 + 3.2619e−003i 2.8666e−002 + 3.2362e−003i 2.9291e−002 + 3.1504e−003i
24.8 26.0 28.5 33.0 42.3 62.4 99.9 172.7
As indicated in Figs. 9–11 eigenvalues lying along the negative real κ-axis may make a significant contribution to the total pressure when the receiver position is at y = 0. Table 8 shows that contour integration along the outermost parabolic contour in Fig. 3 back to Re κ = −2560 results in only about two digit accuracy. The situation improves
ARTICLE IN PRESS 24
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Table 9 f = 100 kHz, y = 0.001 Re κl
p
Time
−20 −40 −80 −160 −320 −640 −1280
2.9157e−002 + 3.5317e−003i 2.8671e−002 + 3.4693e−003i 2.8904e−002 + 3.4328e−003i 2.8918e−002 + 3.4249e−003i 2.8910e−002 + 3.4257e−003i 2.8911e−002 + 3.4256e−003i 2.8911e−002 + 3.4256e−003i
24.7 25.8 28.3 32.7 43.3 60.6 98.1
Table 10 x = 0.03, y = 0.001 Frequency
β
Re κ
| p|
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
−320 −320 −320 −80 −80 −80 −40 −40 −20
2.9113e−002 3.5905e−002 4.6836e−002 6.5416e−003 1.2102e−002 2.6703e−002 9.4016e−004 3.7096e−003 1.5470e−002
43.3 51.6 69.7 40.8 59.2 109.1 65.4 93.0 156.4
250
500
Table 11 x = 0.03, y = 0.001: 0.001: 0.02 Frequency
β
κ
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
−320 −320 −320 −80 −80 −80 −40 −40 −20
48.2 59.7 81.2 57.3 73.1 123.4 122.6 125.0 224.2
Frequency
β
κ
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
−320 −320 −320 −80 −80 −80 −20 −20 −20
59.0 68.3 86.1 95.4 93.8 129.4 164.7 149.7 276.1
250
500
Table 12 x = 0.0, y = 0.001: 0.001: 0.02
250
500
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
25
Table 13 x = −0.01, y = 0.001: 0.001: 0.02 Frequency
β
κ
Time
100
0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95
−320 −320 −320 −80 −80 −40 −20 −20 −20
53.7 60.0 77.4 73.3 75.7 112.5 139.6 130.9 211.4
250
500
Fig. 23. Magnitude of pressure at x = 0.03 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
if y > 0 since the exponential terms in (31) and (32) diminish the contribution of eigenvalues lying near Re κ < 0. Table 9 shows that four digit accuracy is attainable by integrating out to Re κ = −320 when x = 0.03, y = 0.001. Tables 10 and 11 show how far to the left the parabolic contour had to be extended to obtain four digit accuracy when x = 0.03 and y = 0.001 and y = 0.001, 0.002, . . . , 0.02 for several different bone specimens and frequencies. As can be seen the times required to compute the pressure at twenty y-values increased only moderately over the time required for one y-value. Tables 12 and 13 give the corresponding results when the horizontal coordinate of the receiver was x = 0.0 and −0.01 respectively. Figs. 23–28 show the magnitude of pressure for three different bone specimens at different frequencies and receiver positions. As can be seen the magnitude varies more with the bone specimen for receiver positions to the right of the specimen. For receiver positions left of the specimen the pressures were more differentiated for larger y-values.
ARTICLE IN PRESS 26
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 24. Magnitude of pressure at x = 0.0 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 25. Magnitude of pressure at x = −0.01 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
ARTICLE IN PRESS J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
27
Fig. 26. Magnitude of pressure at x = −0.01 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 27. Magnitude of pressure at y = 0.001 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
ARTICLE IN PRESS 28
J.L. Buchanan, R.P. Gilbert / Mathematical and Computer Modelling (
)
–
Fig. 28. Magnitude of pressure at y = 0.001 for bone specimens of porosity β = 0.72 (blue), 0.83 (red) and 0.95 (green) at frequencies 100 (Top), 250 (Middle) and 500 kHz (Bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
References [1] J.G. Berryman, Confirmation of biot’s theory, Appl. Phys. Lett. 37 (1980) 382–384. [2] M.A. Biot, Theory of propogation of elastic waves in a fluid-saturated porous solid. i. lower frequency range, J. Acoust. Soc. Am. 28 (2) (1956) 168–178. [3] M.A. Biot, Theory of propogation of elastic waves in a fluid-saturated porous solid. ii. higher frequency range, J. Acoust. Soc. Am. 28 (2) (1956) 179–191. [4] J.L. Buchanan, R.P. Gilbert, K. Khashanah, Determination of the parameters of cancellous bone using low frequency acoustic measurements, J. Comput. Acoust. 12 (2) (2004). [5] J.L. Buchanan, R.P. Gilbert, K. Khashanah, Recovery of the poroelastic parameters of cancellous bone using low frequency acoustic interrogation, in: A. Wirgin (Ed.), Acoustics, Mechanics, and the Related Topics of Mathematical Analysis, World Scientific, 2002, pp. 41–47. [6] S. Chaffai, F. Padilla, G. Berger, P. Languier, In vitro measurement of the frequency dependent attenuation in cancellous bone between 0.2 and 2 MHz, J. Acoust, Soc. Amer. 108 (2000) 1281–1289. [7] R. Hodgskinson, C.F. Njeh, J.D. Currey, C.M. Langton, The ability of ultrasound velocity to predict the stiffness of cancellous bone in vitro, Bone 21 (1997) 183–190. [8] A. Hosokawa, T. Otani, Ultrasonic wave propagation in bovine cancellous bone, J. Acouist. Soc. Am. 101 (1997) 558–562. [9] M.L. McKelvie, S.B. Palmer, The interaction of ultrasound with cancellous bone, Phys. Med. Biol. 10 (1991) 1331–1340. [10] E.F. Morgon, O.C. Yeh, W.C. Chang, T.M. Keaveny, Nonlinear behavior of trabecular bone at small strains, J. Biomechanical Engineering (2001) 1–9. [11] J.L. Williams, Prediction of some experimental results by biot’s theory, J. Acoust. Soc. Am. 91 (1992) 1106–1112. [12] B. Yavari, A. Bedford, Comparison of numerical calculations of two biot coefficients with analytical solutions, J. Acoust. Soc. Am. 90 (1991) 985–990.
