Ultrasonic wave propagation in heterogeneous elastic and poroelastic media

(1)

Ultrasonic wave propagation in heterogeneous elastic and

poroelastic media

(2)

(3)

Ultrasonic wave propagation in heterogeneous elastic and

poroelastic media

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus Prof.dr.ir. J.T. Fokkema, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op dinsdag 18 januari 2005 om 15:30 uur

door

Jeroen JOCKER

mijnbouwkundig ingenieur geboren te Sliedrecht

(4)

Dit proefschrift is goedgekeurd door de promotor: Prof.dr.ir. J.T. Fokkema

Toegevoegd promotor: Dr.ir. D.M.J. Smeulders

Samenstelling promotiecommissie: Rector Magnificus, voorzitter

Prof.dr.ir. J.T. Fokkema, Technische Universiteit Delft, promotor

Dr.ir. D.M.J. Smeulders, Technische Universiteit Delft, toegevoegd promotor Prof.dr.ir. C.P.A. Wapenaar, Technische Universiteit Delft

Prof.dr. S.M. Luthi, Technische Universiteit Delft

Prof.dr.ir. M.E.H. van Dongen, Technische Universiteit Eindhoven Prof.dr. J.A. Trampert, Universiteit Utrecht

Dr. D.L. Johnson, Schlumberger-Doll Research, USA

Copyright c° 2005 by Jeroen Jocker ISBN 90-9019017-1

Printed by Febodruk

(5)

The answer is hidden in the question

(6)

7 Conclusions 141 A Definition of the temporal and spatial Fourier transforms 145 B High-frequency limits of the Fr´echet kernels 147 C Components of matrices A1 and A2 149 D Components of the single-layer matrix propagator 153 E Kennett’s reflectivity scheme 157 F Conventional poroelastic parameter determination 159 F.1 Porosity φ and grain density ρs . . . 159

F.2 Constrained and shear moduli Kp and µ . . . . 159

F.3 Steady-state permeability k0. . . 160

(8)

CONTENTS iii

Samenvatting 163

Bibliography 165

Dankwoord 173

(9)

(10)

Abstract

The influence of small- and large-scale heterogeneities on an incident wave-field was investigated in ultrasonic experiments. The wave-field generated by ventional piezoelectric transducers was investigated separately. It was con-firmed that the far-field diffraction was correctly predicted by Fraunhofer the-ory, but that the near-field diffraction exhibited significant deviations from the Fresnel theory. These deviations are attributed to the assumption that the source oscillates uniformly over its surface, which seems unrealistic for the actual transducers.

Small rubber and teflon spheres with velocities in the range of the sur-rounding water acted as small-scale heterogeneities in an otherwise homo-geneous medium. The influence on amplitude and traveltime was measured with a needle hydrophone, a focused transducer was utilized for wave gener-ation. These experiments successfully validated first-order scattering theory where a Maslov correction is applied for non-linear effects resulting from a strong contrast between sphere and background medium. Ray theory was found to be unsuitable to properly predict the measured phenomena.

Large-scale heterogeneities were also investigated. Ultrasonic reflection and transmission experiments at different angles of incidence were carried out on plane-parallel layers of aluminum. Here the signal was generated by a conventional plane piezoelectric transducer source and the recording was by means of an identical receiver. Unlike in the conventional bulk-wave transit-time method where successive reverberations within the sample are separated in time, we captured all direct and scattered energy by means of a special recording technique where the individual traces are added and compared with traces recorded in absence of the sample. We found that this recording tech-nique gives much improvement over single-position measurements at normal incidence and is ideally suited for full wavefield analysis. Excellent agreement between experiment and Thomson-Haskell theory was obtained.

(11)

plane-par-2 Abstract

allel layers of poro-elastic materials. The samples were made out of sintered glass powder and had different grain sizes. The Thomson-Haskell theory was extended and new analytical expressions for the reflection and trans-mission coefficients were derived. These expressions were validated through comparison with the Kennett reflectivity method and stability regimes were determined. Effective medium parameters were obtained through inversion of transmission measurements. These effective parameters were found to be in good agreement with the conventional laboratory values although there was a systematic error in the permeabilities. Both the reflection and transmission experiments were generally found to be in good agreement with theoretical predictions although significant deviations were found for samples with larger grains and for larger angles of incidence. In the former case, deviations from continuum Biot theory are observed, whereas in the latter case the effect of reflections from the edges of the samples played an important role.

(12)

List of symbols

Latin

A Frequency-dependent complex amplitude term [N/m]

C Critical exponential instability value [-]

D Directivity term [-]

F Fresnel parameter [-]

R Reflection coefficient [-]

T Transmission coefficient [-]

A Chapter 2: Observation point with respect to source

Chapter 3: Wavefield amplitude [s]

A0 _{Projection of observation point A on the source plane}

a Chapter 2: Source radius [m]

Chapter 5: Characteristic pore scale [m]

A1, A2 Symmetric matrices which yield dispersion relations

for the compressional and shear waves

A, P , Q, R Generalized elastic constants from Biot theory [Pa]

A0_{, Q}0 _{Redefined Biot coefficients} _[Pa]

b Receiver radius [m]

b0 Viscous damping factor [Pa·s/m2]

B Viscous correction factor [-]

B Composition matrix L times the diagonal phase-advancement matrix

C Cross-correlation function

cf Fluid velocity [m/s]

cp, cs Compressional and shear wave velocities (elastic) [m/s] cp1 Fast compressional wave velocity (poroelastic) [m/s] cp2 Slow compressional wave velocity (poroelastic) [m/s]

csh Shear wave velocity (poroelastic) [m/s]

d Layer thickness [m]

(13)

4 List of symbols

eij Strain components for the solid [-] El, G Definitions adopted for notational convenience (Ch.5) [-] Ej Phase-advancement diagonal matrix

ex, ez Unit vectors in the x- and z-direction

f Frequency [Hz]

f First-order array containing six wavefield components

F Definition adopted for notational convenience [Pa2_]

Fl Compressibility factor (Ch.5) [Pa] G Green’s function

G0, G1 Dry, respectively saturated weight of porous sample [N]

g Gravitational constant [m/s2_]

H Total matrix propagator

Chapter 4: Single layer matrix propagator Hj Single layer Haskell matrix propagator (Ch.5) Hij Components of the total matrix propagator hij Components of the single layer matrix propagator H,G,E,F,P,Q Definitions used by Folds and Loggins (1977) (Ch.4)

Hl Slowness (Ch.5) [s/m]

I Identity matrix

Il Intensity [kg/s3m4] K2D

∆t 2D Fr´echet kernel for traveltime residuals

K2D

∆A 2D Fr´echet kernel for amplitude variations

K3D

∆t 3D Fr´echet kernel for traveltime residuals

K3D

∆A 3D Fr´echet kernel for amplitude variations

Ks Bulk modulus of the solid [Pa] Kf Bulk modulus of the fluid [Pa] Kb Bulk modulus of the porous drained matrix [Pa] Kp Constrained modulus [Pa]

k Wavenumber [1/m]

kx Horizontal wavenumber [m−1] k0 Permeability [m2]

L Source-receiver distance [m]

LF Width of the first Fresnel zone [m]

L Matrix to compose the wavefield array f from the plane-wave array ϕ

L1, L2 Submatrices of composition matrix L

M, N Definitions used by Brekhovskikh (1960) (Ch.4)

M Number depending on pore geometry (Ch.5) [-]

M22, M55 Terms introduced for notational convenience (Ch.5) [-]

M52 Term introduced for notational convenience (Ch.5) [Pa·s/m]

