Mathematical Modelling of Ultrasonic Non-Destructive Evaluation

High-frequency asymptotics have been used at our Centre to develop codes for modelling pulse propagation and scattering in the near-field of the ultrasonic transducers used in NDE (Non-Destructive Evaluation), particularly of walls of nuclear reactors. The codes are hundreds of times faster than the direct numerical codes but no less accurate.


Introduction
Our Centre specialises in mathematical modelling of NDE based on high- frequency asymptotics.We have been the first 1-3 to produce the complete asymptotic description of the time-harmonic near field of a circular com- pressional transducer which is directly coupled to isotropic solid, both its geometrical regions and boundary layers in between the geometrical regions (see Fig. 1).Pulse propagation, rectangular transducers and transducers of complex apodisation have been also modelled using this approacha-6.The crux of the method is approximation of integrals containing an ex- ponential factor, such that when observation point moves across the near field, the factor undergoes many oscillations while the amplitudes varies slowly.The main contributions to the integrals of this type come from the critical points: singularities of the amplitude, PSP (the stationary points Due to space limitation this article was omitted from theprevious special issue on McNabb Symposium. Requests for reprints should be sent to Larissa Ju Fradkin, Centre for Waves and Fields, School of EEIE, South Bank University, London SE1 0AA, England b) Figure 1.a) The wave fronts and b) boundary layers underneath a circular compressional transducer.Thick solid linedirect, thin solid lines edge P, dashed lines edge S, and dotted lines edge head wave.P penumbra, AH, AC and zone between lines marked A axial region, C boundary layers around the critical rays.
of the phase functions) and various types of critical boundary points.Such contributions may be evaluated using various formulae of the USPM 9-1   (Uniform Stationary Phase Method), depending on whether the critical points are isolated or coalesce.In the first case the exponential factor undergoes at least several oscillations in between the points, and in the second it does not.The resulting time-harmonic field contains the geo- metrical zones and boundary layers in between.The geometrical zones are described by ray-asymptotic series (in inverse powers of dimensionless wave number) which are contributions of isolated critical points.The boundary layers are described by boundary layer asymptotic series (in other types of functions of dimensionless wave number) which are due to coalescing critical points.We illustrate the approach in the next section.

