Applied Mathematics and Computation 168 (2005) 283–304 www.elsevier.com/locate/amc
Approximate determination of the transmission and reflection coefficients for water-wave flow over a topography M.S. Abou-Dina *, F.M. Hassan Department of Mathematics, Faculty of Science, Cairo University, Giza, Egypt
Abstract The objective of this paper is to present a semi-analytical method for the calculation of the transmitted and the reflected parts of an incident harmonic progressive wave in a homogeneous fluid layer over a general topography. The domain of definition of the problem is properly divided into suitable subdomains. Elementary solutions of the problem, satisfying the field equation and a subset of the boundary conditions as well, are obtained in each subdomain by the combined use of the techniques of separation of variables and finite integral transforms. These solutions are given in the form of infinite series involving an infinite number of unknown coefficients. The series are truncated and only a finite number of the coefficients is determined so as to satisfy approximately the remaining boundary conditions. Numerical illustrations indicate that the method yields highly accurate results and that the accuracy of the procedure is improved by increasing the number of the unknown coefficients taken to account with the presence of certain optimum degrees of freedom. 2004 Elsevier Inc. All rights reserved.
*
Corresponding author. E-mail addresses:
[email protected], moustafa-aboudina@hotmail. com (M.S. AbouDina). 0096-3003/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.08.019
284
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
1. Introduction The problem of determination of the transmitted and the reflected parts of an incident wave inside a homogeneous fluid layer by a topography in the presence of a free surface attracts considerable attention, for both its theoretical and experimental aspects. This type of problems belongs to a class known as the Cauchy–Poisson Problem (CPP). In practice, CPP simulates various problems having engineering and geophysical interests such as flows over weirs, flows under gates, ship motion and fluid flows over a topography. The old work due to Thomson [1,2] contains a comprehensive discussion of this subject. From a mathematical point of view, CPP is a non-linear, mixed free boundary value problem which may be subject to certain specified initial conditions [3–6]. The main difficulties facing the theoretical investigation of the problem are the non-linearity of the system of equations and conditions governing the problem and that the free surface is apriori unknown. The geometrical complexity of the domain occupied by the fluid may also represent a major difficulty for the theoretical investigation of the problem. The difficulty concerning the non-linearity of the problem is remedied by dealing with the problem in the frames of the theory of long waves [7], the shallow water theory [8–11] or if one adopts the hypothesis of the linear theory of motion and theories of higher orders which have a vast range of application [3,7,12–14]. It is usual in the literature to overcome the difficulty arising from the geometrical complexity of the domain of the solution by using a perturbation technique and seeking approximate analytical solutions. The perturbation may be carried out around a vertical line [9–11,15–18] or around a horizontal line [19– 25] according to the geometry of the studied problem. The convergence of this perturbation technique is not proved and the application of which is, therefore, limited to cases with mild slope topographies and cases where the reflected and the transmitted waves are very small compared with the incident wave. The exact analytical solutions, although useful in analyzing the phenomena under consideration, are almost impossible to achieve once the topography deviates from being simple and in addition to the simplifications adopted in the theory, severe limitations are imposed on the geometry of the problem in order to obtain analytical solutions, such as the assumption of infinite depth [15–17] or that of very thin barriers [9–11,18,19] or very short and mild slope topographies [14,20–25]. The numerical methods are the other extreme alternative whenever the analytical approach is inadequate (see [26] and references included therein for an exhaustive review of the numerical techniques applied to CPP). These methods, however, have their own difficulties as well, represented in problems of convergence, stability and error accumulation and require high computational effort.
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
285
The so-called semi-analytical methods, which seem to be the best alternative, lie in between. They provide approximate analytical expressions for the solutions through relatively simple numerical calculations. Among these, we cite Trefftz method [27,28], the methods based on perturbation techniques [14,19,29], the boundary integral methods [30,31] and its variants, e.g. the method of fundamental solutions [32]. Each of these methods has its advantages, but also its drawbacks. The essence of almost all the semi-analytical methods is to construct certain approximating sequences of functions which yield the required solution up to a certain prescribed degree of accuracy. It is the way in which these sequences are formed that usually differentiates one method from another. It is this type of methods which we are presently interested in and the procedure proposed hereafter is inspired from noting the inaptitude encountered in dealing with other analytical or semi-analytical methods and seems to have the advantage to deal with the case of a general topography and also with the other aspects of the CPP. We limit ourselves to deal with topographies with finite support (i.e. bottoms which are horizontal except for an area of a finite horizontal extent). The domain of definition of the solution is properly divided by vertical straight lines into suitable subdomains one of which contains the irregular part of the topography. General solutions satisfying the field equation and the boundary conditions in the subdomains with regular topography were expressed in terms of a set of unknown constant coefficients and elementary functions, obtained by the method of separation of variables. The conditions on the vertical lines separating the subdomains are not considered at this stage. A justified analytical continuation, through the corrugated bottom, for the velocity potential function in the subdomain with irregular topography is carried out and a general solution, in terms of another set of arbitrary constants, is obtained in the extended subdomain using the method of finite cosine Fourier transform. This solution satisfies the field equation and the condition on the free surface and gives rise to a continuous horizontal component of the velocity on the vertical lines bounding this area. Matching the solutions in the different domains, to ensure continuity of other physical quantities (pressure, vertical component of the velocity and free surface elevations) through the vertical lines mentioned earlier, leads to expressing the solution in the whole domain of the problem in terms of one set only of unknown coefficients. This solution satisfies the field equation, everywhere in the physical domain of the problem, and all the boundary conditions except the one on the corrugated part of the bottom. Satisfaction of this latter condition leads to a series equation satisfied on this part of the topography. A truncation process is carried out, the unknown constants are determined as a solution of a finite system of linear equations and satisfactory results satisfying a certain specified error measure are obtained for the special case of a regular topography with a sinusoidal hump.
286
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
The worked example ascertains the efficiency of the proposed procedure and reveals the existence of critical value of the frequency where the reflection coefficient attains its maximum value before decreasing to zero. The reflection coefficient has been found to increase with the increase of the height of the sinusoidal hump.
2. Problem and frame of reference A homogeneous fluid layer, bounded from above by a free surface, occupies an infinite channel of finite depth. The impermeable bottom of the channel, bounding the fluid layer from bellow, is horizontal except for an area of finite horizontal extent where it possesses a general topography (Fig. 1). An incident harmonic wave (with frequency x), propagating inside the fluid layer, is partially obstructed by the bottom topography. Transmitted and reflected waves going towards the downstream and the upstream extremities of the channel, respectively, and local perturbations disappearing far from the corrugated part of the bottom, are therefore generated. The problem is to determine the transmission and the reflection coefficients and also the perturbations occurring in the fluid motion due to the irregularity of the topography. For simplicity of the mathematical treatment, the study is carried out in two dimensions within the frame of the linear theory of motion, the fluid is assumed incompressible and inviscid, the motion is irrotational and steady harmonic case is considered. The rectangular frame of reference O(x, y), shown in Fig. 1, is used for the description of the motion. The origin O is located in the mean level of the free surface at the beginning of the region lying over the non-horizontal part of the bottom, the x-axis is taken in the direction of the incident wave, pointing to the downstream extremity of the channel and the y-axis is directed vertically upwards.
Fig. 1. Problem and frame of reference.
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
287
3. Equations of motion and boundary conditions The upper bound of the domain occupied by the fluid is replaced,within the frame of the linear theory, by the horizontal line y = 0. For a harmonic irrotational fluid motion, the time dependent velocity potential /(x, y, t) is related to a time independent one U(x, y) by /ðx; y; tÞ ¼ Re½Uðx; yÞeixt ;
ð1Þ
where x is the frequency of the motion and U(x, y) is in general a complex function of x and y only, subject to the following system of equations and conditions [5]: (i) In the fluid mass (1 < x < 1, h + k(x) 6 y 6 0): o 2 U o2 U þ 2 ¼ 0; ox2 oy
ð2Þ
where y = h + k(x) is the equation of the bottom. (ii) At the upper bound of the domain of the problem (y = 0): oU x2 U ¼ 0: oy g
ð3Þ
(iii) On the impermeable bottom (y = h + k(x), 1 < x < 1): oU ¼ 0; on
ð4Þ
where ono means differentiation along the direction normal to the bottom at the considered point. (iv) The radiation condition at infinity (jxj ! 1) stating that: no wave is coming from infinity except the prescribed incident one. The time dependent free surface elevation g*(x, t) and pressure P*(x, y, t), take respectively the forms: g ðx; tÞ ¼ Re½gðxÞeixt
ð5Þ
P ðx; y; tÞ ¼ qgy þ Re½P ðx; yÞeixt ;
ð6Þ
and
where g(x) and P(x, y) are, in general, two complex functions given in terms of the time independent velocity potential U(x, y) as: ix Uðx; yÞ at y ¼ 0 ð7Þ gðxÞ ¼ g
288
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
and P ðx; yÞ ¼ iqx Uðx; yÞ: ð8Þ * The time dependent stream function W (x, y, t) is expressed in terms of the time independent one W(x, y) as W ðx; y; tÞ ¼ Re½Wðx; yÞeixt :
ð9Þ
The functions U(x, y) and W(x, y) are related by the Cauchy–Riemann conditions in the form grad U ¼ k ^ grad W:
ð10Þ
4. The perturbation technique The main difficulty facing the analytical treatment of the system of equations and conditions (2)–(4) together with the radiation condition (iv) at infinity (jxj ! 1), is that condition (4) on the irregular bottom of the channel (y = h + k(x)) is not taken along a coordinate curve of the used system of coordinates (i.e. x = const. or y = const.). To overcome this difficulty, it is usual in the literature to assume the function k(x), describing the irregularity of the topography, to be of the order of a small parameter e in the form (see [20–25] and references included therein) kðxÞ ¼ eh1 ðxÞ:
ð11Þ
The solution of the system of equations and conditions (2)–(4) is then sought through the expansion 1 X Uðx; yÞ ¼ U0 ðx; yÞ þ en Un ðx; yÞ; ð12Þ n¼1
where U0(x, y) is the velocity potential of the incident wave and Un(x, y), n = 1, 2, . . . represents the contribution of order n, due to the irregularity of the topography, in the total velocity potential U(x, y). Eqs. (2) and (3) imply the function Un(x, y) to satisfy relations of the form o2 Un o2 Un þ ¼ 0; ox2 oy 2
at 1 < x < 1; h 6 y 6 0
ð13Þ
and oUn x2 Un ¼ 0; at y ¼ 0: ð14Þ oy g The expansion of condition (4), on the irregular bottom about the coordinate curve y = h, replaces this condition by oUn dQn ¼ ; oy dx
at y ¼ h;
ð15aÞ
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
289
where 2 oU1 1 2 o U0 þ ½h1 ðxÞ ; 2 ox oxoy
Q1 ðxÞ ¼ h1 ðxÞ
oU0 ; ox
Q3 ðxÞ ¼ h1 ðxÞ
2 3 oU2 1 1 2 o U1 3 o U0 þ ½h1 ðxÞ þ ½h1 ðxÞ ;... 2 ox oxoy 6 oxo2 y
Q2 ðxÞ ¼ h1 ðxÞ
at y ¼ h: ð15bÞ
The problem reduces thus to find solutions for Eq. (13), with n = 1, 2, 3, . . ., satisfying conditions (14), (15a,b) taken along the coordinate curves y = 0 and y = h in addition to the radiation condition at infinity. Following this procedure, the first order (truncation at n = 1) problem has been solved in [20] for random h1(x) and in [21] for arbitrary h1(x); see also [22–24]. A sequential procedure for the determination of the function Un(x, y), starting from a straight-crested incident wave in water of uniform depth for U0(x, y), has been given in [25]. In this latter publication third order solution (n = 3) was expressed and the results were compared with laboratory measurements due to [21]. The perturbation technique, briefly outlined in the present section, proved efficiency in predicting the phenomenon in certain special cases (see [21,25]). However, the convergence of expression (12) is not studied and the expansion of condition (4) about y = h, leading to formulae (15a) and (15b), is valid only within its radius of convergence, which is also unknown. In addition, as it is clear from expressions (15a), the form of the function Qn for n = 3 is complicated and the complexity increases with the increase of n, hence for cases where the convergence needs to consider large number of functions Un(x, y), it shall not be easy to deal with the expressions of the necessary functions. Further the value of the small parameter e, introduced in (11) and intervening in (12), is not identified with the parameters of the problem and the technique requires both the reflected and the transmitted waves to be very small compared with the incident wave [21]. These drawbacks put major limitations on both the vertical extent and the slope of the irregular topography in the cases studied following this perturbation technique. In the remaining part of the present paper, we shall introduce a procedure, thought to be original, for the study of the problem in the case of a general topography.
5. Procedure for the solution The region with non-horizontal bottom is enclosed between the y-axis and the straight line at x = a. These two vertical lines divide the domain of definition of the problem into three subdomains: V lying on the left to the y-axis,
290
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
Fig. 2. The physical and the extended subdomains.
V+ lying on the right to x = a, and V0 lying in between (Fig. 2). A formal solution of the problem inside each one of these subdomains will be obtained in terms of certain unknown coefficients. This solution is built so as to satisfy the field equation and a subset of the imposed boundary conditions. Matching the solutions in the different domains, to ensure continuity of some of the physical quantities (velocity, pressure, and free surface elevation) in crossing the vertical lines at x = 0 and x = a will serve to get relations between the expressions of the solution in the adjacent subdomains. A complete determination of the unknown coefficients is then performed so as to satisfy the remaining part of the boundary and continuity conditions in an approximate manner. 5.1. Solution in the regions with horizontal bottoms Both of the semi-infinite domains V and V + have the advantage of being bounded by straight lines which are natural boundaries in the used system of coordinates. This permits to apply the usual technique of separation of variables for the solution of the problem in these two regions. Following this procedure and denoting by k0 the only real positive root of the transcendental equation k tanh kh ¼
x2 g
ð16Þ
and by kp (p = 1, 2, 3, . . .) the positive roots of the transcendental equation k tan kh ¼
x2 ; g
ð17Þ
arranged in the order of increasing magnitude, the expressions of the velocity potential function in the regions V and V+ satisfying: (i) the field equation in the fluid mass, (ii) the condition on the free surface, (iii) the condition on the
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
291
bottom and (iv) the radiation condition at infinity, are found to take the forms [3–5]: 1 X Rp cos kp ðy þ hÞekp x U ðx; yÞ ¼ I 0 eik0 x þ R0 eik0 x cosh k0 ðy þ hÞ þ p¼1
ð18Þ and Uþ ðx; yÞ ¼ T 0 cosh k0 ðy þ hÞeik0 ðxaÞ þ
1 X
T p cos kp ðy þ hÞekp ðxaÞ ;
ð19Þ
p¼1
respectively, where I0 is a prescribed complex quantity characterizing the incident wave, R0, T0 are the reflection and transmission coefficients and Rp, Tp with p = 1, 2, . . . are complex constant coefficients. In this solution, one readily recognizes two essential types of the fluid motion possessing distinct features: the first one represents progressive waves propagating towards one of the channel extremities (incident, reflected and transmitted waves) and the other one, written in the form of infinite series, represents local perturbations vanishing far from the corrugated part of the bottom. The reflection and transmission coefficients R0 and T0 and the complex coefficients Rp and Tp are to be determined in the sequel. 5.2. Solution in the region over the topography The intermediate region V0 lying over the topography (Fig. 2) is bounded from its left by the y-axis (x = 0), from its right by the vertical straight line x = a, from below by the corrugated part of the bottom (y = h + k(x), 0 6 x 6 a) and from above by the horizontal line y = 0. The velocity potential function U0(x, y) defined in this region is harmonic in the interior of the domain V0 and extends to be a continuous function on its closure V 0 with vanishing harmonic conjugate on the corrugated bottom (since this latter is a stream line). These conditions are sufficient for the function U0(x, y) to extend to a harmonic function on the rectangular region OABC of Fig. 2, which we denote by V [33]. Without any confusion, the analytic continuation function, defined in the rectangular region V, will be also given the same notation U0(x, y). The equations and conditions satisfied by the extended velocity potential U0(x, y) in its domain of definition V are given as (I) In the interior of the domain V o2 U0 o2 U0 þ ¼ 0: ox2 oy 2
ð20Þ
292
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
(II) On the upper bound of the domain V(y = 0, 0 6 x 6 + a), the free surface condition reads oU0 x2 0 U ¼ 0: oy g
ð21Þ
(III) On the left bound of the domain V(x = 0, h 6 y 6 0), the continuity of the pressure and the velocity of the fluid along this boundary together with expressions (8) and (18) lead to U0 ð0; yÞ ¼ ðI 0 þ R0 Þ cosh k0 ðy þ hÞ þ
1 X
Rp cos kp ðy þ hÞ
ð22Þ
p¼1
and 1 X oU0 ð0; yÞ ¼ ik0 ðI 0 R0 Þ cosh k0 ðy þ hÞ þ kp Rp cos kp ðy þ hÞ: ox p¼1
ð23Þ
(IV) On the right bound of the domain V(x = a, h 6 y 6 0), the continuity of the same physical quantities as in (III), together with expressions (8) and (19) lead to U0 ða; yÞ ¼ T 0 cosh k0 ðy þ hÞ þ
1 X
T p cos kp ðy þ hÞ
ð24Þ
p¼1
and 1 X oU0 ða; yÞ ¼ ik0 T 0 cosh k0 ðy þ hÞ kp T p cos kp ðy þ hÞ: ox p¼1
ð25Þ
The continuity of the free surface elevations at x = 0 and at x = a is guaranteed by the satisfaction of conditions (22) and (24) with the help of expression (7). For the determination of the function U0(x, y) solving the boundary value problem, given by Eqs. (20), (21), (23) and (25), we introduce the finite cosine Fourier transform of this function defined as [34] f0 ðyÞ ¼ U m
Z 0
þa
p U0 ðx; yÞ cos m x dx; a
m ¼ 0; 1; 2; . . .
ð26Þ
Performing this transformation to the system of equations and conditions (20), f0 ðyÞ, m = 0, 1, 2, . . . are found to satisfy a (21), (23) and (25), the functions U m certain system of ordinary differential equations with some specified boundary conditions. The solution of this latter for m = 0, 1, 2, . . . takes the form
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
mp ax2 mp aI 0 f 0 y þ y sinh cosh Um ðyÞ ¼ am a a mpg 2 d0m k0 iðR0 I 0 þ ð1Þm T 0 Þ 2 2 cosh k0 ðy þ hÞ k0 mp a 1 X k p ðRp þ ð1Þm T p Þ 2 2 cos kp ðy þ hÞ; kp þ mp p¼1 a
293
ð27Þ
where am, m = 0, 1, 2, . . . are constant unknown coefficients to be determined. The Inversion formula for the finite Fourier transform (27) reads [34] U0 ðx; yÞ ¼
1 p X 2 d0m f0 Um ðyÞ cos m x : a a m¼0
ð28Þ
Applying this formula to the transformations (27) and using the expressions
1 mp X 2 d0m ð1Þm k0 1 x ¼ cosðk0 xÞ cos ð29Þ
2 2 mp a sinðk0 aÞ a k m¼0 0
a
and
1 mp X 2 d0m ð1Þm kp 1 x ¼ coshðkp xÞ; cos
2 2 mp a sinhðkp aÞ a kp þ a m¼0
ð30Þ
one obtains the function U0(x, y), representing the velocity potential in the region V0 lying over the topography, in the form 1 mp ax2 mp p X y þ y cos m x sinh am cosh U0 ðx; yÞ ¼ I 0 a a a mpg m¼0 ifðR0 I 0 Þ cos k0 ðx aÞ þ T 0 cosðk0 xÞg
1 X p¼1
fRp cosh kp ðx aÞ þ T p coshðkp xÞg
cosh k0 ðy þ hÞ sinðk0 aÞ
cos kp ðy þ hÞ : sinhðkp aÞ
ð31Þ
The set of complex (unknown) coefficients am, m = 0, 1, 2, . . . may now be given its meaning by noting that these constants are related to the Fourier coefficients of the vertical derivative of the extended function U0(x, y) calculated at the lower bound of the rectangle V (i.e. at y = h, 0 6 x 6 +a) by the relation 2 1 X oU0 x mph mp mph sinh ðx; hÞ ¼ I 0 cosh am cosðmpx=aÞ: a a a oy g m¼0 ð32Þ
294
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
Inserting expression (31) for the function U0(x, y) into conditions (22) and (24), confirming the continuity of the pressure, the vertical component of the velocity and the free surface elevations at x = 0 and at x = a respectively, one obtains the two following relations 1 X ðI 0 þ R0 Þ cosh k0 ðy þ hÞ þ Rp cos kp ðy þ hÞ ¼ I0
1 X
am
m¼0
p¼1
mp ax2 mp y þ y sinh cosh a a mpg
ifðR0 I 0 Þ cos k0 a þ T 0 g
cosh k0 ðy þ hÞ sinðk0 aÞ
1 X cos kp ðy þ hÞ fRp cosh kp a þ T p g sinhðkp aÞ p¼1
ð33Þ
and T 0 cosh k0 ðy þ hÞ þ
1 X
T p cos kp ðy þ hÞ
p¼1
¼ I0
1 X
m
am ð1Þ
m¼0
mp ax2 mp y þ y sinh cosh a a mpg
ifðR0 I 0 Þ þ T 0 cos k0 ag
cosh k0 ðy þ hÞ sinðk0 aÞ
1 X cos kp ðy þ hÞ : fRp þ T p cosh kp ag sinhðkp aÞ p¼1
ð34Þ
Using relations (33) and (34) together with the completeness and the orthogonality properties, over the interval h 6 y 6 0, of the system of functions {cosh k0(y + h), cos kp(y + h), p = 1, 2, 3, . . .}, the reflection and transmission coefficients R0 and T0 and also the coefficients of the local perturbations Rp and Tp, p = 1, 2, 3, . . . can be explicitly expressed in terms of the constant coefficients am, m = 0, 1, 2, . . . in the form R0 ¼
1 I0 X m f1 ð1Þ eik0 a ga0;m am ; 2 m¼0
" T 0 ¼ I0 e
Rp ¼
ik0 a
# 1 ik a 1X m e 0 ð1Þ a0;m am ; 2 m¼0
1 I0 X m 1 ð1Þ ekp a ap;m am 2 m¼0
ð35Þ
ð36Þ
ð37Þ
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
295
and Tp ¼
1 kp a I0 X m e ð1Þ ap;m am ; 2 m¼0
ð38Þ
where the coefficients a0,m and ap,m are defined for p = 1, 2, . . . and m = 0, 1, 2, . . . as h mph x2 mph i 4k0 mp sinh cosh a a a g
2 2 a0;m ¼ ð39Þ ½2k0 h þ sinhð2k0 hÞ ma2p k20 and ap;m
h x2 i g cosh mph 4kp mp sinh mph a a a i : ¼ hm2 p2 2kp h þ sinð2kp hÞ a2 þ k2p
ð40Þ
Expressions (35)–(38) enable to express the function U0(x, y) given by (31) in terms of the unknown constant coefficients am, m = 0, 1, 2, 3, . . . in the form " # 1 X U0 ðx; yÞ ¼ I 0 eik0 x cosh k0 ðy þ hÞ þ am U0m ðx; yÞ ; ð41Þ m¼0
U0m ðx; yÞ,
where m = 0, 1, 2,. . . are well determined functions of the variables x and y defined as mp ax2 mp p y þ y cos m x U0m ðx; yÞ ¼ cosh sinh a a a mpg 1 a0;m eik0 x þ ð1Þm eik0 ðxaÞ cosh k0 ðy þ hÞ 2 1 1X m ap;m ekp x þ ð1Þ ekp ðxaÞ cos kp ðy þ hÞ: ð42Þ 2 p¼1 The corresponding stream-functions W(x, y), W0(x, y) and W+(x, y) are given in the regions V, V0 and V+ using relations (10), (18), (19) and (31) as W ðx; yÞ ¼ ifI 0 eik0 x R0 eik0 x g sinh k0 ðy þ hÞ þ
1 X
Rp sin kp ðy þ hÞekp x ;
p¼1
ð43Þ " 0
ik0 x
W ðx; yÞ ¼ I 0 ie
sinh k0 ðy þ hÞ
1 X m¼0
# am W0m ðx; yÞ
ð44Þ
296
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
and Wþ ðx; yÞ ¼ iT 0 sinh k0 ðy þ hÞeik0 ðxaÞ
1 X
T p sin kp ðy þ hÞekp ðxaÞ ;
ð45Þ
p¼1
where R0, T0, Rp and Tp, p = 1, 2, 3, . . . are given in terms of the coefficients am, m = 0, 1, 2, . . . by expressions (35)–(38) and W0m ðx; yÞ, m = 1, 2, . . . are functions of the variables x and y defined as mp ax2 mp p W0m ðx; yÞ ¼ sinh y þ y sin m x cosh a a a mpg i m þ a0;m feik0 x ð1Þ eik0 ðxaÞ g sinh k0 ðy þ hÞ 2 1 1X ap;m fekp x ð1Þm ekp ðxaÞ g sin kp ðy þ hÞ: 2 p¼1
ð46Þ
In the particular case where the corrugation of the bottom is limited to a thin vertical submerged plate at x = 0, it can be easily shown that the conservation of mass requires the satisfaction of a relation of the form T 0 þ R0 ¼ I 0 :
ð47Þ
In fact relation (47), which translates a conservation law, must be satisfied whatever the particular form of the bottom topography. It seems, in view of expressions (35) and (36), that the verification of such a relation for an arbitrary width a of the subdomain V0 of Fig. 2, needs to take into consideration a very large number of the coefficients am, m = 1, 2, 3, . . . However, relation (47) turns to an obvious identity, if we limit ourselves to the choice a ¼ 2kp=k0 ;
ð48Þ
where k is a positive integer. To optimize the numerical computations, the positive integer k of relation (48) can be taken as the smallest one permitting the subdomain V0 to contain all the corrugated part of the topography and, if necessary, another part of the flat bottom. This integral value of k will be denoted hereafter by k0. The coefficient a0,m of expression (39) can be seen by passing to the limit to be equal to unity for m = 2k0. 5.3. Complete determination of the solution The velocity potential function U(x, y) and the stream-function W(x, y), describing the motion of the fluid in the whole domain of definition of the problem are given explicitly, in terms of the set of unknown constant coefficients am, m = 0, 1, 2, 3, . . ., by holding together expressions (18), (19)
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
297
and (41) for U(x, y) and (43)–(45) for W(x, y). These functions give rise to continuous physical quantities in the fluid mass and satisfy the field equation and all the boundary conditions of the problem except the one on the irregular part of the bottom. The satisfaction of this latter condition will be used hereafter for the determination of the set of unknown coefficients. The stream function W(x, y), given by (43)–(45), should have a vanishing value on the irregular part of the topography since it is a continuous function and vanishes on the horizontal part of the bottom as may be verified by expressions (43) and (45). This, together with expression (44) show that the set of unknown coefficients am, m = 0, 1, 2, 3, . . . has to satisfy a series equation of the form 1 X
am W0m ðx; h þ kðxÞÞ ¼ ieik0 x sinhðk0 kðxÞÞ;
0 6 x 6 a;
ð49Þ
m¼0
where the functions W0m ðx; yÞ are given by (46) for m = 1, 2, . . . It is worth noting that the analysis carried out, up to this stage of calculations, is exact and that no approximation process has taken place. For the determination of the unknown coefficients of expression (49), we propose to follow a slightly modified version of the well known collocation technique [35]. We truncate the infinite series on the LHS of (49) by replacing the infinite upper limit of the summation therein by a positive integer M, chosen sufficiently large so as to satisfy the accuracy requirements, and to choose a suitable partition {xn, n = 0, 1, 2, . . ., N} for the interval 0 6 x 6 a, so that 0 = x0 < x1 < x2 < < xN = a, where N (number of collocating points) is a suitably chosen integer with N > M, then we verify the satisfaction of the series equation (49) at the points xn of the partition to obtain an over determined finite system of linear algebraic equations in the unknown coefficients am, m = 0, 1, 2, . . ., M of the form M X
Bn;m aNm ¼ C n ;
n ¼ 1; 2; . . . ; N ;
ð50Þ
m¼0
where the coefficients Bn,m and Cn are given for n, m = 0, 1, 2, . . ., M as Bn;m ¼ W0m ðxn ; h þ kðxn ÞÞ
ð51Þ
C n ¼ ieik0 xn sinhðk0 kðxn ÞÞ:
ð52Þ
and
The vector A of unknown coefficients aNm , m = 0, 1, 2, . . ., M is then obtained as the solution of the square system of linear equations ðBT BÞA ¼ BT C;
ð53Þ
298
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
where BT is the transpose of the matrix B of coefficients Bn,m and C is the vector of constants Cn, for n = 0, 1, 2, . . ., N and m = 0, 1, 2, . . ., M. The number M + 1, of the considered coefficient, is often called degree of freedom. Having solved the square system (53) for the coefficients aNm , the streamfunction W0(x, y) is then approximated by an expression of the form " # M X ik0 x 0ðN ;MÞ N 0 W ðx; yÞ ¼ I 0 ie sinh k0 ðy þ hÞ am Wm ðx; yÞ : ð54Þ m¼0
The values of this stream function on the topography, which have to vanish, represent a sort of local error distribution ER(x), 0 6 x 6 a, for the method ERðxÞ ¼ W0ðN ;MÞ ðx; h þ kðxÞÞ
ð55Þ
and one may test the validity of the chosen values of N and M by checking an inequality for the maximum local error of the form ðN ;MÞ
EL
¼ max jERðxÞj < e 06x6a
ð56Þ
for a certain prefixed e. The validity of the chosen values of N and M may alternatively be tested by checking another inequality, based on the mean global error, presented in the form Z 1 a ðN ;MÞ EG ¼ jERðxÞjdx < e: ð57Þ a 0 The idea of the procedure is to start with a choice of N and M, reasonable with the width a of the subdomain V0, then if (56) (or (57)) is not met, the values of N and M are increased gradually and the procedure is repeated till satisfaction of the accuracy requirement. By the choice N = M, the method reduces to the so called straight forward collocation technique. Approximate expressions for the reflection and transmission coefficients are finally given in terms of the coefficients aNm , m = 0, 1, 2, . . ., M as ðN ;MÞ
R0
¼
M I0 X 1 ð1Þm eik0 a a0;m aNm 2 m¼0
ð58Þ
and " ðN ;MÞ T0
¼ I0
# M 1X m ik0 a N 1 1 ð1Þ e a0;m am ; 2 m¼0
respectively. The coefficients a0,m, m = 1, 2, . . ., M being given by (39).
ð59Þ
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
299
6. Numerical results for a sinusoidal hump The practical use and the efficiency of the procedure explained above are tested for the case of a topography limited to a horizontal bottom with a sinusoidal hump of height H and width W. The function k(x) is taken as 1 2px 1 cos kðxÞ ¼ HU ðW xÞ for 0 6 x 6 a; ð60Þ 2 W where U(z) denotes the Heaviside step function,which equals unity if its variable z is non-negative and vanishes otherwise. According to the above procedure and to optimize the calculations, k0 is chosen to be the smallest integer satisfying W 6 a = 2k0p/k0. 6.1. Accuracy and efficiency of the method ðN ;MÞ
Table 1 exhibits the maximum local errors EL occurring for certain values of the numbers M and N. The results in this table have been obtained for the pffiffiffiffiffiffiffiffi particular choice x ¼ 0:548 g=h for the frequency, H/h = 0.25 and W/h = 10.0 for the height and the width of the hump respectively. For this choice the value of k0h is 0.576 and it is sufficient to take k0 = 1 which yields a/h = 2k0p/(k0h) = 10.907 > W/h as necessary. The collocating points {xn, n = 0, 1, 2, . . ., N} of the partition of the interval 0 6 x 6 a are taken equidistant. The results in Table 1 show that, the accuracy is in general satisfactory and that the maximum local error decreases for the same value of M with the increase of N. The number M of the considered coefficients can not be increased without limit since there is an optimum degree of freedom (M + 1 = 261 for the present application). The best result in the table occurs for M = 260 and N = 300. Table 1 ðN;MÞ The errors EL for the topography of Eq. (60) M 10 20 30 40 50 100 150 200 250 260 270 300
N=M 3
1.101 · 10 2.938 · 104 1.345 · 104 7.578 · 105 4.927 · 105 1.258 · 105 5.521 · 106 3.354 · 106 1.812 · 104 6.484 · 104 0.011 0.022
N = M + 20
N = M + 40
4
6.575 · 104 1.850 · 104 7.927 · 105 4.563 · 105 3.002 · 105 8.117 · 106 3.885 · 106 2.251 · 106 1.489 · 106 1.397 · 106 1.126 · 105 4.591 · 103
6.582 · 10 1.895 · 104 8.246 · 105 4.882 · 105 3.270 · 105 9.198 · 106 4.410 · 106 2.609 · 106 1.642 · 106 1.536 · 106 6.960 · 106 1.051 · 103
300
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
Fig. 3. Absolute error for M = 200 and N = 400.
Parallel computations have shown that the corresponding global error ðN ;MÞ ðN ;MÞ EG is at least one order better than the maximum local error EL given in Table 1. This is because the local error ER(x) is always localized in the neighborhood of the end points of the topography at x = 0 and at x = W. Fig. 3 exhibits the distribution of the absolute local error jER(x)j along the interval 0 6 x 6 a for the choice M = 200 and N = 400. The maximum local ðN ;MÞ ðN ;MÞ error for this choice is EL ¼ 1:710 106 and the global error is EG ¼ 8 3:635 10 . The smoothness of the boundary of the topography plays an important rule in the rapidity of attaining a high and satisfactory accuracy. The topography represented by the sinusoidal hump of Eq. (60) is continuous at the end points x = 0 and x = W with continuous first order derivative. The second order derivative of this topography is discontinuous at both of the points x = 0 and x = W. To detect the influence of the smoothness of the topography on the accuracy of the results, we slightly modify the above sinusoidal hump by replacing it by another one, possessing a continuous second order derivative, of the form 2 1 2px k 1 ðxÞ ¼ HU ðW xÞ 1 cos for 0 6 x 6 a: ð61Þ 2 W ðN ;MÞ
Table 2 shows the the maximum local errors EL corresponding to this latter choice of the topography (61) for the same particular values of the parameters as those considered for Table 1. We see by the aid of Table 2 that for sufficiently smooth topographies the proposed method is a very power tool for the solution of the problem. Using twenty coefficients leads to an accuracy of the order of 106. The optimum
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
301
Table 2 ðN;MÞ The errors EL for the topography of Eq. (61) M
N=M
N = M + 20
N = M + 40
10 20 30 40 50 100 150 200 250 260 270 300
1.122 · 103 7.080 · 106 1.476 · 106 4.625 · 107 1.933 · 107 1.226 · 108 2.415 · 109 1.935 · 108 7.148 · 105 1.413 · 103 1.683 · 103 8.103 · 103
5.953 · 104 4.576 · 106 7.606 · 107 2.430 · 107 1.055 · 107 7.371 · 109 1.605 · 109 5.489 · 1010 2.330 · 1010 2.045 · 1010 3.124 · 1010 8.011 · 109
5.953 · 104 4.540 · 106 7.489 · 107 2.342 · 107 1.000 · 107 6.455 · 109 1.343 · 109 4.630 · 1010 1.896 · 1010 1.691 · 1010 3.718 · 1010 8.871 · 109
choice, occurring here also at M = 260 and N = 300, gives rise to an accuracy of the order of 1010 which is more than necessary in all practical applications. In all cases the results obtained following the straight forward collocation technique, for which N = M, are considerably improved by increasing the value of the number N of the collocating points specially for large numbers of considered coefficients M > 150. 6.2. The effect of the frequency The results tabulated in Table 3 exhibit the influence of the change in the frequency x on the absolute value of the relative reflection coefficient jR0/I0j for the topography limited to a sinusoidal hump given by Eq. (60) with particular
Table 3 The variation of the ratio jR0/I0j with the frequency x x2h/g
jR0/I0j
x2h/g
jR0/I0j
0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.15
0.012 0.021 0.029 0.039 0.050 0.064 0.078 0.093 0.106 0.115 0.119 0.119
1.20 1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00 3.00 4.00 5.00
0.119 0.116 0.112 0.108 0.104 0.101 0.097 0.093 0.088 0.039 0.013 0.004
302
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
Table 4 The variation of the ratio jR0/I0j with the relative height H/h H/h
jR0/I0j
0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
0.022 0.047 0.078 0.115 0.162 0.226 0.319 0.478 0.893
values of the parameters W/h = 4.0 and H/h = 0.25 within a maximum local ðN ;MÞ error EL of the order of 105. The results in Table 3 show, for the considered model, that the absolute value jR0/I0j increases with the increase of the frequency x to attain a maximum value (=0.119) at a critical frequency xc (with x2c h=g ¼ 1:15) before decreasing to zero. Additional calculations show that, here and throughout the paper, the values of the transmission and the reflection coefficients T0 and R0 are independent of the particular choice of the integer k (Pk0) intervening in the definition of the real positive parameter a by relation (48). 6.3. The effect of the height of hump The influence of the variation of the maximum relative height of the topography (H/h) on the absolute value of the relative reflection coefficient jR0/I0j is tabulated in Table 4 for the particular values of the parameters W/h = 10 ðN ;MÞ and x2h/g = 0.1 within a maximum local error EL not exceeding the order 5 of 10 . The results of Table 4 show that the absolute value jR0/I0j increases monotonically with the increase of the relative maximum height of the topography H/h.
7. Conclusions The main purpose of the present work is to present a simple semi-analytical method for the study of the problem of the propagation of waves over a general topography. The obtained solution can be used to check the validity of other numerical or analytical solutions proposed for the same problem. It avoids the drawbacks of the perturbation technique used in [20–25] and puts
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
303
no limitations on the vertical extent or on the slope of the topography. This solution satisfies, identically, the field equation in the bulk and the conditions on free surface, on the flat part of the bottom and at infinity. The condition on the irregular topography is satisfied approximately within an appropriate accuracy. To the best of our knowledge, the present method is original and it can be applied for the study of several boundary value problems with physical interest of the CPP type. The method has shown, through the worked applications, to be easy to handle and needs a few amount of numerical calculations that can be carried out using any available software. The method possesses a high degree of accuracy and is self adjusting, in the sense that the error admitted can be estimated and is controlled by choosing suitable degrees of freedom. Further, the proposed method proves a perfect efficiency in dealing with problems with smooth bottoms. Problems with topographies having corner points are expected to require a larger number of coefficients of the series representing the solution. This may be avoided by smoothing the corner points of the topography to the desired degree. The study of the limiting case where the topography reduces to a vertical partially submerged thin barrier, for which smoothing is meaningless, needs a special care and is in progress. References [1] W. Thomson, On stationary waves in flowing water, Philos. Mag. 5 (XXII) (1886) 353, 454, 516. [2] W. Thomson, On stationary waves in flowing water, Philos. Mag. 5 (XXIII) (1887). [3] H. Lamb, Hydrodynamics, Sixth ed., Dover, New York, 1932. [4] J.J. Stoker, Water Waves, Interscience Pub., New York, 1957. [5] J.V. Wehausen, E.V. Laitone, in: Surface Waves Handbuch der Physik, vol. 9, Springer, Berlin, 1966. [6] P.H. Le Blonde, L.A. Mysak, Waves in Ocean, Elsevier, Amsterdam, 1978. [7] G.B. Whitham, Linear and Nonlinear Waves, Interscience Pub., New York, 1974. [8] J.P. Germain, The´orie ge´ne´rale des mouvements dÕun fluide parfait pe´sant en eau peu profonde de profondeur constante, in: 1000, C.R.A.S. Paris T. 274 (20 mars) (1972) 997. [9] M.S. Abou-Dina, M.A. Helal, The influence of a submerged obstacle on an incident wave in stratified shallow water, Eur. J. Mech. (B) 9 (6) (1990) 545–564. [10] M.S. Abou-Dina, M.A. Helal, The effect of a fixed barrier on an incident progressive wave in shallow water, IL Nuovo Cimento (B) 107 (3) (1992) 331–344. [11] M.S. Abou-Dina, M.A. Helal, The effect of a fixed submerged obstacle on an incident wave in stratified shallow water (Mathematical Aspects.), IL Nuovo Cimento (B) 110 (8) (1995) 927–942. [12] L. Skjelbreja, J.A. Hendrickson, Fifth order gravity wave theory, in: Proc. of 7th Conf. on Coas. Eng., vol. 1, pqrt 10, 1961, pp. 184–196. [13] C.M. Lee, The second order theory of heaving cylinders in a free surface, J. Ship Res. 12 (1968) 313–327.
304
M.S. Abou-Dina, F.M. Hassan / Appl. Math. Comput. 168 (2005) 283–304
[14] M.S. Abou-Dina, Nonlinear transient gravity waves due to an initial free-surface elevation over a topography, J. Comput. Appl. Math. 130 (2001) 173–195. [15] D.V. Evans, A note on the waves produced by small oscillations of a partially immersed vertical plate, J. Inst. Math. Appl. 17 (1976) 135–140. [16] F. Ursell, The effect of a fixed vertical barrier on surface waves in deep water, Proc. Cambridge Philos. Soc. 43 (1947) 374–382. [17] P.F. Rhodes-Robinson, On the forced surface waves due to a partially immersed vertical wavemaker in water of infinite depth, Quart. Appl. Math. LIV (4) (1996) 709–719. [18] E.O. Tuck, Shallow-water flows past slender bodies, J. Fluid Mech. 26 (1966) 81–95. [19] T.F. Ogilvie, Propagation of waves over an obstacle in water of finite depth, Ph.D thesis, Berkeley, California, USA, 1960. [20] R.B. Long, Scattering of surface waves by an irregular bottom, J. Geophys. Res 78 (1973) 7861–7870. [21] A.G. Davies, A.D. Heathershaw, Surface wave propagation over sinusoidally varying topography, J. Fluid. Mech. 144 (1984) 419–443. [22] C.C. Mei, Resonant reflection of surface water waves by periodic sand-bars, J. Fluid. Mech. 152 (1985) 315–335. [23] J.T. Kirby, A general wave equation for waves over rippled beds, J. Fluid. Mech. 162 (1986) 171–186. [24] J. Miles, P. Chamberlain, Topographical scattering of gravity waves, J. Fluid. Mech. 361 (1998) 175–188. [25] J. Miles, On gravity-wave scattering by non-secular changes in depth, J. Fluid. Mech. 376 (1998) 53–60. [26] R.W. Yeung, Numerical methods in free-surface flows, Ann. Rev. Fluid Mech. 14 (1982) 395– 442. [27] I. Herrera, in: Trefftz Method, in: C.A. Brebbia (Ed.), Topics in Boundary Element Research, vol. 1, Springer-Verlag, Berlin, 1984. [28] A.P. Zielinski, I. Herrera, Trefftz Method: fitting boundary conditions, Int. J. Numer. Methods Eng. 24 (1987) 871–891. [29] A.H. Nayfeh, Introduction to Perturbation Techniques, Interscience Pub., New York, 1981. [30] M.S. Abou-Dina, M.A. Helal, Reduction for the nonlinear problem of fluid waves to a system of integro-differential equations with an oceanographical application, J. Comput. Appl. Math. 95 (1998) 65–81. [31] M.S. Abou-Dina, M.A. Helal, Boundary integral method applied to the transient nonlinear wave propagation in fluid with initial free surface elevation, J. Appl. Math. Modell. 24 (8–9) (2000) 535–549. [32] G. Fairweather, A. Karageorghis, The method of fundamental solutions for elliptic boundary value problems, Adv. Comput. Math. 9 (1998) 69–95. [33] S. Lang, Complex Analysis, third ed., Springer-Verlag, New York, 1993. [34] C.J. Tranter, Integral Transforms in Mathematical Physics, Third Ed., Methuen Co. Ltd, London, 1966. [35] J.A. Kolodziej, Review of application of boundary collocation methods in mechanics of continuous media, SM Arch. 12 (4) (1987) 187–231.