Numerical simulation of capillary gravity waves excited by an obstacle in shallow water

Capillary gravity waves excited by an obstacle are investigated by numerical simulations. Under the resonant condition for which large-amplitude solitary waves are generated, solutions of the Euler equations show that the capillary effects induce the generation of short waves both upstream of the solitary waves and downstream of the obstacle. Overall characteristics of these waves agree with the weakly nonlinear theory, although the theory overestimates the wavelength of the upstream short waves.


INTRODUCTION
We consider the capillary gravity waves excited by an obstacle in an open-channel flow of uniform velocity (Fig. 1).In an inviscid fluid, the flow is governed by the Bond number and the Froude number, which are defined by where T is the surface tension, ρ is the density, g is the acceleration due to gravity, D is the fluid depth, and U is the mean flow velocity.The Bond number represents the ratio between the surface tension and the gravity force, and the Froude number represents the mean flow velocity relative to the velocity of the long wave.We should mention here that the Weber number We = ρDU 2 /T , which is related to the Bond number by We = Fr 2 /Bo, can be used instead of Bo.While the small-amplitude waves are described well by the linear theory [1], large-amplitude nonlinear waves are resonantly excited when Fr ≃ 1, for which the mean flow velocity agrees with the long-wave speed.Therefore, several weakly nonlinear theories, including the forced Korteweg-de Vries (fKdV) equation, have been developed to describe the large-amplitude waves.
Nonlinear behaviour of pure gravity waves (Bo = 0) under resonant condition (Fr ≃ 1) has been intensively investigated experimentally [2,3] and numerically [4], and the fKdV equation has been found to give good results, including the generation of unsteady undular bores upstream and downstream of the obstacle [5,6].On the other hand, we have comparatively little knowledge about the waves with capillary effects, and the verification of weakly nonlinear theories such as the 5th-order forced KdV equation has yet to be done.
In this study, we have investigated the unsteady behaviour of capillary gravity waves under the exactly resonant condition of Fr = 1, by numerically solving the Euler equations.At the same time, applicability of the weakly nonlinear theories has been tested by comparing their solutions with the solutions of the Euler equations.

Euler equations
We assume a two-dimensional flow illustrated by Fig. 1.The governing equations are the continuity equation and the non-dimensional Euler equations, which are given by and where u u u = (u, w) is the velocity with u and w being the horizontal and vertical component, respectively; p is the pressure; and ẑ z z is the unit vector in the vertical direction.In this study, we scale the length by D, the velocity by U, and the pressure by ρU 2 .
Two boundary conditions must be imposed on the free surface, namely the kinematic boundary condition and the dynamic boundary condition.The kinematic boundary condition states that the fluid particles on the free surface move always along the surface, i.e., where f (x,t) is the vertical displacement of the free surface.The dynamic boundary condition describes the continuity of the stress at the free surface, i.e., where p 0 is the atmospheric pressure and κ is the curvature of the free surface.Since we assume an inviscid fluid in this study, Eq. ( 5) describes simply the continuity of normal stress.As the initial condition, impulsive start is used, i.e., the velocity (u, w) = (U, 0) is instantly applied at t = 0. We use the body-fitted curvilinear coordinates in which the grid lines are adjacent to the bottom obstacle and also to the free surface.The bottom obstacle has a fixed shape, while the free surface evolves according to Eq. ( 4).We rewrite the governing equations in the transformed curvilinear coordinates and discretize them by the finite difference method, and solve them by the Marker and Cell method.To avoid the computational instability in the numerical solution of inviscid Euler equations, we use a third-order upwind scheme for the convection term, while all the other spatial derivatives are discretized by the central difference of secondorder accuracy.For the time integration, the secondorder Adams-Bashforth method is used for the Euler equations, while the Crank-Nicolson method is used for the kinematic boundary condition at the free surface.
The obstacle shape is given by z = h(x), where the function h(x) is given by An example of the grid for solving the Euler equations, which has the total number of grid points 12 000 (horizontal) × 150 (vertical), is demonstrated in Fig. 2. The minimum horizontal resolution is ∆x = 0.02, corresponding to the shortest wavelength (∼ 1) that appears upstream of the obstacle.The grid is concentrated near the free-surface and the bottom boundary, giving the minimum vertical resolution (∆z = 0.005) there.
Numerical accuracy is assured by checking the dependence of the solution on the computational grid size.
Free-surface displacements for the case of Bo = 0.18, obtained from two computational grids, are presented in Fig. 3, where the 24 000 × 300 grid has twice the resolution of the 12 000 × 150 grid (Fig. 2) everywhere.Despite the different resolution, the difference in the free-surface displacements is negligibly small not only in the overall description of the wave (Fig. 3a) but also in the enlarged figures of upstream short waves (Fig. 3b) and downstream short waves in the depression region (Fig. 3c).Since the wavelength of the short waves is minimum at Bo = 0.18 and increases with Bo, the accuracy at Bo = 0.18 assures the accuracy for larger values of Bo.Then, we have used the 12 000 × 150 grid (Fig. 2) for the computation for Bo < 0.3.For higher Bond numbers (Bo ≥ 0.3) the 6000 × 150 grid, which has minimum intervals ∆x = 0.04 and ∆z = 0.005, was found to be large enough to obtain sufficient accuracy.The time increment is ∆t = 2.0 × 10 −4 , and the computation was continued until t = 400 for Bo > 0 and t = 500 for Bo = 0.
Most of the computational results were obtained using the NEC-SX9 and NEC-SX-ACE computer in the Cybermedia Center of Osaka University.