M25 Term introduced for notational convenience (Ch.5) [m/Pa·s]

(14)

5

O Origin

p Chapter 2: Acoustic pressure [Pa] Chapter 5: Fluid pressure [Pa]

P Wavefield (Ch.3) [-] P Poynting energy flux vector

ql Vertical wave slowness [s/m]

r Positional vector [m]

r Indicating direction normal to source axis (radial)

R Distance between point A and surface element dS [m]

R1, R2 Composite trace respectively with and without

the presence of a sample (reflection configuration)

Rj Total reflection matrix of the stack with

topmost interface zj

R+,−_j Downward and upward reflection matrix

rf Intrinsic pore-fluid resistivity [Ω] rs Intrinsic resistivity of the saturated porous solid [Ω] S Chapter 2: Source surface [m2_]

Chapter 3: Source wavelet

sl Complex phase slowness [s/m] sx Real-valued horizontal wave slowness [s/m]

t Time [s]

tw Arrival time of water signal [s] ts1, ts2 First and second arrival time respectively [s] tm Time between multiples [s] T1, T2 Composite trace with and without

the presence of a sample (transmission configuration)

T Temperature [◦_C] Tj Total transmission matrix of the stack with

topmost interface zj

T_j+,− Downward and upward transmission matrix

u Averaged local displacement in the solid [m] U Averaged local displacement in the fluid [m] v Averaged local velocity in the solid [m/s] w Averaged local velocity in the fluid [m/s]

V0 Velocity of piston surface [m/s]

v0 Reference velocity [m/s]

Vb Bulk volume of porous sample [m3]

Yp1,p2 Impedance (Ch.5) [Pa·s/m]

Zf Acoustic fluid impedance [Pa·s/m]

(15)

6 List of symbols

Greek

α Angle of incidence [deg.]

αc Critical angle of incidence for the compressional wave [deg.] αp Angle of refraction for the compressional wave [deg.]

αs Angle of refraction for the shear wave [deg.]

α∞ Tortuosity [-]

βl Fluid to solid amplitude ratio [-]

γl Definition adopted for notational convenience (Ch.5) [-]

δij Kronecker symbol [-]

δl Goniometric arguments (Ch.5) [-]

εij Strain components for the fluid [-]

ε Dilatation of the fluid in the pores [-]

ζz Vertical volume flux velocity [m/s]

ζ (no subscript) Dimensionless (inverse) distance [-]

η Fluid viscosity [Pa·s]

λ Acoustic wavelength [m]

Chapter 4: Lam´e coefficient [Pa]

λl Adapted Lam´e coefficients (Ch.5) [Pa]

µ Shear modulus [Pa]

ξl Inverse squared wave velocity [s2/m2]

ρ11, ρ12, ρ22 Density terms used in Biot theory [kg/m3]

ρs Solid density [kg/m3]

ρf Fluid density [kg/m3]

σij Intergranular stress [Pa]

σ Maslov index [-]

τij Force per unit bulk area applied

to the solid part of a unit cube [Pa]

τ Total normal tension force per unit bulk area

applied to the fluid part of a unit cube [Pa]

φ Porosity [-]

Appendix F: Core diameter [m]

φeff Effective porosity term used in the definition of

the generalized elastic Biot constants (Ch.5) [-] φl,φ12,φ32,φ13 Porosity terms used for notational convenience (Ch.5) [-]

φ+_{, φ}− _{Partitioned forms of down- and upgoing plane-wave} amplitude arrays

ϕ Chapter 3: Phase [-]

Chapter 5: First-order plane-wave amplitude array

ϕ+,−_l Plane-wave amplitude [m2_]

χ Formation resistivity factor [-]

ω Radial frequency [Hz]

∆s, ∆a Characteristic determinants for (anti-)symmetric guided waves

Λ Weighted pore-volume to surface-area ratio [m]

Φl Short-hand notation for plane-wave amplitudes including phase-advancement term

(16)

7

Superscripts

ˆ Variable transformed to the wavenumber-frequency domain ˜ Variable transformed to the space-frequency domain

+ Downgoing - Upgoing

T Transposed

Subscripts

B Born approximation indicator

e With respect to the source edge

f With respect to the farthest edge of the source

ij Chapter 4: Matrix coefficient indicators

Chapter 5: Directional placeholders in cartesian coordinate system

kk Repeated indices indicating summation is assumed

l Chapter 5: Placeholder for p1, p2, or sh

Poynting energy vector: Placeholder for p1, p2, sh, or f

m Multiple

n With respect to the nearest edge of the source

p1, p2, sh Indicators for the fast, the slow, and the shear wave

pert Perturbed

ray With reference to ray theory ref Reference

R Rytov approximation indicator

r Chapter 2: In radial direction Chapter 3: Receiver property

s Source property

scat With reference to scattering theory

xyz Direction indicators in cartesian coordinate system 0 Chapter 3: Reference wavefield indicator

(17)

8 List of symbols

Other

∆ Difference

∂ Partial derivative

i Imaginary number

< Real part of a complex number

= Imaginary part of a complex number

J0, J1 Bessel functions of the first kind, zeroth and first order

respectively

1D, 2D, 3D Indication of one-, two-, or three dimensions respectively (x, y, z) Cartesian coordinates

(r, φ, z) Cylindrical coordinates (ρ, θ, φ) Spherical coordinates

< > Surface average

∇ Nabla operator

(18)

Chapter 1

Introduction

1.1 Background

The structure of the earth has been studied extensively using seismic waves generated by natural earthquakes, volcanic eruptions and man-made sources. However, the physical mechanisms controlling the intrinsic attenuation of sedimentary rock throughout the seismic band of frequencies (roughly 1 to 104 _{Hz) are still not fully understood. Often more attenuation is measured} than can be explained on the basis of existing theoretical models (Pride et al., 2004). Intrinsic loss is commonly represented by the inverse quality factor Q−1_{which quantifies the portion of wave energy dissipated (i.e., lost to heat)}

in each wave period. For transmission experiments (earthquake recordings, vertical seismic profiling, crosswell tomography, sonic logs), the total attenua-tion inferred from the recordings can be decomposed as Q−1_total= Q−1_scat+ Q−1_,

where Q−1_scat stands for the effective ”scattering attenuation”. In transmis-sion and reflection experiments, multiple scattering transfers energy from first-arrival pulses into the so-called coda (later arrivals) and into directions that are usually not recorded on the seismogram. Observations have shown that these coda are incoherent waves scattered by randomly distributed het-erogeneities having random sizes and contrasts of physical properties. The characteristic scale of the heterogeneity that has the most influence on a given wave is not always much longer than but is sometimes of the same order of magnitude as the wavelength of the seismic wave (Sato and Fehler, 1998).

In this thesis, these different scales of heterogeneity will be investigated in laboratory ultrasonic transmission and reflection experiments. As we use much higher frequencies than in the seismic range (typically 105 Hz and higher), the scatterers have to be scaled accordingly. If we want the scale of

(19)

10 Introduction

the heterogeneity to be much larger than the wavelength, we use sample layers of a few centimeters thickness. For other scattering experiments, spheres with a diameter of a few millimeters are used. All samples have velocity contrasts with the surrounding medium (water).

In all reflection and transmission experiments our aim is to capture all direct and scattered energy by means of a special recording technique where the individual traces are added and compared with traces recorded in ab-sence of a scatterer. Unlike the conventional bulk-wave transit-time method, this technique has the advantage that there is no requirement to separate the signals corresponding to successive reverberations within the sample, or to identify the wave arrivals corresponding to the different bulk waves. There-fore this method can be used on thinner samples and it is more applicable to attenuative (porous) materials where echo separation is difficult. We made no attempt to separate intrinsic from scattering losses, yet aimed at mak-ing a full wavefield comparison between experiment and theory and inferrmak-ing elastic and poro-elastic properties.

