Determination of the parameters of cancellous ... - Semantic Scholar

Determination of the parameters of cancellous bone using high frequency acoustic measurements James L. Buchanan Mathematics Department United States Naval Academy Annapolis MD Robert P. Gilbert Department of Mathematical Sciences University of Delaware Newark DE July 20, 2004 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. The question this article addresses is whether the sort of experiments described by these authors in fact can be used to determine the parameters of the Biot model. A method of computing acoustic pressure in the low 100 kHz range is developed. A parameter recovery algorithm which uses parallel processing is developed and tested. It is found that when it is assumed that the agreement between calculated and measured data is about two digits, porosity can be determined to within about 1-2% and permeability, pore size and the bulk moduli to within about 40%. 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 calci…ed bone matrix with interstial fatty marrow. Hence mathematical models of poroelastic media are applicable. McKelvie and Palmer [7], Williams [8], and Hosokawa and Otani [6] discuss the application of the Biot model for a poroelastic medium to cancellous bone. In the experiments of McKelvie and Palmer and Hosokawa 1

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 …eld 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 an numerical experiment which simulated in two dimensions the physical experiment of Hosokawa and Otani. The …nite element method was used to simulate both the target pressure data and the trial data used in the parameter recovery algorithm. 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 di¤erent 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 in…nite 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 ‡uid. 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 ‡uid 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 ‡uid xx yy xy

= 2 exx + e + Q ; = 2 eyy + e + Q ; = exy ; yx = eyx ; = Qe + R

(1)

where the solid and ‡uid dilatations are given by e=r u=

@U @V @u @v + ; =r U = + : @x @y @x @y 2

(2)

Symbol

Parameter Density of the pore ‡uid Density of frame material Complex frame bulk modulus Complex frame shear modulus Fluid bulk modulus Frame material bulk modulus Porosity Viscosity of pore ‡uid Permeability Structure constant Pore size parameter

f r

Kb Kf Kr

k a

Table 1: Parameters in the Biot model The stress-strain relations are exx =

@u @u @v @v ; exy = eyx = + ; eyy = : @x @y @x @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 = Kb

2 (Kr + 3

Kr2 D Kb Kr ((1 Q= D

2

Kb )

2 Kr (Kr D Kb

Kb ) +

2

Kr2

(4)

2

R=

) Kr Kb

Kb )

:

where D = Kr (1 + (Kr =Kf

1)):

(5)

The bulk and shear moduli Kb and are often given imaginary parts to account for frame inelasticity. Equations (1), (2) and (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 ( @t2 @2 r[Qe + R ] = 2 ( @t

r2 u + r[( + )e + Q ] =

11 u

+

12 U )

12 u

+

22 U )

@ (u @t @ b (u @t

+b

U)

(6)

U ):

Here 11 and 22 are density parameters for the solid and ‡uid, 12 is a density coupling parameter, and b is a dissipation parameter. These are calculated from 3

the inputs of Table 1 using the formulas 11 12 22

= (1 = ( =m

b=

) f

(

r

m )

f

m )

2

p F a !

2

f=

k

where

f

m=

and the multiplicative factor F ( ), which was introduced in [3] to correct for the invalidity of the assumption of Poiseuille ‡ow at high frequencies, is given by 1 T( ) F( ) = (7) 4 1 2T ( )=i where T is de…ned in terms of Kelvin functions T( ) =

ber0 ( ) + ibei0 ( ) : ber( ) + ibei( )

For time-harmonic vibrations u(x; y; t) = u(x; y)e the equations become r2 u + r[( + )e + Q ] = r[Qe + R ] =

i!t

; U (x; y; t) = U (x; y)e

!2 (

11 u

+

12 U )

ib(u

U)

(8)

2

12 u

+

22 U )

+ ib(u

U ):

(9)

