Abstract

A simple radial basis function (RBF) meshless method is used to solve the two-dimensional shallow water equations (SWEs) for simulation of dam break flows over irregular, frictional topography involving wetting and drying. At first, we construct the RBF interpolation corresponding to space derivative operators. Next, we obtain numerical schemes to solve the SWEs, by using the gradient of the interpolant to approximate the spatial derivative of the differential equation and a third-order explicit Runge–Kutta scheme to approximate the temporal derivative of the differential equation. For the problems involving shock or discontinuity solutions, we use an artificial viscosity for shock capturing. Then, we apply our scheme for several theoretical two-dimensional numerical experiments involving dam break flows over nonuniform beds and moving wet-dry fronts over irregular bed topography. Promising results are obtained.

1. Introduction

Development of robust meshless methods for the numerical solution of partial differential equations has attracted considerable interest over the past twenty years, for example, [13]. There are three types of meshless techniques: meshless techniques based on weak forms such as the element-free Galerkin method [4], meshless techniques based on collocation techniques such as the meshless collocation technique based on radial basis functions (RBFs) [5], and meshless techniques based on the combination of weak forms and collocation technique. In the literature, several meshless weak form techniques have been reported; among others, the smooth particle hydrodynamic method [6] and boundary point interpolation method are worth noticing [7]. The weak forms are used to derive a set of algebraic equations through a numerical integration process using a set of quadrature domain that might be constructed globally or locally in the domain of the problem. In this topic, Liu et al. [8] applied the notion of meshless local Petrov–Galerkin and developed the meshless local radial point interpolation method. This method was studied and used on a class of three-dimensional wave equations later by Shivanian [9]. A combination of the natural neighbour finite element method with the radial point interpolation method, using the multiquadric radial basis function, is developed by Dinis et al. [10] and Belinha et al. [11] to analyse a three-dimensional solid.

Initially, the radial basis function meshless method was developed for data surface fitting, and later, with the work developed by Kansa [5], the radial basis function was used for solving partial differential equations. Fedoseyev et al. [12] and Cheng et al. [13] have shown that meshless radial basis functions (RBFs) are attractive options because of the exponential convergence of certain RBFs. Various RBFs have been successfully applied to obtain accurate, efficient solutions of partial differential equations of engineering interest. This method has been applied to solve inviscid compressible flows [14], natural convection [15], heat conduction [16], three-dimensional incompressible viscous flows [17], and long waves in shallow water [18].

The Saint-Venant equations, also called shallow water equations, have often been preferred for the propagation of the flow in the open channel, and they exhibit a simplified mathematical structure, with an ability to take into account the smoothly varying flow conditions and flow discontinuities such as hydraulic jumps, moving bores, and propagation over dry beds, despite it still poses many theoretical and practical problems [19]. The dam break flow problem is an ideal theoretical example that involves these hydraulic challenges. In this paper, we examine the application of the radial basis function meshless method to the numerical solution of the shallow equations for dam break flow problems involving wetting/drying over complicated, frictional bed topography. In the formulation of RBFs, the radial basis functions represented by the multiquadric functions are employed as basis functions to compute weighting coefficients for space differential operators, over a global set of computational collocation points [20]. The friction term is included in the momentum equations and discretized by a splitting implicit scheme [21]. A third-order Runge–Kutta algorithm is used for time integration [22].

It is known that the system of shallow water equations admits nonsmooth solutions that may contain shocks and rarefaction waves and, in the case of nonsmooth bottom topography, also contain contact discontinuities. To perform properly, a numerical method must be nonlinearly stable because linearly stable methods may develop large spurious oscillations and even blow up. To stabilize the proposed RBF model for slowly varying flows, as well as rapidly varying flows involving shocks or discontinuities such as dam breaks and hydraulic jumps, an artificial viscosity technique is used, following [23]. Furthermore, to avoid numerical instability caused by negative water depth near wet/dry fronts, the local flow variables are modified, imposing zero discharge when the water height became very small. As a consequence, the present numerical scheme ensures preservation of nonnegative water depth, and there is no need to limit the fluxes during simulation.