Ultrasonic modelling of an elliptic crack
Let an elliptic crack be irradiated by a plane harmonic P wave u P(i'c) n P(inc) eik'(''c)'x, where (inc) stands for incident; a P for compressional and S for shear; wave number k, w/ca, with w, the frequency and ca, the speed of c-wave.Inside the solid, the total harmonic displacement field u u (i'c) -t-u(scat), where (scat) stands for scattered is described by the reduced elastodynamic equation where the total stress tensor o'ij ,,ij "eii -k-21zeij; $ij is the Kronecker delta; ,k and / are the Lame constants; the total strain tensor is eij 1 -[ui,1 + uj,i], i,j x, y,z whatever the orientation of x-and y-axis; index before the comma denotes a component, and index after the comma denotes differentiation with respect to the corresponding spatial variable.The usual radiation condition at infinity and Meixner's condition at the edge are assumed.We can use the well-known scalar and vector potentials, 0 and (1, 2, 3), respectively and write u e V0, u s 37 .(3) First let us recall the 2.1 Canonical problem Let a plane compressional wave scatter off a semi-infinite crack.Let us introduce the standard Cartesian coordinate system with the x-axis running along the edge of the crack and the y-axis directed into the crack half-plane.
The corresponding spherical coordinates are (r, , 0).Then the unit wave vector n p(inc) is n P(inc) (sin0 (inc) coso(i'C),sinO(inc) sino (in), cos 0(inc)).( 4) Without loss of generality let us assume that the incident wave propagates from above.The problem has been studied by many authors7.We shall employ respective even and odd parts of the potentials of the scattered field as derived in Achenbach et al.r, f A+e () iu+iT()lzld e i" f A[() /""(x) sgn()---- where a total potential component is 'e The descriptor (can) refers to the canonical problem; the projection of the incident wave-vector k P(inc) onto the x-axis is r/= kP(inc)= Equations (5) represent potentials as integrals over plane waves exp(ik x) with complex wave vectors k (,-,gn()()), (O) so that we have ,a() v/k2 /2 2. The functions A() are regular in the upper half -plane and are given in Appendix A. The corresponding scalar and vector potentials involve combinations At() A+()+sgn (z) "Y() AT() .(7) s where et ex,% or ez for 1,2 or 3 respectively.Using (3) the respective displacements involve vectors 3 AP()nP(') A()nP()' AS()n()= Z At()nS() et, (8)  where A s and n are both defined by the above formula.The functions have a pole _inc _kP(inc) ( The integration contour in (5) passes it from above.The other critical points are ff + z" (l) The descriptor (GTD) in ( 10) is justified below.For simplicity of presen- tation, we consider the scattered field only in the upper half-space z > 0.
If the observation point in the upper half-space is such that the pole (9) and the phe stationary points (10) are isolated from each other, we get two contributions: Applying to (5) the standard formulae for the pole contribution (orovikov 9 1994, Sec.1.6) and using (3), the leding order tems are k(re]) /ka(l)'x, (11) where H(x) is a Heaviside nction, the reflected wave vector is Above, the unit reflected wave vector is n a(reI) (sin O a(eI) cos (eI), sin 0a(e) sin (reI), cos 0a(eI)); Rp and Rs are standard reflection coefficients and the respective polarisation vectors pP(ref) llP(ref) and pX(ref) n (ref), with the unit vector of S-displacement being nS (ref) (-cos 0 s(ref) cos o (tel) -cos0 s(ref) sin(rel),sin0S(rel)). Combining ( 4) and (12) it follows that the above angles are defined by ka sin 0 a(ref) kp sin 0(inc), o(ref) Let us now consider a point x a(rel) (xa(reI), ya(reJ:), 0) lying in the plane of the crack, such that we have ko where distance d a(rel) Ix-xa(reI) I.This vector equation results in two scalar equations in two unknowns, X a(re]) and ya(reI).It follows that we have k (tel)   k (ref)   d,(rel) It is clear that if y zk(reI)/k(reI), we get y(reI) 0, that is x lies on the crack edge.Thus, the argument of the Heaviside function in ( 11) is positive when x a(e/) lies on the crack, and negative when it lies in the same plane , but outside, the crack; in both ces the formula applies only when the observation point lies far from the edge.Note that using (12) and ( 14), the total phe may be decomposed into incident and reflected follows: k a(rl) x k P(inc) -x a(rel) + kad a(rel). (16) Thus, the incident ray that hits this point gives rise to the reflected ray which reaches x.In the geometrical acoustics such points are called spec- ular.
It follows that the formula (11) describes the well-known main reflected c-beams, and their geometrical shadows.The beams obey the Shell's law (13).Note also that neither the position of x a(ref) nor the contribution (11) depend on the position or orientation of the crack edge, and (11) constitutes an exact solution of the problem of reflection from an infinite plane.
To simplify the formulae below, let us introduce instead of (r, ,0), the spherical coordinates (r, , ) corresponding to the Cartesian system (e, ez, ex).
Applying to (5) the standard asymptotic formula of the stationary phase method (Borovikov9, 1994, Sec.1.2) and using (3), to the leading order in ksr we obtain Da "pa(GTD)'eika(GTD)'xTir/4 [y-cot Zga(ref where k a(GTD) (k P('c)x ks sin 2 (vTD) cos , ks sin 2 (GTD) sin #) \ ! is the angle between wave vector k and the x-axis; and the gradient of the geometrical a-shadow boundary is given by cot #(re]) k(rel)/kz(el)" It follows that we have ks cos t(GTD) kp cos P(inc) (18) Note that using (10) and ( 12), we also have -k () The conditions in (17) are refined below in the section on the penumbral contribution.To continue, the diffraction coefficients in (17 V/ sin t( GTD coS t( f COS c) the geometrical spreading is ja V/y2 + z 2 sina(GTD); (21) and the respective polarisation vectors are pP(GTD) IIP(GTD) and pS(GTD) nS(reI) (_k(eTD)) _k Let us now consider a point X a(GTD) (Xa(GTD), 0, 0) lying on the edge of the crack, such that we have ka(GTD da(GTD (X Xa(GTD)). (22 The above vector equation results in one scalar equation for X a(GTD).Thus, we have x a(GTD) x-V/y 2 d-z 2 cot a(GTD), da(GTD) V/Y + z2 sina(eTD).(23) Thus, an incident ray that hits this point gives rise to a family of diffracted rays lying on the surface of a cone with the tip at x a(GTD} the axis running along the x-axis and half-angle a(GTD) which are related to P(inc) by the Snell's law of diffraction (18).In GTD such points are known as flash points.Note also that the total phase may be decomposed into incident and diffracted as follows: k (GTD) x k P(inc) x a(GTD) -k kad (GTD). (24) It follows that (17) describes the well-known in GTD edge diffracted waves.
To summarise, in the geometrical regions the high-frequency asymptotics of the scattered field are represented by a sum of the first terms in the GE and GTD ray series uo(x) u(x) + ur(x), where ua(GE)(x), zero in the shadows, is independent of the the location and orientation of the edge and U(GTD)(x) is a term of a higher order in dimensionless frequency.
When the observation point is such that the pole (9) and the phase station- ary point in (10) coalesce, applying to (5) the corresponding asymptotic formula of the stationary phase method (Borovikov9, 1994, Eq.
(2.28)), the leading order contribution is Ua(pen)(x F(Va)Ua(GE)(x + pa(GTD)eik(GTD).x+ir/25)v/ka J where (va)2 _< r, the first term in the right-hand side is expressed as usual via the Fresnel integral f(v and the modified diffraction coefficients /a as in the formulae (20) but with At(a) replaced by At(a) At(-inc).It is clear that v2a is the dif- ference between the phases of the diffracted edge and reflected ray reach- ing the observation point x.The coalescence of the pole (9) with the phase stationary point (10), that is the coalescence of the corresponding specular and flash point occurs when the observation point x crosses a geometrical a-shadow boundary and thus the reflected ray coalesces with the diffracted ray to form the so-called marginal ray, (ref), on the diffraction cone.The transition zone (boundary layer) surrounding this boundary is known as penumbra.As usual, it is defined by (va) 2 _< , so that inside it the phase difference is less than r and the exponential factor in (5) ceases to be rapidly oscillating.To find the boundary of penum- bral region we note that using (18) k a(ref) can be re-written as k a(ref) (kc cos a(GTD), k( sin a(GTD) COS a(ref), k( sin -a(GTD) sin 9((ref)).There- fore, combining the above formula with (12), the equation of the penumbra boundary becomes 2ka v/Y 2 + z 2 sinf a(GTD) sin2(--4a(re/))/2 7r, which is the equation of a parabolic cylindrical surface surrounding the plane ge- ometrical a-shadow boundary.The corresponding conditions in (11), ( 17) and (25) may be refined accordingly to specify that these formulae apply inside or outside the penumbral region respectively.
Evaluating asymptotics of (25) when va --+ c we obtain the behaviour which is entirely consistent with the expected physical picture: the GE reflected and edge diffracted P and S waves in the main beam zone and just edge diffracted P and S waves in the geometrical shadow zone.
To summarise, in the penumbra the high-frequency asymptotics of the scat- tered field may be represented by the first two terms of the penumbral series, u a(pen) (x).
2.2 Scattering of a plane P wave from a planar crack with an elliptic edge Let us now consider scattering from an elliptic crack.Let the new Cartesian coordinate system be (x, y, z), with the origin at the ellipse centre, the x- axis running along the crack major semi-axis and the y-axis, along its minor axis.Let the corresponding spherical polar system be again (r, , 8).The parametric equations of the elliptic crack edge are x'(a) acosa, y'(a) b sin a, with the semi-axes a >_ b and parameter a E [0, 27r].
We proceed by deriving the radiating near field asymptotics (r >> (kp) -1 and r << kpab.)We start by writing the decoupled form of the classical Green's formula / 0 / OG (x,x')/(x')dx'dy' t+(x)-G+(x,x')-z,(X')dx'dy' /(x)-- where the region D is the crack face; O/Oz is a derivative with respect to the outward normal to the crack face; el(x) and G(x, x') are even and isotropic solid Z uP(inc) I x Figure 2. Scattering from an elliptic crack.

First tier
In order to build the first tier asymptotics, we first note that in the canon- ical case, when the pole inc is isolated, and hence the specular points x a(ref) are far from the edge, on substituting ( 5) into (26) the resulting integrals may be approximated by contributions of x a(rel) and the crack edge.These contributions may be obtained using a separation of unity {X1, X2, X3}, where the neutraliser Xx has support in the neighbourhood of xa(el); X2, in the neighbourhood of the crack edge, and X3 1 X1 Therefore, when scattering from the elliptic crack, we can use Assumption 1-When a specular point x a(Spec) is far from the ellipse edge, the potentials may be approximated as follows and, when the specular point is near the ellipse edge, we have cat(x dilI)(x). (32) Since, the expreions in (11) do not depend on the position or location of the crack boundary either, VE) are such that we have (()) u(O)(x) HI1-Rp(rI) ek"('I) .x.(33 We remind the reader that in this paper we present results only for z > 0. On the other hand, the modified diffracted fields in (31) and (32) are to the leading order in kar, contributions to the following integrals D Oa D where D is the elliptic crack.The contributions implied in (al) and are different, since hey depend on the position of he specular point.We evMuate them approimately with the help of often used Assumption (The localisation principle): In the high-frequency approximation we sume that near the crack edge, to the leading order in kPpmin, the leMing terms of the scattered field are the same in the corresponding canonical problem, with the edge of the semi-infinite crack tangent to the elliptic crack edge at this point.
For convenience, in the neighbourhood of the edge we employ local coordi- nates (o, u), where u is the distance Mong the internal normal to the edge along, with the two unit vectors er(o) and e() tgent and normal the ellipse edge correspondingly.Below we do not mention -dependance unless absolutely necessary.Then we can approximate the crack values of the potentials () by suming that at each point x' (, u) the phe of each plane wave in the superposition may be decomposed into incident and scattered, ke( -x'(e)-(u, where x'() is projection of x' onto he crack edge.Then on the irradiated crack side we have 0, (' 1 . e, + ke ike('"'x'() f A[(k)e-i(d Z (,.() Note that in this approximation the components 0t +/Oz'(o, n) and -(a, n) are given when applicable in the local base (er(a), e(a), ez).The non-zero elements of the transition matrix T from this base to {ex, ey, ez } are given Tll T22 er(a) e, T12 -T21 er(a) -e,, T33 (37) By substituting (35) into (34), we obtain the following triple integral 0 0 where he matrix T, is extended to 4x4 matrix, with the only new non- ero element Too 1; ()= (asine) + (bcose); and [1-u/0()] is the Jacobian of the transformation from (,

Second Tier
The critical points in (38) are an amplitude pole, PSP and boundary critical points of the second kindl.Let us consider their contributions in turn.
Co'ibulm o] th isohted bou points: GTD edcje diructed m The PSPs of the double integral in (u, (:) in (38) are where d' V/(X x'(a)) 2 / (y y'(a)) 2 + z2.Let a critical point like this be isolated.Then applying stationary phase method (SPM), we get the Edge Integral 2r Ia(GTD)(x) --e 'k'd'(-) do', where k,d (t) k P(i'c) x'(a) + kad; and the diffracted directivity pat- The boundary critical points, a,, in (38) are PSPs of the Edge Integral which may be found numerically.Their number M may be 2 to 4. Let all a, be isolated.Then applying SPM to the Edge Integral, we have () where, a(TD) are each described by a formula similar to (17), with all () The above formula is in agreement with quantities evaluated at a a m.
GTD. Inside each a-geometrical region, the a-component of the scattered M ,(GTD) field is u (c) u () + m:-m A formula similar to (25) can be re-derived in the framework of the two-tier appraoch, with the result that k (GTD) is replaced by k? (GTD), pa(GTD) by p?(GTD) etc.The coalescence of the specular point x a(rel) with one of the flash points X? (GTD) corresponds to observation point x crossing a geometrical shadow boundary.Inside each c-penumbra, the c-component M (GTD) of the scattered field is u a(sca') u (Pe") + rn-2 Let a and a in Edge Integral coalesce.Applying USPM9', we have ()-{CAi[-k/3]+cAi'[-k/3]}ed (42)   u, [d(O(a)-d(t)(a)]} 2/a.The coalescence of a and a corresponds to the observation point x crossing a surface, such that inside, each x is reached by four edge scattered rays, and outside by two.This type of sur- 3/2 face is known a smooth caustic.Note that kava is proportional to the phe difference between two coalescing edge scattered rays.As usual, the Airy function Ai of-a a describes the smooth caustic contribution.In- side a smooth -caustic boundary layer, the a-component of the scattered field is u a(scae) u a(v) u a(caus)

+ m=3
The boundary layers coalesce may surrounding non-smooth parts of the caustics where three be described via the Pearcey integral or using numerical integration.

Numerical results
The code based on the time-harmonic asymptotic formulae has been fully tested against a numerical code based on the boundary integral equation method (Fig. 3).There is a good agreement for the range of parameters for which the applicability regions of both approaches overlap.
The time-harmonic asymptotics may be used to model the propagation of pulses by means of harmonic synthesis.Let us consider the following model parameters: the wave speeds cp 5840 m/s, cs 3170 m/s, the solid density 7770 kg/m3, and the pressure amplitude P0 1 MPa.For simplicity of presentation, let the pressure input be a narrow band pulse, one cycle of sin(27rft), f 5 MHz.Typical pulse trains are presented in Fig. 4. In Fig. 4a the observation point with lies inside both P and S wave caustic surfaces.In Fig. 4b the observation point lies outside both P and S caustics.

Conclusions
High-frequency asymptotic models for simulating ultrasonic pulse propa- gation and scattering in the transducer near have been developed.The models elucidate the physics of the problem and give explicit dependence on model parameters, thus allowing an easy prediction of pulse amplitudes and shapes.The corresponding asymptotic codes are at least hundred times faster than direct numerical codes, and practically just as accurate.
() ()]. ()with the polarisation vectors pP n(a) and pS n(a).McMaken s derived an analogous edge integral for normal incidence of plane P-wave on a penny-shaped crack.