For the spherical scatterers no intrinsic attenuation was taken into ac-count (chapter 3), whereas the layer scatterers were modelled without at-tenuation (chapter 4), and as Biot media, i.e., a viscous atat-tenuation mecha-nism at the macroscale was assumed (chapter 5). The layer scattering was modelled through the so-called Thomson-Haskell method, where all rever-berations and wave-type conversions are included in a so-called matrix prop-agator. Obviously, this matrix becomes more complicated for porous layers due to the additional wave conversions caused by the presence of the slow compressional wave (chapter 5). For the spherical objects, the scattering was modelled as a first-order Rytov approximation with a Maslov correction for non-linear effects due to the strength of the contrast between object and background. In the ultrasonic experiment the limitations of the so-called ray theory (high-frequency limit) are easily illustrated.

An overview of heterogeneity scales is presented in figure 1.1 (after Aki and Richards, 1980) which graphically classifies the different types of wave propagation problems for a range of anomaly scales (represented by the prod-uct ka) and number of wavelengths travelled (represented by the prodprod-uct kL, where L is the travel distance between source and receiver). The upper left corner of the diagram indicates wave problems in which the correlation distance of inhomogeneity is larger than the travel distance from source to receiver, i.e., a > L. This type of problems belongs by definition to the homo-geneous region. The opposite situation occurs when the inhomogeneity scale is much smaller than the wavelength. Then, scattering may become

(20)

negligi-1.2 Literature review 11 1 10 100 1000 10000 0.01 0.1 1 10 100 1000

Effective medium theory Homogeneous

Full waveform modelling (numerical methods)

a=L

kL (wave number x travel distance)

ka

(wave number x correlation distance)

Figure 1.1: Classification of scattering problems and applicable methods in the ka−kL diagram, where k is the wavenumber, a is the correlation distance of inhomogeneity, and L is the travel distance. After Aki and Richards (1980), p.749.

ble and the medium under observation behaves as a homogeneous body with some averaged properties (effective medium theory). In between, scattering may occur and full waveform modelling is required.

1.2 Literature review

An almost universal feature of sedimentary rocks is that they tend to be lay-ered and heterogeneous at many scales. Heterogeneities such as variations in lithology, porosity, permeability and stiffness impact the amounts and pro-ducibility of reservoir fluids. Consequently, imaging and characterization of these inhomogeneities is the essential problem of reservoir characterization (Marion et al., 1994). In practice, many types of inhomogeneities are en-countered with different geometrical and compositional attributes such as salt domes, turbidites, and igneous intrusions (Press and Siever, 1985). In literature, the concepts of layered and spherical scatterers are often applied to study wave propagation through inhomogeneous media.

There have been numerous theoretical, computational, and experimental studies concerned with wave propagation in stratified media by authors such as Postma (1955), Backus (1962), O’Doherty and Anstey (1971), Berryman

(21)

12 Introduction

(1979) and Marion et al. (1994). Some have dealt with the effective properties of the layered medium, while others have focused on the velocity dispersion caused by ”scattering” effects. Here it is understood that scattering refers to the transfer of energy from first-arrival pulses into coda (later arrivals). Wave propagation in stratified media can be explained by ray theory, effective medium theory, or scattering theory depending on the scales of wavelength and layer spacing. To solve the problem of wave propagation and scatter-ing through stratified media, Thomson (1950) and Haskell (1953) introduced a method to compute a global transfer (”scattering”) matrix that connects stresses and displacements through the layer interfaces. This method was improved and modified by many authors (Nayfeh and Chimenti, 1988, 1991; Kundu, 1988; Datta et al., 1988; Hosten and Castaings, 1993) to deal with wave propagation through composite materials. Measurements were pre-sented by Barnard et al. (1975) and Folds and Loggins (1977). Another technique to describe the multiple scattering of layered media is the reflec-tivity method (Kennett and Kerry, 1979) that builds up the reflection and transmission matrices by simple yet exact intuitive reasoning (Pride et al., 2002). In this thesis the scattering matrices are derived for poro-elastic lay-ered media that satisfy the open-pore boundary conditions, and ultrasonic measurements are used for verifying the expressions.

Scattering from spherical objects is often discussed in terms of Foldy’s ef-fective medium theory (Foldy, 1945) of wave propagation through a random collection of a large number of point scatterers. Scattering of elastic waves by a spherical inhomogeneity has been reviewed by Yamakawa (1962). An exact solution of the scattering by a spherical inhomogeneity in a fluid-saturated porous medium was given by Berryman (1985) and Norris (1985). Gurevich et al. (1992) used the so-called Born approximation to obtain explicit ana-lytical results for moderate frequencies and for inclusions with small relative contrast of properties. Early measurements on spheres in water were reported by Farran (1951) and Hickling (1962, 1964). In this thesis we measure the ultrasonic scattered water signal from spherical inhomogeneities that have a small velocity contrast with the surrounding water. This allows a direct comparison with approximate scattering theory and the possibility to test the deviations from classical ray theory.

In ultrasonic experiments there are mainly two reasons for discrepancies between theoretical predictions and experimental results. The first one is the effect of the limited beam when the waves are generated by conventional ultrasonic transducers of diameters typically 25 mm or less. In this case, plane-wave theory does not provide an accurate model for the experiment

(22)

1.3 Thesis aim and outline 13

when measuring layered media, since multiple reflections at oblique angles displace energy along the layer so that it is not received by the finite re-ceiver. Different ways to overcome this limitation were discussed by Hosten and Castaings (1993) and Cawley and Hosten (1997). One approach which consists of scanning of the wavefield and adding the recorded traces, will be discussed and applied in this thesis.

Intrinsic attenuation of propagating waves in the media establishes the second reason for deviations from theory. For the elastic media in this thesis (made of aluminum), we could not measure any appreciable attenuation so that only scattering was taken into account. For the poroelastic media in this thesis, we simply used Biot’s theory (1956) where the dynamic fluid-solid interaction was described following Johnson et al. (1987). The input parameters for the Biot-Johnson theory were obtained from an inversion pro-cedure on transmission experiments and compared with the ones that were obtained from independent laboratory experiments. In the ultrasonic ex-periments, the measured attenuation was due to a combination of intrinsic attenuation and scattering. This is also the case when sound-absorbing me-dia are investigated, which often consist of high-porosity foams bonded onto a rigid impervious backing (Brouard et al., 1996). Also in this case, matrix propagator methods are used to predict the reflection coefficient and acoustic impedance of porous layers (Lauriks et al., 1990).

1.3 Thesis aim and outline

The aim of this thesis is to develop and validate models for wave propagation in heterogeneous media. Inhomogeneities consist of macro-layering where the layer spacing is much larger than the wavelength, and of spherical anomalies of the order of the wavelength where the velocity contrasts are small so that approximate scattering theory can be used. Ultrasonic reflection and trans-mission measurements are carried out using conventional plane and focused transducers.

Data from ultrasonic experiments will be presented in this thesis. For this reason, we will first investigate and discuss the wavefields emitted by cylindrically shaped, flat-faced piezoelectric transducers. This is the subject of chapter 2 where we model transducers through the concept of a baffled piston. Several wavefield aspects are examined and compared with experi-mental measurements, such as the on-axis and off-axis pressure distributions and the Bass approximation for identical source-receiver combinations.

