Transient Green's tensor for a layered solid half-space with different interface conditions.

Pulsed ultrasonic techniques can be and have been used to examine the interface conditions of a bonded structure. To provide a theoretical basis for such testing techniques we model the structure as a layer on top of a half-space, both of different elastic properties, with various interface bonding conditions. The exact dynamic Green's tensor for such a structure is explicitly derived from the three-dimensional equations of motion. The final solution is a series. Each term of the series corresponds to a successive arrival of a "generalized ray" and each is a definite line integral along a fixed path which can be easily computed numerically. Willis' method is used in the derivation. A new scheme of automatic generation of the arrivals and ray paths using combinatorial analysis, along with the summation of the corresponding products of reflection coefficients is presented. A FORTRAN code is developed for computation of the Green's tensor when both the source and the detector are located on the top surface. The Green's tensor is then used to simulate displacements due to pulsed ultrasonic point sources of known time waveform. Results show that the interface bonding conditions have a great influence on the transient displacements. For example, when the interface bonding conditions vary, some of the first few head waves and regular reflected rays change polarities and amplitudes. This phenomenon can be used to infer the quality of the interface bond of materials in ultrasonic nondestructive evaluation. In addition the results are useful in the study of acoustic microscopy probes, coatings, and geo-exploration.

Key words: bond integrity; dynamic elasticity; fundamental solution; Green's function; interface condition; layered half-space; NDE theoretical mechanics; wave mechanics.

Available online: http://www.nist.gov/jres

1. Introduction

The dynamic Green's tensor of a structure is the fundamental solution to the transient mechanical wave problem of the structure. The theoretical prediction of the behavior of transient waves for arbitrary extended source distributions in space and time can be obtained once the time-domain dynamic Green's tensor is known. The ability to predict the behavior of transient waves in a structure is important in the development and understanding of nondestructive evaluation techniques using ultrasonics or acoustic emission and in other problems of wave propagation in solid and liquid media.

In recent years, the "generalized ray" expansion technique has been applied to compute the Green's tensor in a solid half-space or a solid infinite plate. The governing differential equation of motion is transferred to the Fourier or Laplace domain and the solution in the form of a series is obtained algebraically. Basically, there are two methods to invert the series from the transformed domain to the time-space domain. One is the well known Cagniard-de Hoop inversion method (1); the other is a method developed by Willis in 1973 (2). Willis' method uses the Fourier transform and expresses the resulting transform as a series of "generalized rays" and then inverts the series term by term. For three-dimensional problems, as in the Cagniard-de Hoop method, only one integration remains. The integration path is always around a unit circle and is therefore "fixed" to some extent, but explicit evaluation of the integrand requires the numerical solution of an algebraic equation for each integration variable. Applicat ion of the Cagniard-de Hoop method requires a detailed discussion of the structure of a moderately complicated algebraic function accompanying the transform of the integration path. However, the integration path of the Willis inversion method is "fixed" and succeeds in avoiding the explicit discussion of the structure of the algebraic function and so is applied rather easily even in the anisotropic case. This is the main advantage of the Willis inversion method. The basis for carrying out the Willis inversion is that the solutions of elastodynamic problems are homogeneous functions of time, t, and position, x; that is, the problems are self-similar (2). In addition, the displacements rather than the potentials are used, so the derivations are simplified, and the derivatives of the Green's tensor about spatial coordinates are easily obtained as well. These derivatives can be interpreted as the displacements due to dipole sources. The three-dimensional transient Green's tensor for an isotropic plate and the two -dimensional transient Green's tensor for an anisotropic layer on an isotropic half-space were solved using Willis' inversion method (3,4).

In this paper formulae are derived for computing the three-dimensional transient Green's tensor for an isotropic layer overlay on an isotropic half-space. To understand the influence of the interface bonding condition on the behavior of transient waves, a welded interface, a liquid coupled interface, and a "vacuum" interface are considered. The results for the case of the liquid coupled interface are obtained by artificially casting the boundary conditions into a matrix form similar to that for the case of the welded interface. The results for the "vacuum" interface are obtained by considering the layer with no half-space. The last case has been computed and experimentally confirmed previously, thus it can be checked with independent results.

There are many rays which arrive at the observation point (detector) at the same time owing to the multiple reflection and the mode conversion of the incident P (longitudinal) ray or S (shear) ray emitted from the force source in the layer. These rays are kinematically equivalent and are called "kinematic analogs". Obviously, it is not necessary to separately compute the contribution of each ray to the integration. The question of how many kinematically equivalent rays arrive at the detector at the same time for a given configuration is a problem of combinatorics. It is quite a complicated problem for a multiple layered solid half-space. So in this paper we present a new counting method to deal with this problem. In addition, some new numerical treatments are developed: 1) Automatic generation of the travel paths and the arrival times of various rays, 2) Automatic generation of the products of the reflection coefficients, and 3) An integration method for head wave rays.

The conditions for producing various head waves, surface waves, and interface waves are also examined. These conditions are determined from different singularities in the integrand.

A FORTRAN program has been developed for numerical computations of the response for any choice of materials for the layer and substrate. The computed results for the case of a plexiglass layer and glass substrate show that changes of the interface bonding condition have a great influence on the behavior of transient waves when both the source and detector are located on the top surface of the layer. For example, some of the first few head waves and regular reflected rays change their polarities and amplitudes when the interface bonding conditions change. This phenomenon can be used to infer the quality of the interface bonding of materials in ultrasonic nondestructive evaluation. In addition, results from this fundamental solution are expected to provide insight into the study and optimization of probing tools such as acoustic microscopes and in applications ranging from the study of coatings to geo-exploration.

2. Governing Equations and Boundary Conditions

Consider an elastic structure consisting of an anisotropic homogeneous layer of thickness 2h on a homogeneous half-space as shown in Fig. 1 and suppose there is a point force source of step function time dependence inside the layer. Then the fields in the layer satisfy the equations of motion

[partial][[sigma].sup.I.sub.ij]/[partial][x.sub.j] + [f.sub.i][delta]([x.sub.1])[delta]([x.sub.2])[delta]([x.sub.3] - z)H(t) = [[rho].sup.I] [[partial].sup.2][u.sup.I.sub.i]/[partial][t.sup.2],

- h < [x.sub.3] < h (1.1)

where (0, 0, z) are the position coordinates of the source. The origin of the Cartesian coordinates is at the center of the top layer; [u.sup.I.sub.i], [[sigma].sup.I.sub.ij] and [[rho].sup.I]are the displacement components, the stress components, and the density of layer I, respectively; [f.sub.i], H(t), and [delta](x) are the components of the point force source, Heaviside step function, and Dirac delta function, respectively. The summation convention is used.

The stress and displacement gradient in the layer are related by the generalized Hooke's law 

[[sigma].sup.I.sub.ij] = [C.sup.I.sub.ijkl] [partial] [u.sup.I.sub.k]/[partial] [x.sub.l] (1.2)

where [C.sup.I.sub.ijkl] are the elastic constants in the layer.

Similarly, the field equations and stress and displacement gradient relation in the lower half-space can be written by replacing "I" with "II", thus

[partial] [[sigma].sup.II.sub.ij]/[partial] [x.sub.j] = [[rho].sup.II] [[partial].sup.2][u.sup.II.sub.i]/[partial] [t.sup.2], [x.sub.3] < - h. (1.3)

where [[sigma].sup.II.sub.ij], [u.sup.II.sub.i], and [[rho].sup.II] are the stress components, the displacement componenets, and the density in the half-space.

The stress and displacement gradient in the half-space are related by

[[sigma].sup.II.sub.ij] = [C.sup.II.sub.ijkl] [partial] [u.sup.II.sub.k]/[partial] [x.sub.l], (1.4)

where [C.sup.II.sub.ijkl] are the elastic constants in the half-space.

In addition, the solution should also satisfy the following conditions:

[u.sup.I.sub.ij] = [u.sup.II.sub.k] = 0, t < 0,

[[sigma].sup.I.sub.ij] = [[sigma].sup.II.sub.ij] = 0, t < 0, (1.5)

and

[u.sup.II.sub.k] = 0, [x.sub.3] [right arrow] - [infinity]. (1.6)

In what follows we consider two cases of the interface bonding condition, whereas the top surface boundary conditions remain the same.

Case 1.

Suppose the interface is "welded", whereas the top surface condition is traction free; we have:

[[sigma].sup.I.sub.i3] = 0, [x.sub.3] = h, (1.7)

[[sigma].sup.I.sub.i3] = [[sigma].sup.II.sub.i3], [x.sub.3] = - h, (1.8)

[u.sup.I.sub.i] = [u.sup.II.sub.i], [x.sub.3] = - h. (1.9)

Case 2.

Suppose the interface is intimately "liquid" coupled, whole the top surface condition is again traction free. We have:

[[sigma].sub.i3] = 0, [x.sub.3] = h, (1.10)

[[sigma].sup.I.sub.33] = [[sigma].sup.II.sub.33], [x.sub.3] = - h, (1.11)

[[sigma].sup.I.sub.13] = [[sigma].sup.I.sub.23] = 0, [x.sub.3] = - h, (1.12)

[[sigma].sup.II.sub.13] = [[sigma].sup.II.sub.23] = 0, [x.sub.3] = -h, (1.13)

[u.sup.I.sub.3] = [u.sup.II.sub.3], [x.sub.3] = h. (1.14)

3. Solution Method

The outline of the solution procedure for the problem can be described as follows. First, introduce the Green's tensor and take its Fourier transform in time and space; then, expand the transformed Green's tensor according to the eigenvector of the Christoffel matrix, and decompose the fields into downgoing waves and upgoing waves; third, use boundary conditions to interatively get the solution in the transform domain in a form of "generalized ray" series; finally, use the Willis inversion technique to get the solution.

Defining the "Heaviside Green's tensor", [G.sub.ij], the displacements in the layer can be expressed as

[u.sup.I.sub.i] = [G.sub.ij][f.sub.j]. (1.15)

Substituting Eqs. (1.2) and (1.15) into Eq. (1.1) gives

[C.sup.I.sub.ijkl] [[partial].sup.2][G.sub.kp]/[partial][x.sub.j][partial][x.sub.1] + [[delta].sub.ip][delta]([x.sub.1])[delta]([x.sub.2])[delta]([x.sub.3] - z)H(t) = [[rho].sup.I] [[partial].sup.2][G.sub.ip]/[dt.sup.2], (1.16)

Where [[delta].sub.ip] is the Kronecker delta.

The Green's tensor in the layer may also be expressed as (1.17)

G = [G.sup.[infinity]] + [G.sup.I], (1.17)

where [G.sup.[infinity]] is the infinite body Heaviside Green's tensor and [G.sup.I] is the "image" tensor in the layer formed from the waves reflected on the boundaries [x.sub.3] = [+ or -]h. All matrices are denoted in bold capitals and vectors in bold small letters.

Define matrices [K.sup.L]([omega], [xi]) and [C.sup.L]([xi]) with components

[K.sup.L.sub.ik]([omega], x) = [[rho].sup.L][[omega].sup.2][[delta].sub.ik] - [c.sup.L.sub.ijkl][[xi].sub.j][[xi].sub.1], (1.18)

[C.sup.L.sub.ik](x) = [c.sup.L.sub.i3kl][[xi].sub.l], (1.19)

L = I, II,

