Water waves generated by an impulsive bed upthrust of a rectangular block

Water waves generated by an impulsive bed upthrust of a rectangular block J I I N - J E N LEE and J A W - J O H N C H A N G Department of Civil Engineerin9, University of Southern California, Los Anyeles, California 90007, USA (Received 27 December 1979; revised 11 February 1980) A theoretical and digital simulation on waves generated by an impulsive bed upthrust of a rectangular block with various ratios of length/width is presented. The two-dimensional fast Fourier transform (FFT) algorithm is used to obtain the water surface profile near the generation region. The three-dimensional pictures were constructed from an array of 256 x 256 pixel by using image processing. Comparison of the present three-dimensional results, with previously published two-dimensional results, indicates that for large length/width ratios, the water surface profile is quite similar at certain locations for a small time after the impulsive bed motion is completed. Generally, the water surface profiles deviate significantly from the two-dimensional results.

INTRODUCTION Waves generated by submarine earthquakes, commonly known as tsunamis, have been of interest to ocean and coastal engineers. In many instances, such waves have caused significant damage to the coastal region. Estimation of the tsunami wave form is essential for the prediction of the wave profile in the propagation phase and the final wave form when the wave arrives in the coastal region. Most of the theoretical models of tsunami generation have been based on linearized theory in either a two- or three-dimensional fluid domain of uniform depth. Based on the linear theory, the fundamental solution of the wave forms can be obtained by multiple Fourier transforms such as that used by Driessche and Braddock 3 and Hammack 4 or by the Green's function method used by Kajiura 6,7. The complexity of the integral solution has discouraged many authors from obtaining detailed wave behaviour near the generation region even in the case of a simple bed deformation model. Instead, many authors have been using asymptotic methods to evaluate the integrals. Consequently, only the far-field wave behaviour can be examined. For three-dimensional wave generation models, the computational effort needed to evaluate the integrals by means of conventional numerical integration schemes have been quite laborious. To overcome this difficulty, the modern digital simulation technique can be used to significantly improve the computational efficiency. This paper presents results of a theoretical and digital simulation study on waves generated by an impulsive bed motion. For mathematical simplicity, the extent of the bottom disturbance is assumed to be of rectangular shape although an elliptic shape might be a better tsunami generation model 5. The fundamental solution for the wave amplitude is expressed in terms of multiple integrals which are evaluated by utilizing the two-dimensional fast 0141-1187/80/0401654)6$2.00 © 1980CML Publications

Fourier transform (FFT) algorithm. The threedimensional pictures (based on 256 x 256 samples) showing aspects of wave decay everywhere in the fluid domain at a specified time are obtained by the image processing technique.

PROBLEM FORMULATION Let (x,y,z) constitute a Cartesian coordinate system with z = 0 as the undisturbed water surface as shown, in the definition sketch in Fig. 1. Initially, the fluid is at rest with

A

-~-

DX

A

B

.I.

13

.1

q

z

i ~max~o Figure 1. Definition sketch offluid domain and coordinate system

Applied Ocean Research, 1980, Vol. 2, No. 4 165

Water waves generated by rectangular block upthrust: J. J. Lee and J. J. Chang the free surface and solid boundary defined by z = 0 and z = - h , respectively (h is a constant). For t > 0, the solid boundary is permitted to move as prescribed by ~(x,y;t). The resulting deformation of the free surface is to be determined as z = q(x,y;t). Assuming irrotational flow and an inviscid fluid, the fluid kinematics can be expressed in terms of a velocity potential q~(x,y,z;t). The differential equation and the linearized boundary conditions that ¢p must satisfy can be listed as follows: 0
V2~0 =0

-~
where 4o is the maximum amplitude of the vertical displacement, ~ is the time constant defined as 1.11/6, and 6 is the characteristic time defined by ~/4o = 2/3 at t = 6. The Heavyside step function is defined as:

H(B2-x2)={Io' '

and H(A 1 _y2)=

-h
B2x2<0B2-x2>0

{lo:

A z _y2 < 0

(1) qg, + gq~z= 0

z=0

(2)

~p~= ~t(x,y;t)

z= - h

(3)

In these equations, the subscripts denote partial derivatives and g denotes the gravitational acceleration. The linearized relation between the water surface displacement q and the velocity potential ~0 is:

The transformation of equation (7) yields: ^ sin K1B sin K2A F ~ 1 ~(K1'K2's)=44° KI " K22 ] _ s ( s ~ J

Substituting equation (8) into equation (6) and solving for the surface elevation yields:

rl(x,y,t) =

1

rl(x,y;t ) = -- L~Pt(x,y,O;t) g

(4) f iK x iK ySin K1Bsin K2A e ,e 2 ___

4o ---

The solution for ¢ (and hence, q) is obtained by using the Fourier transform for the spatial variables x, y and the Laplace transform for the time variable t defined by:

n2

JGc

i -- oo - orJ

K1

K2

1 cosh 4 K 1 2 + K 2 2 h

e ~t-cos e~t---sm Ct tot

x,y;t)ddydtx 0

- ~

(8)

(9)

dKldK2

(5)

where

By applying the transformation of equation (5) to the governing equation and boundary conditions with the subsequent inversion, one obtains the water surface elevation as: 1

q(x,y,t) = ~ × ~Znlf

f ( ~ 1. I s~e 2 -iK,x ^ e iK2vest 4(Kx,K2,S).ds']dKldK 2 \ z m J (s z + toZ)coshx/Ka" "-[-K2" ]'/ J

- a, - ~x:,

092 = g~/K 12 _~_K2 2 tanhx/K 12 +

K2 2 h

Equation (9) contains poles at K I = 0 and/or K2=0. These sigularities can be treated easily by L'Hospital's rule. The reader is reminded that equation (9) can be reduced to the two-dimensional solution given by Hammack 4 by integration with respect to K 2 as A--, ~ . (This is easily accomplished by changing the variable of integration to (K2.A)/A and evaluate the integral as A or.) Introducing K 1 = 2nk I and K 2 = 2nk2, equation (9) can be rewritten as:

Br

(6)

rl(x,y,t)=-44offei2"k,xei2"kaf(k1,k2)dkldk2

where

(10)

/a+iF

where f(kl,k2)is a symmetric function defined by:

i =lim~ f Br

IJ - iF

f(kl,k2) = is the Bromwich contour, to is the circular frequency and p is a positive constant. In order to obtain the specific wave profile as defined in equation (6), one must specify the bed defomation history 4(x,y,;t). The bed deformation is assumed to be a rectangular block upthrust with exponential variation in time:

~,(x,y,t)=4o(1 -e-~t)H(B2-x2)H(A2--y 2) for t>~0

166

Applied Ocean Research, 1980, Vol. 2, No. 4

(7)

sin 2nklBsin 2gk2A 1 2nkl 2nk2 cosh x/(2~kO2+(2nk2)2h

and to2 =gx/(2gk02 +(2nk2)~ tanh x/(2nkO2+(2gk2)2h

Water waves generated by rectangular block upthrust. • J. J. Lee and J. J. Chang It should be recognized that equation (10) is a twodimensional Fourier transform offlk~,kz) and can be computed by using the fast Fourier transform (FFT).

cosh x/(2nnl Ak,)2 + (2nnzAkz)2~

NUMERICAL IMPLEMENTATION

I ~ 2 ~ o 2 ( e - S ' - c o s c o t - ~ s i n cot)]

Suppose that f(k~,kz) is defined in the interval of - T/2 T/2. Equation (10) can then be written as:

1

and ~l) 2 "lt- (2nn2Akz) 2 x

G)2 = ( J N / ~ H l ~ k

T T

rl()c,Y,t)=-4~offf(kl,k2)ei2nk'xei2nk~ydkldk2

(11)

T-T

It should be noted that the conventional procedure of computing ~l(x,y,t) in equation (11) by fixing a set of (x,y:t) is quite inefficient. For our purpose, a more advantageous procedure is to fix a time t = t~ and evaluate equation (1 l) for all possible x and y by recognizing equation (11) as a double Fourier transform of f(kj,k2;ti). The integral in equation (11) can be approximated by a Reimann sum to be presented as follows: Choosing s > 0, then there exists c~> 0 such that:

-

~ 111 =

(k., +, - k.,)(k.2 +, - k,)flk',,,k',) x

tanh x/(2nnlAkl) ~ + (2nn2Ak2) 2 h It should be noted that col and u) 2 introduced in equation (13) is a Nlth and N2th root of unity in the complex number field, respectively. Equation (13) is a two-dimensional discrete Fourier transform and can be evaluated by the following two stages of computation: Stage 1: q'(nl,k2,t)=C ~

~

flnl,n2)co2k2,2