This paper is organized as follows: Section 2 outlines the shallow water equations. Section 3 presents the general formulation of partial differential problem interpolation using RBFs. Section 4 describes the application of RBFs to generate the discrete form of shallow water equations. Details are given of the numerical methods used to represent the friction term and carry out time integration. Section 4 presents the way to use the artificial viscosity for shock capturing. Section 6 discusses validation and application of the method; several numerical experiments are carried out for previously published benchmark cases in order to confirm the potential of the proposed scheme. Section 7 summarizes the main findings.

2. Two-Dimensional Shallow Water Equations

In conservation form, the two-dimensional nonlinear shallow water equations are given by (1) below [24]. The depth-averaged continuity and momentum equations in the x- and y-directions arewhere is the total depth from the sea bed to the free surface, and are the depth-averaged velocity components in the Cartesian and directions, is the bed elevation above a fixed horizontal datum, is the acceleration due to gravity, and and are the bed shear stress components, which are defined aswhere is the water density and is the bed friction coefficient, which may be estimated from , where is the Manning coefficient.where is the vector of dependent variables, and are the inviscid flux vectors, and is the vector of source terms. In full, the vectors are

3. The Radial Basis Function Meshless Method

Let be an open domain of . Suppose is a function to be approximated in a set of N pairwise distinct nodes . In the RBF meshless scheme, the approximation of at the node can be written as a linear combination of N RBFs:where are the function values at node , is a RBF centered at , denotes the Euclidean norm between nodes and , and are the coefficients to be determined.

One of the most commonly used RBFs is the multiquadric (MQ) RBFs [25]. Here, we use the infinitely smooth multiquadric radial basis function defined aswhere is a shape parameter that controls the fit of a smooth surface to the data. In the present work, we used the following selection [26]:where denotes the smallest nodal distance. The expansion coefficients in (5) are obtained by solving the following linear system of algebraic equations of :

The expansion coefficients are then calculated bywhere is the vector of approximate solutions and is an matrix of RBFs given as

Space derivatives of the interpolant (5) may be readily calculated, due to its linearity. Generally, the first- and second-order spatial derivatives at the point can be calculated aswhere  = [1, 2] denotes the first- and second-order derivatives. In a compact matrix form, using (9), (11) can be written aswhere

4. Discretized Form of Shallow Water Equations by the RBF Meshless Method

4. 1. Convective Flux and Bottom Topography Term Approx-imations

Let us assume that are known convective fluxes along the axes and at time , where . Using (5), it can be approximated by

In the matrix form, this equation becomes

However, the expansion coefficients are then calculated by

From (14), the partial flux derivative can be written as

Then, the compact matrix form of the partial derivative of the flux vector is

Applying the same procedure described above for the bottom topography function , it is easy to obtain its discrete form. The topography function can be approximated by

Its partial derivative is then expressed as

In the same way, the compact matrix form of the bottom topography term along the axis is given as

4.2. Addition of Friction Effects

To incorporate friction into the present numerical scheme, the friction term is discretized using an operator splitting procedure described by Boushaba et al. [27], which splits the shallow water equations (2) into two equations:where . In the first step of the calculation, the upper ordinary differential equation in (23) is approximated by an implicit method as described in [21]:where the friction term may be expressed using a Taylor series as

Rearranging the above equation leads to the following formula for updating water discharge at the new time step:

In the second step, the value is taken to be the initial condition when solving the second equation in (23).

4.3. Time Integration and Stability Condition

To date, the forward Euler method has been mainly used as the preferred time-stepping scheme for RBF methods. However, the forward Euler method is only the first-order accurate in time and so may introduce excessive numerical dissipation in the computed RBF solutions. To achieve a higher order of accuracy, we use the explicit Runge–Kutta method recommended by [22]. The procedure to advance the solution from the time to the next time is carried out aswhere represents the time level and is the time step. To achieve stability, for this explicit scheme, the time step must meet the following criterion:where is the Courant number, such that 0 < CFL < 1, and denotes the smallest nodal distance between collocation points.