! (

Taking the divergence of these equations yields the system r2 (( + 2 ) e + Q ) + p11 e + p12 = 0

(10)

2

r (Qe + R ) + p12 e + p22 = 0

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 …ve 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: 4

i!t

0:72 0:75 0:81 0:83 0:95

k 5 7 2 3 5

10 10 10 10 10

9 9 8 8 7

a 4:71 8:00 1:20 1:35 2:20

10 10 10 10 10

4

1:10 1:08 1:06 1:05 1:01

4 3 3 3

Re Kb 3:18 2:69 1:80 1:55 2:57

109 109 109 109 108

Re 1:30 1:10 7:38 6:27 1:05

109 109 108 108 108

Table 2: Estimated values of some Biot parameters at di¤erent porosities taken from McKelvie and Palmer or Hosokawa and Otani. The real parts of Kb and

were calculated using the formulas of Williams E Vn 3(1 2 ) f E = Vn 2(1 + ) f

Re Kb = Re

(12)

used by Hosokawa and Otani. Here Vf = 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 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 Kb and were calculated using a log decrement `: Im Kb = ` Re Kb = ; 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 …nite element analysis of Yavari and Bedford [9]. 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. Figure 1 shows that the measurements suggest a linear relationship 5

Figure 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.

between pore size and porosity. We took the uncertainty in this parameter to be the regression line value 0:0002 Permeability is a di¢ cult parameter to estimate. Figure 2 shows the values for the …ve specimens of McKelvie and Palmer and Hosokawa and Otani. The estimates indicate that permeability is approximately a loglinear function of porosity. However McKelvie and Palmer characterize their values without elaboration as ”estimates” and Hosokawa and Otani state that their estimates are based upon those of McKelvie and Palmer. Hence the apparent log-linear relation should be regarded circumspectly. Indeed according the Kozeny-Carmen formula k=

a2 ; 4K

(13)

where K 5 is an empirical constant, the relation is not log-linear. Figure ?? shows that if pore size is indeed a linear function of porosity as indicated by Figure ??, then permeability, as predicted by (13), will deviate signi…cantly 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. 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 6

Figure 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.

Parameter Pore ‡uid density Fluid bulk modulus Pore ‡uid viscosity Frame material density Frame material bulk modulus

Symbol f

Kf r

Kr

Value 950 2:00 109 1.5 1960 2:00 1010

Table 3: Parameters for cancellous bone to be used for all specimens.

7

viscosity . The ‡uid bulk modulus Kf is from Hosokawa and Otani. The frame material densities used by McKelvie and Palmer and Hosokawa and Otani were somewhat di¤erent. Williams reports a range of estimates for bovine 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 Vf = 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 in…nite 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 r2 p + k02 p =

@p + @x

0!

2

(x) (y)

(14)

Uxw = 0

(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; ) = Gx (x; )Gy (y; ) (16) where G00y + k02 Gy = G00x

+

a2w Gx

=

(y)

(17)

(x)

p . with aw := k0 1 p If the branch cut of the square root function is chosen so that Im > 0, then in order for Gy ! 0 as y ! 1 we must have p

Gy (y; ) =

Aeik0 p y ; y 0 : Be ik0 y ; y < 0

To satisfy (17) we require that Gy (0+; ) G0y (0+; ) This gives Gy (y; ) =

(

Gy (0 ; ) = 0 G0y (0 ; ) = 1: p 1p eik0 y ; y 2ik0 p 1p e ik0 y ; y 2ik0

8

0 < 0:

Since pressure is symmetric about y = 0, it su¢ ces to construct the solution for y 0. We henceforth take Gx (x; ) ik0 p p e 2ik0

p(x; y; ) =

y

Assuming that the walls of the tank are rigid, Uxw = 0 at x = wl ; wr . From (15) we than have G0x = 0 at x = wl ; wr . This gives 8 C1 cos aw (x wl ); wl x < 0 < A cos aw x + B sinaaww x ; 0 < x < bl Gx (x; ) = : C8 cos aw (x wr ); bl + w < x wr :

The conditions

Gx (0+; ) G0x (0+; )

Gx (0 ; ) = 0 G0x (0 ; ) = 1:

give A = C1 cos aw wl ; B = C1 aw sin aw wl

1

and hence Gx (x; ) = C1 cos aw wl cos aw x + (C1 aw sin aw wl = C1 cos aw (x

3.2

1)

sin aw x aw

sin aw x ; 0 < x < bl : aw

wl )

Solution of the Biot equations

In (10) We make the change of dependent variables := ( + 2 )e + Q ;

:= Qe + R ;

(18)

a12 + a22

(19)

the inverse transformation for which is e = a11

a12 ;

=

where a11 := R=d;

a12 := Q=d;

a22 := ( + 2 )=d

with d := R( + 2 )

Q2 :

In terms of the new variables (10) becomes r2 + B11 + B12 = 0

r

2

+ B21 + B22 = 0

9

(20)

where B11 := a11 p11 B21 := a11 p12

a12 p12 ; a12 p22 ;

B12 := B22 :=

a12 p11 + a22 p12 ; a12 p12 + a22 p22 :

From (9) the ‡uid displacement vector is U=

1 (r + p12 u) : p22

(21)

Substituting this into (8) gives a partial di¤erential equation for the frame displacement vector r2 u + A31 r + A32 r + B33 u = 0

(22)

where A31 :=

1

a11

;

p12 ; p22

A32 := a12

B33 :=

p11 p22 p212 : 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

p

ux (x; y) = ux (x)eik0

p

y

p

y

; (x; y) = (x)eik0 y

p

y

; uy (x; y) = uy (x)eik0

p

y

and so forth. Substituting this into (20) gives 00 00

where a :=

p

B11

(x) + a2 (x) + B12 (x) = 0 2

(x) + B21 (x) + a

k02 and a := 0000

+ (a2 + a2 )

00

p B22

(x) = 0:

k02 . Elimination of

+ (a2 a2

(23)

in (23) gives

B12 B21 ) = 0

which has general solution (x) = C2 cos m+ (x bl )+C3

sin m+ (x m+

bl )

10

+C4 cos m (x bl )+C5

sin m (x bl ) m (24)

where v u u 2 t a + a2

q 2 (a2 + a2 )

4(a2 a2 B12 B21 ) m := 2 v q u u 2 2 (a2 a2 ) + 4B12 B21 ) t a + a2 = 2 s p B11 + B22 2k02 (B11 B22 )2 + 4B12 B21 = 2 s p B11 + B22 (B11 B22 )2 + 4B12 B21 = k02 : 2 From (24) and (23) we have (x) =

1 B12

00

(x) + a2 (x) :

The horizontal displacements can be calculated from (22) and (21) u00x (x) + A31 0 (x) + A32 0 (x) + a2u ux (x) = 0 where au :=

p

B33

(25)

k02 and Ux (x) =

1 ( 0 (x) + p12 ux (x)) : p22

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 p m = k0 m0 ; au = k0 au0 where

m0 =

B11 + B22

p (B11 B22 )2 + 4B12 B21 B33 ; au0 = 2 : 2k02 k0

p

p

@u

y 0 ik0 y x +uy (x)ik0 Since e(x)eik0 y = r hux ; uy i = @u @x + @y = ux (x)e we have, after a similar computation for the ‡uid dilatation,

1 p (e(x) ik0 1 p ( (x) Uy (x) = ik0 uy (x) =

11

u0x (x)) Ux0 (x)) :

p

eik0

p

y

The frame normal and tangential stress in the bone specimen are given by xx (z)

= e(x) + 2 exx (x) + Q (x) = e(x) + 2 u0x (x) + Q (x)

xy (x)

= =

@uy @ux + @y @x p ik0 ux (x) + u0y (x) :

Finally from (14)2 we have 1 G0 (x; ): 2 x ! 0

Uxw (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 )ux (xd +) Gx (xd ) = xx (xd +) + (xd +) Gx (xd ) = (xd +)= 0 = xy (xd +):

(26)

These conditions derive from continuity of ‡ux, aggregate normal pressure, pore ‡uid 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 2 3 C1 6 7 M 4 ... 5 = v: (27) C8

The general solution has the form p(x; y) =

k02 2 i

I

p(x; y; )d

(28)

C

where C is a suitably chosen contour in the -plane. Substituting this into (14)

12

gives r2 p + k02 p I I k2 k2 = 0 G00x (x; )Gy (y; )d + 0 Gx (x; )G00y (y; )d 2 i C 2 i C I k04 Gx (x; )Gy (y; )d + 2 i C I k02 k02 Gx (x; ) + (x) Gy (y; )d = 2 i C I k02 Gx (x; ) k02 (1 )Gy (y; ) + (y) d 2 i C I k4 + 0 Gx (x; )Gy (y; )d 2 i C I I k2 k2 = (x) 0 Gy (y; )d (y) 0 Gx (x; )d 2 i C 2 i C Since

k02 2 i

I

Gx (x; )d = (x);

C

for a positively oriented contour C which encloses the singularities of Gx , p is a solution of (14) providing C is chosen to exclude the singularities of Gy and enclose those of Gx . Substituting the results of the previous section into (28) gives I Gx (x; ) ik0 p y k2 p d (29) p(x; y) = 0 e 2 i 2ik0 C I k0 Gx (x; ) ik0 p y p = e d 4 C I 8 p C1 ( ) cos aw (x wl ) ik0 k0 y > p > e d ; wl < x < 0 > 4 > > > > C I > > < k0 C1 ( ) cos aw (x wl ) sinaaww x ik0 p y p e d ; 0 < x < bl 4 (30) = > > C I > > p > > C8 ( ) cos aw (x wr ) ik0 k0 y > p e d ; br < x < wr : > 4 > : C

p The only singularity of Gy is the , thus it can be excluded by taking C to be a cut enclosing the positive -axis. Thus the …rst 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 Gx are the solutions f n g 13

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

p(x; y) =

=

3.4

8 > > > >
> > > : 8 < :

k0 4 k0 4

I

IC

1(

8(

) cos aw (x wl ) ik0 p e ( ) ) cos aw (x wr ) ik0 p e ( )

p

p

y

d ; wl < x < bl (31)

y

d ; br < x < wr :

C

ik0 P 1( n) n p n d ( n ) cos aw ( n )(x 2 d ik0 P 8( n) p cos aw ( n )(x d n 2 ( n) n

wl )eik0 ik0

wr )e

p

p

ny ny

; wl < x < b l

; b r < x < wr :

d

(32)

Numerical implementation of the pressure equation

The …rst of the four interface conditions involves displacement, whereas the remaining three involve pressure. Typically the two quantities di¤er 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é coe¢ cient for water. We are interested in determining the pressure …eld 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. Figure 3, which shows the level curves Cn : ( ) = 10n ; n = 0; 50; 100; : : : ; 300 for 100 kHz with the stated parameters, indicates a more formidable di¢ culty in implementing the equations (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 ‡oating point is 10308 . The problem is the presence of terms of the form sin ad and cos ad where a is one of the coe¢ cients aw ; m or au and d is a distance. For instance the outer two columns of the characteristic matrix are of the form 3 2 M11 sin aw dl 0 6 M21 cos aw dl 7 0 6 7 7 6 M31 cos aw dl 0 6 7 6 7 0 0 7: 6 M=6 7 0 M sin a d 58 w r 7 6 6 7 0 M cos a d 68 w r 7 6 4 0 M78 cos aw dr 5 0 0

where dl = bl wl ; dr = wr not grow exponentially as

br and the coe¢ cients Mjk depend upon but do moves away from Re 1. Likewise the interior

14

columns of M contain exponentially growing terms 2 0 M13 M16 6 M22 0 0 6 6 M 0 0 32 6 6 0 M M 43 46 M2:::7 = 6 6 M52 sin m+ w M53 cos m+ w M56 cos au w 6 6 M62 cos m+ w M63 sin m+ w M66 sin au w 6 4 M72 cos m+ w M73 sin m+ w 0 M82 sin m+ w M83 cos m+ w M86 cos au w

p To deal with such terms we let a := k0 a0 as sin ad = e

d

cos ad = e

d

ei d e

= + i;

2 d

e

i d

+e

i d

0 M27 0 0 M57 sin au w M67 cos au w 0 M87 sin au w

2 d

7 7 7 7 7 7: 7 7 7 7 5

0 and rewrite them

(33)

2i ei d e

3

2

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. Figures 4-6 show the level curves of at frequencies 100; 500 and 1000 kHz. The value of j j 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. The right hand side of (27) is 2 cos a b 3

Thus we have i w drx e e w drx e in (31) is 8(

) cos aw (x p ( )

8

= e

bl

2 w drx +e

6 6 6 6 v=6 6 6 6 4

e

w dl

e2

i w drx

2

wr )

eik0

p

y

=e

w l 2 0! sin aw bl 2 0 ! aw sin aw bl 2 0 ! aw

0 .. . 0

+w

7 7 7 7 7: 7 7 7 5

e2 w e2 u w 8 . Finally cos aw (x wr ) = with drx = wr x. Thus the second integrand

w (bl +drx

dr )

ei

w drx

e

2

w drx

2

+e p

i

w drx

8 ik0 p y

e

(34) Since bl + drx 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

15

:

Figure 3: Level curves for

( ) at 100 kHz.

Figure 4: Level curves for ( ) at 100 kHz.

16

Figure 5: Level curves for ( ) at 500 kHz.

Figure 6: Level curves for ( ) at 1000 kHz.

17

gives 1(

) cos aw (x p ( )

wr )

eik0

p

y

=e

w (bl +dlx

dl )

ei

w dlx

e

2

w dlx

2

i

+e p

w dlx

1 ik0 p y

e

(35) We now consider the residue expansion (32). Introducing the operators 1 @ := @ 2

@ @ 1

@ 1 @ @ := +i @ 2 @ 1 @ 2 p we have for a term of the form a = k0 a0 = +i i

@ @ 2

;

@ p k0 a0 @

=

k p 0 2 a0

=

;

=

1

+i

2;

@ @ +i : @ @

Since a is an analytic function of , 0=

@a @ @ = +i @ @ @

)

d @ =i d @

which gives @ @

k02 i : 4a

=

This gives d d

k02 i 4

= e

w dl

dl + dr + 2w aw e

w dr

e2

+w

e2

1 1 1 + + m+ m au w 2

e

uw

+

@ @

:

The term dd is complicated and will be computed numerically. Thus for br < x < wr we have for the terms in (32) p =

p 8( n) cos aw ( n )(x wr )eik0 n y d n d ( n) e w (bl +drx dr ) 8 ei w drx e 2 w drx + e i

p 2

n

h

k0 i 4

dl +dr aw

+ 2w

1 m+

+

1 m

+

w drx

1 au

eik0 +

p

@ @

whereas for wl < x < bl p

p 1( n) cos aw ( n )(x wl )eik0 n y d n d ( n) e w (bl +dlx dl ) 1 ei w dlx e 2 w dlx + e i

= p 2

n

h

k0 i 4

dl +dr aw

+ 2w

1 m+

18

+

1 m

+

w dlx

1 au

eik0 +

ny

i

:

p

ny

@ @

i

:

:

Figure 7: Eigenvalues at 50 kHz for a receiver positioned at x = 0:03; y = 0 m.

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. Figures 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 jpj). 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 = 1. Clearly regions to the left of Re = 10 would have to be included to accurately compute pressure at the frequencies shown. As Figure 12 shows, the eigenvalues near the right most part of Re 1 increasingly dominate the computation of pressure as the receiver position is moved from away from the x-axis. The eigenvalues shown in Figures 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 (Figure 13). The parabolas used approximated the level curves of ( ) = 10 50 and 10500 at 100 kHz 1 (see Figure 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 1 ( ) could be computed up to about 10500 by recaling the terms in the characteristic matrix M.

19

Figure 8: Eigenvalues near Re x = 0:03; y = 0 m.

1 at 50 kHz for a receiver positioned at

Figure 9: Eigenvalues at 100 kHz for a receiver positioned at x = 0:03; y = 0 m.

20

Figure 10: Eigenvalues near Re x = 0:03; y = 0 m.

1 at 100 kHz for a receiver positioned at

Figure 11: Eigenvalues at 250 kHz for a receiver positioned at x = 0:03; y = 0 m.

21

Figure 12: Eigenvalues at 100 kHz for a receiver positioned at x = 0:03; y = 0:01 m.

from consideration. If pq m an attempt was made to …nd 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 subquadrilaterals and the considerations just described were repeated for each subquadrilateral. The procedure is obviously time-consuming, requiring around 40 minutes at 100 kHz and more at higher frequencies2 . It was e¤ective up to around 750 kHz. By 1 MHz it was ine¤ective because the integrands around even a small quadrilateral were sometimes too oscillatory to compute accurately. See Figure 14. The tolerances used for computing the eigenvalues in Figures 7-12 were m = 10 9 ; a = 5 10 9 . As mentioned above …nding the eigenvalues is computationally intensive. Clearly integration along a contour that encloses all of the eigenvalues that make a signi…cant contribution to the pressure will be more e¢ cient, however knowing the location of such eigenvalues is necessary to do this. Figures 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 (Figure 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 …rst quadrant. Figures 16 and 17 however give no indication that 2 Unless

otherwise stated the timed computations were done with a mobile Intel Pentium 4 processor running at 2:4 GHz

22

Figure 13: Subdivision of inner (Top) and region between inner and outer parabolas (Bottom).

Figure 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 :

23

Figure 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 .

subsequent eigenvalues follow these trajectories. This is important since Figures 4-6 indicate that ( ) is di¢ cult to compute far to the right of = 1. Also as frequency increases the trajectories tend to be closer to the axis Re 1:Thus it seems that in the 100 kHz and higher range the parabolic contours shown in the …gures 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 Figure 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 di¢ culty in the contour integration is the integration along the slit between = i and = 1 i. Figure 19 shows that the magnitude of the integrand ‡uctuates 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 (Figure 20) seemed to cause quad no di¢ culty. In the regions where the integrand was oscillatory (Figures 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

24

Figure 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 .

25

Figure 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 .

Figure 18: Typical contour used for computing pressure.

26

accurate to the speci…ed tolerance. The sequence of computations shown here illustrates this: Frequency: 100 kHz, Rel. Tol. = 5.000000e-005, Int. Tol. = 1.466741e-006 jplastj = 2.933482e-002, time = 4.13 N = 1, jpj = 2.929160e-002, rel. err. = 3.028707e-002, time = 4.90 N = 2, jpj = 3.213897e-002, rel. err. = 9.410061e-002, time = 5.36 N = 4, jpj = 2.955957e-002, rel. err. = 1.010306e-001, time = 6.34 N = 8, jpj = 2.913781e-002, rel. err. = 1.531679e-002, time = 7.50 N = 16, jpj = 2.911204e-002, rel. err. = 9.478357e-004, time = 8.65 N = 32, jpj = 2.914829e-002, rel. err. = 1.432803e-003, time = 10.20 N = 64, jpj = 2.915810e-002, rel. err. = 3.482070e-004, time = 12.03 N = 128, jpj = 2.915833e-002, rel. err. = 2.928985e-005, time = 13.91 N = 256, jpj = 2.915802e-002, rel. err. = 2.201273e-005, time = 17.49

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 = jplastj 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 modi…ed to accept vector arguments in the y-variable and require the 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 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

27

Figure 19: Magnitude of integrands along the line

Figure 20: Spike near integrand.

=

+ i to

= 1 + i.

= 0. Real (red) and imaginary (green) parts of the

28

Figure 21: Real (red) and imaginary (green) parts of the integrand.

Figure 22: Real (red) and imaginary (green) parts of the integrand.

29

Frequency 100 250 500 750

N ( = 0:72) 64 128 512 2048

Time 12.03 18.02 47.44 121.13

N ( = 0:83) 32 256 128 512

Time 14.27 34.42 59.10 76.12

N ( = 0:95) 32 32 64 256

Time 21.56 41.27 80.34 110.82

Table 4: Number of intervals required to achieve a tolerance of 5e-5 Frequency 100 250 500 750

N ( = 0:72) 128 64 128 512

Time 14.54 25.41 60.44 211.38

N ( = 0:83) 32 256 128 128

Time 15.35 41.22 55.00 83.59

N ( = 0:95) 32 64 64 64

Table 5: x=0.03; y = 0:.001:.02 N = 512,rel. err. = 8.372496e-006, time = 25.49

As can be seen by comparison to 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. As indicated in Figures 9-11 eigenvalues lying along the negative real Frequency 100

250

500

0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95

# Eigs 894 896 896 2214 2232 2230 4299 4439

2:9696e 2:8956e 2:4116e 7:4316e 1:1638e 9:1995e 1:0134e 4:4549e 3:6933e

p 002 + 2:5700e 002 + 1:7255e 002 + 3:8046e 004 7:6031e 002 + 1:7442e 003 + 2:4580e 003 9:0914e 003 + 3:1661e 003 + 1:5300e

Table 6: Table Caption

30

004i 002i 002i 003i 004i 002i 005i 004i 002i

Time 23.05 50.09 83.32 131.79

Frequency 100

250

500

0.72 0.83 0.95 0.72 0.83 0.95 0.72 0.83 0.95

2:9697e 2:8955e 2:4115e 7:4262e 1:1638e 9:1994e 1:0134e 4:4550e 3:6936e

p 002 + 2:5531e 002 + 1:7256e 002 + 3:8046e 004 7:6030e 002 + 1:7478e 003 + 2:4581e 003 9:0860e 003 + 3:1657e 003 + 1:5300e

004i 002i 002i 003i 004i 002i 005i 004i 002i

Time 24:4 34:1 53:6 39:8 60:8 118:9 81:8 132:3 223:5

Table 7: Table Caption Re l 20 40 80 160 320 640 1280 2560

3:0638e 2:5556e 2:8747e 2:9791e 2:8229e 2:8450e 2:8666e 2:9291e

p 002 + 3:8728e 002 + 3:8674e 002 + 3:4166e 002 + 3:0984e 002 + 3:2654e 002 + 3:2619e 002 + 3:2362e 002 + 3:1504e

003i 003i 003i 003i 003i 003i 003i 003i

Time 24:8 26:0 28:5 33:0 42:3 62:4 99:9 172:7

Table 8: f = 100 kHz,y=0 axis may make a signi…cant 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 Figure 3 back to Re = 2560 results in only about two digit accuracy. The situation improves 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 di¤erent 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. Figures 23-26 show the magnitude of pressure for three di¤erent bone specimens at di¤erent 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 di¤erentiated for larger y-values.

31

Re l 20 40 80 160 320 640 1280

2:9157e 2:8671e 2:8904e 2:8918e 2:8910e 2:8911e 2:8911e

p 002 + 3:5317e 002 + 3:4693e 002 + 3:4328e 002 + 3:4249e 002 + 3:4257e 002 + 3:4256e 002 + 3:4256e

003i 003i 003i 003i 003i 003i 003i

Time 24:7 25:8 28:3 32:7 43:3 60:6 98:1

Table 9: f = 100 kHz,y=0.001

Frequency 100

250

500

0:72 0:83 0:95 0:72 0:83 0:95 0:72 0:83 0:95

Re 320 320 320 80 80 80 40 40 20

jpj 2:9113e 3:5905e 4:6836e 6:5416e 1:2102e 2:6703e 9:4016e 3:7096e 1:5470e

002 002 002 003 002 002 004 003 002

Time 43:3 51:6 69:7 40:8 59:2 109:1 65:4 93:0 156:4

Table 10: x = 0.03, y = 0.001

Frequency 100

250

500

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

Time 48:2 59:7 81:2 57:3 73:1 123:4 122:6 125:0 224:2

Table 11: x = 0.03, y = 0.001:0.001:0.02

32

Frequency 100

250

500

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

Time 59:0 68:3 86:1 95:4 93:8 129:4 164:7 149:7 276:1

Table 12: x = 0.0, y = 0.001:0.001:0.02 Frequency 100

250

500

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

Time 53:7 60:0 77:4 73:3 75:7 112:5 139:6 130:9 211:4

Table 13: x = -0.01, y = 0.001:0.001:0.02

Figure 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).

33

Figure 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).

Figure 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).

