Transient axially asymmetric response of an elastic plate

WAVE MOTION 11 (1989) 271-295 NORTH-HOLLAND

TRANSIENT

AXIALLY ASYMMETRIC

271

RESPONSE

OF AN ELASTIC PLATE

Fadil SANTOSA Department of Mathematical Sciences, Unioersity of Delaware, Newark, DE 19716, U.S.A.

Yih-Hsing PAO Department of Theoretical and Applied Mechanics, Cornell University, Ithaca, N Y 14853, U.S.A.

Received 10 March 1988, Revised 26 July 1988 In Memory of Raymond Mindlin

The problem of finding the response of a thick infinite plate to the application of body forces is considered. Transform techniques are employedto reduce the problem to that of an initial boundary value problem for a one-dimensional hyperbolic system. The latter is solved by the method of eigenfunction expansion. A representation of the displacement field is obtained which may be used to compute the transient response of the plate. As an example we apply this method to obtain the plate response to the application of a point force. We discuss the computerimplementation of this approach. Numerical simulations are presented for the case of an oblique force with a short burst time function.

1. Introduction The purpose of this work is to arrive at a relatively compact representation of transient waves in an elastic plate. The study of wave motions in an infinite linearly elastic extended plate goes back to Rayleigh and Lamb who treated mainly the propagation of steady-state waves. Much of what was known about this subject and extensive analysis of the dispersion relation for Rayleigh-Lamb waves were contained in a monograph by Mindlin [1]. The problem of transient waves excited by an impulsive force was treated in recent years as summarized by Achenbach [2] and Miklowitz [3]. Roughly speaking, there have been three types of approaches to studying this problem: the transform methods, the method of generalized rays, and the method of normal modes. It is not too surprising that all these methods are intimately connected, even though they are different in spirit. For transform methods, the reader is referred to the expositions in Achenbach and in Miklowitz. The solution to the elastodynamic problem is expressed in terms of Laplace or Fourier transform integrals. These integrals, although very complicated, can be evaluated by an asymptotic analysis to obtain approximate solutions, or by applying the calculus of residues. The latter has only been successfully applied recently by Vasuderan and Mal [4] for the case of axisymmetric forces. Ceranoglu and Pao [5] showed that the integral transforms by which the solution of the plate problem is represented has an interpretation in terms of what are called "generalized rays". Indeed, they showed that in problems where the response of the plate near a point source is desired, a highly accurate short-time response can be obtained. However, when the receiver is at a large distance, a very large number of generalized ray integrals are to be evaluated and this method becomes impractical. For the problem of far-field and long-time response, Weaver and Pao [6] applied the method of normal mode or eigenfunction expansion. They considered displacement fields which are axisymmetric, and were 0165-2125/89/$3.50 © 1989, Elsevier Science Publishers B.V. (North-Holland)

272

F. Santosa, Y.H. Pao / Transient plate response

able to arrive at an eigenfunction expansion formula for the solution of the displacement field in a plate. They showed that the resulting expansion is particularly suited for obtaining transient responses for large source-receiver separation, and large times. In arriving at the representation for the response of the infinite plate, Weaver and Pao used the eigenfunctions for a finite circular plate. The results for the infinite plate are then derived by letting the radius of the plate become infinite. This is done formally; the mathematical justification, especially p r o o f of convergence, is very difficult. In this paper, we modified the approach by Weaver and Pao to derive the general solutions for transient waves in an infinite plate due to the excitation of axially asymmetric point forces. The method proceeds with a two-dimensional Fourier transform of the displacements on the plane of the plate, and the resulting transformed displacement problem is solved by resorting to an eigenfunction expansion. Finally, to return to physical displacements, we take inverse Fourier transforms. The general solution thus derived includes the axisymmetric loading as a special case. In fact, the final formulas for the solution are identical for the special case of axisymmetry. What we have avoided is the use of the eigenfunctions for finite plates. This is desirable because the delicate issue of boundary conditions at the edges of the finite plate never arises, especially when the desired displacement fields do not possess symmetries. To obtain the normal modes of the plate, we need to prescribe some boundary conditions at the edge of the plate. These boundary conditions must be such that the desired displacement lies in the span of the eigenmodes. That is, if the eigenmodes are axisymmetric, we can only hope to construct axisymmetric displacement fields. Furthermore, the expansion of the displacement field in an infinite plate has been mathematically justified in the paper by Santosa [7]. In particular, it was shown that the transformed problem has a solution, and that the inverse Fourier transforms are convergent integrals, under some smoothness hypotheses on the excitation. In practice, we will only take a finite number of terms in the eigenfunction expansion, and the inverse Fourer transforms will be truncated. Error estimates for this type of approximations are available in the reference cited. We will show here that, for the case of a point source with smooth time signature, only a finite number of terms are needed to give a good approximation. Moreover, the inverse Fourier transforms mentioned could be replaced by integration over finite domains. Recently, Lu and Felson [8] proposed a hybrid method which combines the method of generalized ray and the method of normal modes to analyze transient waves over a long duration or a large and short distance. So far, they have treated the plane strain problems with a line source in a plate or elastic layer. Extension of the hybrid method to waves without plane or axial symmetry, as treated in this paper, should be of great interest. The p a p e r is organized in six sections. In Section 2, we will give a detailed description of the problem and arrive at a formal representation for the solution. Section 3 is devoted to the calculation of the eigenvalues and eigenfunctions of the Sturm-Liouville system associated with the transformed problem. A specific problem, that of finding the response of a plate to an applied point source, will be solved, and the results presented in Section 4. This is followed in Section 5 by a discussion on the numerical aspect of finding the transient response of the plate. In particular, we will arrive at a rule for specifying the discretization of the problem depending on the input excitation, and the length of the response trace desired. We end in Section 6 with our numerical results and concluding remarks.

2. Problem description and representation of solution We consider the propgation of waves in a plate of finite thickness. The coordinate system is arranged

F. Santosa, Y.H. Pao / Transientplate response

273

SO that the plate occupies the region

