On efficient quantitative analysis in real-time ultrasonic detection of cracks

Ultrasonics 37 (1999) 239–245

On efficient quantitative analysis in real-time ultrasonic detection of cracks Antonio Scalia a, *, Mezhlum A. Sumbatyan b a Dipartimento di Matematica, Universita` di Catania, Viale A. Doria n. 6, Catania 95125, Italy b Research Institute of Mechanics and Applied Mathematics, Stachki Prospect 200/1, Rostov-on-Don 344090, Russia Received 16 January 1998; received in revised form 1 October 1998

Abstract This paper deals with an efficient algorithm applied to a real-time detection of cracks by ultrasonic techniques. Being founded upon dynamic equations of linear isotropic elasticity, it operates with a uniform approximate representation for the Green tensor. The proposed method demonstrates high efficiency for arbitrary shape of the transducer and the crack in plane, as well as for arbitrary distribution of the stress over the probe basis. Some asymptotic estimates are proposed, to test accuracy of the classical ‘engineering’ far-field approximations. The algorithm permits calibration by a bottom echo-signal, as well as by an echo-amplitude reflected from a side-drilled hole. © 1999 Elsevier Science B.V. All rights reserved. Keywords: Crack; Fourier transform; Incidence; Scattering; Transducer

1. Introduction Quantitative ultrasonic methods are very important in detection and evaluation of flaws in solids. A number of analytical and direct numerical techniques have been developed for this aim that are well documented in literature [1–7]. Application of the Fourier transform seems to be very powerful regarding wave fields generated by ultrasonic transducers into elastic half-space, so far as this permits explicit analytical representations for the incident and the scattered wave fields. However, the analytical approach leads to a fourfold or even sixfold multiple integration of highly oscillating functions that requires too long a time when implemented on any modern computer. Thus, one is faced with a very exclusive situation when an explicit analytical result is absolutely useless for real practice. Some possible ways to overcome this obstacle are proposed in a recent study [8]. In the present paper we develop an efficient approach, combining both analytical and numerical treatment, to calculate the echo-response of plane cracks in real time on a PC. Starting from a high-frequency representation for the Green tensor of the point force applied to the * Corresponding author. Fax: +39-095-330094; e-mail: [email protected].

boundary, the problem is then reduced to a multiple integral of slowly varying functions over finite plane areas occupied by the crack and the transducer. For the problem of a round normal piston-type transducer scanning a penny-shaped crack some classical far-field approximations are tested, being compared with results of the computations obtained by the proposed method. It is shown that the proposed method is also quite appropriate for calibration by the echo-impulse reflected from a plane bottom, as well as from a side-drilled, hole.

2. Normal P-probe scanning a horizontal crack The algorithm proposed below is applicable for various problems. To demonstrate its potential, we consider here detection of a plane crack S by a normal transducer c of arbitrary contact area S in the echo-impulse method 0 (see Fig. 1). We assume a harmonic regime with respect to time (the product exp(−ivt) is omitted everywhere below) and accept a linear isotropic elastic model of the medium. Since the longitudinal wave speed is nearly two times greater than the transverse one, in order to simulate impulse propagations we keep a longitudinal component of the wave field only. Thus, the contribution of any normal point force s (x , y ) applied over S to the 0 0 0 0 incident wave incoming to the point (x , y ) µS is c c c

0041-624X/99/$ – see front matter © 1999 Elsevier Science B.V. All rights reserved. PII: S0 0 4 1 -6 2 4 X ( 9 8 ) 0 0 06 0 - 2

240

A. Scalia, M.A. Sumbatyan / Ultrasonics 37 (1999) 239–245

the integral in Eq. (1a):

A B

r eik1R1 ihk 0 1 ∑ ~− 1 s (x , y ) f , 2p 0 0 0 1 R R2 z 1 1 r2 =(x −x )2+( y −y )2, 1 c 0 c 0 (2a2−d2)2 , f (a)= 1 (2a2−d2)2+4a2 앀1−a2 앀d2−a2

(3a)

(3b)

k R2 =r2 +h2, d= 2 >1, 1 1 k 1 as well as for the tangential stress components Fig. 1. Normal ultrasonic transducer above the plane crack parallel to the boundary surface.

given, with the use of the Fourier transform, by [8]

PP

2 (2a2−k2 )2 e−c1h s 0 2 ∑= 0 4p2 (2a2−k2 )2−4a2c c z −2 2 1 2 ×e−[ia1(xc−x0)+ia2(yc−y0)] da da , 1 2 a2=a2 +a2 , c =(a2−k )1/2, j=1,2, 1 2 j j

(1a) (1b)

where k =k and k =k are, respectively, the longitudi1 l 2 t nal and the transverse wave numbers. The total incident wave is then determined as a superposition of the elementary ones:

PP

0 ∑ (x , y , x , y ) dx dy (2) c c 0 0 0 0 S0 z and this is written out only for one component of the incident stress tensor. All other components have a similar form and so are given by some fourfold integrals of highly oscillating functions. It is well known in computational analysis that significant complications arise in numerical calculations of multiple integrals of highly oscillating functions. Analogous difficulties also arise regarding calculations of the scattered wave. So this ‘standard approach’ cannot provide computations in real time – even on modern computers. In recent studies we have developed a different idea. It is well known that in scalar acoustics the incident wave, generated by a point-source applied to the boundary, is trivially given by the exponential function exp(ikR)/R, R=(x2+y2+z2)1/2, unlike the complex integral ( Eq. (1a)). It is tempting to try to propose some elementary representation for the integral (Eq. (1a)), like in the scalar case. It is shown in Ref. [9] that an accurate approximation can be obtained for the integral in Eq. (1a) by using a short-wave asymptotic technique. Actually, the two-dimensional stationary phase method leads to the following high-frequency representation for sinc (x , y )= z c c

ihk eik R T0 ~ 1 s (x , y ) 1 1 f , xz 2p 0 0 0 R2 2 1 ihk eik1R1 T0 ~ 1 s (x , y ) f , yz 2p 0 0 0 R2 3 1 x −x r 0 g 1 , f =2 c 1 2 R R 1 1 y −y r 0 g 1 , f =2 c 3 1 R R 1 1 앀1−a2(2a2−d2) , g (a)= 1 (2a2−d2)2+4a2 앀1−a2 앀d2−a2

A B A B

(4a)

(4b)

(4c)

so the latter involves the total tangential stresses incident to the crack faces, together with the normal stress ( Eq. (2)): tinc (x , y )= xz c c tinc (x , y )= yz c c

P P

S0

T0 (x , y , x , y ) dx dy , xz c c 0 0 0 0

(5a)

T0 (x , y , x , y ) dx dy . (5b) yz c c 0 0 0 0 S0 In practice, the accuracy of Eqs. (3a)–(4c) is uniformly good, out of a small vicinity of the boundary plane xy. It should also be noted that the representation obtained from Eqs. (2)–(5b) for the incident wave field is valid for an arbitrary distribution of the contact pressure s (x , y ), (x , y )µS , as well as for arbitrary 0 0 0 0 0 0 form of the probe basis S in plane. 0 The problem of diffraction of the incident wave by flaws of complex shape can be studied in many different ways, for instance by boundary element methods, Tmatrix methods, as well as by some approximate techniques. A solution for the round crack can be constructed explicitly as an infinite series in Mathieu functions, which seems to become slowly convergent with increasing frequency. We start from the observation that the Kirchhoff ideas, well known in scalar acoustics, are quite appropriate also in elastic dynamics. Namely, we treat this so that the boundary wave field is locally the same as in diffraction by an infinite plane, locally

241

A. Scalia, M.A. Sumbatyan / Ultrasonics 37 (1999) 239–245

tangential to the surface of the flaw. For cracks of arbitrary form in-plane it is a good approximation when its size along any direction is greater than the wave length (that is typically round 2–3 mm). Within this framework, such an approach can give the scattered wave field by using the Fourier transform and the two-dimensional stationary-phase method, as follows [8,9]: ihk ssc (x, y)= 1 z 2p

PP

Sc

[sinc (x , y ) f +tinc (x , y ) f z c c 4 xz c c 5

eik R +tinc (x , y ) f ] 1 2 dx dy , c c yz c c 6 R2 2 x−x r r c g 2 , f =2 2 , f =f 4 1 R 5 2 R R 2 2 2 r y−y 2 , R2 =r2 +h2, c g f =2 2 R 2 2 6 R 2 2 앀d2−a2(2a2−d2) g (a)= , 2 (2a2−d2)2+4a2 앀1−a2 앀d2−a2

A B

(6a)

A B

A B

(6b)

(6c)

K K P∞

S S = 0 c cos q |V (q)|W2(ak sin q)W(2bk sin q) ll 1 1 P l2 R2 0 1 (9)

where a is the radius of the transducer, b is the radius of the crack, S =pa2 and S =pb2 are their respective 0 c squares, R is the distance between their centres, l =2p/k is the longitudinal wave length, cos q=h/R, 1 1 V (q) is the reflection coefficient of the longitudinal ll wave, and the radiation pattern M(x) of the round radiator is known to be M(x)=2J (x)/x, where J (x) is 1 1 the Bessel function. This approximation is worthy of a detailed verification. The far-field limit can be directly extracted from the above representations, if we set there r /R =r /R ¬ 1 1 2 2 sin q, R =R ¬R, (x −x )/R ¬sin q, ( y −y )/R ¬0, 1 2 c 0 1 c 0 1 (x−x )/R ¬−sin q), ( y−y )/R ¬0 everywhere, c 2 c 2 excepting exponential functions. Then one comes to the identity

K K PP P∞

P

r2 =(x−x )2+( y−y )2. 2 c c The integral response is then

PP

transducer is given by [10]

=

0

h2k2 1 [ f 2 (sin q)+4 sin q2 g (sin q) g (sin q)] 1 1 2 4p2R4S 0

×

S0

dx dy 0 0

PP

Sc

dx dy c c

PP

S0

eik1(R1+R2) dx dy

(7)

(10)

S0 which should be estimated relative to the total pressure radiated into the medium (in decibels)

because for a piston-like probe one may put or s (x , y )¬1, and so P =S . Operating with a standard 0 0 0 0 0 approximation for the exponential function in Eq. (10) like in classical acoustics, one then determines that [11]

P∞=

ssc (x, y) dx dy. z

A=20 log(|P∞/P |), 0

P = 0

PP

s (x y ) dx dy . 0 0 0 0 0

(8)

S0 Thus, the proposed algorithm requires a sixfold integration of slowly varying functions over finite domains that are given by Eqs. (2)–(7). The approach, described above as a standard, would require an eightfold integration of highly oscillating functions over infinite regions. Even, when being reduced to a fourfold integral for some particular distribution of the applied stress ( like for a piston-type transducer) and for some canonical forms of the probe basis and the crack (round, elliptic, rectangular, etc.), it still requires multiple integration of strongly oscillating functions in infinite domains.

3. Testing ‘engineering’ models It is known that an ‘engineering’ far-field approximation for the echo-amplitude reflected from a pennyshaped crack, parallel to a non-coaxial piston-like round

R #(x2 +y2 +h2)1/2−(x /R)x 1 c c c 0 =(x2 +y2 +h2)1/2−x sin q, (11a) c c 0 R #(x2 +y2 +h2)1/2−(x /R)x 2 c c c =(x2 +y2 +h2)1/2−x sin q. (11b) c c Taking into account the values arising from the integrals

PP

S0

e−ik1x0 sin q dx dy = 0 0

PP

e−ik1x sin q dx dy

S0 =S W(ak sin q) (12) 0 1 the remaining double integral over the crack surface S c in Eq. (10) can be evaluated by using a change of variables x =R sin q+m cos y, y =m cos y (h=R cos q, c c y /R#0) as follows: c b exp(2ik 앀x2 +y2 +h2) dx dy ~e2ik1R m dm 1 c c c c Sc 0 2p e2imk1 sin q cos y dy=e2ik1R S W(2bk sin q). (13) × c 1 0

PP

P

P

242

A. Scalia, M.A. Sumbatyan / Ultrasonics 37 (1999) 239–245

Note at last that the expression in the square brackets in Eq. (10) is equal to f (sin q), so finally we obtain at 1 R2,

K K

S S = 0 c cos2 q f (sin q)W2(ak sin q)W 1 1 P l2 R2 0 1 ×(2bk sin q). (14) 1 Eq. (14) is different from the one given by Eq. (9), which is based upon some intuitive physical ideas only. It can be shown that, with our notation, the reflection coefficient V in Eq. (9) is ll V (q)=f (sin q), (15) ll ll (2a2−d2)2−4a2 앀1−a2 앀d2−a2 . f (a)=− ll (2a2−d2)2+4a2 앀1−a2 앀d2−a2 P∞

therefore cos q|V (q)|≠cos2 q f (sin q). The reflection ll 1 coefficient would appear in Eq. (14) if a different receiver is placed symmetrically to the right of the main transducer, with the same angle q (see Fig. 1), which is quite natural from the physical point of view. In that case the expression in square brackets in Eq. (10) is equal to f 2 (sin q)−4 sin q2 g (sin q)g (sin q)=−f (sin q)V (q) 1 1 2 1 ll (16) instead of f (sin q). However, even then Eq. (14) differs 1 from Eq. (9).

4. Calibration by a bottom signal and by a side-drilled hole It is proved in recent work [8], for the echo-reflection of the wave beam from the infinite plane, that the wellknown scalar ‘mirror’ ideas about a secondary image of the transducer, symmetric with respect to the bottom, can be directly transferred to the solution of the problem about the normal transducer placed above the elastic layer of a constant thickness. This allows one, recalling Eqs. (2)–(3b), to write out the bottom echo-signal amplitude as a fourfold integral of a slowly varying function P∞=−

PP PP A B

ihk

1 2p

×

S0

f

S0

s (x , y ) dx dy 0 0 0 0 0

r eik1R dx dy, 1 R R2

(17a)

r2=(x−x )2+( y−y )2, R2=r2+(2h)2, (17b) 0 0 with a double integration over the basis of the probe. For a side-drilled hole (see Fig. 2) the calculation of

Fig. 2. Calibration by a side-drilled hole.

the echo-signal amplitude could be reduced to the same sixfold integration, like for the plane crack. However, an infinite interval over the length of the infinite cylinder would arise in this problem (when taking the integral along the variable y ) that, with the oscillating function c exp[ik (R +R )], requires too a long time for computa1 1 2 tion. Fortunately, integration on y µ(−2, 2) permits c explicit evaluation of this integral, so the dimension of the multiple integral can be reduced by a unit, and this is transformed to a fivefold integration of a slowly varying function again. The transformation procedure is not trivial, and should be explained in detail. Let us consider any point (x , y , z )µS that belongs c c c c to an ‘illuminated’ part S of the cylinder boundary c+ surface, with respect to the current point force s (x , y ), (x , y )µS acting over the probe basis. Let 0 0 0 0 0 0 the plane x∞ y∞ be tangential to the surface. Then the c c incident stresses incoming to the plane at the point (x , y , z ) are given as [11] c c c sinc =sinc cos2 Q−tinc sin 2Q+sinc sin2 Q, z∞c z xz x tinc =tinc cos 2Q+1 (sinc −sinc ) sin 2Q, x∞c z∞c xz x 2 z tinc =tinc cos Q−tinc sin Q, y∞c z∞c yz xy z =h−b cos Q, x =b sin Q, c c

(18a) (18b) (18c) (18d )

where h is the distance between the centre of the cylindrical hole and the boundary plane x, y; b is the radius of the cylinder; Q is the angle between the straight lines given by the axes z and z∞. c By using the Fourier transform (marked by an overbar) all components of the stress tensor in Eqs. (18a)– (18d ) can be expressed in explicit form, like in Eqs. (1a)

243

A. Scalia, M.A. Sumbatyan / Ultrasonics 37 (1999) 239–245

and (1b) (see also Ref. [8]): s: inc =s: 0 z

(2a2−k2 )2 2 e−c z , 1c D(a)

the probe basis (19a)

PP P P

ssc = n

D(a)=(2a2−k2 )2−4a2c c , 2 1 2 s: inc =s: 0 x

(2k2 −k2 −2a2 )(2a2−k2 )2 1 2 1 2 e−c z , 1c D(a)

2ia c (2a2−k2 )2 1 1 2 e−c z , t: inc =s: 0 1c xz D(a) t: inc =s: 0 yz

× (19b)

2a a (2a2−k2 )2 1 1 2 e−c z , 1c D(a)

(19d )

written for their two-dimensional Fourier images in the variables (x, y) of the main fixed coordinate system (x, y, z). By analogy, the normal stress component ssc (x, y), n (x, y)µS scattered from the cylinder, in the Kirchhoff 0 approximation, can be determined as a composition of the local contributions (l is the illuminated part of the + generatrix)

PP

(ssc cos2 Q−tsc sin 2Q z∞c x∞c z∞c Sc +ssc sin2 Q) dy dl . (20) x∞c c + Here the quantities at the right-hand side are defined in a sliding coordinate system (x∞ , y∞ , z∞ ) coupled with the c c c sliding tangential plane (see Fig. 2) and have the following form, for the respective Fourier transforms in the variable b=( b , b ): 1 2

ssc (x, y)= n

(2b2−k2 ) e−c1(b)z∞c 2 {(2b2−k2 )sinc (x , y , z ) s: sc =− 2 z∞c c c c z∞c D(b) +2ic [ b tinc (x , y z )+b tinc (x , y , z )]}, 2 y∞c z∞c c c c 2 1 x∞c z∞c c c c t: sc =− x∞c z∞c

(21a)

2ib c e−c1(b)z∞c 1 1 {(2b2−k2 )sinc (x , y , z ) 2 z∞c c c c D(b)

+2ic [ b tinc (x , y z )+b tinc (x , y , z )]}, 2 y∞c z∞c c c c 2 1 x∞c z∞c c c c

(21b)

2ib c e−c1(b)z∞c t: sc =− 2 1 {(2b2−k2 )sinc (x , y , z ) y∞c z∞c 2 z∞c c c c D(b) +2ic [ b tinc (x , y z )+b tinc (x , y , z )]}, 2 y∞z 2 1 x∞c z∞c c c c c ∞c c c c

(21c)

where the incident stresses are given by Eqs. (18a)– (19d ). One then comes to the representation for the normal component of the scattered stress incoming to

S0

−2 2

×

db 2

P

2

−2

P

2

−2

dy c

P P dQ

l+

2

−2

db 1

da 1

A(a , a , b , b , sin Q, cos Q)da 1 2 1 2 2

−2 ×exp −{c (a)z +c (b)z∞+i[a (x −x ) 1 c 1 1 c 0 +a (y −y )+b (x∞−x )+b (y∞−y )]} 2 c 0 1 c 2 c

(19c)

2ia c (2a2−k2 )2 2 1 2 e−c z , 1c D(a)

t: inc =−s: 0 xy

2

s (x , y ) dx dy 0 0 0 0 0

(22a)

where A is a certain algebraic function of its arguments which is not cited here explicitly, so as to reduce the volume of the writing. Note that the values (x∞, y∞, z∞) in Eq. (22a) denote the Cartesian coordinates of the point (x, y)µS in the 0 sliding coordinate system (x∞ , y∞ , z∞ ): c c c x∞=x cos Q−h sin Q,

(22b)

y∞=y,

(22c)

z∞=h cos Q+x sin Q−b.

(22d )

Now integration on y along the directrix of the c cylinder yields the Dirac delta-function

P

2 −2

eiyc(b2−a2) dy =2pd( b −a ) c 2 2

(23)

so integration over a can be omitted, owing to the main 2 property of the delta-function, setting a =b in the 2 2 remaining integrals. Thus one obtains the three-dimensional integral as a kernel in Eqs. (22a)–(22d ),

PPP

2

A(a , b , b , b , sin Q, cos Q) 1 2 1 2 −2 ×eiU(a1,b1,b2) da db db 1 1 2

(24a)

with the phase function U=앀k2 −a2 −b2 z +앀k2 −b2 −b2 z∞−a (x −x ) 1 1 2 c 1 1 2 1 c 0 −b (x∞−x )−b ( y∞−y ) (24b) 1 c 2 0 that should be evaluated by a three-dimensional stationary-phase method, to complete the analogy with the method proposed above for diffraction by plane cracks. It can be seen that the stationary point of the phase, Eq. (24b), defined from the system (see Ref. [12]) ∂U ∂a 1

=0,

∂U ∂b 1

=0,

∂U ∂b 2

=0,

(25)

244

A. Scalia, M.A. Sumbatyan / Ultrasonics 37 (1999) 239–245

is given by a = 1*

(x −x )앀k2 −b2 1 2* , 0 c q

b = 1*

(x −x∞)앀k2 −b2 1 2* , c r

(26a)

y −y∞ 0 b =k , 2* 1 앀(q+r)2+( y −y∞)2 0 q=앀(x −x )2+z2 , r=앀(x −x∞)2+z∞2, 0 c c c

(26b)

with the value of the phase in this stationary point being equal to U =ik 앀(q+r)2+( y −y∞)2. 1 1 0

(27)

We cite also an explicit form for the determinant of the matrix of the second-order derivatives of the function U: det(U◊)=− 1

q3r3[(q+r)2+( y −y∞)2]5/2 0 . k3 z2 z∞2(q+r)4 1 c c

(28)

Finally, the desired approximation for the echo-impulse amplitude, as it follows from Eq. (22a), is P∞=(2p)5/2 ×

P

l+

PP

dx dy S0

PP

S0

s (x , y ) dx dy 0 0 0 0 0 ik U

e 1 1 dQ A(a , b , b , b , sin Q, cos Q) 1* 2* 1* 2* |det(U◊)|1/2 1 (29)

given as a fivefold integral of a slowly varying function. Note that ‘illuminated’ part of the generatrix contains only the points for which the two straight intercepts (x , y )−(x , y , z ) and (x , y , z )−(x, y) do not cross 0 0 c c c c c c the boundary of the cylinder. Many concrete practical problems were solved using the proposed algorithm for various ultrasonic probes, as well as for flaws of various geometries. Fig. 3 demonstrates an example of the scanning of a penny-shaped crack by a round piston-like normal probe, calibrated with respect to a side-drilled hole of the same radius (the hole is in a fixed position just under the transducer). To make computations for any single point in the Fig. 3, it takes nearly 0.12 s when implemented on a PentiumA-120. Note that the echo-signal from a coaxial plane crack is always stronger than from the cylinder of the same diameter placed at the same distance. With increasing distance the echo-amplitude decreases, and the diagram becomes more uniform. This is quite natural from the physical point of view, owing to the divergence of the ultrasonic wave beam.

Fig. 3. Echo-impulse amplitude of a round piston-like transducer (a= 5 mm) scattered from a penny-shaped crack (b=2.5 mm) versus the distance between the probe and the crack centres (in the horizontal direction): c =5.94 km/s, c =3.23 km/s, f= 4 MHz; (1) h=40 mm, 1 2 (2) h=60 mm, (3) h=80 mm. Calibration by a side-drilled hole with b=2.5 mm, h=40 mm.

5. Conclusions (1) The method proposed in the present work permits real-time quantitative analysis in detection of cracks by ultrasonic impulses. The error of the method is less than 1 dB, when the diameters and the distance are greater than a few millimetres. The result obtained can be represented as a multiple integral of a slowly varying function, which takes explicitly a few lines in a written form. (2) ‘Engineering’ models, being based upon intuitive physical ideas, cannot provide correct approximations, and the most reliable way is to extract them as asymptotes from the exact solution. For the problem at hand, no far-field asymptotes (even the one constructed correctly above) give an acceptable approximation, in that the radiation pattern yields the disappearance of the echo-amplitude for a number of values of the angle 0
A. Scalia, M.A. Sumbatyan / Ultrasonics 37 (1999) 239–245

References [1] J. Krautkra¨mer, H. Krautkra¨mer, Ultrasonic Testing of Materials, Springer, Berlin, 1983. [2] J.D. Achenbach, J. Sound Vibr. 159 (1992) 385–401. [3] L.J. Bond, Methods for computer modeling of ultrasonic waves in solids, in: R.S. Sharpe (Ed.), Research Techniques in NDT, Vol. 6 Academic Press, New York, 1982, pp. 107–150. [4] L.W. Schmerr, A. Sedov, J. Acoust. Soc. Am. 86 (1989) 1988–2006. [5] R.K. Chapman, J. Nondestr. Eval. 9 (1990) 197–211.

245

[6 ] K. Harumi, M. Uchida, J. Nondestr. Eval. 9 (1990) 81–99. [7] A. Bostro¨m, H. Wirdelius, J. Acoust. Soc. Am. 97 (1995) 2836–2848. [8] M.A. Sumbatyan, N.V. Boyev, Ultrasonics 32 (1994) 5–11. [9] M.A. Sumbatyan, Akust. Zh. (Sov. Phys. Acoust.) 34 (1) (1988) [10] I.N. Ermolov, Theory and Practice of Ultrasonic Control, Mashinostrojenije, Moscow, 1981 (in Russian). [11] J.D. Achenbach, Wave Propagation in Elastic Solids, North-Holland, Amsterdam, 1973. [12] M.V. Fedorjuk, Zh. Vychisl. Mat. Mat. Fiz. (J. Comp. Math. Math. Phys.) 2 (1) (1962).