34

Figure 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).

Figure 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).

35

Figure 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).

4 4.1

Recovery of the Biot parameters Preliminaries

We now address the question of whether the Biot parameters can be recovered from acoustic measurements in experiments like those conducted by McKelvie and Palmer and Hosokawa and Otani. The general approach to this sort of inverse problem is to assume that a vector of M measurements of acoustic presM sure P = fPm gm=1 is available. It is generally found that inversions are more successful if both magnitude and phase data are available [cite], whence the Pm are complex-valued. The inversion algorithm seeks to recover the parameters by minimizing the di¤erence kP P k in some norm where P is a vector of pressures calculated using some mathematical model. The minimizing procedure we used was the Nelder-Mead simplex method since it does not require the computation of gradients and is known to be good at minimizing functions such as the Rosenbrock function which have contours which deviate considerably from a paraboloid. Indeed Figure 29 which shows the simplex method’s manipulation of porosity and permeability during a minimization is reminiscent of the Rosenbrock function. The objective function used was simply f (P; P ) = kP

P k = kP

P k2 =M

(36)

where M is the number of positions at which the acoustic …eld was measured. Because of the vastly di¤erent magnitudes of the parameters of the Biot model 36

Figure 29: Values of permeability and porosity during a run of the simplex method.

the simplex method manipulated the common logarithms of the parameters, rather than the parameters themselves. Lacking experimental data we were forced to rely upon the approach derived above for both the target data P and the trial data P calculated during the parameter recovery process. In this sense we are committing the so-called inverse crime, since it may be that the Biot model simply produces poor approximations to the pressure …eld generated by cancellous bone, and even if this is not the case, the simpli…cations such as using a two-dimensional model and assuming the bone and tank are in…nite in the y-direction might lead to considerable inaccuracy. This commission of the inverse crime cannot be avoided at this time. We can however avoid one of the more felonious versions of the inverse crime, that of calculating P and P in exactly the same way in which case the minimization problem has a (probably) unique solution P = P . This is of course unrealistic since there will be some modelling and measurement error, and it can generate misleadingly good results. We shall avoid this by reducing the magnitude of three of the parameters upon which the calculation depends: the target accuracy, the number N of subdivisions used integrating along the portion of the branch cut 0 Re 1, and the left hand limit l which controls how far out to the left the contour integration is carried. Table 14 gives the tolerances used for calculating the target data. For the trial data the target tolerance was increased by a factor of 10 and N and l reduced by a factor of 2. Based on the results given above it seems most promising to attempt to recover the Biot parameters using acoustic pressure measurements taken to the right of the specimen where the acoustic pressure varies more with the bone

