Eur. Phys. J. E (2015) 38: 45 DOI 10.1140/epje/i2015-15045-0

THE EUROPEAN PHYSICAL JOURNAL E

Regular Article

Numerical study of laminar, standing hydraulic jumps in a planar geometry Ratul Dasgupta1 , Gaurav Tomar2,a , and Rama Govindarajan3 1 2 3

Department of Chemical Engineering, Indian Institute of Technology, Mumbai 400 076, India Department of Mechanical Engineering, Indian Institute of Science, Bangalore 560 012, India TIFR Centre for Interdisciplinary Sciences, 21 Brundavan Colony, Narsingi, Hyderabad 500 075, India Received 8 June 2014 and Received in final form 28 January 2015 c EDP Sciences / Societ` Published online: 25 May 2015 –  a Italiana di Fisica / Springer-Verlag 2015 Abstract. We solve the two-dimensional, planar Navier-Stokes equations to simulate a laminar, standing hydraulic jump using a Volume-of-Fluid method. The geometry downstream of the jump has been designed to be similar to experimental conditions by including a pit at the edge of the platform over which liquid film flows. We obtain jumps with and without separation. Increasing the inlet Froude number pushes the jump downstream and makes the slope of the jump weaker, consistent with experimental observations of circular jumps, and decreasing the Reynolds number brings the jump upstream while making it steeper. We study the effect of the length of the domain and that of a downstream obstacle on the structure and location of the jump. The transient flow which leads to a final steady jump is described for the first time to our knowledge. In the moderate Reynolds number regime, we obtain steady undular jumps with a separated bubble underneath the first few undulations. Interestingly, surface tension leads to shortening of wavelength of these undulations. We show that the undulations can be explained using the inviscid theory of Benjamin and Lighthill (Proc. R. Soc. London, Ser. A, 1954). We hope this new finding will motivate experimental verification.

1 Introduction A layer of fluid flowing horizontally over a solid boundary is frequently seen to display a sudden increase in height. This phenomenon is known as hydraulic jump. Rivers, canals and estuaries are geophysical examples where such jumps takes place, often being forced by downstream conditions such as presence of obstacles or abrupt constrictions, etc. Hydraulic jumps are also commonly experienced in industrial applications and manufacturing processes. For example, it is used for dissipating energy in water flowing over hydraulic structures. In industries hydraulic jump may be utilized to make mixing of chemicals more efficient. It also finds applications in water purification processes [1]. Such hydraulic jumps in planar geometries were probably well known even before they found mention by Leonardo Da Vinci [2]. They occur in myriad forms —laminar or turbulent, travelling (bores) or standing, undular or regular step, etc. Rayleigh [3] treated the hydraulic jump as a shock and based on continuity of mass and momentum flux under the inviscid approximation, derived the post jump height and velocity. He further showed that loss of energy always a

e-mail: [email protected]

accompanies the jump, ΔE = Q

g(h2 − h1 )3 , 4h2 h1

(1)

where ΔE is the energy loss across the jump, Q is the volume flow rate, g is the gravitational acceleration, and h1 and h2 are the uniform depths upstream and downstream of the jump, respectively. Rayleigh attributed this energy loss to turbulence and other viscous losses. Viscous forces were however not important in the situations such as tidal bores, which Rayleigh intended to study. For hydraulic jumps occurring on thin films (e.g. a circular hydraulic jump), it was recognized by Tani [4] that a finite separated bubble exists underneath the jump, to be consistent with which he included viscosity in the steady shallowwater equation. This boundary layer shallow-water equation (BLSWE) was vertically averaged, and using the assumption of a self-similar velocity profile, an ordinary differential equation was derived. Bohr et al. [5], obtained a scaling relation for the radius of the circular jump from an analysis of the depthaveraged viscous shallow-water equations used earlier by Tani [4]. A similar approach was also used by Singha et al. [6] for a planar geometry. The solution of these

Page 2 of 14

depth-averaged viscous equations with a self-similar velocity profile work quite well upstream and downstream of jumps, but do not resolve the near-jump structure. In this approximation, the solutions display multi-valuedness and the jump is represented as a discontinuity connecting outer and inner solutions. The location of the jump is fixed by using the Rayleigh criterion to obtain the complete downstream and upstream profile connected through a discontinuity. Watanabe et al. [7] proposed an averaged equation model for simulating hydraulic jumps using polynomial forms for velocity profiles thus relaxing the selfsimilarity assumption by introducing a shape factor that allows separation. They demonstrated that the model can be integrated through separation and can describe a typeI jump very well in the upstream and downstream regions while showing small deviations from experimental data in the near jump region. This model was subsequently employed by Bonn et al. [8] to study hydraulic jumps in a channel. They showed that at high Reynolds number the upstream flow does not agree with Watson’s solution, h Re = 1.81 [9], where h is the slope of depth in the flow direction and Re = Q/ν is the Reynolds number with ν as the kinematic viscosity. Suggesting that the effective viscosity increases due to turbulent eddy viscosity at high Reynolds number, they used a mixing-length model of the turbulence to reproduce the experimental results. Bush and Aristoff [10] modified Watson’s analysis [9] for circular hydraulic jumps by including surface tension effects and showed that steady circular jumps do not exist beyond a certain value of surface tension. The downstream effects were ignored by using the assumption of a flat interface profile beyond the jump region. They concluded that although surface tension improves the agreement between the theoretical and experimental results, the persistent difference between the predictions could be due to the effects of downstream boundary conditions. Kasimov [11] built upon the work of Bush and Aristoff [10] with the focus on deriving a base-state for circular jump on plates of finite size which could then be subjected to stability analysis. Laminar standing jumps have also been studied using interactive boundary layer theory. Gajjar and Smith [12] proposed a two layer model with a uniform velocity profile over a viscous sublayer at the wall. They showed that the inviscid-viscous free interaction (or branching) occurs for Froude number (F r) greater than one. This theoretical model was compared with the experiments of Craik et al. [13] by Brotherton-Ratcliff and Smith [14]. They showed that the upstream depth predicted from the inviscid theory is one-fifth of the actual depth measured in the experiments of Craik et al. [13]. Bowles and Smith [15] argued that additional effects should be added to the above theory. Essentially the effect of viscosity and surface tension were included. They showed that the freeinteraction between the surface tension and viscosity upstream of the jump leads to formation of wave-like structures whereas the downstream is determined by the freeinteraction between the gravitational pressure gradient and viscosity. They noted that although the high Reynolds

Eur. Phys. J. E (2015) 38: 45

number theory is successful in describing the separation point in the jump, complete description requires explicit incorporation of the effect of downstream boundary condition. Higuera [16] performed finite difference simulations of boundary layer equations along with the kinematic condition for the evolution of the depth. Considering a twodimensional flow over the flat finite plate, he performed asymptotic expansions around the plate edge to obtain the boundary conditions required for model closure. He also studied the effects of Froude number, surface tension and streamline curvature on the jump location and structure. The analysis was carried out at high Reynolds number (∼ 103 ) assuming the applicability of boundary layer approximation in the whole domain. Direct numerical simulations have been reported in the literature for turbulent planar jumps [17, 18] as well as for laminar jumps in cylindrical configurations [19–21]. Pritchard et al. [22] reported finite element numerical simulations, under the lubrication approximation, of film flow over an inclined surface with perturbations and observed the formation of hydraulic jumps. Given the approximation, their Reynolds number of investigation was restricted to < 25. In an earlier study of film flow on an inclined plane, we predicted the formation of hydraulic jumps below a critical Reynolds number which depends on the angle of inclination [23]. Wols [24] studied the formation of weir-assisted undular jumps for high Reynolds number (∼ 12000) flows. Undulations in height were observed downstream of the weirs. Further, their experimental results were used to evaluate the performance of various one-dimensional models, such as linear and non-linear Boussinesq models and a non-hydrostatic D model, in capturing the undular jumps. Ohtsu et al. [25] experimentally studied the formation of undular hydraulic jumps on horizontal platforms, again at high Reynolds numbers (∼ 6 × 104 ). In particular, they showed the variation in the wavelength and amplitude of the undulations with Froude number, and also identified the critical Froude number below which undular jumps may be expected. They showed that undular jumps with breaking waves are observed when the Reynolds number is decreased. Undular hydraulic jumps, although theoretically possible, have not been reported, to our knowledge, for moderate Reynolds numbers. In the present work, we perform numerical simulations of planar hydraulic jumps using Volume of Fluid method in the moderate Reynolds (∼ O(102 )) and Froude (∼ O(10)) number regime. In particular, we investigate the effect of downstream boundary conditions and flow parameters, namely, Reynolds, Froude and Weber numbers, on the jump location and structure. In the simulations presented here, natural boundary conditions at the platform edge are ensured by allowing the liquid film to fall off from the finite length platform into a pit which is also included in the computational domain. Interestingly, in this regime of Re and F r, we obtain undular hydraulic jumps with wavelengths of the order of depth of the film. Undular hydraulic jumps have been studied earlier using various variants of KdV equation at high Re in the narrow range around the critical F r [26–31]. We show that the hydraulic jumps in the moderate Reynolds

Eur. Phys. J. E (2015) 38: 45

Page 3 of 14

range can also be explained using the inviscid theory of Benjamin and Lighthill [26]. The paper is organized as follows: we lay down the problem formulation and details of numerical method employed in sect. 2, results and discussions for various parametric studies mentioned above are presented in sect. 3 and finally we present conclusions in sect. 4.

be written as [33]   u − un + un+ 12 · ∇un+ 12 = ρn+ 12 Δt   ∇ · μn+ 12 (Dn + D ) + (σκδs n)n+ 12 , cn+ 12 − cn− 12

+ ∇ · (cn un ) = 0,

(7)

Δt ∇pn+ 12 , ρn+ 12

(8)

Δt un+1 = u −

2 Problem formulation

∇ · un+1 = 0.

2.1 Governing equations The governing equation for flow is the two-phase NavierStokes equation given by  ρi

 ∂u + ∇ · uu = −∇p + ∇ · (μi D) + ρi g, ∂t

(2)

where subscript i = l, g stands for liquid and air, respectively. Here u is the velocity field, p is the pressure and D = (∇u+∇uT ) is the deformation rate tensor. The gravitational acceleration is given by g. The fluids are assumed to be incompressible and therefore satisfy the continuity equation ∇ · u = 0. (3) The normal stress balance at the liquid-air interface is given by [pi ] = [n · μi D · n] + [ρi g] · n + σκ,

(4)

where [ · ] denotes the jump in the quantity “ · ” at the liquid-air interface, n is the normal to the interface, σ is the surface tension coefficient and κ is the curvature at the interface. The tangential stress balance yields [t · μi D · n] = [ρi g] · t.

(5)

The liquid-solid interface is considered to be a no-slip surface and an obstacle is placed downstream to assist jump formation. In what follows we discuss the numerical algorithm used to solve the two-phase flow in hydraulic jumps.

2.2 Numerical method Numerical simulations of planar hydraulic jumps are performed using an open source code, GERRIS [32, 33]. Gerris is a Navier-Stokes solver augmented with the volume of fluid algorithm for two phase flows. We describe briefly the algorithm it employs. A second-order accurate staggered time discretisation has been used for the velocity, volume-fraction/density and pressure fields. Using a timesplitting projection method the discretized equations can

(6)

(9)

Here c is the void-fraction field, defined as the ratio of the volume of liquid in a computational cell to the volume of the cell itself. By definition, c is zero for the “gas” or outer fluid, unity for the “liquid” or the fluid of interest, and takes a value between 0 and 1 in the computational cell containing the interface. The subscripts indicate the time step at which the variables have been evaluated. Equation (6) is the discretized Navier-Stokes equation with the terms on the left-hand side representing temporal and convective terms and the terms on the right-hand side are descretized viscous forces and surface tension. The interfacial boundary conditions (eqs. (4) and (5)) are thus implicitly imposed by the two-phase volume-of-fluid algorithm. In the projection scheme used, the pressure terms are not included in the momentum equation but are used for correcting the velocity field, as shown in eq. (8), so that the resultant velocity field is solenoidal, i.e., satisfies eq. (9). Equation (7) is the advection equation for the void-fraction field. The subscript ∗ indicates the auxiliary time step. Convective terms are evaluated at the fractional time step (n + 1/2) using the Godunov procedure [34]. Using eqs. (8) and (9) the Poisson equation governing the pressure field can be written as  Δt ∇pn+ 12 = ∇ · u . (10) ∇· ρn+ 12 This equation for the pressure and the discretized momentum conservation equation (eq. (6)) are solved efficiently using a quad/oct tree-based multigrid solver with an underlying linear (Gauss-Seidel) solver [32]. The viscous terms are discretized using a second-order accurate unconditionally stable Crank-Nicholson scheme. The velocity, pressure and void-fraction are all collocated at the center of the computational finite-volume cell. The volume fraction field is advected using an operator-split algorithm and the velocity field obtained above [35]. The volume flux for advection is computed geometrically to avoid numerical diffusion. The liquid-gas interface is reconstructed at each time step and the interface normal is computed using a mixed-Youngs-Centered (MYC) method [36]. The surface tension force (σκδs n)n+ 12 is calculated using a balanced-force surface tension calculation [37]. A secondorder accurate estimate of the curvature is obtained using the Height-Function technique. The quad/oct tree mesh in Gerris allows efficient mesh refinement and adaptation

Page 4 of 14

Eur. Phys. J. E (2015) 38: 45 S6 Hydraulic Jump S1

Air

Obstacle

g

S3’

S5 Liquid

H

S2 Lo L

Ho

S3 S4

Fig. 1. A schematic of the computational domain. Fluid enters the domain through the side labelled S1 and exits through S4. Note an obstacle of height Ho is placed in the downstream of the jump. By including the “pit” at the edge of the platform (S2), we ensure that downstream boundary conditions for the liquid film flow are appropriately modeled.

Table 1. Boundary conditions for computational domain in fig. 1. Side On velocity S1

On pressure

Free-slip, no penetration Neumann condition except for inlet jet

S2

No-slip, no penetration

Neumann condition

S3

Free-slip, no penetration Neumann condition

S4

Neumann

S5

Free-slip, no penetration Neumann condition

S6

Free-slip, no penetration Neumann condition

Outlet condition, Dirichlet

procedure. We use adaptive mesh refinement in our simulations by defining the vorticity and the gradient of the void-fraction variable as the cost function. A permissible value of 0.01 is chosen, and when the cost function exceeds this limit, the maximum defined refinement is used in that neighborhood. Similarly, where the cost function lies below a lower limit, the mesh of minimum refinement is used. The open source code (Gerris) can be easily modified (see Tomar et al. [38]). Details of the implementation are available in [32, 33]. 2.3 Computational domain Figure 1 shows the computational domain under investigation for simulating hydraulic jumps. The hashed sides represent walls where a no-penetration condition has been imposed. The boundary conditions are given in table 1. In some of our simulations, we place a thin rectangular obstacle of height Ho on the platform at a distance Lo from the inlet (see wall S2 in fig. 1), although this is not a necessary requirement for obtaining jumps as we show later in the next section. Pressure is prescribed at the outflow

boundary where for velocity a Neumann boundary condition is used. In cases involving gravitational acceleration and two-phase flow, the imposition of a constant pressure field at the exit or even a constant hydrostatic field may lead to fictitious upstream flow of the liquid from the outer boundary. Thus, to avoid imposing an artificial condition on S3 , we introduce a “pit” (see fig. 1) and impose the classical outflow boundary condition (Dirichlet in pressure and Neumann for the velocity field) along the horizontal boundary S4 where a constant pressure can be imposed. The introduction of this boundary condition eliminates cases of upstream flow and is able to give grid independent steady-state results. For the simulations shown, the relevant parameters are the height Ho of the obstacle and film thickness H at the inlet. The other dimensional parameters in the problem are the average velocity Uav at the inlet, the length of the domain L, the length L0 from the inlet to the obstacle, kinematic viscosities of the inlet and the ambient fluids, νw and νa respectively, and their respective densities ρw and ρa . The non-dimensional ratios that characterize flow behavior are the Reynolds number Re, the inlet Froude number F ri , the domain length L∗ = L/H, the obstacle location and height Lo /H and Ho /H, respectively, viscosity ratio νa /νw and density ratio ρa /ρw . The Reynolds and inlet Froude √numbers are defined as Re ≡ Uav H/ν and F ri = Uav / gH, respectively. To study the effect of surface tension on hydraulic jump, we define Weber 2 H/γ, where γ is the surface tennumber as W e = ρw Uav sion. In the simulations that are presented here, we use νa /νw = 0.01 and ρa /ρw = 0.01. We argue that due to high density and viscosity ratio for the given geometry and nature of flow, the ambient fluid would have minimal effect on the details of hydraulic jump. To verify this assumption, later we show that the upstream viscous flow in our study is in good agreement with previous free-surface assumption based theoretical results of Watson [9] and Dasgupta and Govindarajan [39] (also see [23]). Grid size used for all the results presented here is L∗ /Δx ∼ 2100. The initial condition for all these simulations comprise of a layer of pre-wetting fluid on the platform. As the simulation is started, the fluid layer due to its initial velocity, flows over the obstacle and into the pit.

Eur. Phys. J. E (2015) 38: 45

Page 5 of 14

6

6

5

y/H

4

4 t/T = 1750 t/T = 1500 t/T = 1000 t/T = 500 t/T = 250 t/T = 125

2 0 0

0.2 0.4 0.6 0.8

x/L

1

Fig. 2. (Color online) Transients showing the formation of an hydraulic jump in time (t/T ; where T = H/Uav ) for Re = 125, F ri = 6 and H0 /H = 1. The length of the platform is L = 300H.

3 Results and discussions The steady-state numerical results presented in this study have been obtained by simulating the complete transients of the hydraulic jump formation on a platform of finite length (see fig. 1). The time scale, T = H/Uav , chosen for non-dimensionalisation is the time that it takes for a fluid parcel to cover the inlet distance H. The transients are depicted in fig. 2. The vertical axis has been scaled with the thickness of the inlet flow (H) and the horizontal axis has been scaled with the distance of the obstacle from the inlet (Lo ∼ 296H). The jump originates from the downstream as an upstream travelling wave-like structure and travels upstream as its front progressively steepens. The undulations start appearing one by one downstream of the jump as the entire jump structure moves towards its upstream steady-state location and reaches its final shape. Finally in the steady state (t/T = 1750), one can see a stationary jump with a series of undulations downstream. All the steady-state results reported in the following sections were obtained by running the simulations for over 10 times the average time it takes for a fluid parcel to traverse from the inlet to the end of the domain (t/T ∼ 3000–4000). As a validation of our numerical procedure, we compare results that we obtain from our numerical solution of the full Navier-Stokes equations to the solutions of the Boundary Layer Shallow-Water Equations (BLSWE) [39]. These equations have been earlier cast in a form suitable for use in this problem in ref. [39]. We found two solutions to these equations, one which displays positive h at any F r, termed the P solution, and, at low Froude numbers, an N solution with negative slope. There is also a similarity solution which is possible at large Froude numbers and is attained far upstream of the jump (h Re = 1.8139) [9]. In addition, another similarity solution is possible at low Froude numbers downstream of the jump [39]. It is shown in fig. 3 that upstream and downstream of the jump, the

y/H 3 2

Navier-Stokes simulations N solution P solution

1 0

0

100

x/H

200

300

Fig. 3. (Color online) Comparison of the simulated height profile with the P and N solution of the BLSWE [39]. Re = 125, F ri = 6, L∗ = 300 and the inlet profile is specified as parabolic. Note that there are no fitting parameters except matching the height at the starting location.

P and N solutions of the BLSWE respectively describe the actual height profiles well. 3.1 Effect of change of geometry In this section, we examine the effect of change of geometrical parameters like the length of the platform over which jump occurs, presence of an obstacle and its height. The structure of a typical jump that we obtain in our simulations is shown in fig. 4. This simulation corresponds to an inlet Froude number of F ri = 6, Re = 125 and without any obstacle. It is interesting to observe that underneath the first few undulations one can locate a steady separated bubble (see zoomed view in fig. 5). To the best of our knowledge, such a flow structure has not been reported before (see [40] for experimental evidence of laboratory scale undular bores). As the amplitude of undulation decays, the size and strength of separation bubbles under undulation crests decrease until attached flow is attained further downstream. The mechanism of jump formation captured through complete transient analysis in our simulations (see fig. 2) suggests that the formation of separation bubbles is preceded by the formation of the undulations and their upstream migration (not shown in the figure). Therefore, the presence of separation bubbles is an effect of the jump formation and not its cause. This mechanism of jump formation has been discussed earlier (see [7, 13]) and mathematical models have been proposed that show formation of jumps without separation [7]. Further, for shorter platforms we obtain jumps without any undulations or separation bubbles. Figure 6(a) shows formation of a weak jump over a platform of length L∗ = 150 for Re = 125 and F ri = 6. Figure 6(b) shows that the velocity profile in the near jump region x/H ∼ 110 (jump region). In figs. 6(a) and (b), we see that (laminar standing) jumps without separation can exist. Our simulations reinforce the idea that separation is an effect and not the cause of hydraulic jumps. Note that the inlet Reynolds

Eur. Phys. J. E (2015) 38: 45

y/H

Page 6 of 14

5 0

0

100

x/H

200

300

Fig. 4. Undular hydraulic jump for Re = 125, F r = 6, L∗ = 290, H0 /H = 0.

y/H

4 2 0 80

90

100

110 x/H 120

130

140

Fig. 5. Zoomed-in view of fig. 4. Note the decrease in the size of the vortices.

y/H

4 2 0 0

50

x/H

100

150

(a)

1

u/umax

0.8 0.6 0.4 0.2 0.2

0.4

0.6

y/ymax

0.8

1

(b) Fig. 6. a) A hydraulic jump without a separated bubble for Re = 125, F r = 6, L∗ = 150, H0 /H = 0 b) The vertical variation of the horizontal velocity profile at x/H = 120 (in the jump region). Note that there is no negative (reverse) velocity (i.e. no flow separation) as is also suggested by the streamlines in fig. 6(a).

and Froude numbers are the same for both the simulations shown in figs. 4 and 6(a): what differs between the two are the downstream boundary conditions viz. the length of the domain. Thus, whether the flow separates underneath a jump seems to be determined by the boundary conditions downstream of the jump rather than those upstream [15]. Interestingly, the jump location for the shorter platform

case (L/H = 150) shifts downstream in comparison to jump on L/H = 300 platform. This is possibly due to the downstream sudden vertical acceleration near the plate edge [16]. To study the effect of an upstream obstacle on the jump structure, we performed simulations with an obstacle of the height of the inlet thickness, H0 /H = 1, at x/H = 96 for platform lengths of L∗ = 150 and 300. Figure 7 shows that the obstacle has little effect on the jump formation for longer platform of L∗ = 300 (also see fig. 8), whereas for L∗ = 150 it forces the formation of a stronger jump located further upstream in comparison to the no-obstacle case shown in fig. 6(a). Placement of an obstacle before the natural (or unforced) jump location leads to rapid loss of energy thus resulting in early formation of the jump. In contrast, if the obstacle is placed substantially downstream of the natural jump location, both momentum and energy have decayed substantially thus suggesting a weaker influence of the obstacle on jump structure. We define a “natural jump” as one where the local Froude number goes through unity purely due to a viscous slowing down mechanism and not because of any downstream forcing. To investigate the above phenomenon further, we study the effect of change of the obstacle location on the jump position. Figure 8 shows the interface profile for Re = 125 and F ri = 6 and a range of obstacle locations ranging from Lo /H = 65 to Lo /H = 296, where the platform length is L∗ = 300. For obstacles located downstream of the natural jump location (see curve in black corresponding to the no obstacle case), the location of the obstacle has a very small effect on the jump location and structure. In fig. 8, such a “natural jump” happens at ≈ x/H = 80 as seen for the case without the obstacle. Placing a downstream obstacle at any location x/H > 80

Eur. Phys. J. E (2015) 38: 45

Page 7 of 14

8

y/H

6 4 Lo /H = 96, L/H = 150 Lo /H = 96, L/H = 300

2 0 0

50

100

150

x/H

200

250

300

Fig. 7. Effect of an obstacle (H0 /H = 1) on hydraulic jump for different lengths of platform (L∗ = 300 and 150) for Re = 125, F r = 6. Obstacle is placed at Lo /H = 96.

8

y/H

6 4

Lo /H = 296 Lo /H = 196 Lo /H = 146 Lo /H = 96 Lo /H = 65 No Obstacle

2 0 0

50

100

150

x/H

200

250

300

Fig. 8. (Color online) Effect of the location of the obstacle (Ho /H = 1) on hydraulic jump for Re = 125, F r = 6 and L/H = 300. The curve in black corresponds to the no-obstacle case.

and varying this location by almost 300% has a small effect of pushing the jump upstream by ≈ 10%. However, once an obstacle is placed upstream of the natural jump length (i.e. Lo /H = 65 in fig. 8), there is an O(1) effect on the location of the jump. Now the jump is pushed upstream by as much as 30% for a 30% change in the obstacle location. In this case, the Froude number is forced to go through unity due to the obstacle and not just due to viscous slowing down. Nevertheless, the downstream structure of the jump, as the fluid film falls off the platform into the pit, is affected by an obstacle placed near the edge in response to the streamline curvature (see the curve in red in fig. 8). For the other obstacle locations, including the case with Lo /H = 65, the negative slope in the post jump region matches well with the no-obstacle case. Figure 9 shows the effect of change in obstacle height on hydraulic jump formation. The height of the obstacle, unlike its location (cf. fig. 8), retains a strong influence on the location of the “natural jump”, being able to push the jump significantly upstream as the height doubles. In addition, this parameter also has a strong influence on the flow upstream of the jump. Figure 10 shows variation in

the scaled slope of the interface (h Re) with local Froude number (F r). Initially the slope varies in agreement with Watson’s solution (h Re = 1.8139) and deviates from it indicating the formation of the jump. Inset in fig. 10 shows the zoom in on the curve corresponding to the no obstacle case. The encircled markers 1 and 2 mark the loops in the curve corresponding to the first and second undulations just downstream of the jump, respectively. Initially, the size of the loops grow suggesting an increase in the steepness of the undulations as the center of the loops travels leftwards towards lower local F r. Subsequently, the loops hit the lowest F r and further downstream the mean profile has a negative slope. With the decrease in the mean height the local F r increases and eventually the curves monotonically decrease (increasingly negative h Re) corresponding to the flow down the edge of the platform. It is seen that when there is no obstacle, the similarity solution of [9] is followed up to a local Froude number of ≈ 1.9 whereas with an increase in the obstacle height, the departure from Watson’s similarity solution is expected to occur earlier. Interestingly, the curves for the no-obstacle (Ho /H = 0) and Ho /H = 1 cases overlap and deviate only

y/H

Page 8 of 14

Eur. Phys. J. E (2015) 38: 45

5

Ho /H = 0 Ho/H = 1 Ho/H = 2

0 0

50

100

150 x/H

200

250

300

Fig. 9. Effect of increasing the height of the obstacle on the jump for Re = 125, F r = 6 and L∗ = 300. The lowest interface has no obstacle H0 /H = 0, the dashed interface is for H0 /H = 1 and the topmost interface is for H0 /H = 2.

20

h’Re

0

10 2

1

0

-20

10

0.5

0.75

Ho /H = 0 Ho /H = 1 Ho /H = 2

-40 0

obstacle height indefinitely. As seen in fig. 11, increasing the obstacle height eventually removes the possibility of reaching a steady-state under the given boundary conditions. The undulations develop increasingly steeper fronts after which they break and this process does not seem to reach any steady state. This is in agreement with the observations in experiments performed by controlling the downstream depth of the film [40]. In what follows, we discuss the effect of flow parameters, namely, Reynolds, Froude and Weber (surface tension) on the structure of the planar hydraulic jump.

Watson’s solution: h’Re = 1.8139

0.5

1

Fr

1.5

2

3.2 Effect of change of flow parameters

2.5

Fig. 10. (Color online) The right end of the curves indicates an upstream location and the down-stream coordinate increases following the spiral from right to left. The horizontal line represents the solution of [9]: h Re = 1.8139. The red curve follows this solution until the local Froude is around 1.5 and then starts differing visibly. For larger obstacle heights, this deviation starts much earlier in F r.

towards the end, whereas, for Ho /H = 2 the curve deviates from Watson’s solution earlier for a higher F r = 2.1. The corresponding effect on the jump location is substantial due to the small slope for Re = 125 (see fig. 9). Due to the presence of the downstream obstacle the flow occurs over the obstacle thus forcing an increase in the downstream depth of the flow √ and consequently resulting in higher phase speed (∼ gh) of gravity waves at the obstacle. For Ho /H = 1, this shift in downstream height is smaller compared to the depth of the hydraulic jump in the no-obstacle case (see fig. 9). Therefore, the effect of the downstream obstacle on the location of the jump is minimal. Whereas, for Ho /H = 2 the downstream depth of the film is increased over that in the no-obstacle case and thus results in upstream shift in the jump location (∼ F r = 2.1). It is not always possible to obtain a steady laminar jump of increasing strength just by increasing the

In various earlier studies employing BLSWE for simulating hydraulic jumps, pressure is assumed to be hydrostatic. Pressure contours from direct numerical simulation for film flow with Re = 125 and F ri = 6 shown in the fig. 12 suggest that the hydrostatic approximation for pressure is even valid for undular jumps. The pressure contours follow the interface undulations and the slope of the increase in pressure across the film is constant and proportional to the gravitational acceleration (p = ρl g(h − y)). Figure 13 shows the effect of Reynolds number on the structure and location of the hydraulic jump. The Froude number, F ri = 6, is same for all the cases. For low Re = 62.5 the jump undergoes steepening of the waves, due to substantial viscous losses, and subsequently resulting into breaking of the waves with the jump location significantly upstream compared to that for the Re = 125 case. Whereas, for high Re = 250, the jump is weak and substantially downstream. We argue, based on our earlier observations of results for Re = 125 and shorter platforms (L∗ = 150; see fig. 7), that even for Re = 250, a stronger jump would have formed had the platform been long enough. The weaker jump for Re = 250 is a result of the vertical acceleration of the fluid into the pit. We suggest that for even higher Re the jump formation would depend on the downstream available length. Nevertheless, with increase in Re a downstream shift in the jump location is expected. Interestingly, there’s a qualitative similarity between our observations for moderate Re and those of [25] at high

y/H

Eur. Phys. J. E (2015) 38: 45

Page 9 of 14

40 20 0

0

50

100

150

200

x/H

250

300

Fig. 11. Breaking undulations due to a tall obstacle (H0 /H = 4) for Re = 125, F ri = 6 and L∗ = 300.

P

y/H

5 0 0

100

x/H

200

300

5.000 4.000 3.000 2.000 1.000 0.100

Fig. 12. (Color online) Pressure contours for no-obstacle case with Re = 125 and F r = 6.

15

Re = 62.5 Re = 125 Re = 250

y/H

10 5 0 0

50

100

150

x/H

200

250

300

Fig. 13. Effect of change of Reynolds number on the structure of the hydraulic jump.

y/H

5 0 0

50

100

Fr = 3 Fr = 6 Fr = 10 Fr = 14 150 x/H 200

250

300

Fig. 14. (Color online) Effect of change of Froude number on the structure of the hydraulic jump.

Page 10 of 14

Eur. Phys. J. E (2015) 38: 45

∞, surface energy has been completely neglected and thus the whole of the excess energy is convected downstream in waves. For finite surface tension values of W e = 50 and 10, some energy is accounted for in the generation of the extra surface of the undulations and thus we observe a rapid decay in the undulation amplitudes (see fig. 15). In the following section, we analyze the undular structure of the jump based on the KdV analysis of Benjamin and Lighthill [26]. 3.3 Analysis of undulations

Fig. 15. (Color online) Effect of surface tension on the structure of the hydraulic jump. Note the decrease in the amplitude and wavelength of the undulations with a decrease in Weber number.

Re ∼ 12000 (also see [24]). As shown in fig. 13, with a change from Re = 125 to Re = 62.5, the non-dimensional amplitude (difference between heights of crest and trough divided by the upstream depth; following the same definition of [25]), increases as Reynolds number goes down; a feature also reported in ref. [25]. We now discuss the role of Froude number on the formation of hydraulic jump. Figure 14 shows interface profiles for Froude numbers, F ri = 3, 6, 10 and 14. The simulations have been performed for the no-obstacle case with Re = 125 and L/H = 300. For lower F ri = 3, jump formation occurs substantially downstream and undulations occur with smaller wavelength in comparison to the F ri = 6 case. We note that for F ri = 3, although the slope of the jump is more in comparison to higher F ri , the overall increase in height of the jump is smaller. This is expected as an increase in F ri can be seen as a relative increase in the strength of gravity, thus a small increase in height corresponds to a relatively higher increase in the potential energy of the flow. For larger F ri = 10, the jump location shifts downstream and a higher jump, but with a weaker slope, is formed. Similarly, for F r = 14, a jump is expected further downstream. Since the length of the platform is not enough for a “natural jump” formation, a rather weak jump forms similar to the case for Re = 250. Similar observations have been reported for circular hydraulic jumps [41]. Surface tension has been shown to play a curious role in the jump structure (see [16]). Figure 15 shows interface 2 h/σ = 10, profiles for three Weber numbers, W e = ρl Uav 50 and ∞. Interestingly, the wavelength of the undulations reduce with decrease in W e (increase in surface tension) contrary to the usual expectation of longer wavelength for surface tension dominant flows. We show in sect. 3.3 that these undulations arise due to the insufficient loss of energy at the jump location which results in dispersion of excess energy in the form of standing wave [26]. For W e =

The formation of undulations downstream of hydraulic jumps is well-known [40] for travelling and standing jumps. The earliest experimental data on travelling jumps or bores [42] suggested that undulations appear only for weak bores. Approximating the jumps as a discontinuity and assuming pressure to be hydrostatic, Rayleigh [3] obtained the following relation between the ratio of (uniform) heights upstream and downstream of the jump H2 /H1 and the upstream Froude number F r1

−1 + 1 + 8F r12 H2 . (11) = H1 2 Experimental data on standing jumps reported that for a predominantly parallel incoming flow [43], post-jump undulations appeared only when the incoming Froude (F r1 ) was not too large compared to unity. The traditional definition of the “strength” of the jump is the ratio H2 /H1 and this in eq. (11) is directly proportional to F r1 . Thus one argues that undulations are formed only for “weak jumps” (i.e. when F r1 or equivalently H2 /H1 is not too large compared to one; empirical ratios of around F r1 ∼ 1.21 are reported in experiments [42, 44]). The observation [42] that the undulations once formed do not change shape led Keulegan and Patterson [45] to suppose that this wave-train was cnoidal as it is the only periodic waveform which could remain stationary in a non-linear dispersive medium. The cnoidal wave forms are stationary solutions of the KdV equation. This idea was furthered in ref. [26] conjecturing that three quantities, namely, the mass, momentum and energy flux determine this cnoidal wave-train uniquely as they appear as coefficients of the KdV equation which governs these undulations. In an inviscid, irrotational flow without any height discontinuities, these three parameters are continuous at any streamwise location. It was shown in ref. [3] that this would not be not true if there was a heightdiscontinuity in the flow. If we insist on making Q and S continuous, then R becomes discontinuous at the location of the discontinuity, the incoming energy-flux being greater than the outgoing when the discontinuity is chosen to be an abrupt increase in the flow direction. One could argue that this discontinuity in energy flux is an artifact of representing the jump as a discontinuity. Benjamin and Lighthill [26] treated the jump as a discontinuity and showed that for an incoming supercritical flow, with no change in mass Q and momentum S and

Eur. Phys. J. E (2015) 38: 45

Page 11 of 14

2

D

0

Re = 125 Re = 180 Re = 250

-2 -4 50

100

150

200

250

x/H

Fig. 16. (Color online) Variation in D with streamwise direction for Re = 125, 180 and 250.

an energy loss ΔR given by eq. (1), no oscillations were possible: there could only be a super-critical to subcritical transition through the discontinuity. However, if there occurred some laminar energy dissipation whose magnitude was less than the value given by eq. (1), oscillations of the cnoidal kind could occur. The amplitude of the resultant wave-train is directly proportional to the strength of the jump [26], thus for strong jumps (H2 /H1 1) large amplitude cnoidal waves appear which are possibly unstable and exhibit breaking. This whole picture agrees remarkably well with the observations in experiments. For F r1 > 1.26, undulations indeed show wave breaking [43] and beyond an incoming Froude of 1.75, the undulations vanish altogether with the flow becoming turbulent. Now for standing jumps which occur on thin films, as in the present study for moderate Reynolds numbers, things can be quite different. Firstly, the assumption of a parallel incoming flow is hardly valid. There are strong variations in the film thickness in the stream wise direction as can be seen from fig. 14. Such linear variations in film thickness have also been reported experimentally by Bonn et al. [8]. A novel feature which emerges and one which has not been reported before, is the presence of separated bubbles underneath the first few undulations. Additionally, it is also clear from figs. 13 and 14 that the presence or absence of undulations is influenced not only by the inlet Froude but also by the inlet Reynolds number. At first sight it appears unlikely that these viscous undulatory jumps can be described at all by the theory of Benjamin and Lighthill [26]. Interestingly, we show here that even for these viscous jumps, a criterion obtained in ref. [26] is useful in deciding whether undulations occur or not. To show this we write down the KdV equation from ref. [26] below. We note here that this is not the KdV equation in its standard form. This form can be obtained by looking for progressive wave solutions and then integrating the resulting equation once (see [46]). ˆ2 Q 3



ˆ dh dˆ x

2 ˆ 3 − 2R ˆ 2 + 2Sˆh ˆ−Q ˆh ˆ 2 = 0, + gh

(12)

ˆ R ˆ where quantities with a hat are dimensional and Q, and Sˆ are the mass, energy and momentum flux discussed earlier. Recognizing that unlike in ref. [26], the streamwise velocity u ˆ in our simulations is not uniform, we generalize

ˆ R ˆ and Sˆ as below the definitions of Q, 

ˆ x) h(ˆ

ˆ= Q

u ˆ dˆ y, 0

1 ˆ2 + Sˆ = g h 2 and ˆ+ 1 ˆ = gh R ˆ 2Q



(13)

ˆ x) h(ˆ

u ˆ2 dˆ y,

(14)

0



ˆ x) h(ˆ

u ˆ3 dˆ y.

(15)

0

It is easily checked that for an uniform u ˆ, one recovers ˆ ˆ the corresponding expressions of Q, R and Sˆ of Benjamin ˆ R ˆ and Sˆ are the coefficients of and Lighthill [26]. Here, Q, the cubic-polynomial in eq. (12). One can solve eq. (12) in terms of elliptic functions and find solutions which would be oscillatory in space. An insightful interpretation due to [26] interprets eq. (12) as the potential plus the kinetic energy equation of an oscillator (not a simple harmonic oscillator else the potential part would be just quadratic) with total energy zero. Oscillatory solutions (in space) become possible only when this cubic equation has three distinct real roots [26] —a single real root (other roots being complex) of the cubic corresponds to non-oscillatory behaviour due to the the total energy of the oscillator ˆ and Sˆ using critbeing zero. Following [26] we rescale R ˆ 2/3 ical energy and momentum values, Rc = (3/2)g 2/3 Q 1/3 ˆ 4/3 and Sc = (3/2)g Q , respectively, taken by a uniform stream with Froude number unity. For three real roots of ˆ c and s ≡ S/S ˆ c need to the cubic in eq. (12), r ≡ R/R satisfy the following inequality [26]: D ≡ 3r2 s2 + 6rs − 1 − 4(r3 + s3 ) > 0.

(16)

Since R and S vary with space due to viscosity, the discriminant D is now a function of streamwise distance in our simulations. Figure 16 shows variation in D, for different Reynolds numbers with streamwise direction. Clearly, D is positive valued for Re = 125 for which we see undulations (fig. 13) whereas for Re = 180 and 250 the discriminant asymptotically goes to nearly zero and corresponding profiles show very weak undulations for Re = 180 (fig. 17) and almost no undulations for Re = 250 (fig. 13). As indicated by fig. 16, there are no noticeable undulations for the Re = 180 case in contrast to the Re = 125 which is

Page 12 of 14

Eur. Phys. J. E (2015) 38: 45

15

Re = 180

y/H

10 5 0 0

50

100

150

x/H

200

250

300

R/Rc

Fig. 17. Hydraulic jump formation for the Re = 180 case showing very weak undulations.

Re = 62.5 Re = 125 Re = 180 Re = 250

5

0

0

50

100

150

200

250 x/H 300

Fig. 18. Variation in r with streamwise direction for Re = 62.5, 125, 180 and 250.

Re = 62.5 Re = 125 Re = 180 Re = 250

S/Sc

4 2 0

0

50

100

150

200

250 x/H 300

Fig. 19. Variation in s with streamwise direction for Re = 62.5, 125, 180 and 250.

shown to allow multiple solutions to the equation U = 0. Figures 18 and 19 show variations of r and s individually for different Reynolds numbers. It was already observed in the inviscid analysis of [26] that the effect of change in r and s on oscillations were in opposite directions. A decrease in r produced oscillations (provided it was less than the value given by eq. (1)) whereas the same effect would be produced by an increase in s. In a given simulation wellupstream of the jump, clearly r and s both decrease with distance and thus oscillations occur downstream depend on their relative rates of decrease as well as their behaviour in the near-jump region. As we decrease Re, both r and s decay rapidly with x due to increased viscous effects but the rate of decrease in r is more than that in s. Note that the energy, r, decreases nearly monotonically whereas the momentum (s), compensated by the potential head, shows undulations at a phase lag of π with the thickness of the

film. For Re = 62.5, both r and s values fall sharply towards one and show breaking of waves (see fig. 13). It is clearly desirable to have differential equations governing the evolution of r and s with streamwise distance. If such equations can be derived, it would become possible to predict purely from the inlet Reynolds and Froude number whether undulations would occur. This work however remains for the future.

4 Conclusions A detailed study of the laminar planar hydraulic jump has been conducted based on 2D numerical simulations of the Navier-Stokes equations. Since there are no stationary solution for viscous flow on an infinitely long plate,

Eur. Phys. J. E (2015) 38: 45

we perform simulations in a computational domain consisted of a platform of finite length followed by a pit which leads to proper boundary conditions for film flow [16]. We have asked whether undular planar jumps can be formed at moderate Reynolds numbers, and answered in the affirmative. We have shown that a boundary layer shallowwater equation is sufficient to describe the flow leading up to the jump, but through most of the region following the jump, these equations are insufficient to describe the dynamics even qualitatively. The undular flow downstream of the jump is actually better described by an inviscid mechanism, and we have provided an argument for this. It is to be remembered of course that viscous effects will always be significant in this flow. One indication that the boundary layer approximation is insufficient to describe this flow is that downstream conditions such as the existence of an obstacle can affect the jump location and height, as well as the post-jump flow. To investigate this, we performed simulations with obstacles of different height placed at the edge of the plate. We showed that obstacles that cause the downstream height to be more than the height of the natural jump (for the no-obstacle case) shift the jump upstream and make the jump steeper. We also performed simulations with an obstacle, of height same as the inlet flow, placed at different locations. In this case, the obstacle has little effect on the jump location and structure, except in the case when the obstacle is placed ahead of the natural jump location. We note that the length of the plate is crucial in determining the structure and location of the jump. For Re = 125, a length of L∗ ∼ 150 gives a weak jump with no flow separation at the jump, whereas, a longer platform with L∗ ∼ 300 yields undular jump with flow separation underneath the undulations. We have also studied the effect of Re, F r and W e on jump location and structure. We showed, as noted earlier in [15, 16], that for high F r the jump location shifts downstream and a relatively weaker jump is obtained. Increasing the Re shows similar effect, with very steep jumps obtained at low Re = 62.5 with undulations undergoing wave breaking. For higher Re, the undulations vanish and a weaker jump is obtained. Surface tension has a curious effect on the jump structure. The undulations now show shorter wavelength with the amplitude decreasing considerably downstream compared to the zero surface tension (W e = ∞) case. We argue, based on the theory of Benjamin and Lighthill [26], that the excess un-dissipated energy at the jump location is radiated downstream in the forms of standing waves. Inclusion of surface tension, provides another mechanism for the consumption of energy loss at the jump by generating larger area and thus shorter waves. By computing, energy and momentum fluxes, and using the numerical values in Benjamin and Lighthill, we show that the undulations are indeed of similar origin as those behind river bores. We hope that the present predictions will motivate carefully designed experiments in the Reynolds and Froude number regimes presented in this study using the geometry we suggest here.

Page 13 of 14

References 1. S. Vigneswaran, C. Visvanathan Water Treatment Processes: Simple Options (CRC Press, 1995). 2. W.H. Hager, Energy Dissipators and Hydraulic Jump (Water Science and Tech. Lib., Kluwer Academic Publishers, 1992). 3. L. Rayleigh, Proc. R. Soc. London Ser. A. 90, 324 (1914). 4. I. Tani, J. Phys. Soc. Japan 4, 212 (1949). 5. T. Bohr, P. Dimon, V. Putkaradze, J. Fluid Mech. 254, 635 (1993). 6. S.B. Singha, J.K. Bhattacharya, A.K. Ray, Eur. Phys. J. B 48, 417 (2005). 7. S. Watanabe, V. Putkaradze, T. Bohr, J. Fluid Mech. 480, 233 (2003). 8. D. Bonn, A. Anderson, T. Bohr, J. Fluid Mech. 618, 71 (2009). 9. E.J. Watson, J. Fluid Mech. 20, 481 (1964). 10. J.W.M. Bush, J. Aristoff, J. Fluid Mech. 489, 229 (2003). 11. A.R. Kasimov, J. Fluid Mech. 601, 189 (2008). 12. J. Gajjar, F.T. Smith, Mathematika 30, 77 (1983). 13. A.D.D. Craik, R.C. Latham, M.J. Fawkes, P.W.F. Gribbon, J. Fluid Mech. 112, 347 (1981). 14. R.V. Brotherton-Ratcliffe, PhD Thesis, University College London (1987). 15. R.I. Bowles, F.T. Smith, J. Fluid Mech. 242, 145 (1992). 16. F.J. Higuera, J. Fluid Mech. 274, 69 (1994). 17. M.M. Rahman, A. Faghri, W.L. Hankey, J. Fluids Eng. 113, 411 (1991). 18. S. Chippada, B. Ramaswamy, M.F. Wheeler, Int. J. Num. Methods Eng. 37, 1381 (1994). 19. K. Yokoi, F. Xiao, Phys. Lett. A 257, 153 (1999). 20. V.G. Ferreira, M.F. Tome, N. Mangiavacchi, A. Castelo, J.A. Cuminato, A.O. Fortuna, S. McKee, Int. J. Num. Methods Fluids 39, 549 (2002). 21. M. Passandideh-Fard, A.R. Teymourtash, M. Khavari, J. Fluids Eng. 133, 014401 (2011). 22. W.G. Pritchard, L.R. Scott, S.J. Tavener, Philos. Trans. R Soc. 340, 1 (1992). 23. S. Ramadurgam, R.V.K. Chakravarthy, G. Tomar, R. Govindarajan, Phys. Fluids 24, 102109 (2012). 24. B.A. Wols, Undular Hydraulic Jumps (Delft University of Technology, 2005). 25. I. Ohtsu, Y. Yasuda, H. Gotoh, J. Hydr. Eng. 129, 948 (2003). 26. T.B. Benjamin, M.J. Lighthill, Proc. R. Soc. London, Ser. A 224, 448 (1954). 27. W. Chester, J. Fluid Mech. 24, 367 (1966). 28. M. David Holland, Rodolfo R. Rosales, Dan Stefanica, Esteban G. Tabak, J. Fluid Mech. 470, 63 (2002). 29. D. Long, P.M. Steffler, N. Rajaratnam, J. Hydr. Res. 29, 293 (1991). 30. L. Quingchao, U. Drewes, J. Hyd. Res. 32, 877 (1994). 31. A.M. Kamchatnov, Y.-H. Kuo, T.-C. Lin, T.-L. Horng, S.C. Gou, R. Clift, R.H.J. Grimshaw, J. Fluid Mech. 736, 495 (2013). 32. S. Popinet, J. Comput. Phys. 190, 572 (2003). 33. S. Popinet, J. Comput. Phys. 228, 5838 (2009). 34. J. Bell, P. Colella, H. Glaz, J. Comput. Phys. 85, 257 (1989). 35. J. Li, C.R. Acad. Sci. Paris, S´er. IIb 320, 391 (1995). 36. E. Aulisa, S. Manservisi, R. Scardovelli, S. Zaleski, J. Comput. Phys. 225, 2301 (2007).

Page 14 of 14 37. M.M. Francois, S.J. Cummins, E.D. Dendy, D.B. Kothe, J.M. Sicilian, M.W. Williams, J. Comput. Phys. 213, 141 (2006). 38. G. Tomar, D. Fuster, S. Popinet, S. Zaleski, Comput. Fluids 39, 1864 (2010). 39. R. Dasgupta, R. Govindarajan, Phys. Fluids 22, 112108 (2010). 40. H. Chanson, Eur. J. Mech. B Fluids 28, 191 (2009). 41. C.T. Avedesian, Z. Zhao, Proc. R. Soc. London, Ser. A 456, 2127 (2000).

Eur. Phys. J. E (2015) 38: 45 42. H. Favre, Etude theorique et experimental des ondes de translation dans les canaux descouverts (Dunod, Paris, 1935). 43. A.M. Binnie, J.C. Orkney, Proc. R. Soc. London Ser. A 230, 237 (1955). 44. B. Sturtevant, Phys. Fluids 8, 1052 (1965). 45. G.H. Keulegan, G.W. Patterson, J. Res. Nat. Bureau Standards 24, 47 (1940). 46. L. Debnath, Nonlinear Water Waves (Academic Press, 1994).