Weakly nonlinear theories
The fKdV equation has been proposed to describe the weakly nonlinear wave near the resonant condition (Fr ≃ 1).It is applicable to the cases with strong capillary effects (Bo > 1/3) [7] or without capillary effects (Bo = 0) [8,9].The equation can be derived if the wave amplitude a and the wavelength L are scaled as where ε is a small parameter.Under the exact resonance condition (Fr = 1), the fKdV equation becomes where the variables in Eq. ( 9) are non-dimensionalized by D and U(= √ gD) in the same way as the Euler equations.
In the limit of Bo → 1/3, the dispersive term f xxx in Eq. ( 9) vanishes so that the fKdV equation no longer contains the dispersion term.This suggests that the inclusion of a higher-order dispersion term becomes necessary to avoid the overturning of the waves.Indeed, in this parameter region (Bo ≃ 1/3), a 5th-order forced KdV (5th-order fKdV) equation with the 5th-order derivative term ( f xxxxx ) has been proposed [10].For the derivation of the 5th-order fKdV equation, different scalings must be used.Under the exact resonant condition (Fr = 1), the equation reduces to We have solved the fKdV equation and the 5th-order fKdV equation numerically.The spectral method and the 4th-order Runge-Kutta method are used for the space and time discretization, respectively.The computed horizontal range is −768 ≤ x ≤ 768, and the computation was continued until t = 400 for Bo > 0, and t = 500 for Bo = 0.
For the test of spatial resolution, we have varied the cut-off wave number N in the spectral method.The solutions of the 5th-order fKdV equation at t = 400 for Bo = 0.3, using N = 1536 and N = 3072, are presented in Fig. 4. It turns out that the two solutions have almost no difference (two lines overlap), thus we use N = 1536 for all the computations of the 5th-order fKdV equation and the fKdV equation, with the time increment of ∆t = 5.0 × 10 −3 .

Numerical results of the Euler equations
We now demonstrate numerical solutions of the Euler equations for different Bond numbers.Without capillary effects (Bo = 0) (Fig. 5), solitary waves are generated by the obstacle and they propagate upstream at a constant speed, while a flat depression region appears and elongates downstream proportional to time, and modulated cnoidal waves emerge downstream of the depression.
When the Bond number becomes non-zero (Bo = 0.18, 0.25; cf.Fig. 6a and 6b), a short-wave train whose wavelength is comparable to the water depth appears both upstream of the solitary waves and downstream of the obstacle in the depression region.Upstream short waves appear as the amplitude of the solitary wave becomes high with time.They have a constant wavelength and their appearing region constantly spreads upstream, although their amplitudes decrease in the upstream direction.On the contrary, downstream short waves show irregular behaviours, both in space and time.
As the Bond number further increases (Bo ≥ 0.30), solitary waves gradually deform into a flat elevated wave, and the downstream cnoidal waves disappear.The counterpart of short waves observed at Bo = 0.18 appears even at large Bond numbers (Bo 0.30), with larger wavelength and amplitude as the Bond number increases.
At the highest Bond number investigated in this study (Bo = 2/3; cf.Fig. 6f), the original short waves become long nonlinear waves with larger wavelengths.Namely, the upstream short waves finally become convex cnoidal waves, while the downstream short waves in the depression region become downstream-advancing solitary waves with negative displacement.This actually reflects the fact that the wave pattern at Bo = 2/3 has a rotation symmetry to the pattern observed at Bo = 0. We can easily verify that the fKdV equation ( 9) is invariant by the simultaneous change of the variables' sign (η → −η ′ , x → −x ′ ) and the parameter (Bo → 2/3 − Bo ′ ).
We now examine the behaviour of short waves upstream of the solitary wave.Figure 7a shows the phase speed c p of the upstream short wave and the upstreamadvancing speed c n of the solitary wave (or that of the flat elevated wave for Bo 0.30).These speeds were measured by tracking one of the wave crests, although we could not measure the values of c n at large Bond numbers (Bo 0.4) because the nearly flat elevated wave does not have a clear peak which can be tracked accurately.We note that the phase speed of the upstream short wave is nearly equal to the speed of the solitary wave (c p ≃ c n ) at least up to Bo ∼ 0.4.We can directly measure the wavelength of upstream short waves λ using the data at t = 400.On the other hand, the wavelength can also be 'estimated' from the phase speed c p of the same wave if we assume that the short wave satisfies the linear dispersion relation The wavelengths estimated from the value of c p in Fig. 7a using the linear dispersion relation (13) are denoted by λ l , and they are given in Fig. 7b for comparison with λ .In this process, we have to solve Eq. ( 13) for λ l when c p is provided.The solution could be obtained by Newton's method.We note that λ l agrees well with λ for Bo ≤ 0.4, while some differences appear for larger Bo.This means that the upstream short waves are almost linear waves, although their amplitude becomes larger for larger Bo, and nonlinear effects should become more significant.
These results show that the upstream short waves are linear waves whose phase speeds are equal to the upstream solitary waves.These characteristics are similar to those discussed by Karpman [11] as 'soliton radiation', in which a soliton emits dispersive waves whose phase speed is equal to that of the soliton.
The motions of downstream short waves appearing in the depression region are much more irregular and complex than the upstream waves.This might be attributable to the wave interactions of short waves radiated by large-amplitude nonlinear waves.However, this is still a conjecture.