where the vector [xi] = ([[xi].sub.1], [[xi].sub.2], [[xi].sub.3]).

In terms of Eqs. (1.15)-(1.19), we obtain the equations involving [G.sup.[infinity]] and [G.sup.I] which satisfy

[K.sup.I]([partial]/[partial]t, [nabla])[G.sup.[infinity]] = I[delta]([x.sub.1])[delta]([x.sub.2])[delta]([x.sub.3] - z)H(t), (1.20)

and

[K.sup.I]([partial]/[partial]t, [nabla])[G.sup.I] = 0, (1.21)

where I is the identity matrix. The Green's tensor, [G.sup.II], in the half-space satisfies

[K.sup.II]([partial]/[partial]t, [nabla])[G.sup.II] = 0. (1.22)

And the boundary conditions corresponding to cases 1 and 2 become:

Case 1.

[C.sup.I]([nabla])([G.sup.[infinity]] + [G.sup.I]) = 0, [x.sub.3] = h, (1.23)

[C.sup.I]([nabla])([G.sup.[infinity]] + [G.sup.I]) = [C.sup.II]([nabla])[G.sup.II], [x.sub.3] = -h, (1.24)

[G.sup.[infinity]] + [G.sup.I] = [G.sup.II], [x.sub.3] = -h, (1.25)

where [C.sup.I]([nabla]) and [C.sup.II]([nabla]) are the matrix operators with components

[C.sup.L.sub.ik]([nabla]) = [c.sup.L.sub.i3kl] [partial]/[partial][x.sub.l], L = I, II. (1.26)

Case 2.

[C.sup.I]([nabla])([G.sup.[infinity]] + [G.sup.I]) = 0, [x.sub.3] = h, (1.27)

[I.sub.1][C.sup.I]([nabla])([G.sup.[infinity]] + [G.sup.I]) = [I.sub.1][C.sup.II]([nabla])[G.sup.I] = 0, [x.sub.3] = -h, (1.28)

[I.sub.2][C.sup.I]([nabla])([G.sup.[infinity]] + [G.sup.I]) = [I.sub.2][C.sup.II]([nabla])[G.sup.I] = 0, [x.sub.3] = -h, (1.29)

[I.sub.3][C.sup.I]([nabla])([G.sup.[infinity]] + [G.sup.I]) = [I.sub.3][C.sup.II]([nabla])[G.sup.I], [x.sub.3] = -h, (1.30)

[I.sub.3]([G.sup.[infinity]] + [G.sup.I]) = [I.sub.3][G.sup.II], [x.sub.3] = -h, (1.31)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

4. Ray Expansion

We follow the method developed in Refs. (2,3,4), but provide only an outline here. Defining the Fourier transform of [G.sup.[infinity]] by

[G.sub.[infinity]]([omega], [xi]) = [[integral].sup.[infinity].sub.0] dt [integral] [[integral].sup.[infinity].sub.-[infinity]] [integral] [dx.sub.1][dx.sub.2][dx.sub.3][G.sup.[infinity]](t, x)exp[i([xi]x - [omega]t)], (2.1)

from Eq. (1.20) we then get

[K.sup.I]([omega], [xi]) [G.sup.[infinity]]([omega], [xi]) = -i/[omega] I exp(i[[xi].sub.3]z), (2.2)

where [omega] is taken to have negative imaginary parts while [xi] = ([[xi].sub.1], [[xi].sub.2], [[xi].sub.3]) has real components. It is easily shown that [G.sub.[infinity]] is an analytic function of [omega] in the lower half of the [omega]-plane, and its inverse transform as a function of time t is equal to zero when t < 0.

In order to express the inverse of the matrix [K.sup.I]([omega], [xi]) in terms of its eigenvectors, we consider the Christoffel equation

[K.sup.I]([omega], [xi])[u.sub.r] = [[LAMBDA].sub.r]([omega], [xi])[u.sub.r], (2.3)

where [[LAMBDA].sub.r]([omega], [xi]) and [u.sub.r] are the eigenvalues and eigenvectors.

For given real [omega], [[xi].sub.1], and [[xi].sub.2], the equation

det [K.sup.I]([omega], [xi]) = 0 (2.4)

has six roots [[xi].sub.3] = [[xi].sup.N.sub.3] ([omega], [[xi].sub.[alpha]]), N = [+ or -] 1, 2, 3, [alpha] = 1, 2, which may either be real or occur in complex conjugate pairs. By means of the concept of Riemann surfaces, the six roots may be considered to define a single-valued algebraic function [[xi].sub.3] ([omega], [[xi].sub.[alpha]]) by Eq. (2.4), if [omega] is allowed to range over the six sheets of its Riemann surface (2). It can be shown by analytic continuation that when Im([omega]) < 0, the algebraic function [[xi].sub.3] has positive imaginary parts on the three sheets of N = -1, -2, -3 and has negative imaginary parts on the other three sheets of N = +1, +2, +3.

Now let us consider the eigenvalues. Since [[LAMBDA].sub.N] = [[rho].sup.I] ([[omega].sup.2] - [[omega].sup.2.sub.N]), where [[omega].sub.N] is inverse to [[xi].sub.3.sup.N] ([omega], [[xi].sub.[alpha]]), let [[xi].sup.N.sub.3] ([[omega].sub.N], [[xi].sub.[alpha]]) or det [K.sup.I]([omega.sup.N], [xi]) = 0, and we obtain the six roots [+ or -] [[omega].sub.N], N = 1, 2, 3 and therefore the three eigenvalues [[LAMBDA].sub.N].

Normalizing the eigenvectors so that [[u.sub.m][u.sub.n].sup.T] = [[delta].sub.mn], we have

I = [SIGMA] [[u.sub.r][u.sub.r].sup.T] = [UU.sup.T] (2.5)

Where the matrix U consists of the three column vectors [u.sub.r], while [U.sup.T] is the transpose of U. As a result we have

[G.sup.[infinity]]([omega], [xi]) = -i/[omega] [summation over (3/r=1)] [[LAMBDA].sup.-1.sub.r] [[u.sub.r][u.sub.r].sup.T] exp(i[[xi].sub.3]z). (2.6)

Using symmetry with respect to the [x.sub.3] axis, and taking the Fourier transform of [G.sup.I] from (t, [x.sub.1], [x.sub.2], [x.sub.3]) to ([omega], [[xi].sub.1], [[xi].sub.2], [[xi].sub.3]) we obtain for the fields in the layer