4.4. Boundary Conditions

In this work, transmissive boundary conditions for open inflow/outflow and reflective boundary conditions for solid walls are applied in the simulations. At a transmissive boundary, the flow variables at a collocation point on the boundary are set to the same values as at the interior point normal to the boundary. At a reflective boundary, the value at a collocation point is simply the mirror image of that at the associated boundary point so that the normal velocity component is zero at the boundary. However, the representation of boundary conditions within an RBF method is less trivial.

4.5. Implementation Algorithms

For a given initial condition , the time integration procedure is described in Algorithm 1. The right-hand side (RHS) used in the algorithm represents the convective flux and bottom topography that are calculated in Algorithm 2.

(1)
(2)   
(3)
(4)
(5)   
(6)   
(7)   
(8)
(9)   
(10)    
(11)  
(12)   
(1)
(2)
(3)
(4)
(5)
(6)

5. Artificial Viscosity for Shock Capturing

For particular hydraulic problems involving shock waves, such as a hydraulic jump or dam break flows, the numerical model is required to represent a steady or unsteady discontinuity, where the presence of oscillations in the solution is expected, and it sometimes may grow over time. However, by introducing a small amount of artificial diffusion, it is possible to damp these oscillations [23, 28]. We therefore augment the right-hand side of the hyperbolic system (2) with the artificial diffusion term:where is a tunable viscosity coefficient. The Peclet number is defined as follows:

It often controls the stability of the numerical solutions [23, 29]. In the case of convection-dominated flow, the Peclet number is large but finite, and the effect of the diffusion term (29) becomes negligible. Therefore, a natural and simple way to stabilize the solution is to reduce the Peclet number.

Using the RBF interpolation of an arbitrary function (12), the augmented term (29) can be obtained as

6. Numerical Results

6.1. One-Dimensional Dam Break over a Wet Bed

A one-dimensional dam break over a wet, horizontal bed is first simulated to demonstrate the shock-capturing capability of the corrected radial basis function method using artificial viscosity. The channel is of length 20 m. The initial water discharge and depth are given by

As boundary conditions, a zero discharge and a free boundary are considered at the left and right ends of the channel. The analytical solution for this simple dam break test consists of a backward-propagating rarefaction and a forward-moving shock wave. Figure 1 shows the numerical results of flow after the dam fails, at three different times in terms of water surface elevation and discharge. Here, the model domain consists of 800 points. An artificial viscosity of 0.001 m2s−1 is applied corresponding to a Peclet number of 300. Figure 2 shows close-up views of the evolving free surface elevation and discharge profiles along the channel at different times. In general, satisfactory agreement is achieved between the numerical and analytical solutions, although a very small amount of numerical diffusion occurs where the surface gradient is steepest.

6.1. Two-Dimensional Partial Dam Break Simulation

For this second test case, we consider a hypothetical two-dimensional dam break problem with nonsymmetric breach that is a typical validation made in many presented papers, for example, [3033]. An illustration of this problem is shown in Figure 3, in which the domain is defined by a 200 m × 200 m channel with the horizontal bed. A dam is located in the middle of the domain, and the nonsymmetric breach is 75 m wide with negligible thickness and is located 95 m from the left side of the domain. The initial water discharge and depth are given by

As boundary conditions, at the left and right ends of the channel, a zero discharge and a free boundary are considered, and a solid wall is considered for others. Figure 4 shows the numerical results of flow after the dam fails, at three different times in terms of water depth and velocity. As we can show, the upstream water is released into the downstream side through the breach, creating the surge wave propagating to downstream and the neglective wave to upstream. Here, the model domain consists of 2000 points. An artificial viscosity of 0.001 m2s−1 is applied corresponding to a Peclet number of 300. No analytical solution is available for this case, but the results can be compared with those of other numerical schemes. Generally, the results are in agreement with the results obtained by other numerical methods in aforementioned literatures.