~ = {(x, y, z) : oo > x > -oo, a 3 > y > - o o ,

Izl<_hL

(2.1)

Thus the plate is of thickness 2h. Let u(x, y, z, t) stand for the particle displacement at (x, y, z) at time t. Then u satisfies the equation of linear isotropic elasticity 02u p - ~ = (A +/z)V(V. u) +/~V2, + f (2.2) where p is the density, A and /~ are the Lain6 moduli of the elastic plate. In (2.2), f = f ( x , y, z, t) is the applied body force, and we assume that f = 0 for t ~<0. Equation (2.2) is satisfied by u for t > 0 and (x, y, z) e D. For t = 0, we have initial conditions which represent a quiescent initial state, given by 0, ot (x,y,z,O)=O

u(x,y,z,O)=O.

and

(2.3)

On the boundary of the plate, we prescribe traction free boundary conditions tr(u), ez =0,

t~>0, Izl= h, oo> x , y > - o o .

(2.4)

In (2.4) e~ is the unit vector in the z-direction. Equations (2.2), (2.3) and (2.4) form an initial boundary value problem (IBVP) which represents a mathematical description of the motion of the plate to the applications of body foces. The objective here is to find the solution u(x, y, z, t) for a prescribed f(x, y, z, t). To this end, we employ the solution method discussed by Santosa in [7]. Define ~ = ~(z, t; ~:, r/) to be the Fourier transform of the displacements u(x, y, z, t) through the formula

d y u ( x , y , z, t ) e -iex-iny.

~(z, t; ~:, 7/)=

(2.5)

The variables ~: and 7/ are real and play the role of wavenumber in our problem. We will adopt the notation that, by ('~), we mean the Fourier transform (2.5), of (.). We denote by (u~, u2, u3) the x, y and z components of u. The inverse Fourier transform of a gives rise to physical displacements,

u(x, y, z, t ) = ( - ~1) z2 f_°°ood~ I_°°oodr/ u(z, A t; sc, r/) e icx+in'. The expression above is a convergent integral, in the sense that it is square-integrable in D, if the forcing term f satisfies the following smoothness hypotheses oods¢

dr/

d z ( s¢ +r/ )If(z, t; g r/)12< oo

(2.6)

for each t < o o and each component f o f f . This fact is established in [7]. Let us suppose for now that (2.6) is satisfied by o u r f . From (2.2), we know that ~(z, t; s~, rt) satisfies

p-~=A~--~z2+A2~+A3a+J;,

(2.7)

t>O, ]zl
where the matrices Az, A2 and A3 are defined to be A~=

0 0

g 0

0 (A+2~,)

,

A 2 = - i ( A +/.Q

0

274

F. Santosa, Y.H. Pao / Transient plate response

I (a + 2/x)~:2+/z~72

(a +/z)~:r/

(a +/z)~r/

(h + 2 / z ) r / 2 + / z ~ 2

0

0

A 3 = --

0

.

° 1

~ ( ~ + n ~)

Initial data for ~ are ~-(z,t;~,~)=0

and

~(z,t;~:,~)=0.

(2.8)

From (2.4), we find the boundary conditions for ~ at [z[ = h, oa, /)"Z- .i~a3 = 0, -

o02" a : - i n a 3 = o,

-

(.~ +2tz) °02" a 3 - a (i~al + int]2) = 0.

(2.9)

The IBVP given by equations (2.7), (2.8) and (2.9) represent the transformed version of our original problem. The aim here is to find d for [z I<~ h and t > 0. Once this is done, we will take the inverse Fourier transform of a to obtain the physical displacements. Observe that while u is real, ~ will be, in general, complex. To arrive at a solution of the transformed problem, we will expand ~ in a linear combination of the eigenfunctions of the operator in (2.7). Thus let the pair (to(~:, ~7), ~(z; ¢, 7)) solve -- pto2~

=

d:t3 + A2 d-7+ dl3 _A3v~ in Izl< h A1 dz---5

(2.10)

and let t3 satisfy boundary conditions (2.9). We will show that there is an infinity of these eigenpairs, and that to depends on (~:2+ ~72). We label each pair by a superscript (n). That is, the set {to ("), ~,)1 v x,=o.1.2.... is the set of eigenpairs for (2.10) with boundary condition (2.9). It is not surprising that the to (") can be interpreted as natural frequencies of the plate, and the ~("), when multiplied by e x p ( - i ~ x - i r / y ) , can be interpreted as the normal modes of the plate. The relation between (~, ~7) and to(n) can be thought of as a dispersion relation. To find the solution a(z, t; ~:, r/), we write fi(z, t; ~, 7 ) = •

a(")(t; ~, n)f3(')(z; ~, 7).

(2.11)

n

The Fourier coefficients a (") are easily found from (2.7) and (2.8) to be a(")(t; ~:, ~7)- pto(,)(~:,

1

,7)11~(.)ll=(~,

7)

Io

d r sin[to(~)(~ :, n)(t - r)]. (], ~(")(~:, r/)),

(2.12)

where (f, 13(")) =~h dzf" 13(")*. We have invoked orthogonality of the eigenfunctions to arrive at the expression above. The physical displacements can then be found by taking the double inverse Fourier transform of (2.11) with (2.12) for the coefficients. In the next section, we will solve the eigenvalue problem in (2.10) and boundary conditions (2.9). Formulas for the eigenfrequencies to(") and eigenfunction ~3(") will be found. It should be noted that, in principle, the coefficients a (") in (2.12) can be calculated for most "reasonable" source functions f Hence, representations of the physical displacements can be written in terms of Fourier integrals. The difficulties in evaluating these integrals in closed form are two-fold: (1) the expressions are very complicated; (2) to(,)(~, ~7) can be found only implicitly. Thus, we are forced to evaluate the solution numerically. Error estimates for the approximation in which the series (2.11) and the Fourier integrals are truncated, can be

F. Santosa, Y.H. Pao / Transient plate response

275

found in [7]. Here we will show that when f is a point source and its time signature is smooth, we can arrive at a rule for finding the n u m b e r of terms in the series needed to give a g o o d a p p r o x i m a t i o n . Moreover, we can replace the Fourier integrals by finite integrals.

3. Eigenvalue problem Let us begin by finding the general solution of the system of differential equations (2.10), rewritten in the form d2~

d~

a~ ~z2 + a2-~z

2 ^ + ( A3 W pto I ) w = 0 .

(3.1)

The matrices A~, A 2 and A 3 have been previously defined. We will try solutions of the form ~ ( z ; 6, 7/) = r(~, 7/) exp(i~/z) where 3' is i n d e p e n d e n t o f z. Inserting this form into (3.1), we arrive at the requirements for r and 7 in pro - ),2/.t - ( k + 2/~)62-/x7/2

- ( k + ~)67/

- ( A +/-Q67/ 3,(3. + / * ) 6

ata 2 - y2/.e - ( X + 2/z)7/2-/z62 3'(3. + p,)7/

7(A + / z ) 6

]

y(X + / z ) 7/ Jr=0. pta 2 - 72(A + 2/*) - / , ( 6 2 + 7/2) (3.2)

For (3.2) to yield nontrivial solutions r, we need to have the determinant of the above matrix being zero. Direct calculation shows that the determinant is zero when

72 = ,f o,2/c -(62+ 7/2), to,2/c - (62 + 7/2). This result is not too surprising because they c o r r e s p o n d exactly to the two speeds o f p r o p a g a t i o n in an isotropic m e d i u m , CL = [(A +21z)/p] ~/2 and Cs = (ix~p) ~/2, which are the longitudinal and shear wave speeds respectively. Let k 2 = 6 2 + 7/2 and a = ( ( ¢ . 0 2 / C L ) - - k2) U2. Then, for 7 = a, the vector rx = (6, 7/, - a ) v satisfies (3.2). For 3' = - a , the corresponding vector is sl = (6, 7/, or) T. Let fl = ((to2/C 2) - k2) 1/2. Then, for 3, =/3, we have two vectors satisfying (3.2): r2 = (6, -7/, 0) v

and

r3 = (6, 7/, k2//3) T.

Similarly, we find the two vectors for 7 = - / 3 . They are s2 = (7/, - 6 , 0) T

and

s3 = (6, 7/, - k 2 / / 3 ) T.

The vectors 1"2 and s 2 have been selected to be p e r p e n d i c u l a r to the vectors r3 and $3 respectively. Having evaluated all the combinations of 7 and r which satisfy the equation, we can p r o c e e d by writing the general solution o f (3.1). We write it as ~ ( z ; ~:, 71) = p~ rt e iaz + qlSl e -iaz +p2r2 e iaz q- q2s2 e -i~z +p3r3 e i~z + q3s3 e -i~z

F. Santosa, Y.H. Pao / Transientplate response

276

where pl, P2, P3 and q~, q2, q3 are arbitrary complex constants. We prefer to work with sines and cosines. So, we rearrange the expression above to get =

i cos z] [ cosz

~ c o s f l z p~+ 71 cosl3z P2+ - i a sin flz -ik2/B sin/3z

+

O ~ P3 cos flz

[

~ sin/3z ,sinflz][~7} rt sin flz q~ + [ r1 sin flz q2 + [ (-ik2//3) cos flz ia cos/3z

( q3 sin flz.

(3.3)

O

The constants, although redefined, are still arbitrary. The splitting in (3.3) is quite natural and physically meaningful. Let us call a transformed displacement vector ~ symmetric (with respect to the z-axis) if ~(z) x e z = ~ ( - z ) X e z ,

~(z).e~=-~(-z).ez;

and antisymmetric if # ( z ) x e, = - ~ ( - z ) x e~,

~ ( z ) ' e z = ~ ( - z ) . ez.

Then any displacement vector can be separated into its symmetric and anti-symmetric parts. Observe that the first three terms in (3.3) correspond to a general symmetric solution of (3.1), and the last three to the antisymmetric solution. We will henceforth use ~3 to denote the symmetric solution and ~ to denote the antisymmetric solution. Thus, let t3 and ~ be cos az P l + rl cos/3z P2 + - i a sin az -ikg/fl sin/3z

if,=

Isinz] I 77 s i n a z q l + i~ cos az

]

rlsin~z q2+ -ikg/fl cos flz

P3 cos flz,

O•

0~: q3 sin flz.

(3.4a)

(3.4b)

So ~ and ~ represents the two types of solutions to our differential equation. We will force each of them to satisfy the prescribed traction-free boundary conditions (2.9). Although the boundary conditions for (2.9) are for z - - + h , we see that by this splitting we need only have that ~ and ~ satisfy (2.9) for z -- h. This is because the conditions at z = - h are satisfied when the conditions at z = h are met. These boundary conditions will determine the eigenvalues {to (")} and the eigenfunctions. First consider the symmetric solution 13. Substituting 13 into (2.9), for z = h, and collecting terms, we get d~----!1- i~v3 =

dz

dr----!-i~/~33 =-2~Tadz

--2~:~ sin ahp~ +

- fl ~ sin flhp2 - "Off sin ¢thp3 = 0,

sin ahpl+(k-~fl-fl)rt sin ~hp2+ ~ sin/3hp3 = 0.

d~3+ ( h + 2 / z ) dz A ( - i ~ v ' - i r l v 2 ) = i [ - ( A + 2 t z ) a 2 - A k 2 ] c ° s a h p ' + i [ ( h + 2 1 z ) k 2 - h k 2 ] c ° s f l J p 2 = O "

F. Santosa, Y.H. Pao / Transient plate response

277

We put the above into a matrix form as a system of equations for p~, P2 and P3:

] - 2 ~ a sin ~h _(fiE _ k 2) cos t~h

~--/3

(3.5)

~7 sin/3h

2k 2 cos/3h

0

Jl. PaJ

For (3.5) to yield nontrivial solutions, i.e., Pl, P2, and P3 not all zeroes, the determinant of the matrix must vanish. The determinant is det[. ] = (~2 + ~q2) sin flh[(fl 2 - k2) 2 cos t~h sin flh + 4c~flk 2 sin ah cos flh]. This is the equation by which the eigenvalues {to<")} are determined for each ~: and r/. The determinant is zero when sin/3h = 0

(SH-frequency equation),

(3.6a)

(f12_ k 2) cos t~h sin flh + 4clflk 2 sin ah cos flh = 0

(Rayleigh-Lamb frequency equation).

(3.6b)

From these equations, we seek, for each ~ and ~, a set of values {to~n)} for which they are satisfied. The first of the equations above involves only the shear wave speed. It will be shown that the corresponding displacements have no z-components. That is, these eigenfunctions correspond to horizontally (the xy-plane) polarized shear waves. The second is the well known Rayleigh-Lamb frequency equation. Similarly, we obtain the matrix equation satisfied by ql, q2 and q3 for the antisymmetric solution,

r2cosh - 2 r / a cos ah

[_(f12 _ k 2) sin ah

k2

coshllq1

2 k 2 sin flh

q2

~o

"= 0

(3.7)

q3

The determinant of the matrix is det[. ] = (s~2+ ~72) cos/3h[(fl 2 - k2) 2 sin c~h cos flh +4aflk 2 cos ah sin flh]. The two frequency equations are cos flh = 0

(SH-frequency equation),

(f12_ k 2) sin ah cos flh + 4otflk 2 cos oth sin flh = 0

(3.8a) (Rayleigh-Lamb frequency equation).

(3.8b)

These correspond to the antisymmetric displacements. Notice that because t~ and/3 depend on k, the resulting eigenvalues will be functions of k only. The variable k can be thought of as a radial wavenumber; ~ and ~ are the x and y wavenumbers. And since the eigenvalues {to <")} are frequencies, we can treat their dependence on k as dispersion relations. There will be four sets of dispersion relations. For every fixed k, we will find four sequences {to(")(k)} which satisfy (3.6a, b) and (3.8a, b). In addition, we will also have four sets of eigenfunctions corresponding to the four sequences. These dispersion relations have been studied extensively by Mindlin [1] (see also a review by Pao and Kaul, [9]). Here we record some of the results which are needed for the evaluation of the transient responses.

F. Santosa, Y . H . Pao / Transient plate response

278

Before doing that, we first normalize the variables w, c~,/3, ~: and ~7:

if-

2h oJ

2h

2h

g~=--o~,

11" C S '

2h

fi=--/3,

Ti"

2h ~=--*/.

~ = - - s ~, and

"IT

~

"iT

Naturally, it follows that /~= (2h/'rr)k. Let us also normalize the physical coordinates for convenience. We use the notation ff

¢rx 2h

~ry )7=

rrz i-2h

and

We can now label our eigenpairs. Let /~ be fixed; we will use two subscripts on &, fi, and o3. The first subscript is either S or A, depending on whether they correspond to the symmetric or antisymmetric displacements. The second subscript is either S or R depending on whether they satisfy the SH-frequency equation or the Rayleigh-Lamb frequency equation. We will continue to use 13 to denote symmetric displacements, and ~ for antisymmetric displacements. Hence, we will only have one subscript on the displacements--either S or R. Under this notation, we list all the variables involved in Table 1. Thus for each fixed /~, displacement and wave type, we can in principle find the set of eigenvalues a5("~, and the associated eigenfunctions. Table 1 List of dispersion relation relations, equations for frequencies, and eigenfunctions for the plate Dispersion relation symmetric

SH waves

-(n) ( k-) toss

displacements

RL waves

o5¢s'~)(/~)

Frequency equation

Eigenfunction

sin(½wfi(s~ ') = 0 ( t~(n)2

E'_2 ", 2

--

lI

-(tl)x

"

1

--(tT)

PSR -- ~ / COS~.~'n'aSR) Sln(~TrflSR) -1 -- .°tO~-(,,I,~O,)Z2 sin(½1r8~('~) COS(~r/3SR~ -i,,,) SR/J SR K

antisymmetric

SH waves

o3~ ( / ~ )

displacements

RL waves

a3~'R)(/C)

cos(~-d~)

-(n)

--(n) -2

+ 4Or ARflARk ~ t 2 = ~ 2 / r 2 - £2 ' fi2 = ~ 2 -

F£2; r = C p / C S.

Table 2

List of normalized variables Natural

Normalized

coordinate

x,y,z

~x "try X=~,fi=~,Z=2h

time

t

? = ~rCs t 2h

radial wavenumber

k

/~= 2 h k 'IT

frequency

¢.o

o~ =

2h

"rrC s

o9

~"'(z's, L ~)

= o

"rr2

I

-(n)

-

I

-In)

COS(~Ti'a AR) sln(~Ttfl mR)

~'~(z; L ~)

F. Santosa, Y.H. Pao / Transientplate response

279

For the SH waves, we can explicitly calculate the dispersion relations, and the eigenfunctions associated. For the Rayleigh-Lamb waves, the dispersion relation can be found numerically only; hence, we do not have closed-form formulas for the eigenfunctions either. The symmetric SH waves satisfy (3.6a). Hence, we must have /3~)=2n

for n=O, 1 , 2 , . . . ,

from which we can find o5~) using ¢i(s~)2 = OS(s~)2-/~2. We will take OS(s~)(/~)=(4n2+/¢2) '/2

for n =0, 1 , 2 , . . . ,

(3.9a)

with the understanding that -OS(s~)(/~) is also an eigenvalue. Inserting this result in the matrix equation (3.5), we find that P3 =-i('rr/2h) 3, hence ~(s")(z; ~ ~ ) =

~ ) [ - i 0 ; ] cos(2nS),

p, =Pz = O. We

choose the normalization

n =0, 1 , 2 , . . . .

(3.9b)

For the antisymmetric SH waves, we find

tfik")s(k)=((2n+l)Z+k2) '/2, n =0,

1,2,...,

(3.10a)

and

,,

~/s")(z; (~ ~)=t-2~) [ i,Jsin(2n+l)~,

(3.10b)

n:0,1,2,...

Observe that tS~s~)(/~) and 05~A~(/~)are simply hyperbolae. We display them in Figs. l(a) and l(b).

symmetric SH dispersion

antisymmetric

12

SH

dispersion

12

9

b

9

U

6

6 N

~

3

~

o

....

0

I,,

,I

,,,I

....

3

o

.....

L,

,,I

....

I ....

b

3 6 9 12 0 3 6 9 12 normalized wave number normalized wave number Fig. 1. (a) The dispersion relation for SH-modeswhich are z-symmetric.These are hyperbolaedefinedin (3.9a); (b) the dispersion relation for SH-modeswhich are z-antisymmetric.These are hyperbolaedefined in (3.10a).

K Santosa, Y.H. Pao / Transient plate response

280

Next, we look at the eigenvalues corresponding to the Rayleigh-Lamb waves. There are no closed form solutions O3~R}(/~) and O3(A~(/~) for every /~ We can, at best, hope for accurate numerical values for them. However, f o r / ~ = 0, these eigenvalues can be found. They are appropriately called cut-off frequencies in the elastic waves literature. Let us look at the Rayleigh-Lamb frequency equation for the symmetric waves. From Table 1, for/~ = 0, we have (~(n)'i4

wsR)

1

-(n)

•

--(n)

1

COS~qTCeSR s m ~ q ' r f l s R = 0 .

The solutions to the equation are 6(S"R)(0) = 1, 3, 5 , . . .

and

/3~R)(0) = 0, 2, 4 . . . . .

This means that

) [ raVd(o), =

Recall that r is the ratio o f the P-wave speed to S-wave speed ratio. We call the dispersion relation (oSCs~ as a function of/~) for a fixed n, the nth branch o f the dispersion curves. Moreover, for each fixed /~, O3~s~ must be an increasing sequence in increasing n's. From the numerical observation that the Rayleigh-Lamb frequency equation does not yield double roots for a fixed/~, we will by the nth branch henceforth mean the dispersion relation whose cut-off frequency is the nth smallest value (with the zeroth branch having the smallest cut-off value). We have evaluated o3~(/~) using a N e w t o n method for finding the root. The ratio was selected for r = 1.985, corresponding to Poisson ratio v = 0.333. The dispersion relations are shown in Fig. 2(a).

symmetric 12

t

f

I

J

i

antisymmetric

RL d i s p e r s i o n i

i

i

i

i

RL d i s p e r s i o n

12

i

/ i

/ / / / / / /

G

0 0

3 normalized

6 wave

9 number

12

0

/ / /

.... 3 normalized

I,

9 6 wave number

b 12

Fig. 2. (a) The dispersion relation for the Rayleigh-Lamb modes which are z-symmetric (see Table 1); (b) the dispersion relation for the Rayleigh-Lamb modes which are z-antisymmetric (see Table 1).

F. Santosa, Y.H. Pao / Transient plate response

281

Although we cannot find a3~)(/~) explicitly, it is still possible to arrive at relatively simple equations for the eigenfunctions. The first two rows of (3.5) work out to the conclusion that P3 = 0. If we choose p, = i0r/2h) sin(½"rr/3~)), we get ,~ ~ ( n ) ~ ( n )

~r ~=S._______R__~S.__R_R . /! _(~)a P2 = - i ~ (/~R) 2 _/~2) smt2~r~sR,. Whence, 'iT

i~E~")(g;/~) J , i~E~")(g;/~) _.:(") ~:(")r ~./~) t,~ S R ~-, 2 ~.¢-,

(3.11a)

with ,~x,(n)H(n)

E~,)(e;/~) = sin(½~r/~R)) cos(e6~))

,-~,SR.SR (~,~~) sin(½,tr6~) ) COS(g/~R))' 2/~ 2

(3.11b)

E~2")(e;/~) = sin(½~r/3~R)) sin(e6 ~R)) ~ (/~R)2 _/~2) sin(½~rti~R)) sin(g/~)) • Similarly, we can compute the cut-off frequencies for the antisymmetric Rayliegh-Lamb waves. This time 6~"R)(0) = 0, 2, 4 , . . . and /3~A~(0)= 1, 3, 5 , . . .

/

ra ~ ( 0 ) , ~<£d(E = o) = [ ~i~(o),

We will again use the same labelling convention discussed earlier. These dispersion curves are shown in Fig. 2(b). The eigenfunctions are i~F~")(~;/~) ] ~n)(g; ~ ¢/) =2"h

i~F~n)(z;/~) J, 8 (;R)F~2")(g;/~)

(3.12a)

with ,).,~(n) ~(n) F~"~(~;/~) =

COS(~rflAR) "ix AR/JAR COS(~WOLAR) , -~.~ sin(;~6
(3.12b)

F~:)(~; ~) = cos(~$~) cos(~ak~)-~ ($k~_ ¢:) cos(½~a~) cos(~$k~). Finally, we need to evaluate the norms of the calculated eigenfunctions. These quantities, defined using the inner product in Section 2, are needed to find the transformed displacements due to arbitrary forces (cf. (2.12)). The calculations are quite straightforward although extremely tedious. The results are f 2 h ( "rr ~6/¢ 2 for n =0, 1[~3~s")112= [ (\2i]6 h ~ /~2

(3.13a) for n>O.

F. Santosa, Y.H. Pao / Transient plate response

282

h('rr/2h)2

{(fi~R~2_£2)2sin2(½,rrfits~)[(&~s~2+£2)

•

2 1 - ( n )

+4ff'~R~2£2 sin (i'rrasa)

sin('rrt~
£2)]

[--(n)2 sin('rr/3~ )) ] (J~SR + £ 2 ) + Tr~(sn)- (~SR--(")2--£2)

8 tc2~m.)2 \ ~L,,SR _ ~72~(°) iv I t.t SR sin(,rra ~R~) sinZ(½w/~s~)}

(3.13b)

,.ff

II,~s")ll~= h

£~ for all n, h(w/2h)2

Jtff(.)2_£2)2

II,~"ql~-(t~k~_£~)~[,,_,A.

.-(.)2,--2

(3.13c) ~,

-{., [ _ ( . , 2

cos (~'~t~.) (~A. +£~)~ 2,-(.,[(-(.)2

-2

+~t~fAR K COS (~'Ti'O/AR) ~AR + k )

sin&rG~A~)t~(.)z_/~2)]~r~(;R )

,"AR

sin(cr~(A~),r~{.)2_/~2)] r~(n) \/.-~AR "l'rp A R