37

Frequency 100 250 500

tolerance 5 10 5 5 10 5 5 10 5

N 128 256 512

l

320 80 40

Table 14: Table Caption

Figure 30: Relative error as a function of receiver position x due to reducing the accuaracy of the computation of pressure. The vertical coordinate was y = 0:001: Top: 100 kHz; Middle: 250 kHz; Bottom: 500 kHz. The porosities of the specimens were Blue: 0:72; Red: 0:83; Green: 0:95.

specimen. Figure 30 shows that for receiver positions 0:03 x 0:05; y = 0:001 the relative discrepancy j(Pm Pm )=Pm j was about 0:005 for bone specimens of porosity 0:72; 0:83 and 0:95 at 100 kHz. At frequencies 250 and 500 kHz the discrepancy for = 0:95 was considerably lower. This might distort results and as indicated above the computations are quicker for lower frequencies. For this reason we chose to test our parameter recovery algorithm at 100 kHz. The simplex method requires an initial guess. The topography of the objective function space is far too complicated for one guess to su¢ ce. It is expected that there will be numerous local minima. Either a means of generating a good initial guess has to be found or multiple guesses must be tried. We have yet to …nd a way to do the former, and so the latter approach seems more promising. This sort of problem is thus amenable to parallel processing. The results presented below were obtained using a group of nine 1:5 GHz PC’s whose only connection was common access to a set of …les on a network drive.

38

4.2

An algorithm for the recovery of the Biot parameters: Phase 1

In devising a means to calculate pressure we used Biot parameters of the …ve specimens of Table 2. In developing a parameter recovery algorithm it is desirable to have a larger set of test problems. This was accomplished by using a pseudo-random number generator to generate parameter values based on the uncertainty ranges given in Section 2. As in [4] the initial stage of the algorithm is to reduce the problem to one of minimizing the objective function with respect to porosity alone. To do this relations which estimate the value of a parameter at a given porosity are required: Frame moduli Re Kb and Re : To generate values of these parameters for the test problems n was chosen randomly from the interval [1:1; 1:6] and the real parts of the moduli were calculated (12). The imaginary parts were calculated using a log decrement of 0:1. For the univariate stage of the algorithm the median value n = 1:35 is used in (12). Since the imaginary parts were not believed to be important, the values used were those for the 0:81 specimen of Table 2. Pore size a: For a given value of porosity the value assigned to pore size in a test problem was chosen randomly from the interval indicated in Figure 1. The univariate phase of the algorithm used the regression line value shown in the Figure. Permeability k: For a given value of porosity the value assigned to permeability for a test problem was chosen randomly from the interval indicated in Figure 2. The univariate phase of the algorithm used the regression line value shown in the Figure. The structure constant was calculated for a test problem from = 1 s(1 1= ) with s chosen randomly from [0:2; 0:3] for a test problem. For the univariate phase of the algorithm s = 0:25 was used. The frame material density r was chosen randomly from the interval [1930; 2000] for a test problem. For the univariate phase r = 1965 was used. The values of all other Biot parameters were taken from Table 3 and were not treated as being subject to uncertainty. For the univariate phase of the algorithm the minimizing value of porosity was sought in the range 0:67 0:96. Twenty-…ve target parameter sets were chosen randomly for each of the seven porosities = 0:71; 0:75; : : : ; 0:95. It was assumed that measurements were taken at eight evenly spaced positions between x = 0:03 and x = 0:05 with y = 0:001. Figure 31 shows that, for the uncertainty bounds given above, this phase of the algorithm was not successful, the value of being in error by as much as 20%. In particular a returned value of porosity 39

Target Univariate Error Target Univariate Error Target Univariate Error Target Univariate Error Target Univariate Error Target Univariate Error Target Univariate Error

0.71 0.846 19.14% 0.75 0.868 15.76% 0.79 0.903 14.26% 0.83 0.745 10.22% 0.87 0.953 9.58% 0.91 0.857 5.87% 0.95 0.670 29.47%

k 1.32E-08 5.05E-08 282.55% 3.16E-08 7.96E-08 151.99% 8.64E-08 1.60E-07 85.46% 7.93E-09 6.52E-09 17.79% 3.43E-07 4.49E-07 31.10% 3.50E-08 6.28E-08 79.34% 1.79E-06 1.41E-09 99.92%