Comparison between the Euler equations and the weakly nonlinear theories
We first compare the solution of the Euler equations with that of the fKdV equation (Fig. 8).Since the fKdV equation would be applicable to the cases with large Bo (≫ 1/3) and very small Bo as discussed in Section 3, we compare the results for the case of Bo = 2/3 (and Fr = 1) in Fig. 6.Two solutions are in good agreement, showing that the fKdV equation is a good approximation to the Euler equation at the resonant condition.
We next compare the solution of the Euler equations with that of the 5th-order fKdV equation (Fig. 9).Since the 5th-order fKdV equation would be applicable to the cases of Bo ∼ 1/3 as discussed in Section 3, we compare the results for the case of Bo = 0.3 (and Fr = 1).Two solutions give qualitative agreement in the generation of the flat elevated wave and the short waves upstream of the obstacle and the downstream short waves in the depression region.A major quantitative difference is in that the 5th-order fKdV equation predicts longer wavelengths of the upstream and downstream short waves than the Euler equations.Since the 5th-order fKdV equation assumes a long wave, it is natural that the nearly flat elevated wave is predicted well including its propagation speed, while the short waves, especially their 'short' wavelengths and the corresponding 'slow' propagation speed, are not predicted well.

CONCLUSIONS
We have investigated the capillary gravity waves excited by an obstacle by solving the Euler equations.It was found that the capillary effects excite short waves further upstream of the upstream solitary waves.Irregular short waves also appear downstream of the obstacle in the depression region.The observed phase speed of the upstream waves is equal to the speed of the solitary wave, in agreement with the soliton radiation theory.
Comparison of results with the weakly nonlinear theories reveals that the fKdV equation is applicable to the cases with strong capillary effects (Bo ≫ 1/3), while the 5th-order fKdV equation with a 5th-order dispersion term is qualitatively applicable to the case of Bo ≃ 1/3, for which the fKdV equation is not applicable.We note, however, that the 5th-order fKdV equation is derived under the assumption of a long wave, and it generally overestimates the wavelength and the propagation speed of the upstream short waves.

Fig. 2 .
Fig. 2. Grid used for the computation.Total number of grid points is 12 000 (horizontal) × 150 (vertical).(a) Whole region of the computation.Only every 100 grid points are depicted in the horizontal direction and only every two points are depicted in the vertical direction.(b) Enlarged image of the grid near the obstacle.Only every 20 grid points are depicted in the horizontal direction.

Fig. 5 .
Fig. 5. Time development of the free-surface displacement a(x,t) obtained by the numerical solution of the Euler equations without capillary effects (Bo = 0).

Fig. 6 .
Fig. 6.Time development of the free-surface displacement a(x,t) obtained by the numerical solution of the Euler equations.

Fig. 7 .
Fig. 7. (a) Phase speed c p of the upstream short wave and the propagation speed c n of the solitary wave or the flat elevated wave.(b) Wavelength λ of the short wave measured at t = 400 and the wavelength λ l estimated from c p given in Fig. 4a using the linear dispersion relation (13).