Edge Diffraction Coefficients around Critical Rays

The classical GTD (Geometrical Theory of Diffraction) gives a recipe, based on high-frequency asymptotics, for calculating edge diffraction coefficients in the geometrical regions where only diffracted waves propagate. The Uniform GTD extends this recipe to transition zones between irradiated and silent regions, known as penumbra. For many industrial materials, e.g. steels, and frequencies utlized in industrial ultrasonic transducers, that is, around 5 MHz, asymptotics suggested for description of geometrical regions supporting the head waves or transition regions surrounding their boundaries, known as critical rays, prove unsatisfactory. We present a numerical extension of GTD, which is based on a regularized, variable step Simpson's method for evaluating the edge diffraction coefficients in the regions of interference between head waves, diffracted waves and/or reflected waves. In mathematical terms, these are the regions of coalescence of three critical points - a branch point, stationary point and/or pole, respectively. We show that away from the shadow boundaries, near the critical rays the GTD still produces correct values of the edge diffraction coefficients.


Introduction
The geometrical theory of diffraction (GTD) is an extension of the standard ray theory. According to GTD, in the high frequency approximation, i.e. under the assumption that the distance from the observation point to the edge diffrcation point as well as both the curvature of the wave front and crack edge at this point are larger than the wavelength, a wave incident on a sharp edge produces a cone of diffracted rays. Diffraction coefficients describe the directivity patterns of the resulting diffracted fields.
GTD is used widely in a variety of applications and in particular, in ultrasonic NDE (Non-Destructive Evaluation). Formulae for evaluating elastodynamic diffraction coefficients for a semi-infinite planar crack in an isotropic homogeneous solid have been presented e.g. in [1] and [2]. A code based on these formulae was incorporated into CIVA, an NDE simulation software developed by CEA [3]. A new simpler version of these formulae has been presented recently in Gautesen [4] and also incorporated into CIVA [5].
The aim of this paper is to describe a method for evaluating the diffraction coefficients of the half-plane crack in an isotropic homogeneous solid -in the regions of interference between head waves, diffracted waves and/or reflected waves, where the standard GTD fails. In mathematical terms, the regions under consideration are regions of coalescence of three critical points -a pole, stationary point and branch point, respectively. The evaluation is carried out numerically, using a regularised, variable step version of Simpson's method.

Statement of the canonical edge diffraction problem
Let us assume that the three-dimensional space is filled with an isotropic homogeneous solid supporting two speeds, c L , the speed of longitudinal waves and c T , the speed of transverse waves. In the presence of a scatterer, u, the total displacement field in this solid can be decomposed into the sum of the incident field u i and scattered field, u sc . In the high frequency approximation, the field scattered by a scatterer with an edge is the sum of the reflected field u r and the diffracted field u d , which describes the effect of the edge, so we can write u = u i + u sc = u i + u r + u d . Each of the components of the geometrical elastodynamic field u i + u r vanishes in the respective shadow (silent) zone. A transition zone in between an irradiated region and the corresponding silent zone is called a penumbra.
Let us now assume that the scatterer is a semi-infinite planar crack. Let us introduce the Cartesian coordinate system such that it lies in the (x 1 , x 2 )plane. Let the crack edge coincide with the x 2 -axis and let the crack contain the positive x 1 -axis (see Fig. 1). Let it be irradiated by a time-harmonic wave represented by U(x, t) = u(x) exp (iωt), where i is the imaginary unity; x = (x 1 , x 2 , x 3 ) is an observation point; t the time; and ω is the circular frequency. Below only time-harmonic waves are considered and the factor exp (iωt) is implied but omitted everywhere. The following conventions are followed above and everywhere below: The bold font is used to denote vector quantities. Whenever appropriate, these symbols are supplied with superscripts to describe the nature of the wave. In particular, the descriptors i and d introduced above (and sometimes put in curly brackets below) denote, respectively, the incident and diffracted waves. The descriptors α, β = L, T H, T V or R introduced below refer to longitudinal, transverse horizontal, transverse vertical or Rayleigh wave, respectively, with α normally used to denote the incident wave mode and β, the scattered one. Subscripts are used to refer to the corresponding vector components or, when dealing with scalars, to the wave mode.
To continue, the equation of motion for this problem can be written as where ∇ = (∂/∂x 1 , ∂/∂x 2 , ∂/∂x 3 ) is the nabla operator with ∂ -the symbol of partial differentiation; σ(x) is the stress tensor; and ρ is the density of the medium. Using Hooke's law, the Cartesian components of the stress tensor are where λ and µ are the Lamé parameters and δ ij is the Kronecker delta.
A realistic boundary condition [2] prescribes traction-free crack faces, where t(x) = ν · σ(x) and ν is the inner normal to the upper face pointing into the solid. It is further assumed that the scattered waves satisfy the standard radiation condition at infinity demanding that they be outgoing.