+8~r £:(/3{A~2-- £2)~A~ sin(~rG{g~) COS:(~r/3Aa)}.~ -(")

(3.13d)

Notice.that all the norms depend only on /~ Let us summarize what we have done in this section. We solved the eigenvalue problem associated with the transformed equations of motion. The solutions were separated into four distinct sets. We have symmetric SH waves symmetric RL waves

{,~d(£), ,%")(e; L ,~)},

antisymmetric SH waves antisymmetric RL waves Together, the eigenfunctions form a complete orthogonal set for each ~ and O in [L2(-~r,' ~r)]~3. With the eigenvalues and eigenfunctions evaluated, we can now use the expansion formula in (2.11)-(2.12) to find the transformed displacements due to applied forces. Indeed, in the next section, we will calculate the response of a plate to a single point force with arbitrary time function.

4. Displacements due to a point force

In this section, we will obtain formulas describing the response of an infinite plate to a point force. The excitation is assumed to have a smooth time signature, and has arbitrary orientation. We choose f ( x , y, z, t) = f ( t)~(x)5(y)5(z - ~)d

where f ( t ) is a smooth function of time with compact support contained in [0, ts]. The constant vector d is a direction vector with components d = (cos y, 0, sin y)W, 0 ~ y ~ 2~r. Thus, we have located a point force at (0, 0, ~') oriented in the xz-plane at angle y. The objective is to find u(x, y, z, t) for t > 0.

F. Santosa, Y.H. Pao / Transient plate response

283

We proceed by computing the Fourier transform o f f ( x , y, z, t) in x and y, f ( z , t; ~ 4 ) = f ( t ) ~ ( z - ~ ) d .

In order to find the transformed displacements ~(z, t ; ~ 4 ) (see (2.11)-(2.12)), we need to compute the inner products o f f with the eigenfunctions ~(s")(z; ~, 4), ~")(z; ~, ~), ~s~)(z; ~, 4) and ~")(z; ~ 4). These eigenfunctions are known to us from our calculations in the previous section. The functions ~(s") and ~(s") are known explicitly; see (3.9b) and (3.10b). Their inner products with .l~ are easily found: (f, ~(s")) =

~

( f ~3(s"))= -

i4 cos y c o s ( 2 n ~ ) f ( t ) ,

(4.1)

~-~ i 4 cos y s i n ( 2 n + l ) ( f ( t ) ,

(4.2)

where ( = ~ff/(2h). The inner products with 13~") and ~3(R "), which are known only implicitly (see (3.11)-(3.12)), can also be calculated. They are found to be ( f z3~")) = ~ (f, ~ " ) ) =

[ - c o s yi~E~")(ff; k ) - s i n y~(s~(k)E(f)(ff; k ) ] f ( t ) ,

(4.3)

~-~ [ - c o s yi~F~)(ff; k ) - s i n yd~(k)F,(")((;/~)]f(t).

(4.4)

Observe that even though (4.3)-(4.4) are known only implicitly, the inner products have a very simple structure. Their dependence on ~ is that they are linear functions of ~ with coefficients which depend on ~ but keep in mind that ~2 = ~2+ 42. Before obtaining the Fourier coefficients a (~) described in Section 2, let us separate the coefficients into four groups. The Fourier coefficients multiplying the eigenfunctions ~3(s") and ~(s~) are denoted by a(s~) and b(s~) respectively. The Fourier coefficients multiplying the eigenfunctions t3~") and ~ " ) are denoted by a~~) and b~") respectively. They will be functions of t, ~ and ~. Proceeding with the calculations, we obtain from (2.12), (4.1) and (4.2) i4('n'/2h) 3 cos y ,~ ~-, f ~ a~s")(t") = o ~ ~ 2 costzn~) Jo d r f ( r ) sin(o3(s~)(/~)(f - r)), /-, s ss~ 111 s ii b(s~)(t-) -

i4(~/2h)32 -~,) _ cos y sin(2n + 1)g

OCso~As(k)ll ~s ~)11=

I0

d r f ( ¢ ) sin(o3(k~(k)(~ "- r)).

(4.5a) (4.5b)

We have introduced another normalization ?= ~ r C s t / ( 2 h ) . Here all the terms going into (4.5) have been explicitly solved. They can be looked up in (3.9b), (3.10b), (3.13a) and (3.13c). For the Rayleigh-Lamb modes, we get from (4.3) and (4.4) a~.) -

2 h i - c o s yi~E~")((; k ) - s i n ya~s~(/¢)E~2"~((;/~)] f o 2 -(.) - ^(.) 2 d r f ( r ) sin(a3SR(/C)(/---¢)), ~rpCst°sR(k)llvR II

b ~')-2h[-c°sy~iF'~)(~;k)+snifl

--

~" £)] y a-(") Aa( k- ) F2(") (~,

-2

-(n)

--

pCs~A.(k) II

II =

r

fo drf(r)sin(o3Aa(/~)(/--r)).

(4.6a) (4.6b)

For the terms going into (4.6), they are given implicitly in Table 1 (for the frequency equation to find o3~g) and WAR~, .~(.)~ (3.11b), (3.12b), (3.13b) and (3.13d).

F. Santosa, Y..H. Pao / Transient plate response

284

We assume henceforth that a(s")(?; ~ ~), a(R~)(?;~ ~), b(s")(?; ~ ~) and b~R")(?;~ r~) are known. This means that we know the transformed displacement field given by

a(z, t; ~ ~7)=E "S"(n)'~(~a-a(R~)V~R ) - a--- h ( " ) ' : (--"us ~ a -~'s h ( ~__) 'u sR

(4.7)

n,R~(").

n

The eigenfunctions have been determined earlier. They are displayed in (3.9b), (3.10b), (3.11) and (3.12). We can consider ~(L ?; ~ r~) as known too. The only task left is to invert (4.7) to find the physical displacement. In the subsequent section, we will discuss the numerical considerations in finding the inverse Fourier transforms to obtain the physical response of the plate.

5. Numerical considerations

The inverse transform formula (from (~ fl) to (~, 37)), together with the representation in (4.7) allows for a natural physical interpretation. To obtain the symmetric SH-response, we would take Fourier coefficients a~s")(?; ~ ~), and the eigenfunctions ~3(s")(~;~ ~), and insert them into the inverse Fourier transform formula. We would arrive at an expression such as Us(g, 37, L t-) = ~

d(

-co d~ ~ " ~ ( L ~; :?) ei(¢~+'~Y)"

d r f ( r ) sin(o3("~(? - ~-)).

(5.1)

Here, ~ " ) is called the modal amplitude function, and is a combination of contributions from the eigenfunctions, their norms and the nature of the applied force. The explicit form of • for our problem will be given later. The issue at hand is the numerical determination of (5.1). We can therefore view u as the convolution of the Green's function with the source function f(t). That integrals such as (5.1) make sense is guaranteed by smoothness hypotheses on f(t) (see [7]). We expect the Green's function itself to be a distribution. In [7], estimates were given on the accuracy of representing the solution with truncated integrhls and only summing a few terms. It is inevitable that one has to approximate the true solution this way, but one must go about it in a rational manner! To understand (5.1), let us take the Fourier transform in time. We use the variable ~o to denote the transformed variable. Thus, lis(~, 37, ~; q~) =

dt e

u(x, y, :g, ?).

(5.2)

Using the transform in (5.1), we arrive at . j?(~) ~i(~, 37, ~; ~o) = 52

d~

00

d~ qCn~(~ ~/; if) ei(g~+~f~;~" 1i { ~ ( ~ - ° ~ ° ~ ( ~ ) ) - ~(~° + ° ~ " ~ ( ~ ) ) } '

(5.3)

where g(. ) denotes the Dirac delta. Equation (5.3) reveals a simple rule by which we should truncate the integrals in ( and ~. Suppose we have a source f(t) whose temporal Fourier transforms are such that

If(,~)l~
,1"°<~ I'PI~< ~

for[~0+<~[~[
(5.4)

where e is small, q~_/> 0 and ~+ << 0o. Thus the excitation wavelet has its energy mostly within the passband [~0_, q~+]. The second derivative of a Gaussian wavelet falls into this class.

285

F. Santosa, Y.H. Pao / Transient plate response

With a source having characteristics as (5.4), it is clear that ~(£, )7, ~, ¢) will be small outside the passband ~o~ [q~_, ~o+]. Since we will only be concerned with the passband response, we replace the integration in (5.3) by the truncated version

(5.5)

where B, = {(~ ~ ) : ~ _ ~< [o5(")(/c)[~< q~+}. The absolute value signs around o3(")(/~) reflect the fact that if o3(")(/~) > 0 is an eigenvalue, so is -o5(")(/~). Notice that we now have a finite sum; this is natural because there are only a finite number of branches (n's) that satisfy ~_ ~< ~+. Recall that for a fixed/~, o5(")(/~) is an increasing sequence. Moreover, o5(")(/~) are increasing functions of/~ for/~ large. Under this rule, we see the following picture in the o3 - / ~ plane (Fig. 3). Integrations are performed for the part of the branches of the dispersion relations which are "trapped" within the band [q~_, q~+]. Thus the cut-off radius B, varies with respect to each n.

12

....

I . . . . . . .

I'

''

9

"e.

6

¢=¢+

o

I .... 0

3

I

I .... 9

B' 0

-

12

k

Fig. 3. Graphical description of finding the integration domains B ' . On this illustration, we find B~ for the symmetric SH-modes. The excitation is supposed to be bandlimited to [¢_, ~+].

Observe however that (5.5) may be tricky to evaluate because of the delta functions. Thus we revert to the original temporal representation. Only now, we have the finite sums and the truncated integrals. The approximation to (5.2) is

u(g, fi, ~., t-)=

drf([-r)

~,

n~N

n

d~Td~ q}(")(~, ¢/; £) e~(¢~+'~) sin(o3(")(/~)r).

(5.6a)

If we look back to ~ ( " ) ( ~ 7}; z) for all the modes of vibration, we can see that all the ~(")'s have a factor

F. Santosa, Y.H. Pao / Transient plate response

286

for [o5(")(/~)] -~. This means that we can integrate over positive o5(")'s if we multiply by a factor of 2 the expression in (5.6a). Thus change B, to

B.I __ - {(~ ~): 0 <~ ~ _ ~ oS("~(/~) <~ ~+}. We have u(~, y, ~, ?) = 2

L

d ~ ' f ( ? - ~) Y~

d(d¢/q~(")(~ ~; ~) e~((~+~Ylsin(o3("~(/¢) r).

(5.6b)

In applications, we normally fix the coordinate point (~, y, ~) and seek a time trace representing the particle displacement. Looking at (5.6) above, we would first evaluate the integrals in ~ and ~, then perform the sum. To obtain the desired output trace, we convolute the result with the source. Of course, for each f(t), we have a different set of B " s and hence different integrals. We can, however, set [q>_, ¢+] to be sufficiently large. Then the inner expression which convolves with f ( ? - r ) in (5.6) may be thought of as a bandlimited Green's function. Convolution of this bandlimited Green's function with sources f ( ? ) whose passband are contained in [q~_, ¢+] gives the transient response trace at the chosen coordinate point. As far as sampling the time trace is concerned, any sampling rate for which the Nyquist frequency is greater than the upper end of the passband ¢p+ will do. Thus, the rule is to choose A[ such that A?<~ -rr/¢+. The only issue left to address is the numerical integration in the ~- and ~-variables. In the applications that we have in mind, we seek response traces for Ix2+y21 large. By "large", we mean 1 2+y21 '/2 i.e., many times greater than the plate thickness. This mean that the exponential term will be very oscillatory. Notice also that for relatively large r, because o3(")(/~)~constant /~ for large /~, we expect the term sin(oS(")(/~)r) to be very oscillatory. On the other hand, the ~ ( ' ) ( ( , ~; z) are nonoscillatory. For the z-symmetric shear modes, the modal amplitude function can be found from the Fourier coefficients a(s")(?; ~ ~) and the eigenfunction t3(s")(~; ~ ~). It is given by q~(,)(~, ~; ~)

C , cos(2n~) (~2, - oS~s~(~)ll ~11~ -~,

o) w

where C, = ('rr/2h) 3 cos Y cos(2nO/(16h2pC~). Now inserting (3.9) and (3.13a), we arrive at _ff~ (~2, _s¢-~, O)X, Z)

n=O

/

=]

(5.7) C" cos(2nz) ,_2

.>0

where C~ = (CJh)(~/2h) 6. The modal amplitude functions decay like 0(/~-~), with the one for n = 0 singular at/~ = O. This poses no difficulty as we are only interested in the response to a bandlimited pulse, and will only integrate ~¢"~(~ ~; z) against the oscillatory terms over a closed set B', not including zero. Moreover, it is an integrable singularity because of the sin(~3°(/~)~-) term multiplying ¢P. The modal amplitude functions for the z-antisymmetric shear waves are ~(-)t~ D" sin(2n + 1)~ ~s~s, ~; z) = ((2n + 1)2+/~2),/2/~ ( ~ , -~-~, 0) x, where D" = (~r/2h) 3 cos 3' sin[(2n +

(5.8)

1)(]/(16h~pC~). For the Rayleigh-Lamb modes, we only know the

F. Santosa, Y.H. Pao / Transient plate response

287

modal amplitude implicitly. These are given by ~E~")(~,/~) '~R)(~

,~; ~ ) -

2 _(,,)

pcs,,.,s~1 II~k"ll: I

cos ~ E ? ) ( ~ ,

~)

]

~)

~E?)(~,

~-~-~-~(~)~(")~ ~ /~) Jt~ o t S R -L, 2

- s i n y E(2")(~/~)

~-%

i:Z(-) r:(,,)~ ~/t.,t S R .L., 1 \ Z~, , , k)

]1

(5.9)

for the symmetric RL-modes. All the terms entering the expression above may be looked up in Section 3. For the antisymmetric RL-modes, we have

1

{

=,°,,

~s~,~,-~

F ~2F~[)(~'/~) ]

cos

E)

---() () ['~aARF2 (Z, k)

,,

-

-

[ i~--~(A~F~")($,/~) ]) .--(n) +sin yF2 (K, k)[ ,*;a kRF,( n ) (z,- k) ]~. (n)

-

-

(5.10)

These two functions can be evaluated numerically once o~(")(/~) has been found. Notice that all the modal amplitude functions depend on ff and @ in a simple way. All the terms entering into (5.7)-(5.10) are functions of/~; the dependence on ff and @ are explicitly displayed. From (5.9) and (5.10), we can clearly see the contributions to displacements by each component of the single force. Having evaluated the four types of modal amplitudes, the remaining task is to insert them in (5.5b) and evaluate the integrals to obtain the displacements. We shall go a bit further and find B'~ for the case of the shear modes. The purpose is to indicate a way of along the same for the RL modes. Suppose we have chosen ~o_ and c~+. We insert it in the equation defining B" when o5(')(/~) corresponds to oS(')t/~Xss ~ J or OS(A~(/~).Suppose further for the sake of illustration that 1 < ~p_ <2, and 5 < ¢+ <6. For the z-symmetric case, we solve for ((, @) such that ~o_<~(4n2+/~2)~/2<~0+,

n =0, 1 , 2 , . . . .

(5.11)

We get (see Fig. 3)

But it is only for n<<-2 that we can solve (5.11) for the chosen [~0_,¢+]. Therefore, for this case, N={n:n=O, 1,2}. Similarly, for the z-antisymmetric modes, we have O3(A~(/~)=((2n + 1)2+/~2)~/2, and we get

B~ = {(~ @) : I/~[e [x/~2_--~1, x/~02- 1]},

B'={(~):l~le[O,,/¢~+-(2n+l)2]},

N={n:n=O, 1,2}. Naturally, this process is carried out for the RL modes as well. Only now, because 03~R)(/~) and O~(A~(/~) are only known implicitly, we determine B" and N numerically. So, for each n, for a particular mode of vibration, we need to evaluate / a : / d ~ d @ ~(")(~ @; ~) e'(~+~') sin(o3(")(/~) ~')

(5.12)

F. Santosa, Y.H. Pao/ Transientplate response

288

where the B " s disks (or anuli) in the ( - ~ plane. In view of the simple dependence on ~ and ~ of the modal amplitude functions, we shall integrate in polar coordinates and integrate out the angular variable to obtain an integral in the radial (/~) variable. The change of variables is given below. (~ ¢/) -~ (/~ cos q~,/~-sin 0),

0 ~ [-'rr, "rr],

d ~ d ~ = k dk dO,

(~,)7)=(RcosO, RsinO),

0 c [-'rr, "rr],

f~zf d~d~'"-" f~, r(d~f_\dO'" .

The integrals we will encounter can be integrated d ( d ( / e i((~+~y).

i~7

in 0. We list them below.

• .~ fB k dkJo(kR)" • ",

d ~ d ~ i~: e i(~+'~).

d(d~

exactly

• • ~ fB k dkJ'o(kR) ,;

e i(~+'~y) • • "~

d~?dr~ ~2 e i ( f ~ + o ~ ) . . . ~

cos 0. • • ,

fa k dkJ~(kR) sin

0. • • ,

fB k dk[-J~(kR)(1

-2

sin 2 0) +Jo(kR) sin 2 0]. • . ,

(E)

fB fd~-d(l~'12ei((~+"Y'.'.-)fa kdk[J~(kR)(1-2sin20)+Jo(kR)coseo]'"

(F)

fa f d(dfl-~ei(¢~+~)...-, fB kdk[-2J'~(kR)-Jo(kR)]sinOcosO..

,

.

Here the dots refer to the components of the integrand that depend only on /~ The kernel Jo(kR) is Bessel's function of order zero. Thus, integrals like (5.12), which contribute to the displacement in (5.6b) are now single integrals over the variable /~ Our question remains: How finely should we sample the integrand in /~ in the numerical quadrature? This is determined by the rate of oscillation of the exponential and sine terms. In the application we have in mind, (~,)7, ~) = (R cos 0, R sin 0, ~) is fixed with g >>½"rrand )7 >>½~r. Suppose the time trace is the interval z e [0, T]. We assume that T is several multiples of R. This means that the sine term will oscillate very rapidly at r = T, more so than the exponential (Bessel function) term. Looking back to Fig. 1 and Fig. 2 (the dispersion relation) we see that ~

03(k)(/~) <~C~

where Cp is the P-wabe speed. Thus the shortest wavelength in the ~ - ~ plane of the oscillatory term is 27r

This is because o3(")(/~) ~

(Cp/Cs)k+

constant, and the Bessel functions have arguments (kR). This shortest

F. Santosa, Y.H. Pao / Transient plate response

289

wavelength is adequately represented by about six sample points. Therefore, 2

6[CPT+R) "

(5.13)

\Cs

We describe our numerical implementation below. (0) Problem parameters: coordinate point (x, y, z), angle of force y, plate thickness h, wavespeeds Cp and Cs, and source wavelet f(t-). (1) Compute dispersion relations (z-symmetric SH and RL, z-antisymmetric SH and RL), cf. Table 1. (2) From f(t-), determine passband [tp_, ~o+]. (3) Find the largest set B" = {(~, ~) :0 <~~p_~
for the contributions which are regular at/~ = 0, and contributions which are singular at/~ = 0, (5) From the maximum time of interest T, determine the sampling rate A/~ in (5.13). Determine desired temporal sampling rate using the criterion Ar <~~-/~0+. Perform numerical quadratures in/~. (6) Increment r<- r + At. (7) Add both contributions. Perform numerical convolution with f(t-) to obtain the transient response at the designated coordinate point (cf. (5.6)). Even though we determined the regions of integrations B ' , which may be different for each mode and index n, we choose to integrate from/~= 0 (/~_ for the singular case) to the radius of the set covering all of B" for all n's and all modes. This is because we want to do the least number of integrations, which may be costly. This means that the Green's function is computed for some branches to higher frequencies than ~0+. But this presents no distortion to the output course because the wavelet is bandlimited to [~0_, ~p+]. However, we cannot overlook the importance of step (3); in this step, the number of branches to be summed to give the needed approximation is done. From Figs. 1 and 2, we note that the only instances for which B" could be a disk are for the branches o3
290

F. Santosa, Y.H. Pao / Transient plate response

In our implementation, we chose a trapezoidal rule with sampling rate A/~ found using (5.13). The calculation can be made more efficient by making A/~ depend on T, accomplished by replacing T by ?. However, this mean we need to resample the integrand for each evaluation, which could be costly. We also note that since r is typically large, we could possibly use the asymptotic forms of the Bessel functions, which would lend itself to using techniques of integration specifically designed to integrate oscillatory functions such as FFT or those described by Davis and Rabinowitz [10]. Since computer time was not an issue for us, we chose a scheme which was easiest to implement, though clearly not the most efficient.

6.

Numerical

example

For a demonstration, we selected a plate of high-strength aluminium. The specifications are density

p = 2800 k g / m 3,

Young's modulus

E = 72 x 109 N / m 2,

shear modulus

# = 27 x 109 N / m 2,

thickness

h = 0.02 m.

From these data, we calculated Poisson's ratio

v = 0.33,

shear wave speed

Cs = 3100 m/s,

P-to-S speeds ratio

r = 1.985.

We shall calculate the transient responses at receiver locations separated from the source by (A) 20h = 0.4 m (R = 10T), and (B) 40h =0.8 m (R = 2 0 7 ) . To excite the medium, we chose the time function displayed in Fig. 4. The high end of the frequency spectrum o f f ( t ) is about 0.5 x 106 Hz. This normalizes to q~+= 12.90. The low end is about 0.8 x 105 Hz, or ~_ = 2.06.

excitation 80

I

. . . .

power spectrum I

'

1000

'

8O0

40 "~ 600 0

~400 -40

2OO ,

-80 0

L

. . . .

I

50 100 t i m e in m i c r o s e c o n d s

,

a

0 150

t

100

200 frequency

300 400 in kilo H e r t z

Fig. 4. (a) The excitation wavelet used in our numerical simulation; (b) The power spectrum of the wavelet in (a).

i00

F. Santosa, Y.H. Pao / Transientplate response

291

The number cp+ allows us to find B" and N, which tells us which branches need to be taken into consideration (refer to Figs. 2(a) and 2(b)). We find that for both the z-symmetric and z-antisymmetric displacements, we need the index set

N={O,

1,2,3,4,5,6,7,8,9}.

The largest of the B" being the one corresponding to n = 0. For B~, we take/~+ = 13.0. In Table 3 a list is shown of the cut-off frequencies of the branches of the dispersion curves that enter into our calculations. Table 3 Index

0 1 2 3 4 5 6 7 8 9

z- s y m m e t r i c

z-antisymmetric

cut-off frequency ~(~)(o)

cut-off frequency ~k~(o)

0 1.985 2 4 5.955 6 8 9.925 10 12

0 1 3 3.970 5 7 7.940 9 11 11.910

Both the receiver and the excitation will be placed on the surface z--h. Thus, g = ½"rr and ~ = l'tr. We will be interested in the response up to time 6 x 10 -4 s. This means that T = 146.089. The sampling rate in time is A f = 0.109956, corresponding to 4.5 x 10 -7 s. With these specifications, we use rule (5.13) in the previous section to arrive at the wave number sampling rate of A/~= 3.0 × 10 -3. We shall only consider z-displacements. Hence, we can completely ignore the SH-displacement fields. The integrals we need to evaluate are (from (5.9), (5.10) and (5.6)) ___ gs(t)-pC~

d/~/~ oEN

cos Y cos 0

+sin'),

fo ~+

d/~/~ ~N

[

~sR~

t ~ , k)E~")(l~; k) sin(o3(s~T ) J~(kR) o3 dllahO)ll 2

{ ~ s ' ; ~ E r ' ( ~ ; .k)} . . _(.)_. 1 - ] o~dli,V)ll~ s,n(tosRt)JJo(kR)~, -(~) (~), -

(6.1)

and

1

gA(t)=-'~(COSyCOSOpL's t + s i n ~,

dkk

IOl~+ [

d/~/~ ,,ZN

~ Ln~N

g'(n)r~(n)'¢l "

~ARr, t~r,/¢)F~n)(l~r;/¢) sin(o3k~?) ak~ll a~.)ll2

-(~) ~(.)

O.)ARt/R

2

] J~(kR)

•

In the expressions above, notice that the contributions by the vertical and horizontal components of the force separated in the two sums. To obtain the response due to the chosen time function for a given force orientation y and receiver location (R, 0), we simply compute uz(T) = f * [gs + gA](t').

F. Santosa, Y.H. Pao

292

/ Transient plate response

In each of the integrals above, the oscillatory terms in the integrands are the Bessel functions Jo(kR) and J~(/~R), and the sin(o3(")?) terms. We removed these oscillatory components from the integrans and computed the resulting expression as a function o f / 7 for several n's. We shall refer to them as the modal amplitudes (in cylindrical coordinates). These modal amplitudes are displayed in Figs. 5 and 6. In Fig. 5(a), we showed the modal amplitudes for z-symmetric displacements due to the vertical component of the force. In Fig. 5(b), we display the corresponding amplitudes due to the horizontal component of the force. Figures 6(a) and 6(b) show the same functions for the z-antisymmetric displacements. The results of our computations are displayed in Figs. 7-9. In the first set of simulations, we simply looked at the separate responses due to the horizontal and vertical components of the force. When y = 0

symm. 0.04

modal

. . . .

amp.

I . . . .

T

horizontal

. . . .

i

symm.

, ~-,a

o.2o

,,-,,

modal '-~'

amp'

'~

I

vertical

. . . .

T''

b

n = 2

0.00

0.15

.g

n

g

n=l -0.04

0.10

-

E -0.08

-

n=O

0.05

-0.12

0.00

0

3 normalized

6 wave

9

12

t_~__~

0

,

3 normalized

number

~

k

6 wave

,

,

t_

~ L

9

12

number

Fig. 5. ( a ) T h e mode amplitudes for the z-displacement due to the horizontal component of the applied force. The displacements are z-symmetric, both source and receiver are on the surface (refer to (6.1) with 3' = 0 ) ; (b) same as in (a), except that y =½~r here, thus these mode amplitudes are due to the vertical component of the force.

antisym, 0.04

....

modal

I ....

antisym,

amp. :horizontal

I ....

i,i1~4

-1

0.00 :D

modal

amp.:

vertical

0.20

~

-0.04

0.15

0.i0

E

:e

-0.08

0.05

n=O

n-3

-0.12

....

L . . . . . . . 3 normalized

6 wave

2 ..... 9 number

,

0.00 12

0

, in-i 3

normalized

i i i ,n 6 wave

. 9

12

number

Fig. 6. ( a ) T h e mode amplitudes for the z-displacement due to the horizontal components of the applied force. The displacements are z-antisymmetric here (refer to (6.2) with 7 = 0 ) ; ( b ) same as (a), except that 3' =½1r here--vertical force.

F. Santosa, Y.H. Pao / Transient plate response 50

'1

....

J ....

I ....

50

I ....

. . . .

293

I ....

I ....

I ....

I ....

I,,

I .....

[ ....

I,,,

Q) © o

25

m "-U

0

Q

-

25

0

-

"u N

-25

-25 c~

,I

-50

....

90

J ....

I ....

t ....

270 360 t i m e in m i c r o s e c o n d s 180

a

-50 450

90

180 270 360 t i m e in m i c r o s e c o n d s

b

450

Fig. 7. (a) The z-displacement response of the plate to a horizontal force. The receiver location is directly ahead of the applied force at a distance of 10 plate thicknesses; (b) the z-displacement response of the plate to a vertical force. The displacement is radially symmetric, the receiver is 10 plate thicknesees from the force.

50

, ,

~ - '

[ ....

I ....

50

I ....

0)

25

-

25 o

m

0 ~D bl

-25

~-25

-

8

I .... I .... I, , , i! -50 180 270 360 45 0 90 180 270 360 450 t i m e in m i c r o s e c o n d s t i m e in m i c r o s e c o n d s Fig. 8. (a) The z-displacement response of the plate to a horizontal force. The receiver location is directly ahead of the applied force at a distance of 20 plate thicknesses; (b) the z-displacement response of the plate to a vertical force. The receiver is 20 plate thicknesses from the force. -50

~ d

0

,,I

....

90

in (6.1) a n d (6.2), w e o b t a i n t h e r e s p o n s e s d u e to a h o r i z o n t a l force. T h e d e p e n d e n c e o f t h e r e s p o n s e o n t h e a n g l e is r a t h e r t r i v i a l ; it is s o m e t - d e p e n d e n t f u n c t i o n t i m e s cos 0. H e n c e , f o r 0 = 0, w e w o u l d o b t a i n t h e r e s p o n s e o f t h e p l a t e at a p o i n t d i r e c t l y a h e a d o f t h e a p p l i e d h o r i z o n t a l force. T h i s is p r e c i s e l y w h a t we p l o t t e d in Figs. 7 ( a ) a n d 8(a). I n Fig. 7(a), R = 10~r (10 p l a t e t h i c k n e s s e s ) , a n d in Fig. 8(a), R = 20~r. W e h a v e s c a l e d the figure b y ( p C ~ ) . In Figs. 7(b) a n d 8(b), w e s h o w t h e r e s p o n s e o f t h e p l a t e at R = 10T a n d 2 0 7 r e s p e c t i v e l y , d u e to a v e r t i c a l l o a d . T h e s e are o b t a i n e d b y setting y = ½~r in (6.1) a n d (6.2). T h e s e d i s p l a c e m e n t s a r e 0 - i n d e p e n d e n t . T h e final set o f s i m u l a t i o n s c o n c e r n t h e r e s p o n s e o f t h e p l a t e to a f o r c e w h i c h is i n c l i n e d at y = ~r. W e h a v e c o m p u t e d t h e r e c e i v e r r e s p o n s e at R = 10~r f o r 0 = 0, z~r, ~ ~xr, ~ z~r, 3 a n o 7r. T h e s e r e s p o n s e s a r e d i s p l a y e d in Fig. 9.

294

K Santosa, Y.H. Pao / Transient plate response 60

I

. . . .

I

. . . .

I

. . . .

I

. . . . [

30

8

-

o

? -30

-

-60 90 60

-

. . . .

r

] . . . .

i

~

I

n

I

180 270 360 t i m e in m i c r o s e c o n d s

I . . . .

-

6e

r,,

450

i

''I

....

I ....

I ....

0=45 °

0=90 °

30

O3 -,.j

30

-30

-30

-60

-6o

90

.-

~

450

180 270 360 time in m i c r o s e c o n d s

.....

90

60

,

,

I, .

I . . . . .

60

i

''''I''

n

r

I

]

r

I

I

i

'I''''

0=135 ° 30

450

180 270 360 time in m i c r o s e c o n d s

0=180 ° - -

.~

30

!o -30

-60

-

d ~ L

0

J

90

,~

I .... ~

....

180 270 360 t i m e in m i c r o s e c o n d s

~

-3o

-60 450

0

90

180 270 360 t i m e in m i c r o s e c o n d s

450

Fig. 9. The z-displacement responses of the plate to a force at angle of ½n from horizontal (y =½~ in (6.1) and (6.2)). The receiver is 10 plate thicknesses away. Shown are the responses at various polar angles.

It s h o u l d b e r e m a r k e d t h a t t h e s e s i m u l a t i o n s a r e v e r y c o m p u t a t i o n a l l y

intensive. In view of the remarks

m a d e a t t h e e n d o f t h e l a s t s e c t i o n , it is c l e a r t h a t a m o r e e f f i c i e n t s c h e m e , a l l o w i n g f o r n o l o w - e n d c u t - o f f , can be developed.

Conclusions

We now summarize our findings, We have used the method of normal-mode

expansion

[7] t o f i n d t h e

transient response of an elastic plate to applied forces. The method makes use of integral transforms on

F. Santosa, Y.H. Pao / Transient plate response

295

the plane of the plate to reduce the equations of elasticity to a system of one-dimensional hyperbolic equations. These equations described the transformed displacement variables, and are solved using normal-mode methods. Our results are perfectly general and may be used to describe the response of a plate to very arbitrary disturbances. As the first application, we used our results to find the displacements caused by the application of a single force. We further showed that the representation for the displacement can be used to compute transient responses of the plate. Certain numerical consideration must be taken into account to assure the desired level of accuracy as well as computational feasibility. Finally, we illustrated the use of this computational technique in numerical simulations.

Acknowledgment We are grateful to Richard Weaver for many hours of helpful discussion, and for providing computational data (axisymmetirc waves) for independent check of our code. Thanks are also due to Willaim Symes for encouragement and useful suggestions. The research was supported by a grant from the National Science Foundation (MEA 81-17988) to Cornell University.

References [1] R.D. Mindlin, An Introduction to the Mathematical Theory of Vibrations of Elastic Plates, U.S. Army Signal Corps Engineering Laboratories, Fort Monmouth, NJ (1955). [2] J.D. Achenbach, Wave Propagation in Elastic Solids, North-Holland, Amsterdam (1973). [3] J. Miklowitz, The Theory of Elastic Waves and Wave Guides, North-Holland, Amsterdam (1978). [4] N. Vasuvedan and A.K. Mal, "Response of an elastic plate to localized transient sources", ASME J. Appl, Mech. 52, 356-362 (1985). [5] A. Ceranoglu and Y.H. Pao, "Propagation of elastic pulses and acoustic emission in a plate, Parts I-III", ASMEJ. Appl. Mech. 48, 125-147 (1981). [6l R. Weaver and Y.H. Pao, "Axisymmetric elastic waves excited by a point source in a plate", ASMEJ. Appl. Mech. 49, 821-836 (1982). [7] F. Santosa, "On elastodynamic problems in infinite cylinders", Q. J. Mech. Appl. Math. (1989), to appear. [8] I.T. Lu and L.B. Felsen, "Ray, mode and hybrid options for time-dependent source excited propagation in an elastic layer," Geophys. J. Roy. Astr. Soc. 86, 177-201 (1986). [9] Y.H. Pao and R.K. Kaul, "Waves and vibrations in isotropic and anisotropic plates", in: G. Herrmann, eds., R.D. Mindlin and Applied Mechanics, Pergamon Press, New York (1974) 149-195. [10] P.J. Davis and P. Rabinowitz, Methods of Numerical Integration, Academic Press, Orlando, FL (1976).