G([omega], [xi]) = [[integral].sup.[infinity].sub.0] dt[integral] [[integral].sup.[infinity].sub.-[infinity]] d[x.sub.1]d[x.sub.2][G.sup.I](t, x)exp[i([[xi].sub.[alpha]] [[x.sub.[alpha]] - [[omega].sub.t])], (2.7)

where

[[xi].sub.[alpha]][[x.sub.[alpha]] = [[xi].sub.1][x.sub.1] + [[xi].sub.2][x.sub.2]. (2.8)

In what follows we consider the cases of downgoing waves and upgoing waves in the layer respectively.

1. Downgoing waves.

Consider the case [x.sub.3] < z. Here z is the position of the point force source. Taking the inverse Fourier transform of Eq. (2.6) about [[xi].sub.3] gives

[G.sup.[infinity]] ([omega], [[xi].sub.1], [[xi].sub.2], [x.sub.3]) = -i/2[pi][omega] [[integral].sup.[infinity].sub.-[infinity]] d[[xi].sub.3] [summation over (3/r=1)] [[LAMBDA].sup.-1.sub.r] [[u.sub.r][u.sub.r].sup.T] exp[-i[[xi].sub.3] ([x.sub.3] - z)]. (2.9)

When [x.sub.3] < z, [x.sub.3] - z is negative. Then the integral can be evaluated by closing the contour in the upper half of the [[xi].sub.3]-plane and by using Cauchy's residue theorem. This gives

[G.sup.[infinity]] ([omega], [[xi].sub.1], [[xi].sub.2], [x.sub.3]) = 1/[omega] [summation over (3/M=1)] [[u.sup.-.sub.M][u.sup.-.sub.M].sup.T] exp[-i[[xi].sup.M-.sub.3] ([x.sub.3] - z)]/[([partial] [LAMBDA] ([omega], [xi])/[partial][[xi].sub.3]).sub.[[xi].sub.3]=[[xi].sup.M-.sub.3]] (2.10)

where [[xi].sup.M-.sub.3] = [[xi].sup.M-.sub.3] ([omega], [[xi].sub.1], [[xi].sub.2]) above, which are located in the upper half of the [[xi].sub.3]-plane when Im([omega]) < 0. The subscript or superscript "-" denotes downgoing waves.

2. Upgoing waves.

In the case [x.sub.3] > z, the contour of integral Eq. (2.9) can be closed in the lower half-plane. Similarly we obtain

[G.sup.[infinity].sub.+] ([omega], [[xi].sub.1], [[xi].sub.2], [x.sub.3]) = 1/[omega] [summation over (3/N=1)] [[u.sup.+.sub.N][u.sup.+.sub.N].sup.T] exp[-i[[xi].sup.N+.sub.3] ([x.sub.3] - z)]/[([partial][LAMBDA]([omega],[xi])/[partial][[xi].sub.3]).sub.[[xi ].sub.3]=[[xi].sup.N+.sub.3] (2.11)

where [[xi].sup.N+.sub.3] = [[xi].sup.N+.sub.3] ([omega], [[xi].sub.1], [[xi].sub.2]) are the three roots of Eq. (2.4) which are located in the lower half of the [[xi].sub.3]-plane, when Im ([omega]) < 0. The subscript or superscript "+" denotes upgoing waves.

The general solution of [G.sup.[infinity]] ([omega], [[xi].sub.1], [[xi].sub.2], [x.sub.3]) should include three downgoing plane waves and three upgoing plane waves. The reflected waves in the layer should also include the downgoing waves and the upgoing waves, but there exist only the downgoing waves in the lower half-space. Then we respectively obtain

[G.sup.[infinity]] ([omega], [[xi].sub.1], [[xi].sub.2], [x.sub.3]) = [G.sup.[infinity].sub.-] ([omega], [[xi].sub.1], [[xi].sub.2], [x.sub.3]) + [G.sup.[infinity].sub.+]([omega], [[xi].sub.1], [[xi].sub.2], [x.sub.3]), (2.12)

[G.sup.I] = [summation over (3/M=1)] [u.sup.-.sub.M][v.sup.-T.sub.M] exp[-i[[xi].sup.M-.sub.3] ([x.sub.3] - h)] + [summation over (3/N=1)] [u.sup.+.sub.N] [v.sup.+.sub.N] exp[-i[[xi].sup.N+.sub.3] ([x.sub.3] - h)], -h < [x.sub.3] < h, (2.13)

[G.sup.II] = [summation over (3/P=1)] [q.sup.-.sub.P] [w.sup.-T.sub.P] exp[-i[[xi].sup.P-.sub.3] ([x.sub.3] - h)], [x.sub.3] < -h, (2.14)

where [q.sup.-.sub.P] represents the eigenvectors of an infinite body in the lower half-space while [[xi].sup.P-.sub.3] represent the roots of Eq. (2.4) when the superscript "I" is replaced by "II", for the same reasoning as mentioned above, which roots locate on the upper half-plane of [[xi].sub.3]-plane when Im ([omega]) < 0, and [v.sup.-.sub.M], [v.sup.+.sub.N], and [W.sup.-.sub.P] are the coefficients to be determined.

Further matrix notation is introduced as follows,

[U.sup.-] = [[u.sup.-.sub.1], [u.sup.-.sub.2], [u.sup.-.sub.3]], (2.15)

[V.sup.-] = [[v.sup.-.sub.1], [v.sup.-.sub.2], [v.sup.-.sub.3]], (2.16)

with corresponding definitions when superscript "-" is replaced by "+" and

[Q.sup.-] = [[q.sup.-.sub.1], [q.sup.-.sub.2], [q.sup.-.sub.3]], (2.17)

[W.sup.-] = [[w.sup.-.sub.1], [w.sup.-.sub.2], [w.sup.-.sub.3]], (2.18)

Then the previous equations can be written in more compact form as

[G.sup.[infinity]] ([omega], [[xi].sub.1], [[xi].sub.2], [[xi].sub.3]) = H(z-[x.sub.3])[U.sup.-][P.sup.-](z-[x.sub.3]-h)[D.sup.-][U.sup.-T]

+ H([x.sub.3]-z)[U.sup.+][P.sup.+](z-[x.sub.3]+h)[D.sup.+][U.sup.+T], (2.19)

where [P.sup.-], [P.sup.+], [D.sup.-], and [D.sup.+] are diagonal matrices with components

[([P.sup.-]).sup.l.sub.k]([x.sub.3]) = [[delta].sup.l.sub.k]exp[i[[xi].sup.l-.sub.3]([x.sub.3] + h)], (2.20)

where [[delta].sup.l.sub.k] is the Kronicker delta. Also,

[([P.sup.+]).sup.l.sub.k]([x.sub.3] = [[delta].sup.l.sub.k]exp[i[[xi].sup.l+.sub.3]([x.sub.3] - h)], (2.21)

[([D.sup.-]).sup.l.sub.k] = [[delta].sup.l.sub.k][[[omega][([partial][LAMBDA]([omega], [xi])/[partial][[xi].sub.3]).sub.[[xi].sub.3]=[[xi].sup.l-.sub.3]]).s up.-1], (2.22)

[([D.sup.+]).sup.l.sub.k] = [[delta].sup.l.sub.k][[[omega][([partial][LAMBDA]([omega], [xi])/[partial][[xi].sub.3]).sub.[[xi].sub.3]=[[xi].sup.l+.sub.3]]).s up.-1], (2.23)

and

[G.sup.I] = [U.sup.-][P.sup.-]([-x.sub.3])[V.sup.-T] + [U.sup.+][P.sup.+]([-x.sub.3])[V.sup.+T], (2.24)

[G.sup.II] = [Q.sup.-][S.sup.-]([-x.sub.3])[W.sup.-T], (2.25)

where [S.sup.-] is the diagonal matrix with components

[([S.sup.-]).sup.l.sub.k]([x.sub.3]) = [[delta].sup.l.sub.k]exp [i[[xi].sup.l-.sub.3]([x.sub.3] - h)]. (2.26)

Case 1.

If we are only interested in the fields in the layer, substituting Eqs. (2.19), (2.24), and (2.25) into the Fourier transformed equations corresponding to the boundary conditions Eqs. (1.23)-(1.25) and eliminating the coefficient matrix [W.sup.-], we get

[B.sup.-.sub.1][V.sup.-T] - [B.sup.+.sub.1][P.sup.+][V.sup.+T] = - [B.sup.+.sub.1][P.sup.+](z)[D.sup.+][U.sup.+T] (2.27)

and

[B.sup.+.sub.2][V.sup.+T] - [B.sup.-.sub.2][P.sup.-][V.sup.-T] = - [B.sup.-.sub.2][P.sup.-](z)[D.sup.-][U.sup.-T] (2.28)

where

[B.sup.+.sub.1] = [C.sup.I]([xi])[U.sup.+], (2.29)

[B.sup.+.sub.2] = [[C.sup.I]([xi]) - [C.sup.II]([xi])[Q.sup.-][([Q.sup.-]).sup.-1]] [U.sup.+], (2.30)

[B.sup.-.sub.1] = [C.sup.I]([xi])[U.sup.-], (2.31)

[B.sup.-.sub.2] = [[C.sup.I]([xi]) - [C.sup.II]([xi])[Q.sup.-][([Q.sup.-]).sup.-1]] [U.sup.-], (2.32)

[P.sup.+] = -[P.sup.+](-h), (2.33)

[P.sup.-] = - [P.sup.-](h). (2.34)

As h approaches infinity the exponential functions [P.sup.+] and [P.sup.-] on the left-hand sides of Eqs. (2.27) and (2.28) approach zero because the exponential factors have negative and positive imaginary parts, respectively. Therefore, Eqs. (2.27) and (2.28) uncouple into two independent equations and the problem becomes two separate problems for half-spaces [x.sub.3] < h and [x.sub.3] > -h. In addition, the larger Im ([omega]) becomes the smaller exponential function [P.sup.+] and [P.sup.-], and so the equations can be solved iteratively to give a uniformly convergent series of "generalized rays" for given [[xi].sub.1], [[xi].sub.2], and the real part of [omega]. Finding [V.sup.-T] and [V.sup.+T] from Eqs. (2.27) and (2.28) by means of an iteration technique similar to Ref. [5] and substituting the results obtained into the previous relevant equations we get the general field in the layer

G = [G.sup.-] + [G.sup.+], (2.35)

where

[G.sup.-] = H(z - [x.sub.3])[U.sup.-][P.sup.-](z - [x.sub.3] - h)[S.sup.-]

- [U.sup.-][P.sup.-]([-x.sub.3])[R.sup.+][P.sup.+](z)[S.sup.+] - [U.sup.-][P.sup.-]([-x.sub.3])[R.sup.+][P.sup.+][R.sup.-][P.sup.-](z) [S.sup.-]

- [U.sup.-][P.sup.-]([-x.sub.3])[R.sup.+][P.sup.+][R.sub.-][P.sub.-][R. sup.+][P.sub.+](z)[S.sub.+] - ... (2.36)

and

[G.sup.+] = H([x.sub.3] - z)[U.sup.+][P.sup.+](z - [x.sub.3] + h)[S.sup.+]

- [U.sup.+][P.sup.+](-[x.sub.3])[R.sup.-][P.sup.-](z)[S.sup.-] -[U.sup.+][P.sup.+](-[x.sub.3])[R.sup.-][P.sup.-][R.sup.+][P.sup.+](z )[S.sup.+]

- [U.sup.+][P.sup.+](-[x.sub.3])[R.sup.-][P.sup.-][R.sup.+][P.sup.+][R. sup.-][P.sup.-](z)[S.sup.-]-... (2.37)

which are the downgoing waves and the upgoing waves in the layer, respectively. Their physical meanings are obvious. If we notice that [R.sup.+] and [R.sup.-] are respectively the matrices of the reflection coefficients at the top surface and the lower surface of the layer, and the diagonal matrices [P.sup.+] and [P.sup.-] represent phase delays, then the first term of Eq. (2.36) represents the downgoing direct rays from the source to the observation point (detector), the second term describes the rays which reflect once at the top surface of the layer and so on. Similarly, we can explain each term of Eq. (2.37). Using the concept of an infinite linear array of image sources, then, with the exception of the first (direct wave) term, all the other terms of Eqs. (2.36) or (2.37) can be considered to be produced from the corresponding image sources located above or below the layer. The definitions of the notation of Eqs. (2.36) and (2.37) are as follows:

[L.sup.-.sub.1] = [([B.sup.-.sub.1]).sup.-1], [L.sup.+.sub.2] = [([B.sup.+.sub.2]).sup.-1], (2.38)

[R.sup.-] = [L.sup.+.sub.2][B.sup.-.sub.2], [R.sup.+] = [L.sup.-.sub.1][B.sup.+.sub.1], (2.39)

[S.sup.-] = [D.sup.-][U.sup.-T], [S.sup.+] = [D.sup.+][U.sup.+T]. (2.40)

A typical term of Eq. (2.36) except for the first term may be written as

[U.sup.-][P.sup.-](-[x.sub.3])[R.sup.+][P.sup.+][R.sup.-][P.sup.-]... [R.sup.[(-1).sup.k]][P.sup.[(-1).sup.k]](z)[S.sup.[(-1).sup.k]], (2.41)

where k counts the number of P's, R's or S's in the term. Interchanging superscripts "+" and "-" in Eq. (2.41) gives the typical term of Eq. (2.37).

Suppose an element of the matrix of a typical term of Eq. (2.41) is written in the form

F([omega], [xi])exp[-I([[xi].sub.[alpha]][x.sub.[alpha]] - [omega]t + [[lambda].sub.1][[xi].sup.M.sub.3] + [[lambda].sub.2][[xi].sup.N.sub.3])].

Then the inverse Fourier transform of Eq. (2.36) gives

[([G.sup.-]).sup.ml] 1/8[[pi].sup.3] [[integral].sup.[infinity]-0i.sub.-[infinity]-0i] d[omega] [integral][integral] d[[xi].sub.1]d[[xi].sub.2] [summation over (F([omega], [xi]))]

exp[-I([[xi].sub.[alpha]][x.sub.[alpha]] - [omega]t + [[lambda].sub.1][[xi].sup.N.sub.3])], (2.42)

where [([G.sup.-]).sup.ml] is the element of the matrix [G.sup.-] and F([omega], [xi]) is a homogeneous function of degree -2, because the elements of the matrices [U.sup.[+ or -]] and [R.sup.[+ or -]] are homogeneous functions of degree zero, while the elements of the matrices [S.sup.[+ or -]] are homogeneous functions of degree -2. In addition, since the series of "generalized rays" is uniformly convergent, we can invert the integration term by term using the Willis inversion method (Appendix A) to obtain

[([G.sup.-]).sup.ml] = -1/4[[pi].sup.2] [summation over ()] ** ds F([OMEGA], [[eta].sub.[alpha]])/-t + [[lambda].sub.1][[xi].sup.M.sub.3,[OMEGA]] + [[lambda].sub.2][[xi].sup.M.sub.3,[OMEGA]], (2.43)

M=1/[empty set]=0

Where [([G.sup.-]).sup.ml] is the element of [G.sup.-], while the subscript [OMEGA] denotes [partial]/[partial][OMEGA]

[OMEGA]=[omega]/\[xi]\, [[eta].sub.[alpha]]=[[xi].sub.[alpha]]/\[xi]\. (2.44)

[OMEGA] in the integrand is now taken as the root in the lower half-plane of the equation

-[OMEGA]t + [[eta].sub.[alpha]][x.sub.[alpha]] + [[lambda].sub.1][[xi].sup.M.sub.3]([OMEGA], [[eta].sub.[alpha]]) + [[lambda].sub.2][[xi].sup.N.sub.3]([OMEGA], [[eta].sub.[alpha]]) = [epsilon]i, (2.45)

where [epsilon] is an arbitrary infinitesimal number.

For the elements of [G.sup.+], we have similar results.

Case 2.

In this case, substituting Eqs. (2.19), (2.24), and (2.25) into the Fourier transformed equations corresponding to the boundary conditions Eqs. (1.27)-(1.29), we get the same Eqs. (2.27) and (2.28) when [B.sup.+]2 and [B.sup.-]2 are replaced by [B.sup.*+]2 and [B.sup.*-]2, respectively, and

[B.sup.*[+ or -]] = {[C.sup.I]([xi]) - ([I.sub.3][C.sup.II]([xi])[Q.sup.-])

[[([I.sub.1]+[I.sub.2])[C.sup.II]([xi])[Q.sup.-] + [I.sub.3][Q.sup.-]].sup.-1]}[U.sup.[+ or -]]. (2.46)

Defining

[L.sup.*+.sub.2] = [([B.sup.*+.sub.2]).sup.-1], (2.47)

[R.sup.*-] = [L.sup.*+.sub.2][B.sup.*-.sub.2], (2.48)

and using the matrices with "*" instead of the corresponding matrices without "*" in Eqs. (2.35)-(2.36), we immediately obtain the solution of case 2. The inversion of the transform is the same.

Up to now we have not considered the properties of the materials, so all the previous results can be applied to the general case of arbitrary anisotropic materials.

5. Special Case of Isotropic Materials

When the layer and the lower half-space consist of two different isotropic materials, then the components of [K.sub.L]([omega], [xi]) become

[K.sup.L.sub.ik]([omega], [xi])

= ([[rho].sub.L][[omega].sup.2] - [[mu].sub.L][[xi].sub.j][[xi].sub.j])[[delta].sub.ik] - ([[lambda].sub.L] + [[mu].sub.L])[[xi].sub.i][[xi].sub.k], L = I, II. (3.1)

where [[lambda].sub.L] and [[mu].sub.L] are Lame's elastic constants of the materials in the layer (L = I) and the lower half-space (L = II).

The roots of the equations

det [K.sup.L]([omega], [xi]) = 0, L = I, II (3.2)

are

[[xi].sup.[+ or -]aL.sub.3] = [+ or -][p.sub.L] (once), (3.3)

[[xi].sup.[+ or -]cL.sub.3] = [+ or -][q.sub.L] (twice), (3.4)

where

[p.sub.L] = i[([[xi].sup.2] - [[omega].sup.2]/[a.sup.2.sub.L]).sup.1/2], L = I, II, (3.5)

[q.sub.L] = i[([[xi].sup.2] - [[omega].sup.2]/[c.sup.2.sub.L]).sup.1/2], L = I, II, (3.6)

[a.sup.2.sub.L] = ([[lambda].sub.L] + [[mu].sub.L])/[[rho].sub.L], L = I, II, (3.7)

[c.sup.2.sub.L] = [[mu].sub.L]/[[rho].sub.L], L = I, II, (3.8)

[[xi].sup.2] = [[xi].sup.2.sub.1] + [[xi].sup.2.sup.2]. (3.9)

Making a comparsion with the previous sections gives

[[xi].sup.M-.sub.3] = [p.sub.I], [q.sub.I], [q.sub.I] (3.10)

[[xi].sup.N+.sub.3] = -[p.sub.I], -[q.sub.I], -[q.sub.I].

Three roots, [p.sub.I], [q.sub.I], and [q.sub.I], of the six roots in Eq. (3.10) are in the upper half of the [xi]-plane, whereas the other three roots, -[p.sub.I], -[q.sub.I], and -[q.sub.I], are in the lower half of the [[xi].sup.3]-plane. For the details of the distribution of these roots, refer to Ref. (2)

The normalized eigenvectors corresponding to the roots of Eq. (3.10) are

[u.sup.[+ or -].sub.1] = [[[xi].sub.1] [[xi].sub.2] [- or +][p.sub.I]][([[xi].sup.2] + [([p.sub.I).sup.2]).sup.-1/2], (3.11)

[u.sup.[+ or -].sub.2] = [[+ or -][[xi].sub.1][q.sub.I] [+ or -][[xi].sub.2][q.sub.I] -[[xi].sup.2]] [([[xi].sup.2] + [([q.sub.I]).sup.2]).sup.-1/2] [[xi].sup.-1], (3.12)

and

[u.sup.[+ or -].sub.3] = [-[[xi].sub.2][[xi].sup.-1] [[xi].sub.1][[xi].sup.-1] 0]. (3.13)

For the lower half-space there exist only the downgoing waves; the following three roots should be chosen as

[[xi].sup.p-.sub.3] = [p.sub.II], [q.sub.II], [q.sub.II]. (3.14)

The corresponding eigenvectors are

[q.sup.-.sub.1] = [[[xi].sub.1] [[xi].sub.2] [p.sub.II]] [[[xi].sub.2] + [([p.sub.II].sup.2]].sup.-1/2], (3.15)

[q.sup.-.sub.2] = [[[xi].sub.1][q.sub.II] [[xi].sub.2][q.sub.II] -[[xi].sup.2]] [[[xi].sub.2] + [([q.sub.II]).sup.2]].sup.-1/2]/[xi], (3.16)

and

[q.sup.-.sub.3] = [-[[xi].sub.2][[xi].sup.-1] [[xi].sub.1][[xi].sup.-1] 0]. (3.17)

The matrix operators [C.sup.I]([xi]) and [C.sup.II]([xi]) should always operate on the respective regions of eigenvectors first to obtain the proper stress components. We have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.18)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.19)

Substituting the relevant matrices into Eqs. (2.39) and (2.48) we obtain the reflection coefficient matrices [R.sup.+], [R.sup.-], and [R.sup.*-] respectively.

[R.sup.+] = [[R.sup.+.sub.ij]], (3,20)

[R.sup.+.sub.11] = [R.sup.+.sub.22] = [[([q.sup.2.sub.1] - [[xi].sup.2]).sup.2] - 4[p.sub.1][q.sub.I][[xi].sup.2]]/[[([q.sup.2.sub.1] - [[xi].sup.2]).sup.2] + 4[p.sub.1][q.sub.I][[xi].sup.2]],

[R.sup.+.sub.12] = [4[q.sub.1][xi]([q.sup.2.sub.1] - [[xi].sup.2])[([p.sup.2.sub.1] + [[xi].sup.2]).sup.1/2]]/{[([q.sup.2.sub.1] + [[xi].sup.2]).sup.1/2][[([q.sup.2.sub.1] - [[xi].sup.2]).sup.2] + 4[p.sub.1][q.sub.1][[xi].sup.2]},

[R.sup.+.sub.21] = - [4[p.sub.1][xi]([q.sub.2.sub.1] - [[xi].sup.2])[([q.sup.2.sub.1] + [[xi].sup.2]).sup.1/2]]/{[([p.sup.2.sub.I] + [[xi].sup.2]).sup.1/2][[([q.sup.2.sub.1] - [[xi].sup.2]).sup.2] + 4[p.sub.1][q.sub.I][[xi].sup.2]]},

[R.sup.+.sub.33]= -1,

[R.sup.+.sub.13] = [R.sup.+.sub.23] = [R.sup.+.sub.32] = 0,

[R.sup.-] = [[R.sup.+.sub.ij]], (3.21)

[R.sup.-.sub.11] = [-[[DELTA].sup.+.sub.1][[DELTA].sup.-.sub.2] + [[DELTA].sup.-.sub.3][[DELTA].sup.+.sub.4]]/[[DELTA].sub.II],

[R.sup.-.sub.12] = [-[[DELTA].sup.+.sub.1][[DELTA].sup.-.sub.4] + [[DELTA].sup.+.sub.4][[DELTA].sup.-.sub.1]]/ [[DELTA].sub.II],

[R.sup.-.sub.21] = [[[DELTA].sup.+.sub.3][[DELTA].sup.-.sub.2] + [[DELTA].sup.+.sub.2][[DELTA].sup.-.sub.3]]/ [[DELTA].sub.II],

[R.sup.-.sub.22] = [[DELTA].sup.+.sub.3][[DELTA].sup.-.sub.4] - [[DELTA].sup.+.sub.2][[DELTA].sup.-.sub.1]]/ [[DELTA].sub.II],

[R.sup.-.sub.33] = ([[mu].sub.1][q.sub.I] - [[mu].sub.II][q.sub.II])/([[mu].sub.1][q.sub.I] + [[mu].sub.II][q.sub.II]),

[R.sup.-.sub.13] = [R.sup.-.sub.23] = [R.sup.-.sub.31] = [R.sup.-.sub.32] = 0,

where

[[DELTA].sup.[+ or -].sub.1] = {[+ or 1]2[[mu].sub.1][q.sub.1]] + [[mu].sub.II][[q.sub.2]([q.sup.2.sub.2] + [[xi].sup.2]) [+ or -][q.sub.1](2[p.sub.2][q.sub.2]] + [[xi].sup.2] - [q.sup.2.sub.2])]/([p.sub.2][q.sub.2] + [[xi].sup.2])}/[([q.sup.2.sub.1] + [[xi].sup.2]).sup.1/2],