Solution of the canonical edge diffraction problem
Solving the above problem shows that the isotropic medium supports two types of plane waves u = Ad exp (ikp · x), L (longitudinal), also known as P (compressional) waves and T (transverse), also known as S (shear) waves. The polarization vector (unit displacement vector) d L is parallel to propagation vector p. The T wave is decomposed into a TV vertical component and a TH horizontal component with respect to a plane. Choosing this plane to be a crack plane and using the decomposition in [4], the corresponding polarization vectors are It is shown in [4] that the single Fourier transform of the scattered field u sc (ξ 1 , x 2 , x 3 ) can be represented aŝ where the single Fourier transform of f (x 1 ) iŝ the expressions for vectors d and v can be found in [4] and are not reproduced here. We just emphasize that they depend on The scattered field can be evaluated using the Inverse Fourier Transform In the high-frequency approximation the above integrand is a product of a fast-oscillating exponential factor and often a slowly varying amplitude. Such integrals can often be evaluated asymptotically and when the phase stationary point is isolated each diffracted field can be expressed as where exp(ik L ξ α{i} · x β{d} ) is the incident wave, the diffraction point , sγ β (ξ β{d} 1 )); r β is the distance from the diffraction point to the observation point x, so that we If the stationary point is not isolated and is coalescing with the branch point or pole or both, then comparing (7) and (8) and projecting (8) onto the displacement vector d β (ξ β ) we can still represent the scattered filed in the form (8) but in this case the scattering coefficient S α β can be calculated using where the far field parameter a = k β r β and the phase is given by for the incident pole and ξ br = −γ L (0) = − 1 − ξ 2 2 for the branch point. We denote the angles of the incoming wave by φ i and θ i , and the angles for the scattered wave by φ d and θ d , using the spherical coordinate system shown in Fig. 1.

A new numerical integration technique
The integral in (10) can be evaluated using Simpson's rule adjusted to accommodate the presence of the fast oscillating exponent, a pole and branch point. We describe each of the aspects of the adjusted Simpson's method below.

Integration interval
Since we are interested in head waves and their interaction with other waves we concentrate on the T scattered field only. The integration in (10) is carried over the entire real axis and a typical integrand is presented in Fig. 2. We note that the phase f the sense that the derivative of the phase (11) increases without bound as we approach the endpoints, see Fig. 2, and when these end points are included numerical integration becomes unstable. Another reason to avoid the endpoints is the presence of γ T (ξ 1 ) in the denominator of the integrand in (10). Therefore, when the diffraction angle is away from grazing we carry out numerical integration over the interval where the offsets δ ± = δ · min (|γ T (0) ∓ ξ st |, |γ T (0) ∓ ξ br |) and δ is a prescribed small number. The choice assures that the main critical points ξ st and ξ br lie inside the integration interval, the excluded regions are small. The real and imaginary parts of the integrand oscillate around zero fast, thus making the neglected contribution very small. When the incident pole lies close to the lower/upper bound −γ T /γ T , we redefine δ ± as δ ± = 0.001δ · (|γ T (0) ∓ ξ p |) . When the diffraction angles are near grazing the interval contains no phase stationary point and we choose δ ± to be δ ± = δ · |γ T (0) ∓ ξ p |. Using numerical experimentation we have established that in both cases the most robust choice for δ is 0.1. When the interval I contains neither stationary point nor pole we choose δ ± = ±10 −9 . When all three critical points coalesce, numerical experimentation suggests the choice δ ± = ±10 −6 .
• The integrands decay exponentially outside the interval I. We make allowance for this by following advice given in [6], that is, integrating

Standard Simpson's rule
The standard Simpson's rule for a function H (ξ) and the integration interval specified above is with the step size h = (ξ + − ξ − ) /n and nodes ξ j = ξ − + j · h.

Variable step Simpson's rule
As mentioned above, as we approach the accumulation points, the number of cycles in our highly oscillating integrands increases without bound. To accommodate this behaviour we implement a variable step version of Simpson's rule. Namely, we halve I δ , and then every new subinterval, repeating the process until 2(n+1) subintervals have been generated, with n specified by the user. For every subinterval, we estimate the length of the smallest cycle of the integrand there. The estimate is given by 2π/f (ξ 1 ) as calculated at the point of the highest oscillation. In accordance with the usual numerical practice, each cycle should be covered by at least n eval 10 nodes. When dealing with subintervals, which lie close to the accumulation points we further multiply the number of nodes there by 1/γ T (ξ 1 ).