6.2. Circular Dam Break

This second test case consists of the instantaneous breaking of a cylindrical tank of diameter 20 m (Figure 5), initially filled with 2 m of water at rest. The wave generated by the breaking of the tank propagates into still water with an initial depth of 0.5 m. This classic problem is widely used to test the shock capture capability of numerical schemes [34, 35].

Figure 6 illustrates the wave propagation on computational 8000 collocation points. An artificial viscosity of 0.04 m2s−1 is applied corresponding to a Peclet number of 100. At the beginning, that is, at t = 0 s, the dam is broken instantaneously and the column of water is released, and the shock wave results in the increase of water depth in the lower depth region, propagating in the radial direction. Our state solutions obtained at time t = 1 and 2.5 s by different methods, that is, the RBF meshless method and the triangular finite volume method, as a based mesh method and a 1D reference solution are in an agreement (Figure 7). No singular corner effects on smoothness of the solution can be noticed.

6.3. Dam Break over Three Humps

The nonlinear shallow water equations are not too many analytical solutions for the problem of the flow over wetting and drying topography. In this context, Antuono et al. [36] studied an analytical solution of the shallow water flow over uneven topography, and this is for analysis of wave run-up on a beach. The current shallow water model is applied to simulate a dam break over an initially dry floodplain with three humps, recommended by Kawahara and Umetsu [37] for this challenging problem that involves complex flow hydrodynamics, bottom topography, and wetting and drying. The simulated setup is sketched in Figure 8, where the dam break occurs in a 75 × 30 m rectangular domain with the dam located 16 m away from the upstream end. A reservoir with a still water surface elevation of 1.875 m is assumed to be upstream of the dam. The bed topography of the domain is defined as

The floodplain is initially dry, and a constant Manning coefficient of 0.018 s/m1/3 is used throughout the domain. The domain walls are assumed to be solid. The dam collapses instantly at t = 0 s, and simulation is run for 300 s on 2000 collocation points.

After the dam fails, the initial still water in the reservoir rushes onto the downstream floodplain. In Figure 9, we present the profile of water depth and contours at different output times t = 1, 6, 12, 30, and 300 s. As can be observed from these results, after about 1 s, the wet-dry front reaches the two small humps and begins to climb over them. At t = 6.0 s, the two small humps are entirely submerged and the wave front hits the large hump. At t = 12 s, the wave front passes through both sides of the large hump and begins to submerge the lee of the hump. Finally, the flow becomes steady due to the dissipation caused by bed friction, as shown at t = 300 s when the flow is nearly motionless and the peaks of the humps are no longer submerged.

The overall flow pattern for this example is preserved with no spurious oscillations appearing in the results by the RBF meshless method. Obviously, the obtained results verify the stability and the shock-capturing properties of the proposed meshless method. The RBF meshless method performs well for this test problem since it does not diffuse the moving fronts, and no spurious oscillations have been observed when the water flows over the humps.

7. Conclusions

This paper has investigated a simple RBF meshless method for numerical simulation of dam break flows by solving the shallow water equations including bed friction and nonuniform bed elevation terms. The space derivative has been approximated by radial basis functions on global collocation points. An explicit third-order Runge–Kutta scheme was used for time integration. The RBF meshless method has been found to be flexible and straightforward to implement. The method is particularly attractive for the shallow water equations in comparison with classical mesh-based methods because it is inherently mesh-free and requires no special treatment of the wet-dry interface. The method was tested against several theoretical dam break problems, including a dam break with a nonsymmetric breach, a circular dam break, and a dam break on an irregular bed involving the wet/dry interface. The numerical model gives promising predictions in comparison with previously references and published results.

Conflicts of Interest

The author declares that there are no conflicts of interest regarding the publication of this paper.