(23)

14 Introduction

on wave propagation. An ultrasonic experiment was designed to record local changes in a wavefield which was scattered by a small sphere. The sphere represents either a positive or a negative velocity anomaly compared with the surrounding water. Using needle hydrophones and a focused transducer source, we are able to determine the effects of the anomaly on wavefield attributes such as traveltime and amplitude. The experimental results are compared with results from finite-frequency wave diffraction theory based on the first-order Rytov approximation.

An experimental procedure to isolate the normal incidence plane-wave component from an incident field is tested through a series of normal and oblique reflection and transmission measurements on a plate of aluminum in chapter 4. The thickness of the plate equals several wavelengths, i.e., the plate represents a large-scale heterogeneity to the propagating wavefield. It is found that experimental plane-wave reflection and transmission coefficients can successfully be measured when the so-called plane-scan technique is ap-plied. In this technique, traces recorded in a plane normal to the source axis at different equidistant receiver positions are added in time to simulate an infinitely large planar receiver.

In chapter 5, Biot theory and the Thomson-Haskell propagator matrix approach are combined to yield closed-form expressions for the reflection and transmission coefficients of a stratified poroelastic medium. In order to assess the numerical stability of the resulting expressions, the expressions are compared with results from the inherently stable Kennett reflectivity method. It is shown that propagator matrix approaches remain stable for a large and relevant spectrum of frequencies, layer thicknesses and permeabilities.

In chapter 6, data from ultrasonic experiments on porous slabs will be presented and a data inversion process to obtain poroelastic parameter values such as wave velocities, thickness of the slab, permeability, and tortuosity will be discussed. The acoustically determined values are compared with values obtained through independent laboratory measurements. Additionally, we will discuss a single case of acoustic measurements performed on a two-layered medium.

(24)

Chapter 2

Wavefield characterization

2.1 Introduction

An important concept in ultrasonics is the circular piston radiator. It is rep-resentative for many ultrasonic transducers and the concept can be applied to study problems such as the far-/nearfield transition and the directivity pattern. In this chapter the wavefield of a conventional ultrasonic trans-ducer is described and compared with experimental results. Transtrans-ducers and hydrophones are used for the recordings.

2.2 Theory

The pressure ˆp in point A as a result of an oscillating surface S is described using the Rayleigh integral of the first kind (Berkhout, 1987):

ˆ p = iωρf 2π Z S V0 e−ikR R dS, (2.1)

where R is the distance between the observation point A and the surface element dS (see figure 2.1). The projection of A on the plane defined by the source surface S is indicated by A0, and z is the distance between A and A0_{. S is considered to be the surface of an uniformly oscillating piston source}

(velocity amplitude V0, radius a) within an infinitely rigid plane (baffle). The fluid in front of the plane has density ρf and velocity cf. Integral (2.1) is

(25)

16 Wavefield characterization

s

A’

a) Side View b) Top view

S O j O A z A’ S V₀ dS R a dS se R e( )j

Figure 2.1: Geometrical variables used for the description of the pressure field of a circular planar piston source with surface S and radius a. The projection A0 _{of the observation point A on the source plane lies inside the} source surface S.

(26)

2.2 Theory 17

Region 1: A0 _{∈ S}

Using projection A0 _{as the origin, equation (2.1) is rewritten in cylindrical}

coordinates (see figure 2.1):

ˆ p = iωρfV0 2π Z _2π 0 Z _σ_e(ϕ) 0 e−ikR R σdσdϕ, (2.2)

where we have used dS = σdσdϕ. The distance from the origin A0 to the edge is defined σe(ϕ). Recognizing that σdσ = RdR yields for equation (2.2):

ˆ p = iωρfV0 2π Z _2π 0 Z _R_e(ϕ) z e−ikRdRdϕ, (2.3)

where Re(ϕ) is the distance from A to the edge. Evaluating the integral over

R yields that ˆ

p = ρfcfV0e−ikz−ρf_2πcfV0 I

e−ikRe(ϕ)dϕ. (2.4) Equation (2.4) states that for region 1, the pressure field is the superposition of a plane-wave and a wave originating from the edge of the source.

Region 2: A0 _{∈ S}_/

If A0 _{lies outside surface S, then equation (2.2) becomes (see figure 2.2):}

ˆ p = iωρfV0 2π Z _ϕ_C ϕB Z _σ f(ϕ) σn(ϕ) e−ikR R σdσdϕ, (2.5)

where the distance to the nearest circle element BC is defined σn(ϕ) and the

distance to the farthest circle element BC is defined σf(ϕ). The points B

and C are defined by circle tangents A0_{B and A}0_{C and ϕ}

B and ϕC are the

corresponding angles (see figure 2.2). Substituting σdσ = RdR into equation (2.5) yields ˆ p = iωρfV0 2π Z _ϕ_C ϕB Z _R_f_(ϕ) Rn(ϕ) e−ikRdRdϕ, (2.6)

where Rn(ϕ) is the distance from A to the nearest edge and Rf(ϕ) is the

distance from A to the farthest edge. Evaluating the integral over R yields that ˆ p = −ρfcfV0 2π Z _ϕ_C ϕB [e−ikRf(ϕ)_{− e}−ikRn(ϕ)_]dϕ. _(2.7)

(27)

s

A’

S

O

j

O

A

z

A’

S

R

_{f( )}

_j

sf

jC

jB

B

C

R

_{n( )}

_j

dS

V

₀

Figure 2.2: Geometrical variables used for the description of the pressure field of a circular planar piston source with surface S and radius a. The projection A0 _{of the observation point A on the source plane lies outside the} source surface S.

Equation (2.7) no longer contains a plane-wave contribution while the near-and far-edge contributions are subtracted. According to equation (2.7), a certain amount of acoustic energy is observed even at observation points for which the projection lies outside the source surface S. This effect is due to edge diffraction; the corresponding waves are called edge waves.

2.2.1 On-axis pressure

The pressure field observed when point A is located on the symmetry axis of the cylindrical source is found by recognizing that, in this case, distance R(ϕ) is always equal to: R =√z2_{+ a}2_{. With this, equation (2.4) can now} be rewritten as ˆ p(z) ρfcfV0 = e −ikz_{− e}−ik√z2_+a2 , (2.8)

which can be evaluated analytically. Setting β =√z2_{+ a}2_{, we obtain} ˆ p(z) ρfcfV0 = e −ik1₂(β+z)h_eik1₂(β−z)_{− e}−ik1₂(β−z)i = 2ie−ik2(β+z)sin · k 2(β − z) ¸ . (2.9)

(28)

2.2 Theory 19

Now, it can easily be seen that ¯ ¯ ¯ ¯_ρp(z)ˆ fcfV0 ¯ ¯ ¯ ¯ = ¯ ¯ ¯ ¯2 sin µ k 2( p z2_{+ a}2_{− z)} ¶¯_¯ ¯ ¯ , (2.10)

so that the maxima are located at z a = a λ(2m + 1) − (2m + 1) 4 λ a, for m = 0, 1, 2, · · · , (2.11) where λ is the wavelength. For relatively small values of λ/a, this may be approximated by z a ≈ 1 2m + 1 a λ. (2.12)

We find that for increasing z, the last maximum is found for m = 0, or: zλ

a2 = F ≈ 1, (2.13)

where we have introduced the Fresnel parameter F. In literature, at distances for which F < 1 the field is often called the nearfield or Fresnel zone. Here, the wavefield is constructed from local spherical waves and has a more or less constant diameter. The resulting wavefront is approximately plane. For distances F > 1, the field is indicated farfield or Fraunhofer zone. Here, the wavefield is constructed from local plane-waves and is divergent. The resulting wavefront is no longer plane.