a 3.27E-04 1.45E-03 344.10% 7.02E-04 1.62E-03 130.19% 1.02E-03 1.87E-03 83.63% 1.46E-03 7.13E-04 51.18% 0.0016575 2.24E-03 35.26% 2.00E-03 1.53E-03 23.61% 2.07E-03 1.62E-04 92.16%

1.113 1.046 6.08% 1.100 1.038 5.60% 1.070 1.027 3.98% 1.053 1.086 3.08% 1.0319 1.012 1.91% 1.029 1.042 1.23% 1.014 1.123 10.81%

Re Kb 3.80E+09 1.63E+09 57.12% 4.25E+09 1.32E+09 68.93% 2.45E+09 8.77E+08 64.25% 2.39E+09 3.22E+09 34.46% 2.09E+09 3.25E+08 84.46% 7.01E+08 1.48E+09 111.07% 1.81E+08 4.56E+09 2417.97%

Re 1.56E+09 6.67E+08 57.12% 1.74E+09 5.40E+08 68.93% 1.00E+09 3.59E+08 64.25% 9.79E+08 1.32E+09 34.46% 8.56E+08 1.33E+08 84.46% 2.87E+08 6.06E+08 111.07% 7.41E+07 1.87E+09 2418.01%

Table 15: Result of the univariate phase for Problems 71w,...,95w. of around = 0:85 was compatible with practically any porosity between 0:71 and 0:95. It was decided to use the specimens at each of the seven porosities with the greatest error in porosity in the univariate phase as target specimens. These will be designated as Problems 71w; : : : ; 95w respectively. Table 15 gives the values of the target parameters and the results of the univariate phase. Observe that the results for Problem 95w are not consistent with Figure 31. The Figure indicates that the worst outcome of the univariate algorithm returned an estimate of about = 0:85 rather than the endpoint minimum of = 0:67 indicated in the Table. The explanation appears to be that in generating the Figure the algorithm was applied to the target parameters at full IEEE double precision (about 15 digits), but the parameter values were stored to six digits. In generating the results shown in the Table these stored values were used by the univariate algorithm and that was su¢ cient to drastically change the outcome. A more serious problem emerged when the simplex algorithm was applied to trial specimens generated by the algorithm which is described below. The values of the parameters Re Kb and Re tended to drift to huge numbers, on the order of 1013 or 1014 or to very small numbers on the order of 105 and this apparently caused MATLAB’s routine quad, which was used to generate the trial pressure …elds in the manner indicated in the previous section, to exceed its permitted number of function evaluations. This phenomenon was con…ned to the = 0:95

40

r

1999.9 1965 1.75% 1952.8 1965 0.62% 1967.1 1965 0.11% 1977.1 1965 0.61% 1981.1 1965 0.81% 1954.8 1965 0.52% 1930.1 1965 1.81%

Figure 31: Outcome of applying the univariate algorithm to randomly selected parameter sets.

specimen and since it is not clear how important this extreme case is, it was dropped from consideration.

4.3

An algorithm for the recovery of the Biot parameters: Phase 2

The next phase of the algorithm was done in parallel using a group of 1:5 GHz PC’s whose only connection was common access to a set of …les on a network drive. This phase consisted of three stages. For the …rst stage the interval [(1 0 ) 0 ; (1 + 0 ) 0 ] was subdivided into Np equally spaced nodes, where Np was the number of processors, 0 is the result of the univariate phase, and 0 = 0:2=(1 + 1=Np ). The value of 0 was chosen so that the range of values tried over the …rst two stages could be 0 0:2 if necessary. For each n ; n = 1; : : : ; Np , the remaining parameters were calculated from n as they were for the univariate phase. The Nelder-Mead simplex algorithm was then applied by one of the processors to minimize the objective function (36) with this set of parameters as the initial guess. The number of iterations permitted to the simplex method was 64; 64; 256 for stages 1; 2; 3 respectively. When a processor …nished a stage it took its result and used it to compute Np candidate parameter s 1 for stage s. The sets for the next stage. This was done setting = 0 =Np interval [(1 ) ; (1 + ) ] was again subdivided into Np equally spaced nodes n as was the interval [(1 ) ; (1 + ) ] for the structure factor. For permeability the interval [(1 2 ) log k; (1+2 ) log k] was subdivided into Np equally spaced nodes ln to produce Np guesses for permeability kn = 10ln . The same was done for pore size and the two frame moduli. Preliminary exploration 41

indicated that the simplex method often returned a value of frame density r which was well outside of its interval of uncertainty [1930; 2000] and so pursuit of this parameters value was discontinued. It was set at the median value r = 1965 at all stages. The guesses for the parameters kn ; an ; n ; Re Kbn ; Re n were put in ascending or descending order according to whether they are expected to increase or decrease with porosity and an objective function value was computed for each of the Np assemblages ( n ; kn ; an ; n ; Re Kbn ; Re n ). The results were sorted according the objective function value and interleaved with current result stack with the Np best objective function values being retained. As an illustration of how Phase 2 worked we present the cases of Problems 71w and 83w. Np = 9 computers were used. For 71w the nine guesses at the start of Stage 1 were 3.47E-03 3.54E-03 4.34E-03 4.80E-03 5.41E-03 6.41E-03 6.57E-03 8.28E-03 9.65E-03

0.860 0.827 0.893 0.794 0.927 0.960 0.760 0.727 0.694

6.75E-08 3.43E-08 1.33E-07 1.74E-08 2.61E-07 5.14E-07 8.85E-09 4.50E-09 2.29E-09

1.56E-03 1.31E-03 1.80E-03 1.07E-03 2.05E-03 2.29E-03 8.24E-04 5.80E-04 3.35E-04

1.041 1.052 1.030 1.065 1.020 1.010 1.079 1.094 1.110

1.43E+09 1.91E+09 9.92E+08 2.42E+09 5.98E+08 2.64E+08 2.96E+09 3.53E+09 4.12E+09

5.86E+08 7.81E+08 4.06E+08 9.91E+08 2.45E+08 1.08E+08 1.21E+09 1.44E+09 1.69E+09

The …rst number is the value of the objective function and remaining numbers are n ; kn ; an ; n ; Re Kbn ; Re n . Observe that the two guesses for porosity which are in fact closest to correct have the highest objective function values at this point. At the end of Stage 1 after applying 64 iterations of the simplex method to each guess and and using the interleaving process described above the nine guesses for Stage 2 were 4.54E-05 5.49E-05 1.01E-04 1.43E-04 1.57E-04 1.99E-04 2.10E-04 2.28E-04 2.64E-04

0.687 0.721 0.754 0.693 0.690 0.717 0.725 0.821 0.683

9.74E-09 1.57E-08 3.80E-08 1.41E-08 1.17E-08 1.31E-08 1.88E-08 1.41E-08 8.10E-09

2.52E-04 4.01E-04 7.08E-04 2.97E-04 2.73E-04 3.71E-04 4.33E-04 4.23E-04 2.32E-04

1.111 1.094 1.083 1.099 1.105 1.099 1.088 1.068 1.116

5.96E+09 2.07E+09 1.20E+09 3.80E+09 4.76E+09 2.57E+09 1.67E+09 1.01E+07 7.47E+09

1.76E+09 2.55E+09 1.37E+09 1.15E+09 1.42E+09 3.16E+09 2.05E+09 1.96E+11 2.18E+09

Now eight of the nine guesses for porosity are close to the correct value. At the end of Stage 2 after another 64 iterations of the simplex method the guesses for Stage 3 were 1.93E-05 1.97E-05 2.25E-05

0.709 0.709 0.695

1.01E-08 9.86E-09 1.34E-08

2.46E-04 2.44E-04 3.21E-04

42

1.112 1.112 1.099

4.77E+09 4.89E+09 3.84E+09

1.37E+09 1.41E+09 1.64E+09

2.38E-05 2.54E-05 2.70E-05 2.78E-05 2.94E-05 3.03E-05

0.695 0.695 0.717 0.694 0.720 0.710

1.32E-08 1.37E-08 1.33E-08 1.29E-08 1.57E-08 1.03E-08

3.18E-04 3.24E-04 3.63E-04 3.15E-04 4.03E-04 2.48E-04

1.100 1.098 1.100 1.100 1.094 1.111

3.94E+09 3.75E+09 2.63E+09 4.03E+09 2.15E+09 4.66E+09

1.68E+09 1.61E+09 3.13E+09 1.72E+09 2.62E+09 1.34E+09

Now all of the guesses for porosity are close to correct and the ranges for the other parameters are much narrower. The result at the end of Stage 3 was 7.69E-06 7.72E-06 7.77E-06 7.77E-06 7.79E-06 7.80E-06 7.81E-06 7.82E-06 7.82E-06

0.710 0.710 0.710 0.710 0.710 0.710 0.710 0.710 0.710

1.28E-08 1.28E-08 1.30E-08 1.31E-08 1.28E-08 1.30E-08 1.15E-08 1.15E-08 1.31E-08