nz = - Nz/2 + 1

for - ~N- -2

+ 1 ~k2 ~9,

where

C 1 = -4(oAklAk

(14a)

n 2 =

-N~/2+ 1 -Nff2+l

Stage 2: q(kl,k2,t)= ei2rtk'nIXei2nk'n2Y < ,E

2

q'(nl,k2,t)col k....

n t = -N#'2 + 1

provided that:

N 2 for -~--&2+ 1 ~
(14b)

k,,<.Nk',
q(x,y,t)= - 4 ~ o

2

Na/2

~

~

hi=

n2=

f(nlAkl;n2Ak2) ei2. . . . ak, x

-Nt/2+l -NJ2+ 1

ei2~°'n'ak2"AklAk2

(12)

where Ak i = 2T/N~ for i = 1,2. This definition of Ak, implies that A x = A y = l / 2 ( l O . By introducing X = k l . A X , y = k 2 • Ay, co~ = e z"~/N,and co2 = e2'~/N~we can rewrite equation (12) as follows: N ff2

r/(k 1,k23) -~ - 4~oAk 1Ak2

N2/2

E

E

n I ~

n 2 =

-Nj2+l

f(nl,n2)colGn'co2

- N2/2+ I

for - N ff2 + 1 <~k1 <~N1/2 , - N 2 / 2 + 1 ~ k 2 <~N ff2 where

•n 2 n n l A k ~ B sin 2nn2AkzA f (nl,n2) = sl2nn~l Ak l - - . 2nnzAkE

k2n2

(13)

We observe that stage 1 and stage 2 are two onedimensional discrete Fourier transforms of](n 1,n2). If we let N~ = N 2 = 2 " where m is an integer, then the F F T algorithm can be used to compute rf(nl,k2,t ) and q(kl,k2,t) s. The number of real multiplications and additions needed to compute equations (14a) and (14b) is 4N21og2 N and 6N21ogzN, respectively. An array of 256 x 256 points was used to compute equations (14a) and (14b) (N1 = N 2 = 2 5 6 ). To avoid the spurious short-period oscillations, k should be kept small. However, a small Ak implies that the region of computation is also correspondingly limited. All numerical results presented in this paper are obtained for Ak = 0.03. This value was arrived at through a series of numerical experiments. RESULTS AND DISCUSSION Numerical results obtained from equations (13) and (14) for an impulsive bed upthrust of a rectangular block with various aspect ratios (A/B) have been obtained. The twodimensional case of the same bed motion characteristics, i.e., for A--,oc, has been studied by Hammack ¢ both theoretically and experimentally in which it was found that the theoretical results agreed well with the experimental results. As the three-dimensional experimental results are not available, it is advantageous for us to study such a case with various aspect ratios. It is reasonable to expect that for large values of A/B, our results should approach that reported by Hammack 4. The water surface profiles as a

Applied Ocean Research, 1980, Vol. 2, No. 4

167

Water waves generated by rectangular block upthrust." J. J. Lee and J. J. Chang x=O v-O

x-B v=o

1.0- . ~ y

eHAMMACK& A/B =1052

• HAMMACK& A/B = 1052

05 eHAMMACK& IB=10,5

• HAMMACK& / ~A/B=IO5 q/~o o

20 / V ~40 / A/B =2~ ' 1

1o

,

/x,

r•

A/B =105,2

110U ° ~ . , ~ ° ~ . . I 0

....

V • Hammack(1972) correspondsto A/B ~

~./;7-

Figure 2. Water surjdce profilesjbr diJJerentlength~width ratios function of the dimensionless time parameter (tx/9/h) at three different locations are shown in Fig. 2. The dimensionless parameters for this case are: texf~/B = 0.069, B/h = 12.2 (B=2.0 ft), ~0/h=0.2. The time history of bed motion is an exponential form as seen in equation (7) with a characteristic time parameter, ~ = 18.46. In Fig. 2, the time history of the water surface profiles for A/B = 10, 5, and 2 are shown for three locations (x = y = 0), (x = B, y =0), and ( x = 0 , y=A). It is seen that the curves for A/B = 10, 5 and 2 are very close to that of H a m m a c k ' s at the first two locations. For A/B--2, the three-dimensional effect becomes quite significant for tx/gfh > 20 as shown in Fig. 2. Based on the physical dimensions given, the time associated with tx/~-h = 20 is t = 1.426 sec. This is smaller than the time required for disturbance to propagate from the edge of the rectangular block (y = A) to the centre (y = 0) if one assumes the propagation speed to be xflgh = 2.3 fps. (This latter time is computed to be 1.734 sec.) It should be noted that the early arrival of end effects (with relation to the propagation speed of xflgh) as A is reduced suggests that some of the Fourier components composing the corner region of the initial disturbance have phase shifts so that they have distances less than A to travel before reaching the location (B,0). The phase shifts will be relatively smaller as A is increased. The third graph shows the water surface profile as a function of time at the edge of the major axis (x = 0, y = A). This provides another view of the three-dimensional effect. We note that for the twodimensional case, the water surface profile at x = 0, y = A, would be the same as that shown for x = 0, y = 0. However, for a finite ratio of A/B, the edge effect of the generation region is being felt instantly after the impulsive bed motion is produced. The water surface profiles as a function of time at other locations away from the region of rectangular block are shown in Fig. 3 for A/B = 5 and A/B = 2. The three locations chosen are all along the minor axis (y =0) and at x = 1.5B, 2.0B, and 2.5B. The three-dimensional effect is also quite obvious; for A/B = 2, a negative wave occurs at a smaller value of t x / g ~ . It is interesting to see some three-dimensional pictures of the water surface profile. A series of three-dimensional pictures for the case A/B--:2 is present in Figs. 4 and 5. These three-dimensional pictures were constructed from

168

Applied Ocean Research, 1980, Vol. 2, No. 4

256 x 256 samples using image processing. (Readers interested in image processing are referred to Andrew1.) They represent an overall pictorial view of the wave profile at t'x/g/h=4.20, 12.61, 21.02, and 29.43. The vertical viewing angle for these pictures is 6ff ~ and the horizontal angle is 60:. These pictures allow us to visualize the manner by which the waves are transformed from the original form of a parallelepiped into the complicated wave form at some later time. In order to demonstrate further the transformation of wave profile as a function of time, the water surface profile along a certain axis is also obtained. Figure 6 shows the water surface profile along x axis (y = 0, for x > 0) for two length/width ratios: A/B=5, 2. They are obtained by cutting the three-dimensional pictures such as that shown in Figs. 4 and 5. This series of water surface profiles shows how waves are propagated from the origin of disturbance towards the surrounding region. Using these wave profiles, we can compute the propagation speed as 2.2 ft/sec. This propagation speed is very close to the long wave celerity C = x f g h = 2 . 3 fps. From Fig. 6 we can see that the water surface profile is the same for A/B = 5 and 2 for t x / ~ - < 16.81. At a later time, the water surface profile differs significantly, especially at the tail region. However, the leading wave portion is quite similar in both cases for the time variables presented. Water surface profiles along the y axis are presented in Fig. 7 for two different values of A/B for various dimensionless time parameters, tx/g/h. The region where y = A is marked for reference. By examining this series of wave profiles, a number of observations can be made. (a) The leading wave is similar for the case of A/B = 5 and 2 at almost all the time parameters presented. (b) A reverse symmetric wave profile with respect to y = A is found at almost all the time parameters shown. (y = A corresponds to the edge of the rectangular block). (c) As waves propagate away from the generation region, interesting wave profile histories in the region y < A can be found. The leading negative wave which exists within the region y

0"5

0

t

~ ' ~ 0l ~ ~

10

x : 15B, y ~0 A/B=5

20

0"5 -

- - ~ A/

x = 2 0 B , y=O 0-5q/(0

0

-

0 5

-

A/B= 2 x = 2.5B, y=0

Q.5~

0 05-

10

20

30

40

~

/

50

%

~

A/B= 2

Figure 3. Watersurface profilesfor A/B = 5 and A / B - - 2

Water waves yenerated by rectangular block upthrust: J. J. Lee and J. J. Chan9

4.20

12.61

< A appears to be converging towards the centre first and then propagating outward (this can be clearly seen from the curves A/B=2). This further demonstrates that the waves are indeed propagating in every direction resulting in complicated wave profiles, especially when A/B is small. (d) The wave profile along the major axis (Fig. 7) and that along the minor axis differs considerably for tx/9~h greater than a certain value. (In the example of Figs. 6 and 7, t x / 9 ~ > 12.61). (e) Due to radiation of waves in every direction, significant modulation of wave profile has already taken place (as evidenced from the curves presented) even though the time parameter t x / 9 ~ is still small. This is one of the aspects that represents a significant departure from the two-dimensional results. The numerical method used for computing the waves generated through the impulsive bed upthrust of a rectangular block is much more efficient than the direct integration method. For a particular time parameter, to obtain the wave profiles for an array of 256 x 256 data points (in the x - y plane) would require only 37 sec of C P U time on a DEC-10 computer. Using these data, one can then construct the three-dimensional picture of the water surface elevation such as those presented in Figs. 4 or 5. CONCLUSIONS

Figure 4. Three-dimensional pictures showing the wave amplitude q/~o near the generation region for specified time parameter (viewing angle 0 = 60 °, &= 60°); aspect ratio A/B =2

An effective numerical method has been presented in this paper. The method is used to compute the water waves generated through the impulsive bed upthrust of a rectangular block of varying ratio of major axis/minor axis (A/B). Three-dimensional pictures showing water

= 21.02

Along x axis

x>O y=O t9~=

4"20

1.0- " ~ 5 . 2

0

I

I

~

2

3

I

BIh=12"2

[••

4

1

t~Zh= 8'41 5,2

2

I 3

I 4

t g,/~-= 12"61 1-0-

A/B= 5.2

I 1

ql{o o

2

1

B

2

4

t 9 4 ~ - 21.02

= 29.43

1,0-

t g~h" = 25'22

A/B= 5 ~ ~l~

A 1

/

B

= 2

2 3

A/B= 5

0

=

I 4

. . . . ,.I

A/B= 2

,, .......

1

2

3

i

4

t ,/~i" = 33'63

t g , ~ = 29"43

1"0-

Figure 5. Three-dimensional pictures showin9 the wave amplitude q/~o near the generation region for specified time parameter (viewin9 angle 0 = 60 °, 6 = 60°); aspect ratio A/B =2

t g , ~ = 16"81

A/B= 5.2

A/8= 5

A/B=2 i

. . . . . .

2

3

4

x/e

Figure 6. Water surface profiles for two length~width ratios (A/B = 5 and A/B = 2) at x > O, y = 0

Applied Ocean Research, 1980, Vol. 2, No. 4 169

Water waves generated by rectangular block upthrust." J. J. Lee and J. J. Chang Along y axis

x=O

B/h=122

y>O

A =5B

A-5B ,

t g~=4"20

._

.~

"..:--'~

I""..

[ % d

2

4

6

I,

A/B=2

~'~A/B~5

I

8

2

4

6

t g~g~= 12.61

A/B=2

2f: 2

4

6

"

•"-.. I 4

V

AIB=2......~

~[~

ACKNOWLEDGEMENTS t 8

t ~ [ h = 25"22

A/B=5

."/"4_

16"81

A/B- 5

6

t g,/~-= 21"02

A/B=2..~E.

8

t~ -

A/8 =5

q/{o o

,~=8.41

,/'-.

°

/ ......

t ~g~g~g~g~g~g~g=~g~29"43 g~h

The work reported herein is supported by the N a t i o n a l Science F o u n d a t i o n u n d e r G r a n t No. E N V 77-01599. The numerical c o m p u t a t i o n is d o n e at Engineering C o m p u t e r Laboratory, University of S o u t h e r n California.

A]B=5

............ ~....

I

As can be inferred from the results presented in this paper, the three-dimensional effect is quite significant, even with the assumed simple bed deformation. Therefore, in order to more fully u n d e r s t a n d the n a t u r e of t s u n a m i s due to a s u b m a r i n e earthquake, one must take into account the three-dimensional effect.

I

°

REFERENCES

t g,/~-= 33"63

1 AIB=2-~..., AIB=5~t, ~ 2 . ~4~

..:...,,.,- ~

,

AIB=2 ~.~...~:.. AIB=5.~ " '

8 ..../....... ~ , - , ~

,

2 8

3 4

rio

Figure 7. Water surface profiles for two length/width ratios (A/B = 5 and A/B = 2) at x = O, y > 0 surface elevation in the x y plane for a certain specified time parameter have been shown, F o r the cases computed, along with the time parameter presented, the leading wave shape was similar for the different ratios of A/B. However, the water surface profile of the tail region differs considerably as A / B is reduced. The wave profile along the major axis a n d that along the m i n o r axis differs considerably for tx/,q~h greater than a certain value in all cases.

170

Applied Ocean Research, 1980, Vol. 2, No. 4

5 6 7 8 9 10

Andrews,H. C. Computer Techniques in Image Processing, Academic Press, New York, 1970 Brigham, E. O. 7he Fast Fourier Transform Prentice-Hall, Englewood Cliffs, New Jersey, 1974 Driessche, P. V. D. and Braddock, R. D. On the elliptic generating region of a tsunami, J. Marine Res. 1972, p. 217 Hammack, J. L. A note on tsunamis: their generation and propagation in an ocean of uniform depth, Journal of Fluid Mechanics 1973,60, 769 {For more detail, see Report No. KH-R28, W. M. Keck Laboratory of Hydraulics and Water Resources, California Institute of Technology, Pasadena, California, 1972.) Horikawa,K. Coastal Engineering An Introduction to Ocean Engineering John Wiley, New York, 1978, Chapter 3 Kajiura, K. The leading wave of a tsunami, Bull. Earthquake Research Institute, 1963, 41, 535 Kajiura, K. Tsunami source, energy and the directivity of wave radiation, Bull. Earthquake Research Institute Tokyo University 1970, 48, 835 Oppenheim,A. V. and Schafer, R. W. Digital Signal Processing, Prentice-Hall, Englewood Cliffs, New Jersey, 1975 Reed,I. S., Truong, T. K., Kwoh, Y. S. and Hall, E. L. Image processing by transformsover a finitedfield,IEEE Transactionson Computers, 1977, C-20, (9) Van Dorn, W. G. Tsunamis, Adv. Hydrosci. 1965, 2, 1