The transition distance at z = a2_{/λ between the Fresnel and Fraunhofer} zone can also be derived in a simpler way when it is realized that the Fresnel zone indicates the region close to the source where both constructive and destructive interference can occur between different elements on the surface of the radiator. For an on-axis observation point, the maximum travel difference between a wave coming from the edge of the piston and a wave coming from the center of the source equals ∆l =√z2_{+ a}2_{− z. The point beyond} which interference no longer occurs is defined by ∆l ≤ λ/2 which yields that z ≥ a2/λ for z >> a.

In figure 2.3, the absolute value of the normalized pressure is shown as a function of z for a 500 kHz source with a radius of 12.7 mm (solid line). Additionally, three values for the non-dimensional inverse distance parame-ter ζ defined by Bass (1958) as ζ = (k/2)[(z2_{+ 4a}2₎12 − z] are indicated on the horizontal axis. Note that ζ decreases for increasing z. The transition from Fresnel to Fraunhofer zone takes place at z ∼= 53 mm. Also shown in the figure are the experimental results of an on-axis needle measurement.

(29)

20 Wavefield characterization 0 50 100 150 200 250 300 350 400 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 on−axis distance (mm) |p/( ρ f c f V 0 )| ζ=6.459 ζ=2.268 ζ=1.367

Figure 2.3: Theoretical (solid line) and experimental (◦ ◦ ◦) results from on-axis needle hydrophone measurements of a 500 kHz transducer source.

Dz

z r source needle receiver

a) on-axis needle measurement

Dr

z r source needle receiver

b) off-axis needle measurement

Figure 2.4: Schematic overview of the experimental determination of the on-axis and off-on-axis pressure field of a cylindrical source using a needle receiver.

(30)

2.2 Theory 21 0 25 50 75 100 125 150 175 200 0 2 4 6 8 10 12 z ref (mm) Σ (| p(z)/p(z ref )| th −| p(z)/p(z ref )| exp ) 2

Figure 2.5: Summed squared differences between theory and experiment for different values of z_ref.

The wavefield of a 500 kHz flat, piezoelectric transducer (Panametrics, type V301) was scanned using a needle hydrophone (Specialty Engineering As-sociates, type PZT-Z44-1000) which has a frontal diameter of 1 mm. The needle hydrophone is moved along the axis of the source in the z-direction with increments of ∆z = 2 mm (see figure 2.4a). At each needle position, the pressure was recorded and subsequently transformed to the frequency domain. We determined an appropriate normalization pressure in order to be able to compare the data with theory by eliminating the transfer function of the measurement equipment. From equation (2.8) it follows that for the theoretical value ¯ ¯ ¯ ¯_p(z_ˆp(z)ˆ ref) ¯ ¯ ¯ ¯ th = ¯ ¯ ¯ ¯ ¯ e−ikz_{− e}−ik√z2_+a2 e−ikzref− e−ik

√ z2 ref+a2 ¯ ¯ ¯ ¯ ¯. (2.14)

We may also determine |ˆp(z)/ˆp(z_ref)|_exp directly from measurements. We used a least-squares method to determine the appropriate zref (see figure 2.5). We found zref = 66 mm, the corresponding results for this zref are shown in figure 2.3. The pattern of consecutive maxima and minima, characteristic for the Fresnel zone, is clearly visible in both theory and experiment. This pattern is followed by the steady pressure decline in the Fraunhofer zone.

(31)

22 Wavefield characterization 0 2 4 6 8 −10 0 10 20 30 40 50 ζ/π attenuation (dB) z=52mm z=21mm z=9mm

Figure 2.6: Apparent attenuation due to interference and divergence as a function of ζ/π. Apparent attenuation is determined as −20 log [ˆp(z)/ˆp(zref)], where the normalization pressure ˆp(z_ref) is determined using a least-squares method. On-axis needle hydrophone measurements of a 500 kHz transducer source. Solid line: theory; (◦ ◦ ◦): On-axis needle hydrophone measurements.

(32)

2.2 Theory 23

The same data are plotted as apparent attenuation versus ζ/π in fig-ure 2.6. Here, the term ”apparent” relates to the decrease of signal amplitude due to interference in the nearfield, and divergence in the farfield. The atten-uation scale is defined as −20 log [ˆp(z)/ˆp(zref)]. For illustrational purposes, some z-values are indicated in figure 2.6.

At large source-receiver distances, theory is in good agreement with the experimental results. At smaller distances however, deviations from theory are found. The deviations are explained by the fact that we assumed an uniform oscillation of the source, whereas this is not realistic for the ac-tual transducer. The influence of different amplitude distributions on the on-axis pressure profile is investigated in figure 2.8, the tested amplitude dis-tributions are shown in figure 2.7. The analytic, axisymmetric functions for the nonuniform amplitude distributions referred to as ”simply supported” or ”Fermi”, have been used extensively by many authors such as Harris (1981). The expressions for the ”simply supported” and ”Fermi” distributions are given by

A(r) = 1 − (r/a)2, (simply supported) A(r) = [1 + exp(b1r/a − b2)]−1, (Fermi)

where the constant parameters b1 and b2 control the shape of the amplitude distribution. In figure 2.7, b₁ = 25 and b₂ = 20. The pressure profiles in figure 2.8 result from a numerical evaluation of the time domain variant of the Rayleigh integral (2.1) where we assumed V0 = V (r).

From figure 2.8 it becomes clear that in case of continuous wave exci-tation the amplitude distribution of the source plays an important role in the establishment of the on- (and off-)axis pressure profile. This is explained by the observation that interference between two spherical waves emanating from two different surface elements of the source only results in zero total amplitude if the two sources have the same amplitude while oscillating out of phase. If however the amplitude distribution depends on location, then the total on-axis pressure resulting from interference may never reach zero amplitude.

2.2.2 Off-axis pressure

We also performed off-axis measurements with the needle hydrophones. An example is given in figure 2.9 where we show the needle-scan in the (r,t)-domain at an axial source-receiver distance of 20 cm. The wavefield recorded in time and space as shown in figure 2.9 is decomposed into its

(33)

frequency-24 Wavefield characterization 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (r/a) (−)

Normalized amplitude distribution

Fermi Simply supported

Uniform

Figure 2.7: Three different amplitude distributions over the source surface. Solid line: Uniform; Dashed: Fermi (b1 = 25, b2 = 20); Dash-dot: Simply supported.

and plane-wave components using the Fourier transforms (A.1) and (A.3). The resulting amplitude spectrum is shown in figure 2.10. We note that next to the 500 kHz center-frequency component many other frequencies appear in the spectrum. The 1.5 MHz component even appears higher in amplitude than the 500 kHz component which indicates that the needle receiver is more sensitive in the higher frequency range. Also a significant amount of non-normal radial wavenumbers kr appear in the recorded spectrum.

For six different axial positions z (solid lines), the pressure profiles are shown in figure 2.11. These curves result from a numerical evaluation of equa-tions (2.4) and (2.7), where Simpson’s rule is used to compute the integrals. The experimental data result from a needle-scan in the radial direction where ∆r = 0.5 mm (see figure 2.4b). Here, no fitting procedure is applied. For the normalization pressure we used the on-axis value (at r/a = 0). In figure 2.11 we notice that for z = 20, 50, and 100 mm the theory accurately predicts the measurements. For z = 5 mm the theory predicts more pronounced minima than actually recorded. For z = 10 mm the experimental values are gener-ally underestimated by theory for r/a < 1, and for z = 30 mm there is a serious underestimation by theory at r/a < 0.35. The discrepancies between