[[DELTA].sup.[+ or -].sub.2] = {[+ or -]2[[mu].sub.1][p.sub.1]] + [[mu].sub.II][ -[p.sub.2]([q.sup.2.sub.2] + [[xi].sup.2])

[+ or 1][p.sub.1](2[p.sub.2][q.sub.2] + [[xi].sup.2] - [q.sup.2.sub.2])]/([p.sub.2][q.sub.2] + [[xi].sup.2])}/[([p.sup.2.sub.1] + [[xi].sup.2]).sup.1/2],

[[DELTA].sup.[+ or -].sub.3] = {[[mu].sub.1]([q.sup.2.sub.1] - [[xi].sup.2]) + [[mu].sub.II][(2[p.sub.2][q.sub.2] + [[xi].sup.2] - [q.sup.2.sub.2])

[+ or -][p.sub.1][q.sub.2]([q.sup.2.sub.2] + 1)]/([p.sub.2][q.sub.2] + [[xi].sup.2])}/[([p.sup.2.sub.1] + [[xi].sup.2]).sup.1/2],

[[DELTA].sup.[+ or -].sub.4] = {[[mu].sub.1]([q.sup.2.sub.1] - [[xi].sup.2]) + [[mu].sub.II][([+ or -][p.sub.2][q.sub.I] - 1)([q.sup.2.sub.2] + 1) +(2[p.sub.2][q.sub.2] + 2[xi] - [q.sup.2.sub.2])]/([p.sub.2][q.sub.2] + [[xi].sup.2])}/[([q.sup.2.sub.2] + [[xi].sup.2]).sup.1/2],

[[DELTA].sub.II] = [[DELTA].sup.+.sub.3][[DELTA].sup.+.sub.4] - [[DELTA].sup.+.sub.1][[DELTA].sup.+.sub.2],

and

[R.sup.*-] = [[R.sup.*-.sub.ij]], (3.22)

[R.sup.*-.sub.11] = [- [[DELTA].sup.*+.sub.1] [[DELTA].sup.*-.sub.2] + [[DELTA].sup.*-.sub.3] [[DELTA].sup.*+.sub.4]]/[[DELTA].sup.*.sub.II],

[R.sup.*-.sub.12] = (- [[DELTA].sup.*+.sub.1] [[DELTA].sup.*-.sub.4] + [[DELTA].sup.*+.sub.4] [[DELTA].sup.*-.sub.1])/[[DELTA].sup.*.sub.II],

[R.sup.*-.sub.21] = ([[DELTA].sup.*+.sub.3] [[DELTA].sup.*-.sub.2] + [[DELTA].sup.*+.sub.2] [[DELTA].sup.*-.sub.3])/[[DELTA].sup.*.sub.II],