Regularized Simpson's rule
When an integrand H(ξ) contains a pole (say, the geometric pole ξ p = −ξ α 1 , see (10) we can write it in the form H(ξ) = G(ξ)/(ξ − ξ p ) and regularise the integral as follows: This recipe fails when the pole is close to a node. When it is we shift the whole integration interval so that the pole lies half-way between two nodes.

Contribution of the pole
The integral (13) can also be seen as that is, the sum of the Cauchy Principal Value integral (evaluated using the regularised Simpson's rule) and πiG (−ξ α 1 ), half residue, the integral over a small half-circle passing under the pole −ξ α 1 . Let us justify this choice of the interation path: in the vicinity of the pole the phase can be represented as f (ξ) ≈ f (ξ p ) + f (ξ p )(ξ − ξ p ). This approximation works when the integrand is fast-oscillating, leading to cancellations outside the neighbourhood of the stationary point. Therefore we can extend the limits of integration beyond the limits of validity of the above Taylor expansion. Taking into account that the phase contains one maximum point (f (x st ) < 0), it can be justified more rigorously by introducing a new variable that allows to represent the phase as a quadratic [7]. One way to evaluate the resulting integral is by closing the integration contour in the upper ξ 1 half-plane. This choice is dictated by the radiation condition at infinity: when the pole ξ p is isolated and f (ξ p ) > 0 the Jordan Lemma can be applied to evaluate its contribution to give the outgoing reflected wave. When f (ξ p ) < 0 the Jordan Lemma is not applicable and the pole makes no contribution. Since the stationary point is a maximum, f (ξ p ) > 0 when ξ p < ξ st .
When wishing to describe the diffracted field only and ξ p < ξ st the full residue, 2πiG (−ξ α 1 ) should be subtracted from the result, thus removing the asymptotic limit to the reflected field [6], see Fig. 3.

Contribution of the branch point
When evaluating the transverse scattered field the third component of each v β in (10) contains the factor γ L (ξ 1 ) and its domain contains a branch cut that runs along the integration contour from the branch point ξ br = −γ L to the left. When the branch point ξ br is isolated or coalescing with the stationary point and f (ξ br ) ≥ 0 the Jordan Lemma can be applied to establish that the branch point makes no contribution. When f (ξ br ) < 0 the Jordan Lemma is inapplicable. Thus, when f (ξ br ) < 0, even when integrating over the upper side of the cut, the neighbourhood of the branch point makes a contribution. This contribution is known as the head wave. Since the stationary point is a maximum, f (ξ br ) < 0 (and therefore, the head wave contribution) is non-zero when ξ br > ξ st . The situation becomes more complicated when the branch point coalecses not only with the stationary point but with the pole as well. When f (ξ p ) < 0 and f (ξ br ) > 0 and integration is carried out over the upper side of the cut (equivalent to closing the contour in the upper half-plane) neither the pole no branch point make a contribution. It follows that when the integral (10) is evaluated along the upper side of the cut the results can be contaminated with the contribution of the head wave.

Numerical Results and Discussion
Our code is specifically designed to calculate the scattering coefficient when a stationary point coalesces with the branch point ξ st ξ br or else when all three critical points ξ st , ξ br and ξ p coalesce. For this reason it employs the modified Simpson's method only when the stationary and branch point are coalescing, ξ st ξ br , relying on the classical GTD otherwise. We refer to it as the mixed GTD/numerical code. We say that the stationary point and branch point are coalescing when the phase (11) satisfies condition |a · f (ξ st ) − a · f (ξ br )| < 6π. For validation purposes, below we use only the full numerical model, which utilises the Simpson's method at all observation angles. At observation angles near grazing, where we have 0 • , 180 • , the calculations are always carried out using GTD. Everywhere we choose δ = 0.1, n eval = 10 and 2n = 10. Let us illustrate the performance the code in the region of interference of the reflected, diffracted and head waves. In all examples below c P = 5890 m · s −1 and c S = 3210 m · s −1 as in mild steel.
In Fig. 4 we compare performance of the full numerical model with GTD in the transition region where the head wave interacts with both a diffracted wave and reflected wave. The figure shows that at critical angles, at smaller distances from the diffraction point the numerical method predicts a smaller value of the scattering coefficient than GTD. Other authors, see e.g. [8] and [9], report similar behaviour for diffraction by a finite strip. As explained in [9], for the finite strips the smoothing of the critical spikes in the diffracted field is due to interference between the head waves and diffracted waves generated by the strip edges as well as the reflected Rayleigh waves. We believe that a similar phenomenon explains our results: at smaller distances where the head wave is still significant there is destructive interference between it, diffracted wave and maybe the reflected wave. This conclusion is supported further by Fig. 4b) which demonstrates that for larger values of the far field parameter, the full numerical method produces the scattering coefficient with the peaks at critical angles, which are higher than for a = 2π × 10. They are smaller than those in GTD coefficients, since the latter fails in penumbras and in the reported configuration shadow boundaries and critical rays coalesce. Other numerical experiments carried out where GTD applies, away from the shadow boundary, confirm that the GTD peaks at critical angles are correct: they are the same as the peaks obtained with the full numerical model.