(34)

2.2 Theory 25 0 10 20 30 40 50 60 70 0 0.2 0.4 0.6 0.8 1 on−axis distance (mm) normalized amplitude (−)

Figure 2.8: Normalized on-axis pressure for three different amplitude distri-butions for a 500 kHz transducer source. Solid line: Uniform; Dashed: Fermi; Dash-dot: Simply supported.

(35)

Figure 2.9: Radial scan of the wavefield emitted by a 500 kHz flat, cylindrical transducer. Wavefield scanned using a needle hydrophone at an axial source receiver distance of 20 cm. The unit of the numbers indicated on the grayscale is milliVolt (mV).

(36)

2.2 Theory 27

Figure 2.10: 2D Fourier transform of the wavefield emitted by a 500 kHz flat, cylindrical transducer and recorded with a needle hydrophone at an axial source receiver distance of 20 cm (see also figure 2.9).

(37)

28 Wavefield characterization 0 0.5 1 1.5 2 0 0.5 1 1.5 2 |p/( ρ f c f V 0 )| z = 5 mm 0 0.5 1 1.5 2 0 0.5 1 1.5 2 z = 10 mm 0 0.5 1 1.5 2 0 0.5 1 1.5 2 |p/( ρ f c f V 0 )| z = 20 mm 0 0.5 1 1.5 2 0 0.5 1 1.5 2 z = 30 mm 0 0.5 1 1.5 2 0 0.5 1 1.5 2 r/a |p/( ρ f c f V 0 )| z = 50 mm 0 0.5 1 1.5 2 0 0.5 1 1.5 2 r/a z = 100 mm

Figure 2.11: Dimensionless pressure distribution in radial direction, plotted versus the dimensionless parameter r/a. Solid line: theory; o: experiment.

(38)

2.3 Directivity 29

s

A’

S O j dS j’ O A z q A’ r V 0 dS R

Figure 2.12: Geometrical variables used in the derivation of the directivity characteristics of a circular planar piston source.

theory and experiment are again thought to be caused by non-uniformity of the piston radiator.

2.3 Directivity

An interesting result is obtained from the Rayleigh integral (2.1) if the fol-lowing coordinate relations are applied (see figure 2.12). Now we define the surface element dS with respect to the origin O by σ and ϕ. We need addi-tional coordinates (ρ, θ, ϕ0) to define A (and consequently A0). Defining l as the distance between A0 _{and dS, we obtain}

OA0 = ρ sin θ, (2.15)

ρ2 = z2+ (OA0)2, (2.16)

l2 = σ2+ (OA0)2− 2σ(OA0) cos(ϕ0− ϕ), (2.17)

dS = σdσdϕ, (2.18)

R2 = z2+ l2, (2.19)

where R is again the distance between the observation point A and dS, and ρ is the distance between A and origin O (see figure 2.12). Combining equations

(39)

(2.15)-(2.19) yields

R =pρ2_{+ σ}2_{− 2σρ sin θ cos(ϕ}0_{− ϕ).} _(2.20)

Assuming ρ À a ≥ σ (farfield approach) and neglecting all but the first order term of the Taylor expansion of (2.20) yields:

R ≈ ρ − σ sin θ cos(ϕ0− ϕ). (2.21) Substituting equation (2.21) into the exponential term of equation (2.1) yields

ˆ

p = iωρfV0 2πρ

Z

S

e−ik(ρ−σ sin θ cos(ϕ0−ϕ))dS,

= iωρfV0e −ikρ ρ Z _a σ=0 · 1 2π Z _2π ϕ=0

eikσ sin θ cos(ϕ0−ϕ)dϕ ¸

σdσ, (2.22)

where we have used R ≈ ρ for the denominator of (2.1). In equation (2.22), the part between brackets is recognized as a Bessel identity:

ˆ p = iωρ_fV₀e−ikρ ρ Z _a σ=0 J₀(kσ sin θ)σdσ. (2.23) Finally, using the substitution ξ = kσ sin θ we arrive at (Hall, 1987, p. 117)

ˆ p = iωρfV0 e −ikρ ρk2_sin2_θ Z _{ka sin θ} ξ=0 ξJ0(ξ)dξ, = iωρ_fV₀a2e−ikρ ρ · J₁(ka sin θ) ka sin θ ¸ , (2.24)

which is the farfield- or Fraunhofer approximation of the wavefield of a cir-cular piston source with radius a and wavenumber k. According to (2.24), the farfield can be approximated as the product of a point source and a directivity function:

ˆ

p = A(ω)e

−ikρ

ρ D(θ), (2.25)

where A(ω) is the frequency dependent complex amplitude term, while the directivity-term D(θ, ω) is defined as (Hall, 1987)

D(θ) = J1(ka sin θ)

ka sin θ . (2.26)

(40)

dif-2.3 Directivity 31

Figure 2.13: Radiation patterns for different values of ka. Sources become more directive as the product of wavenumber and source radius increases. A: ka = 0.1, B: ka = 5, C: ka = 10, D: ka = 25. The grayscales indicate amplitude levels. Black: -1; Gray: 0; White: +1 (in arbitrary units).

(41)

Freq θ−20dB Fresnel λ in water

(degr) distance (mm) (mm)

500 kHz 7.29 54.55 2.96

1 MHz 3.64 108.98 1.48

2.25 MHz 1.62 245.20 0.66

Table 2.1: Divergence angle, Fresnel distance a2/λ in water (cf = 1480 m/s), and wavelength for three different frequencies. A source radius of a = 1.27 cm was used for the calculations.

ferent values of ka. For ka ¿ 1, the radiation pattern is essentially omni-directional. For larger values of ka, the source becomes increasingly directive with constructive interference in the vicinity of the axis of the source and destructive interference further away and more off-axis. Clearly, at higher frequencies energy stays more confined to the central beam. Representative values for ka of the sources used in our reflection and transmission experi-ments vary from 25 to 60.

For a −20 dB drop in sound pressure it can be shown that J1(ka sin θ−20dB) ka sin θ−20dB = 1 10, (2.27) so that sin θ−20dB∼= 1.089λ/2a, (2.28)

where sin θ−20dB is the divergence angle. Table 2.1 displays the divergence

angle, wavelength in water, and the Fresnel distance a2_{/λ for three different} frequencies. The divergence angle decreases with increasing frequency, again emphasizing that at higher frequencies energy stays more confined to the central axis.

2.4 Identical source-receiver pairs

We also used receivers identical to the source to perform measurements of the emitted field. Figures 2.14 and 2.15 show the field-scans in the time and fre-quency domain respectively for an identical source-receiver pair. We notice that with respect to figure 2.10 the spectrum has become much smoother. The piston receiver spectrum is confined to a much smaller bandwidth in both frequency and wavenumber direction. In this sense, a large planar receiver

(42)

2.4 Identical source-receiver pairs 33

Figure 2.14: Radial scan of the wavefield emitted by a 500 kHz flat, cylindrical transducer. A receiver identical to the source was used. The axial source receiver distance equals 20 cm. The unit of the numbers indicated on the colorbar is milliVolt (mV).

acts as a filter for both frequency and wavenumber components. Obliquely incident plane-wave components are suppressed in favor of plane-waves in-cident at normal angle, and higher frequency components are suppressed in favor of the center-frequency component. Apparently summation of record-ings from different receiver positions results in the enhanced isolation of the normal incidence plane-wave component and the center-frequency compo-nent. We also notice that the spectrum amplitudes are much larger for the identical receiver case than for the needle results. Obviously this is due to the fact that because of its larger surface, a larger receiver records more acoustic energy at its receiver position.