[R.sup.*-.sub.22] = [[[DELTA].sup.*+.sub.3] [[DELTA].sup.*-.sub.4] - [[DELTA].sup.*+.sub.2] [[DELTA].sup.*-.sub.1]]/[[DELTA].sup.*.sub.II]]

[R.sup.*-.sub.33] = ([[mu].sub.I][q.sub.I] - [[mu].sub.II][q.sub.II]),

[R.sup.*-.sub.13] = [R.sup.*-.sub.23] = [R.sup.*-.sub.31] = [R.sup.*-.sub.32] = 0,

where

[[DELTA].sup.*[+ or -].sub.1] = {[+ or -]2[[mu].sub.1][q.sub.1] + [[mu].sub.2][[([q.sup.2.sub.2] - 1).sup.2] + 4[p.sub.2][q.sub.2]]/[p.sub.2]([q.sup.2.sub.2] + 1)}/[([q.sup.2.sub.1] + 1).sup.1/2],

[[DELTA].sup.*[+ or -].sub.2] = [+ or -]2[p.sub.1][[mu].sub.1]/[([p.sup.2.sub.1] + 1).sup.1/2],

[[DELTA].sup.*[+ or -].sub.3] = {[[mu].sub.1]([q.sup.2.sub.1] - 1) [+ or -] [[mu].sub.2][p.sub.1][[([q.sup.2.sub.2] - 1).sup.2] + 4[p.sub.2][q.sub.2]]/[p.sub.2]([q.sup.2.sub.2] + 1)}/[([[p.sup.2.sub.1] + 1).sup.1/2],

[[DELTA}.sup.-[+ or -].sub.4] = [[mu].sub.2] 4[q.sub.2]([p.sub.2] - [q.sub.2])/[p.sub.2]([q.sup.2.sub.2] + 1) + [[rho].sub.2][[omega].sup.2]/[p.sub.2].

In what follows let us first consider the interface condition for case 1, i.e., the welded interface. In the case of isotropic materials the inverse transform Eq. (2.43) takes the following form:

[([G.sup.-]).sup.ml] = [summation over ([infinity]/k=1])][(-1).sup.k][1 + [[delta].sup.l.sub.k](H(z - [x.sub.3]) - 1)] 1/[4[pi].sup.2]

[summation over (g=1,2,3 j=1,2,3)] [[**].sub.\[eta]\=1 [phi]=0] [([U.sup.-]).sup.m.sub.j] [R.sup.-] (j, g, k) [[([U.sup.[(-1).sup.k]]).sup.T]].sup.g.sub.l] [([D.sup.[(-1).sup.k]]).sup.g.sub.g]/[partial] [[phi].sup.-] ([OMEGA], [eta])/[partial][OMEGA] ds, (3.29)

where j and g take 1, 2, and 3 and respectively represented the type of the final trip and the initial trip of the ray path to be a P ray, SV ray or SH ray in the layer, and

[[phi].sup.-]([OMEGA], [eta]) = -[OMEGA]t + [[eta].sub.[alpha]][x.sub.[alpha]] - 2h[k.sup.-.sub.1][p.sub.I], - 2[hk.sup.-.sub.2][q.sub.I], (3.30)

where [k.sup.-.sub.1] and [k.sup.-.sub.2] are the numbers of times the layer is traversed by a P ray and S ray, respectively. They may take fractional values when the source or the detector is not on the top surface or the interface. Also,

[([D.sup.-]).sup.l.sub.k]] = -[(2[[rho].sub.I][OMEGA][b.sub.k]).sup.-1][[delta].sup.l.sub.k], (3.31)

where

[b.sub.1] = [a.sup.2.sub.I][p.sub.I], [b.sub.2] = [b.sub.3] = [c.sup.2.sub.I][q.sub.I]. (3.32)

R (j,g,k) represents the sum of the products of all the reflection coefficients produced by the rays which travel from the source to the observation point (detector) and touch the top surface and the bottom surface of the layer. These rays have the same final trip j and the same initial trip g as well as the same total number of layer traverses k of the layer as shown in Fig. 2. I.e., they have the same arrival time. The concrete expression of [R.sup.-](j,g,k) will be discussed in the next section.

Similarly, we have

[([G.sup.+]).sup.ml] = [summation over ([infinity]/k=1)][(-1).sup.k] {1 + [[delta].sup.l.sub.k] [H ([x.sub.3] - z) - 1]} 1/4[[pi].sup.2]

[[SIGMA].sub.g=1,2,3 j=1,2,3] [[**].sub.\[eta]\=1 [phi]=0] [([U.sup.+]).sup.m.sub.j][R.sup.+](j,g,k)[[[([U.sup.[(-1).sup.k+1]]). sup.T].sup.g.sub.l][([D.sup.[(-1).sup.k+1]])].sup.g.sub.g]/[partial][ [phi].sup.+]([OMEGA], [eta])/[partial][OMEGA] ds, (3.33)

where

[[phi].sup.+]([OMEGA], [eta]) = - [OMEGA]t + [[eta].sub.[alpha]][x.sub.[alpha]] - 2[hk.sup.+.sub.1][p.sub.I] - 2[hk.sup.+.sub.2][q.sub.I], (3.34)

[([D.sup.+]).sup.l.sub.k] = -[(2[[rho].sub.I][OMEGA][b.sub.k]).sup.-1] - [[delta].sup.l.sub.k]. (3.35)

The meanings of the various quantities in Eqs. (3.33)-(3.35) are the same as before or are similar to that of the quantities in Eqs. (3.29)-(3.32).

Using

[phi]([OMEGA], [eta]) = 0, (3.36)

then Eqs. (3.30) and (3.34) can be written as

[partial][[phi].sup.[+ or -]]/[partial] [OMEGA][\.sub.[phi]=0] = 1/[OMEGA](-[[eta].sub.[alpha]][x.sub.[alpha]] - 2[hk.sup.[+ or -].sub.1]/[p.sub.I] - 2[hk.sup.[+ or -].sub.2]/[q.sub.I]). (3.37)

5.1 Detector on the Top Surface

Usually, the detector (observation point) is put on the top surface of the layer; then Eqs. (3.29) and (3.33) can be further combined. If the detector was considered to be located an infinitesimal amount below the upper surface, we should have to take into account at almost the same instant the wave going upward before reflection and the reflected wave going downward which come from the same ray and have almost the same phase function. Substituting into Eqs. (3.29) and (3.33) with k = n + 1 and k = n + 2, respectively, and defining

[GAMMA] = [U.sup.+] - [U.sub.-][R.sup.+], (3.38)

and substituting the previous relevant expressions into Eq. (3.38) we obtain

[GAMMA] = [[[gamma].sub.ij]]. (3.39)

Here

[[gamma].sub.11] = [[[eta].sub.1]4[p.sub.I][q.sub.I]([q.sup.2.sub.I] + 1)]/[[[p.sup.2.sub.I] + 1].sup.1/2 [[DELTA].sub.I],

[[gamma].sub.12] = [[-[eta].sub.1]2[q.sub.I]([q.sup.2.sub.I] - 1)[([q.sup.2.sub.I] + 1)].sup.1/2]]/ [[DELTA].sub.I],

[[gamma].sub.13] = -2[[eta].sub.2]

[[gamma].sub.21] = [[[eta].sub.2]4[p.sub.I][q.sub.I]([q.sup.2.sub.I] + 1)]/[[[p.sup.2.sub.I] + 1].sup.1/2] [[DELTA].sub.I],

[[gamma].sub.22] = [-[eta].sub.2]2[q.sub.I]([q.sup.2.sub.I] - 1)[([q.sup.2.sub.I] + 1).sup.1/2]/[[DELTA].sub.I],

[[gamma].sub.23] = 2[[eta].sub.1],

[[gamma].sub.31] = [-2[p.sub.I]([q.sup.2.sub.I] - 1)([q.sup.2.sub.I] + 1)]/[[[p.sup.2.sub.I] + 1].sup.1/2] [[DELTA].sub.I],

[[gamma].sub.32] = [-4[p.sub.I][q.sub.I][([q.sup.2.sub.I] + 1).sup.1/2]]/[[DELTA].sub.I],

[[gamma].sub.33] = 0,

[[DELTA].sub.I] = [([q.sup.2.sub.I] - 1).sup.2] + 4[p.sub.I][q.sub.I].

Then Eqs. (329) and (3.33) can be combined and written as one formula:

[G.sup.ml] = [summation over ([infinity]/n=0)] [[SIGMA].sub.g=1,2,3 j=1,2,3] [[**].sub.\[eta]\=1 [phi]=0]

[(-1).sup.n][[gamma].sup.m.sub.j]R(j,g,n)[([U.sup.[(-1).sup.n]]).sup. l.sub.g] D [phi] ([k.sub.1], [k.sub.2])P(g)ds, (3.40)

where

[phi] = -[tau]y + d - [k.sub.1][p.sub.I] - [k.sub.2][q.sub.I], (3.41)

[tau] = [c.sub.I]t/2h, y = [OMEGA]/[c.sub.I] [[xi].sub.[alpha]]/2h, d = [[xi].sub.[alpha]][[eta].sub.[alpha]], [p.sub.I] = [([[alpha].sup.2.sub.1][y.sup.2] - 1).sup.1/2], [q.sub.I] = [([y.sup.2] - 1).sup.1/2], [zeta] = z/2h, [[alpha].sub.1] = [c.sub.I]/[a.sub.I], [k.sub.[beta]] = [[delta].sup.g.sub.[beta]] (1/2 - [zeta]), [beta] = 1,2.} (3.42)

D[phi]([k.sub.1],[k.sub.2]) = 1/16[[pi].sup.2][[rho].sub.I]h(d + [k.sub.1]/[p.sub.I] + [k.sub.2]/[q.sub.I]), (3.43)

P(g) = 1/[b.sub.g]. (3.44)

5.2 Source and Detector on the Top Surface

Similarly, when the source is also located on the top surface of the layer then Eq. (3.40) may be written as

[G.sup.ml] = [summation over ([infinity]/n=0)] [[sigma].sub.g=1,2,3 j=1,2,3] [[**].sub.\[eta]\=I [phi]=0]

[(-1).sup.n] [[omega].sup.m.sub.j] R(j, g, n)[([psi]).sup.l.sub.g] D[phi]([k.sub.1], [k.sub.2])P(g)ds, (3.45)

where

[phi] = -[tau]y + d - [k.sub.1][p.sub.I] - [k.sub.2][q.sub.I]. (3.46)

[([psi]).sup.l.sub.g] is a component of the following matrix,

[PSI] = [R.sup.+][([U.sup.+]).sup.T] - [([U.sup.-]).sup.T]. (3.47)

Obviously, substituting the corresponding quantities with "*" into Eqs. (3.29), (3.33), (3.40), and (3.45) gives the formulas corresponding to case 2, i.e., the liquid coupled interface.

If we are only interested in the early time or short time arrival of the signal received by the detector, only the finite terms of all the previous integrals need be computed, because the reflected rays with more reflections do not have enough time to arrive at the detector.

Following Ref. (3), we introduce the following transform

[[eta].sub.1] = cos[theta] = d/x,

[[eta].sub.2] = sin[theta] = [+ or -][([x.sup.2] - [d.sup.2]).sup.1/2]/x