3.11E-04 3.11E-04 3.19E-04 3.20E-04 3.12E-04 3.19E-04 2.69E-04 2.69E-04 3.20E-04

1.110 1.110 1.112 1.112 1.110 1.112 1.101 1.101 1.112

3.80E+09 3.79E+09 3.80E+09 3.79E+09 3.78E+09 3.81E+09 3.85E+09 3.86E+09 3.78E+09

1.53E+09 1.53E+09 1.52E+09 1.52E+09 1.52E+09 1.53E+09 1.54E+09 1.54E+09 1.52E+09

Thus starting from the guesses at the end of Stage 2, all nine processors arrived at nearly the same answer. As stated above the trial pressures used by the simplex method were computed to less accuracy than the target pressures. Thus the objective function value for the target parameters will not be zero. The value of the objective function for the target specimen was in this case 8:06e 06. Since the answers obtained all have smaller objective function values, there is no reason why further reduction of the objective function value should produce more accurate parameter values. We shall refer to such answers as subcorrect. When an algorithm can attain subcorrect answers it must judged a success irrespective of the accuracy of the answers. Table 16 shows the progress that was made in reducing the percentage error over the three stages for the result that had the lowest objective function value at each stage. While the solution of Problem 71w proceeded as might be hoped, this was not always the case. Here is how the solution of Problem 83w evolved. The guesses for Stage 1 were: 1.46E-03 1.96E-03 2.52E-03 3.31E-03 4.19E-03 5.48E-03 5.94E-03 7.59E-03 9.08E-03

0.748 0.722 0.775 0.696 0.801 0.670 0.827 0.853 0.879

6.97E-09 4.10E-09 1.19E-08 2.41E-09 2.02E-08 1.41E-09 3.44E-08 5.85E-08 9.96E-08

7.38E-04 5.46E-04 9.30E-04 3.54E-04 1.12E-03 1.62E-04 1.31E-03 1.51E-03 1.70E-03

43

1.084 1.096 1.073 1.109 1.062 1.123 1.052 1.043 1.034

3.16E+09 3.61E+09 2.73E+09 4.08E+09 2.31E+09 4.56E+09 1.91E+09 1.53E+09 1.17E+09

1.29E+09 1.48E+09 1.11E+09 1.67E+09 9.44E+08 1.87E+09 7.80E+08 6.26E+08 4.80E+08

Problem 71w Univariate Stage 1 Stage 2 Stage 3 83w Univariate Stage 1 Stage 2 Stage 3

f (P; P )

k

a

Re Kb

Re

3.47E-03 4.54E-05 1.93E-05 7.69E-06

19.14% 3.31% 0.10% 0.00%

282.55% 26.29% 23.78% 3.28%

344.10% 23.09% 24.82% 4.94%

6.08% 0.24% 0.13% 0.27%

57.12% 56.75% 25.52% 0.20%

57.12% 13.22% 11.73% 1.66%

3.34E-03 1.28E-05 9.62E-06 8.20E-06

10.22% 6.96% 0.32% 1.13%

17.79% 22.39% 26.03% 9.03%

51.18% 27.19% 25.15% 9.89%

3.08% 1.84% 0.15% 2.02%

34.46% 2.96% 0.91% 0.75%

34.46% 0.07% 0.38% 0.48%

Table 16: Percentage errors after each stage of the algorithm for Problems 71w and 83w.. The guesses at the start of Stage 2 were 1.28E-05 1.31E-05 1.52E-05 2.99E-05 6.98E-05 8.78E-05 1.77E-04 1.80E-04 1.87E-04

0.772 0.745 0.799 0.694 0.825 0.730 0.849 0.664 0.873

6.15E-09 4.84E-09 7.59E-09 2.02E-09 1.03E-08 3.29E-09 1.50E-08 1.19E-09 2.55E-08

1.06E-03 8.12E-04 1.35E-03 3.08E-04 1.88E-03 5.19E-04 2.09E-03 1.55E-04 3.38E-03

1.073 1.084 1.062 1.110 1.052 1.100 1.042 1.124 1.034

2.46E+09 2.45E+09 2.44E+09 2.20E+09 2.28E+09 2.26E+09 1.33E+09 2.06E+09 5.84E+08

9.78E+08 9.91E+08 9.76E+08 1.06E+09 1.12E+09 1.19E+09 1.93E+09 8.74E+08 2.37E+09

The guesses at the start of Stage 3 were 9.62E-06 1.05E-05 1.15E-05 1.24E-05 1.87E-05 2.86E-05 7.26E-05 1.26E-04 1.37E-04

0.827 0.799 0.773 0.747 0.726 0.691 0.679 0.680 0.692

9.99E-09 7.60E-09 6.11E-09 4.82E-09 2.99E-09 2.07E-09 1.11E-09 1.13E-09 2.12E-09

1.83E-03 1.35E-03 1.06E-03 8.10E-04 4.89E-04 3.15E-04 1.66E-04 1.67E-04 3.18E-04

1.055 1.062 1.073 1.085 1.098 1.109 1.128 1.127 1.108

2.41E+09 2.45E+09 2.45E+09 2.44E+09 2.31E+09 2.21E+09 1.95E+09 1.90E+09 2.16E+09

9.75E+08 9.78E+08 9.88E+08 9.97E+08 1.05E+09 1.06E+09 1.10E+09 1.07E+09 1.04E+09

7.21E-09 9.81E-09 4.43E-09 5.67E-09 5.05E-09

1.32E-03 1.84E-03 7.55E-04 1.02E-03 8.84E-04

1.032 1.068 0.900 1.050 1.013

2.37E+09 2.37E+09 2.37E+09 2.37E+09 2.40E+09

9.74E+08 9.75E+08 9.64E+08 9.83E+08 9.80E+08

The …nal result was 8.20E-06 8.21E-06 8.25E-06 8.66E-06 8.71E-06

0.839 0.839 0.859 0.818 0.815

44

Problem 71w 75w 79w 83w 87w 91w Average Worst

f (P; P ) 7.69E-06 4.35E-06 5.07E-06 8.20E-06 5.96E-06 9.71E-06

f (Pt ; P ) 8.06E-06 8.52E-06 5.23E-06 8.63E-06 6.25E-06 4.74E-06

0.00% 0.20% 0.04% 1.13% 0.02% 0.74% 0.36% 1.13%

k 3.28% 18.60% 25.75% 9.03% 0.78% 4.53% 10.33% 25.75%

a 4.94% 25.22% 28.94% 9.89% 3.76% 0.02% 12.13% 28.94%

0.27% 1.13% 0.40% 2.02% 0.21% 1.08% 0.85% 2.02%

Re Kb 0.20% 7.54% 0.48% 0.75% 5.89% 17.59% 5.41% 17.59%

Table 17: Phase 2 percentage errors for Problems 71w, ..., 91w. 9.48E-06 1.16E-05 1.74E-05 1.75E-05

0.796 0.762 0.695 0.839

5.50E-09 3.56E-09 4.10E-09 7.23E-09

9.73E-04 6.00E-04 6.42E-04 1.32E-03

1.062 1.058 1.104 1.032

2.40E+09 2.37E+09 2.42E+09 2.37E+09

9.91E+08 1.00E+09 1.01E+09 9.72E+08

The objective function value for the target was 8:63E 06 and thus the best three answers are subcorrect and the next two best answers nearly so, yet the porosities of the …rst three answers spanned a range of more than 0:02 and the best …ve a range of more than 0:04. The estimates of permeability and pore size also varied considerably. Moreover since the application of the simplex method in the third stage was limited to 256 iterations it is possible that fewer or more iterations could have changed the outcome signi…cantly. Table 17 shows the …nal results of Phase 2 of the algorithm for the six Problems 71w; : : : ; 91w. The third column gives the value of the objective function for the target parameters. The algorithm was successful in obtaining subcorrect answers except for Problem 91w.

4.4

Further explorations

In view of the fact that allowing for a possible 20% error in the outcome of Phase 1 was sometimes tantamount to examining the entire assumed range of porosities we decided to test an alternative Phase 1 in which the assumed range of possible porosities [0:67; 0:96] was simply subdivided into Np + 1 intervals and the Np interior nodes used as the porosities for the initial guesses for Phase 2. The other parameters were calculated from porosity as in the univariate algorithm. Phase 2 was the same except that 0 was changed to 0:27=Np . Table 18 shows the result. The most noticeable di¤erences are in the answers to Problems 83w where the answers were worse than when the univariate Phase 1 was used and 91w where the answers for the alternative Phase 1 were better. For Problem 83w the alternative Phase 1 value for porosity was almost 2:4% in error and the error for the structure constant was very large, in fact being non-physical by virtue of being less than 1. The nine answers at the end of Phase 2 for Problem 83w for the alternative Phase 1 algorithm were