In case of a receiver with finite size, the incident wavefield is essentially averaged over the receiver surface:

hˆpi = R

ˆ pdS

πb2 , (2.29)

where b is the radius of the receiver. Knowing the exact pressure field (see figure 2.11), hˆpi can simply be computed numerically. In the experimental

(43)

Figure 2.15: 2D Fourier transform of the wavefield emitted by a 500 kHz flat, cylindrical transducer. A receiver identical to the source was used. The axial source receiver distance equals 20 cm (see also figure 2.14).

(44)

2.4 Identical source-receiver pairs 35

determination of reflection and transmission coefficients of elastic and poroe-lastic media, we will use two identical transducers (b = a): one as the source, the other as the receiver. For this case, Williams (1951) derived an expression for the acoustic pressure:

hˆpi = ρfcfV0exp[−ikz] −

4ρ_fc_fV0 π

Z _π/2 0

exp[−ikpz2_{+ 4a}2_cos2_{θ] sin}2_θdθ, (2.30) where θ is supposed to be an integration variable without any geometrical significance. The first term on the right hand-side represents the plane-wave solution and the second term represents the influence of edge diffraction. Unfortunately, the integral in the second term of equation (2.30) cannot be evaluated in a closed form. Probably the best approximation of equation (2.30) was given by Bass (1958) who arrived at