ds = d[theta]

y([theta]) = y(-[theta])} (3.48)

where x = \x\ and [theta] is the angle measured around \[eta]\ = 1 from the direction of the vector x. Note that if y is a root of

-[tau]y + d - [k.sub.1]i[(1 - [[alpha].sup.2.sub.1][y.sup.2]).sup.1/2] - [[k.sub.2]i(1 - [y.sup.2]).sup.1/2] = 0, (3.49)

for a given value of d, then -y will be a root of the following equation when d is replaced by -d:

[tau]y - d - [k.sub.1]i[(1 - [[alpha].sup.2.sub.1][y.sub.2]).sup.1/2] - [k.sub.2]i[(1-[y.sub.2]).sup.1/2] = 0. (3.50)

Thus, if d and y form a solution to Eq. (3.49) or (3.50), then d, and y (they are the complex conjugates of d, y) are solutions to Eqs. (3.50) or (3.49), while -d, and -y are solutions to Eqs. (3.49) or (3.50). When all the branch cuts for the square root functions in the above equations were taken along the real axis between each pair of branch points, using the above relations Eqs. (3.48)-(3.50) and, finally, the alternate form [phi] = [epsilon]i for the integration locus, then every term of the integrals Eqs. (3.40) and (3.45) can be simplified to the form

2[[integral].sup.x-i0.sub.-x-i0] I[y(d),d] dd/[([x.sup.2] - [d.sup.2]).sup.1/2], (3.51)

where the integrand I[y(d), d] represents the integrand of each term, including all the factors of Eq. (3.40) or (3.45). The denominator [([x.sup.2] -[d.sup.2]).sup.1/2] and the constant factor 2 are introduced by the variable transform Eq. (3.48) and the change of the integration path. All the singularities of the integrand of Eq. (3.51) are located in the upper half or on the real axis of the d-plane. For a detailed derivation of Eq. (3.51), refer to Ref. (3).

6. Automatic Generation of Ray Paths and of Corresponding Products of Reflection Coefficients

In general, there are many rays which arrive at the detector at the same time owing to the multiple reflection and the mode conversion of the incident P ray or SV ray emitted from the force source in the layer, and therefore it is not necessary to separately compute each term in the integrals mentioned above. The question of how many rays will arrive at the detector at the same time for a given configuration is a problem of combinatorics. Although explicit formulae for summing up the rays with the same arrival time to form a single integral can be derived for the case of a plate because the reflection coefficients at the top surface and bottom surface are the same (3,12), similar formulae are not easily derived for the present case of a layer on a substrate. Even though these formulae can be obtained, the actual program would also be rather complicated to carry out for a summation with reflection coefficients at the top surface which are different from those at the lower surface. In the following we will give a special counting scheme for this problem. First, let us change the notation in Fig. 2 and we obtain Fig. 3. Notation "1" and "0" in Fig. 3 respectively corresponds to notation "s" and "p" in Fig. 2. As a result it is not difficult to understand the meaning of R(0,1,1), R(1,1,2), etc. The first argument denotes the type of incident ray while the second one the reflected ray. The third arguments, "1" and "2", respectively, represent the lower surface and the top surface of the layer. Obviously, the configuration of the rays in Fig. 3 may be expressed in a sequence consisting of "0" and "1" as shown in Fig. 4. On the other hand, the sequence in Fig. 4 may also be considered as a binary notation of some integer.

To find how many rays travel through the thickness of the layer three times with p-wave speed and three times with s-wave speed, therefore, is the same problem as to find how many different six-bit binary integers can be constructed by three 1's and three 0's. Using a bit-test program it is easy to get the count and find all those binary numbers that correspond to all the ray paths with the proper p-wave and s-wave sequences. An outline of the procedures for assembling the proper sequence of the reflection coefficients for each ray path is given in the following for the case of the source and the detector located on the top surface of the layer.

1. For given [k.sub.1] and [k.sub.2], the number of layer traverses with p-wave speed and s-wave speed respectively, there is in general more than one ray path which will have the same arrival time. But each ray may have a different reflection sequence. From Figs. 3 and 4, the sequence for each ray path has a corresponding N-bit binary number representation of [k.sub.1] 1's and [k.sub.2] 0's; N= [k.sub.1] + [k.sub.2].

2. Using the bit-test program, we can determine IP, the number of distinct binary sequence with [k.sub.1] 1's and [k.sub.2] 0's, store all such N-bit binary numbers in an array, IA, of dimension IP. Each binary number will have [k.sub.1] 1's and [k.sub.2] 0's, and they are all different.

3. For each binary number in IA we can construct a product of the reflection coefficients according to Fig. 3.

4. Summing all the products thus obtained gives all the rays which arrive at the detector at the same time.

The procedure outlined above is rather efficient in terms of computation speed as long as N is limited to a small number, say less than 16. The results in the case of a plate without the lower half space have been checked against the results computed using explicit formulae previously derived (12). These same formulae for a plate cannot be used for the case of a layer on a half-space.

7. Singularities and Wave Front Arrivals

The integrand of Eqs. (3.40) or (3.45) may contain singularities, and they can be divided into three categories that respectively correspond to three different types of wave front arrivals: 1) The body waves relevant to the regular reflected rays which are determined from the singularities of [partial][phi]/[partial]y = 0,2) The interface waves, including Rayleigh surface waves, Stoneley interface waves, and the other possible leaky waves, which are determined from the singularities of [[DELTA].sub.I] = 0 or [[DELTA].sub.II] = 0 in the denominators of the reflection coefficients, and 3) The head waves determined from the singularities of the branch points of the square root functions.

Following the analysis in Ref. (3) we will show that all the above mentioned singularities of the integrands in the y-plane will be located in the upper half or on the real axis of the y-plane and their mapping onto the d-plane, which is determined from Eq. (3.48), will locate them in the upper half or on the real axis of the d-plane. So the integration path will never touch the singularities of the d-plane when we use formula (3.51) to compute the Green's functions and choose the integration path slightly below the real axis. This is indeed a main advantage of the Willis method.

Let us first consider the head wave arrivals which are related to the singularities of the branch points.

8. Head Waves

Case 1. Welded interface.

In this case we have the following four square root functions:

[p.sub.I] = [([[alpha].sup.2.sub.1][y.sup.2] - 1).sup.1/2], [q.sub.I] = [([[alpha].sup.2.sub.2][y.sup.2] - 1).sup.1/2],

(4.1)

[p.sub.II] = [([[alpha].sup.2.sub.3][y.sup.2] - 1).sup.1/2], [q.sub.II] = [([[alpha].sup.2.sub.4][y.sup.2] - 1).sup.1/2],

where

[[alpha].sub.1] = [c.sub.I]/[a.sub.I], [[alpha].sub.2] = 1, [[alpha].sub.3] = [c.sub.I]/[a.sub.II], [[alpha].sub.4] = [C.sub.I]/[C.sub.II]. (4.2)

Their branch points are respectively

[y.sub.1] = [+ or -] [[alpha].sup.-1.sub.1], [y.sub.2] = [+ or -] 1, [y.sub.3] = [+ or -] [[alpha].sup.-1.sub.3], [y.sub.4] = [+ or -] [[alpha].sup.-1.sub.4]. (4.3)

For convenience of discussion, without losing generality we assume the source and detector are located on the top surface of the layer and

[a.sub.II] > [c.sub.II] > [a.sub.I] > [c.sub.I]. (4.4)

Then from Eq. (4.2) we have

[[alpha].sub.3] < [[alpha].sub.4] < [[alpha].sub.1] < 1. (4.5)

The relation between y and d satisfies the equation

d = [tau]y + [k.sub.1]i[(1 - [[alpha].sup.2.sub.1][y.sup.2]).sup.1/2] + [k.sub.2]i[(1 - [y.sup.2]).sup.1/2], (4.6)

where

[tau] = [tc.sub.1]/2h, d = [x.sub.[alpha]][[eta].sub.[alpha]]/2h. (4.7)

A simple analysis of Eq. (4.6) shows that when all the branch cuts were taken on the real axis of the y-plane between each pair of branch points, the mapping of the branch cuts in the d-plane will be located in the upper half and on the real axis of the d-plane.

The following gives some examples of head wave arrivals corresponding to the branch points.

1. The head wave arrivals corresponding to branch point [y.sub.1].

From Eq. (4.3) we have

[y.sub.1] = 1/[[alpha].sub.1] = [a.sub.I]/[c.sub.I], (4.8)

and

y = csc[theta], (4.9)

where [theta] is the incident angle. Substituting Eq. (4.9) into Eq. (4.8) and writing [theta] as [[theta].sub.[alpha]1] we obtain

sin[theta].sub.[alpha]1] = [c.sub.I]/[a.sub.I]. (4.10)

Equation (4.10) is the critical condition satisfied by the critical angle [[theta].sub.[alpha]1] which leads to the generation of one type of head wave. The head waves of this type is named SP*S head waves. Here S and P respectively denote an SV ray and P ray in the layer. For a detailed physical mechanisms of the generation of head waves, refer to Refs. (6) and (7).

Substituting Eq. (4.8) into Eq. (4.6) and letting [k.sub.1] = 0 we obtain the normalized arrival time of the SP*S head waves,

[tau] = d [c.sub.I]/[a.sub.I] + [k.sub.2] [c.sub.I]/[a.sub.I] [[[([a.sub.I]/[c.sub.I]).sup.2] - 1].sup.1/2]. (4.11)

For example, when [k.sub.2] = 2 it is just the arrival time of the SP*S head wave, where SP*S also represents the propagation path of the head wave ray. It consists of three sections, the first section of the notation, S, represents a critically incident S ray in the layer, the second, P*, represents a P ray propagating along the interface but on the side of the layer, and the final S is a critically reflected S ray in the layer. Later on, the lower case letter p or s with a star will represent a P ray or S ray propagating along the interface on the side of the half-space.

2. The wave arrival corresponding to branch point [y.sub.2]. From Eq. (4.3) we have

[y.sub.2] = 1. (4.12)

The corresponding critical angle satisfies

sin[[theta].sub.[alpha]2] = 1, or [[theta].sub.[alpha]2] = [pi]/2; (4.13)

this condition seems as if the S ray emitted from the source propagates along the top surface of the layer. We denote this ray by S*. It is possible only if [k.sub.1] = 0 and [k.sub.2] = 0. Letting [y.sub.2] = 1 in Eq. (4.6), we obtain the normalized arrival time of this ray

[tau] = d. (4.14)

3. The head wave arrivals corresponding to branch point [y.sub.3].

From Eq. (4.3) we have

[y.sub.3] = [a.sub.II]/[c.sub.I]. (4.15)

The corresponding critical angle is

[[theta].sub.[alpha]3] = [sin.sup.-1] ([c.sub.I]/[a.sub.II]). (4.16)

In this case we may have two kinds of head waves that are excited by the same critically incident S ray but have different critically reflected S and P rays. With [k.sub.1] = 0 and using Eq. (4.15) in Eq. (4.6) we obtain the normalized arrival time of one kind of head wave

[tau] = d [c.sub.I]/[a.sub.II] + [k.sub.2] [c.sub.I]/[a.sub.II] [[[([a.sub.II]/[c.sub.I]).sup.2] - 1].sup.1/2]. (4.17)

When [k.sub.2] = 2 we obtain the arrival time of the Sp*S head wave.