45

Re 1.66% 9.29% 0.50% 0.48% 5.19% 28.59% 7.62% 28.59%

Problem 71w 75w 79w 83w 87w 91w Average Worst

f (P; P ) 7.79E-06 4.36E-06 5.07E-06 8.14E-06 5.88E-06 5.17E-06

f (Pt ; P ) 8.06E-06 8.52E-06 5.23E-06 8.63E-06 6.25E-06 4.74E-06

0.01% 0.21% 0.04% 2.41% 0.10% 0.73% 0.58% 2.41%

k 12.66% 8.94% 26.09% 21.13% 12.63% 11.37% 15.47% 26.09%

a 17.42% 14.06% 29.46% 23.09% 13.61% 7.57% 17.54% 29.46%

1.11% 0.72% 0.43% 5.87% 0.05% 1.05% 1.54% 5.87%

Re Kb 1.00% 6.92% 0.72% 1.08% 4.43% 6.89% 3.50% 6.92%

Table 18: Phase 2 percentage errors for Problems 71w, ..., 91 with alternative Phase 1 8.14E-06 8.22E-06 8.23E-06 8.27E-06 8.28E-06 8.37E-06 8.50E-06 8.60E-06 8.62E-06

0.850 0.838 0.837 0.834 0.833 0.828 0.850 0.850 0.838

6.25E-09 8.93E-09 8.98E-09 8.11E-09 7.63E-09 5.00E-09 6.25E-09 6.25E-09 8.93E-09

1.12E-03 1.66E-03 1.67E-03 1.50E-03 1.40E-03 8.75E-04 1.12E-03 1.12E-03 1.66E-03

0.991 1.063 1.064 1.060 1.053 0.990 0.991 0.991 1.063

2.37E+09 2.37E+09 2.37E+09 2.37E+09 2.37E+09 2.39E+09 2.37E+09 2.37E+09 2.37E+09

9.72E+08 9.76E+08 9.76E+08 9.77E+08 9.77E+08 9.76E+08 9.71E+08 9.72E+08 9.75E+08

All were subcorrect and yet there was a good bit of variation in the values for porosity, permeability, pore size and structure factor. This illustrates the danger in simply taking as the answer the result with the lowest objective function value. For Problem 91w the alternative Phase 1 values for the moduli Re Kb and Re were more accurate, presumably because the objective function value was closer to subcorrect. The …nal question we addressed was whether the availability of multiple answers permitted estimation of the errors made when the actual answer was not known. Unfortunately the nine answers at the end of Phase 2 were not always suitable for this purpose. Sometimes the agreement among the answers was much better than the agreement with the actual answer. For Phase 3 the average of all answers at the end of Phase 2 that were within a factor of 2 of the lowest objective function value was computed. The Np guesses were computed as in Phase 2, that is for an average value v of porosity or the structure factor the interval [(1 )v; (1 + )v] was subdivided into Np points and for permeability, pore size, and the two moduli the interval [(1 2 ) log v; (1+2 ) log v] was subdivided. For each parameter the guesses were ordered according to whether they were expected to increase or decrease with respect to increasing porosity. The parameter was taken to be 0:02. To increase the chance that process would …nd a local minimum the simplex method was allowed to run until the simplex size was 10 4 , MATLAB’s default convergence criterion. For Problem 83w this led to the guesses

46

Re 0.34% 9.52% 0.37% 0.74% 5.23% 12.42% 4.77% 12.42%

8.08E-03 5.64E-03 2.86E-03 1.19E-03 2.30E-05 1.36E-03 2.76E-03 4.12E-03 5.37E-03

0.802 0.806 0.810 0.814 0.818 0.822 0.826 0.830 0.835

2.76E-09 3.34E-09 4.03E-09 4.87E-09 5.89E-09 7.12E-09 8.60E-09 1.04E-08 1.26E-08

8.02E-04 8.59E-04 9.20E-04 9.85E-04 1.06E-03 1.13E-03 1.21E-03 1.30E-03 1.39E-03

1.047 1.042 1.036 1.031 1.026 1.021 1.016 1.011 1.006

5.65E+09 4.55E+09 3.67E+09 2.95E+09 2.38E+09 1.92E+09 1.55E+09 1.25E+09 1.00E+09

2.25E+09 1.83E+09 1.49E+09 1.21E+09 9.82E+08 7.98E+08 6.49E+08 5.28E+08 4.29E+08

The result after the application of the simplex method was 8.83E-06 9.15E-06 8.99E-06 1.01E-05 9.28E-06 1.07E-05 8.34E-06 9.39E-05 8.15E-06

0.814 0.823 0.805 0.800 0.828 0.832 0.829 0.772 0.847

4.81E-09 6.65E-09 4.72E-09 4.04E-09 7.12E-09 8.82E-09 5.85E-09 5.85E-09 6.61E-09

8.53E-04 1.18E-03 8.32E-04 7.11E-04 1.27E-03 1.57E-03 1.05E-03 9.74E-04 1.20E-03

1.031 1.021 1.036 1.043 1.016 1.011 1.027 1.047 1.006

2.37E+09 2.42E+09 2.38E+09 2.35E+09 2.42E+09 2.45E+09 2.37E+09 4.13E+09 2.37E+09

9.86E+08 9.74E+08 9.87E+08 9.94E+08 9.71E+08 9.67E+08 9.79E+08 7.90E+07 9.72E+08

For each parameter the smallest and largest values among those whose objective function value was within a factor of 2 of the lowest obtained objective function value were taken as the con…dence interval for the parameter and the midpoint of this interval was taken as the parameter estimate. The results are shown in Tables 19-24. The con…dence interval encompassed the target values in all but three instances (marked with * in the Tables)3 . In the case of porosity and the real part of the shear modulus for Problem 79w the underestimation of the error was very slight and would likely be of no practical consequence. The underestimation of the possible error for the structure factor for Problem 83w was somewhat larger. A more signi…cant di¢ culty is that at the two highest porosities, = 0:87 and 0:91, the error estimates were sometimes much higher than the actual errors. It is also noteworthy that despite starting from nine guesses that were fairly near the target answer, the best objective function value for Problem 91w was slightly higher than that for the Phase 2 answer when the univariate Phase 1 was used (Table 17), and much higher than the Phase 2 objective function value when the alternative Phase 1 was used (Table 18). Thus this is a particularly di¢ cult set of target parameters and the success of the alternative Phase 1 algorithm is likely fortuitous. For purposes of comparison with the results of Tables 19-24 we also computed the mean and standard deviation of all Phase 3 answers whose objective 3 The slight underestimation of the error of Re K in Problem 79w was due to the use of b the midpoint value in normalizing the estimated error vice the use of the smaller target value in the computation of the actual percentage error.

47

Problem 71w Low Target High Midpoint Error Est. Error f (P; P ) f (Pt ; P )

0.709 0.71 0.716 0.712 0.33% 0.50% 7.69E-06 8.06E-06

k 1.04E-08 1.32E-08 1.52E-08 1.28E-08 2.96% 18.73%

a 2.29E-04 3.27E-04 3.58E-04 2.94E-04 10.31% 22.04%

1.073 1.113 1.125 1.099 1.27% 2.35%

Re Kb 3.16E+09 3.80E+09 4.96E+09 4.06E+09 6.73% 22.09%

Re 1.38E+09 1.56E+09 1.59E+09 1.49E+09 4.57% 7.04%

Table 19: Phase 3 results for Problem 71w, The objective function value is lowest of the nine results.

Problem 75w Low Target High Midpoint Error Est. Error f (P; P ) f (Pt ; P )

0.748 0.75 0.752 0.750 0.01% 0.22% 4.35E-06 8.52E-06

k 2.62E-08 3.16E-08 5.17E-08 3.90E-08 23.41% 32.79%

a 5.39E-04 7.02E-04 1.25E-03 8.95E-04 27.44% 39.84%

1.082 1.100 1.112 1.097 0.19% 1.37%

Re Kb 3.34E+09 4.25E+09 4.48E+09 3.91E+09 8.13% 14.60%

Re 1.57E+09 1.74E+09 2.24E+09 1.91E+09 9.61% 17.40%

Table 20: Phase 3 results for Problem 75w, The objective function value is lowest of the nine results.

Problem 79w Low Target High Midpoint Error Est. Error f (P; P ) f (Pt ; P )

0.788 0.79 0.790 0.789 0.17% 0.12%* 5.07E-06 5.23E-06

k 5.85E-08 8.64E-08 1.24E-07 9.12E-08 5.60% 35.85%

a 6.52E-04 1.02E-03 1.65E-03 1.15E-03 13.18% 43.43%

1.064 1.070 1.086 1.075 0.53% 1.04%

Re Kb 2.44E+09 2.45E+09 2.94E+09 2.69E+09 9.52% 9.35%