hˆpi ≈ ρ_fc_fV0e−ikz[1 − a1(ζ) − a2(ζ)], (2.31) where ζ was previously defined as ζ = (k/2)[(z2_{+ 4a}2₎1₂ _{− z], and}

a₁ = µ 1 − ζ2 2k2_a2 ¶ [J₀(ζ) + iJ₁(ζ)]e−iζ, (2.32) a2 = ζ 2 k2_a2 · iJ1(ζ) ζ ¸ e−iζ. (2.33)

Here, J0 and J1 are Bessel functions of the first kind and zeroth and first order respectively. The approximation (2.31) is within 2.8% accuracy for z > a and within 0.6% for z > 2a (Bass, 1958).

The average pressure is plotted against source-receiver distance z in fig-ure 2.16. In figfig-ure 2.17, this information is plotted as apparent attenuation with respect to a reference pressure, or −20 log |< ˆp(z) > / < ˆp(zref) >| ver-sus ζ/π. Some values for ζ are indicated in figure 2.16, and some z-values are indicated in figure 2.17. Also shown in both figures are the experimen-tal results from an axial scan with an identical receiver. The setup of the experiment is similar to the on-axis needle measurement (see figure 2.4a), ex-cept for the receiver which in this case is identical to the source. Again, the data from the experiment are normalized following the procedure described in section 2.2.1. In this case, the average pressure at z = 52 mm was found appropriate for normalization. In the original paper by Bass (1958) it was suggested to use the value at minimum source-receiver distance for normal-ization. We found that this hardly influences the results. Finally, the dashed

(45)

36 Wavefield characterization 0 100 200 300 400 500 0 0.25 0.5 0.75 1 on−axis distance (mm) |<p>/( ρ f c f V 0 )| ζ=3.375 ζ=1.368 ζ=0.856

Figure 2.16: Theoretical on-axis dimensionless pressure according to Bass (solid line), numerical evaluation of equation (2.29) (dashed line), and on-axis identical receiver measurement (◦ ◦ ◦). The reference distance used to normalize the data is z_ref = 52 mm (500 kHz source-receiver combination).

(46)

2.4 Identical source-receiver pairs 37 0.5 1 1.5 2 2.5 3 3.5 4 −1 0 1 2 3 4 ζ/π attenuation (dB) z=686mm z=171mm z=96mm

Figure 2.17: Apparent attenuation caused by edge wave diffraction as a func-tion of ζ/π (500 kHz source-receiver combinafunc-tion). Apparent attenuafunc-tion is determined as −20 log [ˆp(z)/ˆp(zref)], where the normalization pressure ˆp(zref) is determined using a least-squares method. Solid line: theory; Dashed line: numerical evaluation of equation (2.29); (◦ ◦ ◦): On-axis identical receiver measurement.

(47)

lines in figures 2.16 and 2.17 indicate the numerical evaluation of equation (2.29) where the local pressure is known through equation (2.4), and where Simpson’s rule is used to compute the integral. For the most part, Bass the-ory and the numerical result coincide, thus confirming the accuracy of the Bass theory.

Typical values of ζ for the ultrasonic measurements presented in this re-port are equal to or larger than ζ ≥ 0.7 corresponding to source-receiver distances smaller than 50 cm for a 500 kHz source. From figure 2.17 it ap-pears that here, the apparent attenuation of the field due to diffraction effects from the edge of the source is small. In general, theory is confirmed by the ul-trasonic measurements. The characteristic minima and maxima encountered previously in figures 2.3 and 2.6 are still visible in the theoretical curves though not as pronounced compared with the on-axis case. This is a result of the averaging effect of the finite size of the receiving transducer. The ex-perimental data show even less interference effects. It is believed that this is due to the fact that the assumption that the receiver can be represented as an ideal piston source may not be adequately fulfilled in practice (Bass, 1958).

(48)

Chapter 3

Experimental validation of

first-order wave diffraction

theory

3.1 Introduction

The characteristics of finite-frequency waves propagating in heterogeneous media are fundamentally different from what is predicted by standard ray theory. Ray theory is a high-frequency solution to the wave equation which physically means that energy travels along lines in space (called rays) joining the source and receiver. It is valid to use ray theory to model high-frequency wave propagation in simple media consisting of large-scale velocity struc-tures. On the other hand, the high-frequency approximation in ray theory breaks down for finite-frequency waves propagating in complex media with small-scale velocity structures. In such heterogeneous media, the effect of wavefield scattering becomes important so that the phase and amplitude variation of propagating waves behaves in a way that can not be explained by the geometrical ray approximation.

We present a theory for wavefield propagation in heterogeneous media with sizes of anomalies so small that the conditions for ray theory are not satisfied. The theory for the traveltime and amplitude variations of finite-frequency wavefields takes the single-scattering effect of waves due to small-scale velocity contrasts into account. Practically, we add an extra term for the single-scattering process to the standard point source solution for a trans-mitted wavefield. The finite-frequency wave theory for the traveltime and

(49)

40 Experimental validation of first-order wave diffraction theory

amplitude is then linear with respect to velocity perturbations. In case of triplications in the propagating wavefield, non-linear effects get increasingly dominating. Here, the term ”non-linear” refers to the dependence of the scattered field on the strength of the contrast between the inclusion and the background. By combining the linear wave scattering theory with the Maslov theory for the formation of caustics, it is possible to account for the finite-frequency effect on wavefronts, including triplications to leading order. The wave scattering theory is validated in a 3D ultrasonic wave experi-ment where a spherical ball of teflon or rubber is placed in a reference medium consisting of water. There is a sharp difference between the velocity of the ball and the speed of sound in water. The diameter of the sphere is chosen such that the regime of wave scattering theory is significant. From the ultrasonic data it becomes evident that the properties of the finite-frequency wave the-ory cannot be predicted by ray thethe-ory. For instance, in the single-scattering approach, the maximum traveltime sensitivity to velocity anomalies appears off the geometrical ray, while the amplitude has a significant sensitivity along the ray path. In ray theory, the traveltime shift is only sensitive to velocity contrasts on the geometrical ray path while there is no amplitude variation due to wave scattering effects. Hung et al. (2001) performed a similar 3D numerical wave experiment with a smooth anomaly. In their numerical simu-lation, the acoustic 3D wave equation was solved by a pseudo-spectral method and the smoothness was a numerical requirement. Our experimental setup allowed us to investigate the effect of sharp discontinuities.

The theoretical work for the finite-frequency wave theory is well docu-mented by Aki and Richards (1980), Woodward (1992), Snieder and Lomax (1996) and Spetzler and Snieder (2001a) who applied the Rytov approxima-tion to include the effect of heterogeneous small-scale structures on the phase of propagating waves. In a similar vein, Marquering et al. (1999), Zhao et al. (2000) and Hung et al. (2001) chose to work with the Born approximation in order to account for single-scattering wave effects on the traveltime of propagating waves. Aki and Richards (1980), Woodward (1992) and Dahlen and Baig (2002) focused on the finite-frequency wave theory for amplitude variations. Spetzler et al. (2002a) show an example of testing the wave scat-tering theory using a stochastic approach. On the contrary, the experiments validating finite-frequency wave theory as presented in this chapter are de-terministic which has the advantage of clearly illustrating the properties of the finite-frequency wave theory.

This chapter consists of four parts. First, the theory for traveltime and amplitude variations of finite-frequency wavefields is discussed both in

(50)

two-3.2 Theory 41

and three dimensions. Next, the ultrasonic wave experiment to validate the 3D theory is introduced followed by the section with theory and data com-parison. Finally, conclusions are drawn.

3.2 Theory

In this section, the theory for wave propagation in heterogeneous media with anomalies smaller in size than the Fresnel zone is presented. The wave mod-elling theory is based on the first order Rytov approximation which has a larger validity than the first order Born approximation (Snieder and Lomax, 1996). The Rytov approximation considered here is limited to contrasts in the speed of sound only (i.e., not contrasts in density), and it is limited to fluid inclusions only. As the experiments described in section 3.3 were performed using solid inclusions, this means that density contrasts as well as shear wave effects will be neglected. We will first derive expressions for traveltime shifts and amplitude variations for the 2D case followed by a generalization to three dimensions.

The Rytov wavefield PR(rr, rs, ω) at the angular frequency ω = 2πf

emit-ted from the source position rs and recorded at the receiver position rr is

given by P_R(r_r, r_s, ω) = P₀(r_r, r_s, ω) exp³ PB P0 (r_r, r_s, ω) ´ , (3.1)

which accounts for the single-scattering process of a propagating wavefield in heterogeneous media (Aki and Richards, 1980; Snieder and Lomax, 1996). The reference wavefield inherent to the reference velocity model v0(r) is de-noted P0(rr, rs, ω), while PB(rr, rs, ω) denotes the first order Born wavefield.

The phase and amplitude perturbations of the scattered field with respect to the reference field are derived by separation of the real and imaginary part of the exponential function in equation (3.1). Hence, the phase delay due to the presence of a velocity anomaly is given by

∆ϕ(rr, rs, ω) = = µ PB P0(rr, rs, ω) ¶ . (3.2)

which yields for the traveltime residual ∆t

∆t(rr, rs, ω) = _ω1= µ PB P0(rr, rs, ω) ¶ , (3.3)

(51)

42 Experimental validation of first-order wave diffraction theory

where we have used the linear expression ∆ϕ = ω∆t (e.g. Spetzler and Snieder (2001a)).

In a similar vein, the amplitude AR(rr, rs, ω) of the Rytov wavefield with

respect to the amplitude A0(rr, rs, ω) of the reference wavefield is given by

ln Ã AR A0(rr, rs, ω) ! = < Ã PB P0(rr, rs, ω) ! . (3.4)

Let the amplitude difference be defined as ∆A = A_R− A₀. For small per-turbations of the amplitude of the Rytov wavefield, the logarithm to leading order is ln(1 + x) ≈ x for x ¿ 1. With this in mind, the amplitude variation with respect to the reference wavefield is approximated as

∆A A0 (rr, rs, ω) ≈ < Ã P_B P0(rr, rs, ω) ! , (3.5) for ∆A/A0 ¿ 1.

Aki and Richards (1980) and Snieder and Lomax (1996) derived the first Born wavefield (from now on simply called the Born wavefield) which is a solution of the acoustic wave equation for a constant density and in the presence of a velocity perturbation. The Born wavefield is given by

P_B(r_r, r_s, ω) = Z V 2∆v(r)ω2 v3 0(r) P₀(r, r_s, ω)G(r_r, r, ω)dV, (3.6)

where the velocity perturbation with respect to a homogeneous reference velocity v0(r) = v0 is denoted ∆v(r), and the 2D volume integration is indi-cated by V. The 2D far-field Green’s function for the homogeneous reference medium is given by Snieder and Lomax (1996):

G(r_r, r, ω) = − r v0 8πω | rr− r | exp i(ω | rr− r | v0 + π 4). (3.7)

As was shown previously by Spetzler and Snieder (2001a), equation (3.7) can be used to express the Born wavefield in equation (3.6) explicitly in the system coordinates. For a point source this yields

PB(rr, rs, ω) = _4πvω₂ 0 exp i µ ωL v0 + π 4 ¶ (3.8) × Z _∞ −∞ Z _L 0 ∆v(r)exp i( ωLr2 2v0z(L−z) +π4) p z(L − z) dzdr.

(52)

3.2 Theory 43

Figure 3.1: Illustration of the geometric variables to describe the single-scattering process for a 2D geometry.

A sketch of the system coordinates is shown in figure 3.1. Let the point source be located at the origin, the receiver position is placed on the z-axis at z = L, and the position of the scatterer point is (z, r). The Born wavefield in equation (3.8) is then divided with the point source solution

equation (3.3) and the monochromatic amplitude variation in equation (3.5). Recorded waves are never monochromatic but have broadband frequency properties. For this reason, the monochromatic traveltime delay and ampli-tude variation are integrated over the frequency band [f₀− ∆f ; f₀+ ∆f ]. In turn, the broadband frequency characteristics of recorded wavefields with the normalized amplitude spectrum, that is, Rf0+∆f

f0−∆f A(f )df = 1 are taken into account. Hence, the broadband traveltime shift for source position rs and

receiver position rr is given by

∆t(rr, rs) = Z _∞ −∞ Z _L 0 ∆v(z, r)K_∆t2D(z, r)dzdr, (3.10)

with the frequency averaged sensitivity function (also known as the Fr´echet kernel) for traveltime residuals defined as

K_∆t2D(z, r) = − s L v5 0z(L − z) (3.11) × Z _f_0+∆f f0−∆f A(f )pf sin µ f πLr2 v0z(L − z)+ π 4 ¶ df,

Ultrasonic wave propagation in heterogeneous elastic and poroelastic media

Ultrasonic wave propagation in heterogeneous elastic and

poroelastic media

Ultrasonic wave propagation in heterogeneous elastic and

poroelastic media

Contents

Abstract

List of symbols

Chapter 1

Introduction

1.1

Background

1.2

Literature review

1.3

Thesis aim and outline

Chapter 2

Wavefield characterization

2.1

Introduction

2.2

Theory

s

A’

S

O

j

O

A

z

A’

S

R

j

sf

jC

jB

B

C

R

j

dS

V

2.3

Directivity

2.4

Identical source-receiver pairs

Chapter 3

Experimental validation of

first-order wave diffraction

theory

3.1

Introduction

3.2

Theory

_j

_j