Using Eq. (4.15) in Eq. (4.6) we obtain the normalized arrival time of the other kind of head wave

[tau] = d [c.sub.I]/[a.sub.II] + [k.sub.1] [c.sub.I]/[a.sub.II] [[[([a.sub.II]/[a.sub.I]).sup.2] - 1].sup.1/2] + [k.sub.2] [C.sub.I]/[a.sub.II] [[[([a.sub.II]/[a.sub.I]).sup.2] - 1].sup.1/2]. (4.18)

When [k.sub.1] = 1 and [k.sub.2] = 1. It is the arrival time of the Sp*P head wave.

4. The head wave arrivals corresponding to branch point [y.sub.4].

From Eq. (4.3) we have

[y.sub.4] = [c.sub.II]/[c.sub.I]. (4.19)

The corresponding critical angle is

[[theta].sub.[alpha]4] - [sin.sup.-1]([c.sub.I]/[c.sub.II]). (4.20)

In this case we also have two kinds of head waves. With [k.sub.1] = 0 and using Eq. (4.19) in Eq. (4.6) we obtain the normalized arrival time of one kind of head wave

[tau] = d [c.sub.I]/[a.sub.II] + [k.sub.2] [c.sub.I]/[c.sub.II] [[[([c.sub.II]/[c.sub.I]).sup.2] - 1].sup.1/2]. (4.21)

When [k.sub.2] = 2 it is the arrival time of the Ss*S head wave.

Using Eq. (4.19) in Eq. (4.6), we obtain the normalized arrival time of the other kind of head wave

[tau] = d [c.sub.I]/[c.sub.II] + [k.sub.1] [c.sub.I]/[c.sub.II] [[[([c.sub.II]/[a.sub.I]).sup.2] - 1].sup.1/2] + [k.sub.2] [c.sub.I]/[a.sub.II] [[[([c.sub.II]/[c.sub.I]).sup.2] - 1].sup.1/2]. (4.22)

When [k.sub.1] = 1 and [k.sub.2] = 1, it is the arrival time of the Ss*P head wave.

It may be seen that the branch points [y.sub.1], [y.sub.3], and [y.sub.4] correspond to the arrivals of various head wave rays, while the branch point [y.sub.2] seems to correspond to the arrival of the S ray propagating along the top surface of the layer. These head waves all are excited by the S ray emitted from the source under the conditions of Eq. (4.4). A picture of the wave fronts of the head waves in the layer generated by reflection and refraction of a spherical P wave and a spherical S wave impinging on the interface is shown in Fig. 5.

In order to describe the arrivals of various head wave rays and the P ray propagating along the top surface of the layer excited by the P ray emitted from the source, we need to change the form of the relations mentioned above. We write

y = [OMEGA]/[c.sub.I] = [a.sub.I][OMEGA]/[a.sub.I][c.sub.I] = [a.sub.I][OMEGA]/[c.sub.I][a.sub.I] = [a.sub.I]/[c.sub.I] [y.sub.p], (4.23)

where

[y.sub.p] = [OMEGA]/[a.sub.I]. (4.24)

Substituting Eq. (4.23) into Eq. (4.1) gives

[p.sub.I] = [([[beta].sup.2.sub.1][y.sup.2.sub.p] - 1).sup.1/2], [q.sub.I] = [([[beta].sup.2.sub.2][y.sup.2.sub.p] - 1).sup.1/2],

[p.sub.II] = [([[beta].sup.2.sub.3][y.sup.2.sub.p] - 1).sup.1/2], [q.sub.II] = [([[beta].sup.2.sub.4][y.sup.2.sub.p] - 1).sup.1/2], (4.25)

where

[[beta].sub.1] = 1, [[beta].sub.2] = [a.sub.I]/[c.sub.I], [[beta].sub.3] = [a.sub.I]/[a.sub.II], [[beta].sub.4] = [a.sub.I]/[c.sub.II]. (4.26)

Substituting Eq. (4.23) into Eq. (4.6) we obtain

d = [tau] [a.sub.I]/[c.sub.I] [y.sub.b] + [k.sub.1]I[(1 - [y.sup.2.sub.p]).sup.1/2] + [k.sub.2]I[[1 - [([a.sub.I]/[c.sub.I] [y.sub.p]).sup.2]].sup.1/2]. (4.27)

Obviously, the transform of Eq. (4.23) maps the original y-plane onto the [y.sub.p]-plane, whereas the original branch points of the y-plane become those of the [y.sub.p]-plane as follows:

[y.sub.p1] = [+ or -]1, [y.sub.p2] = [+ or -][[beta].sup.-1.sub.2], [y.sub.p3] = [+ or -][[beta].sup.-1.sub.3], [y.sub.p4] = [+ or -][[beta].sup.-1.sub.4], (4.28)

Similarly, we may discuss the arrivals corresponding to the branch points of Eq. (4.28)

5. The head wave arrivals corresponding to the SH ray.

Since the SH rays do not produce mode conversion and their speed is equal to that of the SV rays, the arrivals of the head waves excited by the incident SH rays emitted from the source are the same as those of the

Ss*S head waves. These head waves may be called the Hh*H head waves.

Case 2. Liquid coupled interface.

In this case the number and the distribution of the branch points are the same as those in the case of the welded interface. The arrival time of various head waves and surface P and S rays are the same as those of the welded interface. (The polarities and amplitudes have important diagnostic differences; see Computed Results and Discussion.)

9. Interface Waves

The interface waves include the Stoneley waves and the leaky waves which propagate along the interface of two different solid media (8). When one of the media is vacuum we obtain the Rayleigh waves propagating along the free surface of the solid. The Stoneley waves do not attenuate and the leaky waves attenuate when propagating along the interface. Unlike the Rayleigh waves, the Stoneley waves can exist only if some parameters of the media satisfy certain conditions (9). The singularities corresponding to the Stoneley waves come from the real roots of the so called Stoneley equation, i.e., [[DELTA].sub.II] = 0 on the top Riemman sheet of the y-plane, while the singularities corresponding to the leaky waves come from the complex roots of the Stoneley equation on the other Riemman sheets of the y-plane and therefore have some attenuation (10).

Suppose the source and the detector are located on the top surface; in this case it can be shown from Eq. (3.49) that the mapping of the real root singularities corresponding to the Stoneley waves onto the d-plane will locate them in the upper half of the d-plane, while the mapping of the complex root singularities corresponding to the leaky waves onto the d-plane will never occur on the top Riemman sheet of the d-plane. Thus only the roots of the Rayleigh equation [[DELTA].sub.t] = 0 will be considered. Letting the denominator of the reflection coefficients of the top surface be equal to zero, we obtain the following Rayleigh equation:

[([q.sup.2.sub.1] - 1).sup.2] + [4p.sub.I][q.sub.I] = 0. (4.29)

Using the principle of the argument, it can be proved that Eq. (4.29) has only two real roots, [+y.sub.R] and [-y.sub.R]. Obviously, only the positive real root is of interest. The speed of the Rayleigh waves is less than that of the S waves, and 50 [y.sub.R] < 1.

Substituting [y.sub.R] into Eq. (4.6) and letting [k.sub.1] = 0 and [k.sub.2] = 0, because the source and the detector are located on the top surface, we have

d = t[y.sub.R]. (4.30)

Since t and [y.sub.R] are real numbers, the d corresponding to the singularity [y.sub.R] should be located on the real axis of the d-plane. The arrival time of the Rayleigh waves is determined from Eq. (4.30) if d and [y.sub.R] are known.

10. Regular Reflected Rays

Analysis shows that the singularities of [partial][theta][phi]/[partial]y = 0 which relate to the regular reflected rays are located on the real axis of the y-plane and so the mapping of these singularities onto the d-plane will locate them on the real axis of the d-plane from Eq. (3.49). It is not difficult to understand this consequence if we note that y = csc [theta], where [theta] is the incident angle of the ray, and d is the wave front distance satisfying Snell's law. When d, [[alpha].sub.1], [k.sub.1], and [k.sub.2] are given, we have a group of specified rays that have the same arrival time t satisfying Eq. (4.6) and the following equation:

d = [k.sub.1]/[([[alpha].sup.2.sub.1][y.sub.2] - 1).sup.1/2] + [k.sub.2]/[([y.sup.2] - 1).sup.1/2]. (4.31)

Equation (4.31) may be obtained using Snell's law. Obviously, we must first know the value of y in order to get the arrival time t. It is seen from Eqs. (4.6) and (4.31) that this is a problem of solving a nonlinear algebraic equation in y. Following Ref. (12), the solution for y is obtained using a simple, intuitive, iterative technique. The arrival time t is determined by substituting the solution obtained into Eq. (4.6).

11. Integration Technique and Computation Procedures

The previous analysis shows that when we use Eq. (3.51) and choose the path of integration below the real axis of the d-plane, it will not touch the singularities of the integrand, which include Rayleigh poles and all the branch points of the square root functions. When all the branch cuts are taken along the real axis of the y-plane, then their mapping in the d-plane will locate them in the upper half-plane and on the real axis. This is indeed a main advantage of the Willis method.

The actual procedure for numerically computing the Green's tensor can be summarized as follows:

1. For given material properties and test configuration represented by x, the distance between source and detector, a time of arrival table is computed first. The arrivals include all the regular reflected rays and all possible head waves and Rayleigh wave; each has an associated pair of number [k.sub.1] and [k.sub.2], denoting the group of the ray paths. The arrivals are sorted according to successive time sequences.

2. To compute a particular component of the Green's tensor for a specific time, t, the number of possible arrivals, N, can be determined by a comparison with the time of arrival table computed in step 1. N is also the number of terms in Eq. (3.45) that need to be computed; each term consists of one definite integral.

3. For each integral to be numerically computed, the integrand consists of IP terms and each is a product of a unique sequence of reflection coefficients and some other factors. Both IP and the sequence series are computed by calling a bit test program as explained in Sec. 6.

4. The numerical integration is done by providing a function that computes the integrand for a given integration variable, d, along with the integration limits. The function is computed first by solving the equation [empty set] = 0 for y for a given d, then computing the IP terms one at a time and summing; each term has a factor which corresponds to the product of all the sequences of reflection coefficients mapped by a binary number. A numerical integration subroutine is applied which handles integrands with removable singularities particularly well and also provides error estimates.

5. The current program was tested by considering the case when the lower half-space is a vacuum; i.e., the case when the structure is a plate. Results obtained are identical to the results obtained by the program to compute the Green's tensor of a plate developed previously (12) which had also been checked experimentally.

12. Computed Results and Discussion

A FORTRAN program has been developed to numerically compute the Green's tensor for a layer on a half-space with three different bonding conditions between the layer and the half-space. The program is written in such a way that for given isotropic material parameters, maximum observation time, subscript of the component of the Green's tensor, number of sampling points, and distance between the source and detector, the program will compute the displacement at each sampling time. Furthermore, the arrivals of various rays are also computed and identified. The current limitation of the program is that the positions of the source and detector must be located on the top surface of the layer. However, the program will be modified to be applicable to the case of the source and detector located at arbitrary positions in the layer.