Re 1.01E+09 1.00E+09 1.04E+09 1.02E+09 1.91% 1.52%*

Table 21: Phase 3 results for Problem 79w, The objective function value is lowest of the nine results.

48

Problem 83w Low Target High Midpoint Error Est. Error f (P; P ) f (Pt ; P )

0.800 0.83 0.847 0.824 0.77% 2.81% 8.15E-06 8.63E-06

k 4.04E-09 7.93E-09 8.82E-09 6.43E-09 18.89% 37.22%

a 7.11E-04 1.46E-03 1.57E-03 1.14E-03 22.09% 37.60%

1.006 1.053 1.043 1.024 2.76% 1.80%*

Re Kb 2.35E+09 2.39E+09 2.45E+09 2.40E+09 0.27% 2.05%

Re 9.67E+08 9.79E+08 9.94E+08 9.80E+08 0.13% 1.38%

Table 22: Phase 3 results for Problem 83w, The objective function value is lowest of the nine results.

Problem 87w Low Target High Midpoint Error Est. Error f (P; P ) f (Pt ; P )

0.864 0.87 0.876 0.870 0.03% 0.70% 5.94E-06 6.25E-06

k 7.84E-08 3.43E-07 6.55E-07 3.67E-07 7.02% 78.61%

a 3.44E-04 1.66E-03 2.61E-03 1.48E-03 10.85% 76.71%

1.009 1.032 1.052 1.030 0.15% 2.10%

Re Kb 7.12E+08 2.09E+09 4.44E+09 2.57E+09 23.06% 72.33%

Re 8.14E+08 8.56E+08 1.27E+09 1.04E+09 21.86% 21.95%

Table 23: Phase 3 results for Problem 87w, The objective function value is lowest of the nine results.

Problem 91w Low Target High Midpoint Error Est. Error f (P; P ) f (Pt ; P )

0.886 0.91 0.925 0.905 0.52% 2.14% 9.86E-06 4.74E-06

k 2.70E-08 3.50E-08 5.51E-08 4.10E-08 17.23% 34.31%

a 1.50E-03 2.00E-03 2.96E-03 2.23E-03 11.36% 32.78%

0.999 1.029 1.034 1.016 1.24% 1.76%

Re Kb 2.89E+08 7.01E+08 9.78E+08 6.33E+08 9.71% 54.40%

Re 1.42E+08 2.87E+08 5.43E+08 3.43E+08 19.50% 58.48%

Table 24: Phase 3 results for Problem 91w, The objective function value is lowest of the nine results.

49

Problem 71w 75w 79w 83w 87w 91w

Error Est. Err Error Est. Err Error Est. Err Error Est. Err Error Est. Err Error Est. Err

0.12% 0.22% 0.00% 0.19% 0.14% 0.12%* 0.94% 1.27% 0.03% 0.30% 0.56% 0.99%

k 1.83% 8.34% 21.20% 30.41% 1.17% 36.14% 23.30% 17.71%* 2.57% 32.76% 23.52% 17.14%*

a 5.41% 9.49% 22.86% 38.40% 4.10% 44.77% 25.94% 17.87%* 4.15% 27.61% 17.41% 16.44%*

0.55% 1.10% 0.34% 1.19% 0.32% 0.99% 2.78% 0.87%* 0.33% 0.90% 1.24% 0.86%*

Re Kb 0.71% 11.27% 9.19% 12.24% 7.67% 9.00% 0.02% 0.99% 6.86% 39.87% 21.12% 26.49%

Re 3.44% 3.61% 6.34% 15.15% 1.64% 1.44%* 0.02% 0.66% 7.10% 13.19% 34.73% 22.22%*

Average Worst Table 25: Phase 3 percentage errors when using mean values for Problems 71w, ..., 91w. Estimated errors are calculated from 95 percent con…dence intervals. function value was within a factor of 2 of the lowest value and used these to …nd a 95% con…dence interval for the mean. The result is shown in Table 25. Instances of underestimation, indicated by *, were more common, however only the underestimation of the error for the structure factor in Problem 83w was severe. On the other hand the overestimations of the error were less severe than with minimum/maximum/midpoint approach and on the whole better characterize the actual errors.

5

Summary and conclusions

Computation of the acoustic pressure …eld by contour integration of the Green’s function appears feasible up to about 500 kHz, but the times required increase with frequency. Beyond 500 kHz the integrands are so oscillatory that computation times are prohibitively long. That acoustic pressure can be computed this way requires that all "strongly" propagating eigenvalues lie within the contour used. We have made an empirical case for this. Table 26 summarizes the results of the various inversion schemes considered. The time to calculate a single objective function value varied with the trial porosity, but was generally on the order of 30 seconds. Thus the inversion times were long, ranging from six to nine hours for the Phase 2 results. The Phase 3 inversion times were quite long because the simplex method was allowed to run to convergence in part in order to con…rm the Phase 2 results; it probably could be limited to a …xed number of iterations, but this was not tested. As noted earlier the Phase 3 results indicate that the di¤erences between the Phase 2 and Phase 2 Alt outcomes were fortuitous in the sense that the more thorough Phase 3 Obj had a higher error for porosity than Phase 2, and the low Phase 2 50

Variant Phase 2 Phase 2 Alt Phase 3 Obj Phase 3 Mid Phase 3 Mean

Error Avg. Worst Avg. Worst Avg. Worst Avg. Worst Avg. Worst

0.36% 1.13% 0.58% 2.41% 0.53% 2.02% 0.30% 0.77% 0.30% 0.94%

k 10.33% 25.75% 15.47% 26.09% 13.99% 30.24% 12.52% 23.41% 12.27% 23.52%

a 12.13% 28.94% 17.54% 29.46% 15.50% 34.10% 15.87% 27.44% 13.31% 25.94%

0.85% 2.02% 1.54% 5.87% 1.31% 4.51% 1.02% 2.76% 0.93% 2.78%

Re Kb 5.41% 17.59% 3.50% 6.92% 6.11% 17.59% 9.57% 23.06% 7.60% 21.12%

Re 7.62% 28.59% 4.77% 12.42% 7.58% 28.59% 9.60% 21.86% 8.88% 34.73%

Table 26: Summary of the average and worst errors made by …ve variants of the inversion algorithm for Problems 71w,...,91w.: Two phase algorithm with univariate Phase 1; Two phase algorithm with alternative Phase 1; Three phase algorithm using the result with the lowest objective function value; Three phase algorithm using the midpoint of the lowest and highest values; Three phase algorithm using the mean. See the text for greater detail. Alt errors for the moduli could not be replicated in Phase 3. The group of nine computers used to conduct these simulations was available only for a limited time and thus it was not possible to conduct further simulations. Perhaps the most interesting question is whether the Phase 2 algorithm can be modi…ed so that it will give useful error estimates, thereby eliminating the need for Phase 3, which requires several more hours of computation time. Another question is whether three or …ve computers could produce comparable results. It should be pointed out that the number of logical computers need not be the same as the number of actual computers, but for the sake of e¢ ciency the number of logical computers should be an integral multiple of the number of actual computers. Thus three computers should be able to simulate nine logical computers with about a trebling of the runtime.

References [1] J. G. Berryman. Con…rmation of biot’s theory. Appl. Phys. Lett., 37:382– 384, 1980. [2] M.A. Biot. Theory of propogation of elastic waves in a ‡uid-saturated porous solid. i. lower frequency range. J. Acoust. Soc Am., 28(2):168–178, 1956. [3] M.A. Biot. Theory of propogation of elastic waves in a ‡uid-saturated porous solid. ii. higher frequency range. J. Acoust. Soc. Am., 28(2):179–191, 1956. [4] J. L. Buchanan, R. P. Gilbert, and K. Khashanah. Determination of the parameters of cancellous bone using low frequency acoustic measurements. J. of Computational Acoustics, 12(2), 2004. 51

[5] J.L. Buchanan, R.P. Gilbert, and K. Khashanah. Recovery of the poroelastic parameters of cancellous bone using low frequency acoustic interrogation. In A. Wirgin, editor, Acoustics, Mechanics, and the Related Topics of Mathematical Analysis, pages 41–47. World Scienti…c, 2002. [6] A. Hosokawa and T. Otani. Ultrasonic wave propagation in bovine cancellous bone. J. Acouist. Soc. Am., 101:558–562, 1997. [7] M. L. McKelvie and S. B. Palmer. The interaction of ultrasound with cancellous bone. Phys. Med. Biol., 10:1331–1340, 1991. [8] J.L. Williams. Prediction of some experimental results by biot’s theory. J. Acoust. Soc. Am., 91:1106–1112, 1992. [9] B. Yavari and A. Bedford. Comparison of numerical calculations of two biot coe¢ cients with analytical solutions. J. Acoust. Soc. Am., 90:985–990, 1991.

52