Figures 6-27 show the computed results of the components of the Green's tensor and their spatial derivatives for a plexiglass layer on a glass substrate. They were carried out on an IBM compatible personal computer. The abscissa coordinate, time, is normalized by the time required for a shear wave to vertically travel the thickness of the layer. The solid curves are the results for a welded interface condition, while the dashed curves are for a liquid coupled interface condition. The subscripts ij (11, 12,13, ..etc.) or ijk (111...etc.) indicate the response in a specified coordinate direction, i, of a point detector located on the top surface of the layer to a point force with step function time dependence exerted on the top surface in a specified coordinate direction, j. The index k denotes the specific spatial derivative direction. Physically, the spatial derivative function Gij, k can be considered as the displacement in the i-direction due to a differentiated force which is equivalent to a couple or a di pole. The distance between the source and detector is three times the thickness of the layer (except for Figs. 7 and 8) and is chosen in such a way that the very large Rayleigh wave arrivals are avoided within the given observation time, and therefore the details of the early time arrivals can be examined. The other reason for choosing this distance is to provide a basis of comparison for our experiments which we have recently conducted (15) with a geometry that is convenient to arrange and carry out. It is seen from these figures that the differences in the results with different interface conditions are very significant. The responses for a delta function time dependence can be obtained by numerical differentiation of the responses mentioned above. The results for different distances are shown in Figs. 6-8. If we compare the early time arrivals in Fig. 6 (x = 3) and Fig. 7 (x = 6), they are very different only because of the change of the distance between the source and detector. This occurs because when th e distance increases, more head waves arrive, and different arrivals change their order of occurrence. The Rayleigh wave arrivals with welded and liquid coupled interface conditions can be clearly seen from Fig. 8 and they are coincident. The early time arrivals and their differences between the two different interface conditions are totally masked by the amplitude of the Rayleigh wave arrivals. Little information about the interface can be derived.

Determination of the Green's function for a structural configuration results in a method for determining the response of the structure to a temporally varying and spatially distributed input loading. The motion of a point on the structure can be computed due to, for example, a transient contact input force. The transient waveform input is the force that would be produced by a point source transducer. For an input waveform of interest, the motional response is calculated by a point-by-point convolution of the time differential of the source waveform with the Green's function component that represents the step function response in the motional direction of interest.

The response of a structure due to a damped sinusoid input can be easily determined and may be of interest in analyzing the response to pulses used for probing materials and structures. The response of a layered structure interrogated by a pulsed laser beam can also be determined by convolving an input waveform with the derivative of the components of the Green's tensor that represents a dipole source. By comparing the response in the case of a vacuum lower half-space with the results of a welded solid half-space, the behavior of a large debond can be calculated. The large differences in amplitude and polarity of the first few head waves and regular reflected rays are strong indicators that there is no bond. These results and others will be described in a future paper. Many other responses of interest for particular applications, such as nondestructive evaluation, are envisioned.

13. Appendix A. The Willis Inversion Method

This appendix directly quotes from Ref. [3]. Consider the inversion of transforms of the form

f(t, x) = 1/8[[pi].sup.3] [[integral].sup.[infinity]-0i.sub.-[infinity]-0i] d[omega] [integral] [integral] d[[xi].sub.1]d[[xi].sub.2] [SIGMA] F ([omega], [xi])

exp[-i([[xi].sub.[alpha]][x.sub.[alpha]] - [omega]t + [[lambda].sub.1][[xi].sup.M.sub.3] + [[lambda].sub.2][[xi].sup.N.sub.3])] (A.1)

where F([omega],[xi]) is a homogeneous function of degree -2. Each term of the integrand, which represents a generalized ray in the text, is analytic for Im([omega]) <0 and the integration of each term vanishes for t < 0. So f(t, x), i.e., the integration of the sum of the integrands should vanish for t < 0. The functions [[xi].sup.M.sub.3] and [[xi].sup.N.sub.3] homogeneous of degree 1 in [omega], [[xi].sub.1], and [[xi].sub.x]; [[lambda].sub.1][[xi].sup.M.sub.3] and [[lambda].sub.1][[xi].sup.N.sub.3] have imaginary parts.

Setting

[OMEGA] = [omega]/\[xi]\, [[eta].sub.[alpha]] = [[xi].sub.[alpha]]/\[xi]\, (A.2)

where \[xi]\ = [([[xi].sup.2.sub.1] + [[xi].sup.2.sub.2]).sup.1/2], gives

f(t, x) = 1/8[[pi].sup.2] [lim.sub.[epsilon][right arrow]0] [[integral].sup.[infinity]-0i.sub.-[infinity]-0i] d[OMEGA] [[**].sub.\[eta]\=1 ds [SIGMA] F ([OMEGA], [eta]) [[integral].sup.[infinity].sub.0] d\[xi]

exp[--i\[xi]\(-- [OMEGA]t + [[eta].sub.[alpha]][x.sub.[alpha]] + [[lambda].sub.1][[xi].sup.M.sub.3] + [[lambda].sub.2][[xi].sup.N.sub.3] - i[epsilon])] (A.3)

where [[xi].sup.M.sub.3] = [[xi].sup.M.sub.3]([OMEGA], [eta]). The "convergence factor" [e.sup.-\[xi]\] is inserted to facilitate evaluation of the integrals in Eq. (A.3) by simple quadrature; generalized functions are continuous with respect to limiting operations of the type introduced, and the limit [epsilon] [right arrow] 0 can be taken at any convenient stage. Evaluating the integral with respect to \[xi]\ gives

f(t, x) = 1/8[[pi].sup.2]i [lim.sub.[epsilon][right arrow]0] [SIGMA] [[**].sub.\[eta]\]ds [[integral].sup.[infinity]-0i.sub.-[infinity]-0i]

d[OMEGA](F ([OMEGA], [eta])/-- [OMEGA]T + [[eta].sub.[alpha]][x.sub.[alpha]] + [[lambda].sub.1][[xi].sup.M.sub.3] + [[lambda].sub.2][[xi].sup.N.sub.3] - i0). (A.4)

Setting

[PHI] = - [OMEGA]t + [[eta].sub.[alpha]][x.sub.[alpha]] + [[lambda].sub.1][[xi].sup.M.sub.3] + [[lambda].sub.2][[xi].sup.N.sub.3] - i[epsilon] (A.5)

the integral with respect to [OMEGA] is now evaluated by closing the contour in the lower half-plane because the integrand has behavior as O([[OMEGA].sup.-2]) for large [OMEGA] and using Cauchy's theorem. This gives

f(t, x) = -1/4[[pi].sup.2] [SIGMA] [[**].sub.\[eta]\=1 [phi]=0] ds [F(OMEGA], [eta])/-t + [[lambda].sub.1][[xi].sup.M.sub.3],[OMEGA] + [[lambda].sub.2][[xi].sup.N.sub.3],[OMEGA]\.sub.[phi]=0], (A.6)

where [OMEGA] is now taken as the root in the lower half-plane of the equation

-- [OMEGA]t + [[eta].sub.[alpha]][x.sub.[alpha]] + [[lambda].sub.1][[xi].sup.M.sub.3]([OMEGA], [eta]) + [[lambda].sub.2][[xi].sup.N.sub.3]([OMEGA], [eta]) = 0i, (A.7)

if there is such a root. If there is more than one root, then the contribution from each root should be included; if none, then that term of the integrand is replaced by zero.

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

[FIGURE 14 OMITTED]

[FIGURE 15 OMITTED]

[FIGURE 16 OMITTED]

[FIGURE 17 OMITTED]

[FIGURE 18 OMITTED]

[FIGURE 19 OMITTED]

[FIGURE 20 OMITTED]

[FIGURE 21 OMITTED]

[FIGURE 22 OMITTED]

[FIGURE 23 OMITTED]

[FIGURE 24 OMITTED]

[FIGURE 25 OMITTED]

[FIGURE 26 OMITTED]

[FIGURE 27 OMITTED]

Accepted: September 20, 2002

14. References

(1.) Y. H. Pao and R. R. Gajewski, Chap. 6 in The generalized ray theory and transient responses of layered elastic solids, in Physical Acoustics, Vol. 13, W. P. Mason and R. N. Thurston, eds., Academic, New York (1977).

(2.) J. R. Willis, Self-similar problems in elastodynamics, Phil. Trans. Roy. Soc. A247, 435-491 (1973).

(3.) J. A. Simmons, J. R. Willis, and N. N. Hsu, The dynamic Green's tensor for an elastic plate, Research Report (Draft), National Bureau of Standards, Washington, DC 20234 (1981). Sec also D. G. Eitzen, F. R. Breckenridge, R. B. Clough, E. R. Fuller, N. N. Hsu, and J. A. Simmons, Fundamental Development for Quantitative Acoustic Emission Measurements, EPRI NP-2089, Project 608-1 Interim Report, Prepared by National Bureau of Standards, Washington, D.C., Prepared for Electrical Power Research Institute, 3412 Hillview Avenue, Palo Alto, California 94304, October 1981.

(4.) R. J. Bedding and J. R. Willis, The elastodynamic green's tensor for a half-space with an embedded anisotropic layer, Wave Motion 2, 51-62 (1980).

(5.) M. Shmuely, Response of plates to transient sources, J. Sound Vibration 32, 491-506 (1974).

(6.) L. M. Brekhovskikh, Waves in Layered Media, Academic Press, New York (1960).

(7.) V. Cerveny and R. Ravindra, Theory of Seismic Head Waves, University of Toronto Press, Toronto (1971).

(8.) W. L. Roever, T. P. Vining, and E. Strick, Propagation of elastic wave motion from an impulsive source along a fluid/solid interface, Phil. Trans. Roy. Soc. A251, 255-523 (1959).

(9.) A. S. Ginzbarg and E. Strick, Stoneley-wave velocities for a solid- solid interface, Bull. Seism. Soc. Am. 48, 51-63 (1958).

(10.) W. L. Pilant, Complex roots of the Stoneley-wave equation, Bull Seism. Soc. Am. 62, 285-299 (1972).

(11.) J. D. Achenbach, Wave Propagation in Elastic Solids, North-Holland Publishing Company, Amsterdam (1973) p. 189.

(12.) N. N. Hsu, Dynamic Green's Functions of an Infinite Plate--A Computer Program, NBSIR 85-3234, U.S. Department of Commerce, National Bureau of Standards, Gaithersburg (1985). Available through NTIS PB 86-143856/AS.

(13.) D. A. Hutchins and L. F. Bresse, Laser generation of ultrasound, Ultrasonics International 87 Conference Proceedings (1987) p. 135.

(14.) S. I. Rokhlin and D. Marom, Study of adhesive bonds using low-frequency obliquely incident ultrasonic waves, J. Acoust. Soc. Am. 80, 585-590 (1986).

(15.) F. R. Breckenridge, N. N. Hsu, and D. G. Eitzen, Experimental waveforms for a layered solid half-space with different interface conditions, NISTIR, to be published.

Shu-Chu Ren (1)

(1.) Deceased; was on leave from Institute of Acoustics, Academia Sinica, Beijing, People's Republic of China.

About the authors: Shu-Chu Ren was a physicist in the NIST Manufacturing Engineering Laboratory, on leave from the Institute of Acoustics, Academia Sinica, Beijing. Nelson N. Hsu is a guest researcher and Donald G. Eitzen is a mechanical engineer in the Mechanical Metrology Division of the NIST Manufacturing Engineering Laboratory. The National Institute of Standards and Technology is an agency of the Technology Administration, U.S. Department of Commerce.

*[Text unreadable in original source]