Numerical Simulation of Dam Break Flows Using a Radial Basis Function Meshless Method with Artificial Viscosity

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. 'en, we apply our scheme for several theoretical twodimensional numerical experiments involving dam break flows over nonuniform beds and moving wet-dry fronts over irregular bed topography. Promising results are obtained.


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, [1][2][3].
ere 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].e 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.is 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.is 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].
e 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].e 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].e 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.
is 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.

Two-Dimensional Shallow Water Equations
In conservation form, the two-dimensional nonlinear shallow water equations are given by (1) below [24] where h is the total depth from the sea bed to the free surface, u and v are the depth-averaged velocity components in the Cartesian x and y directions, z b is the bed elevation above a fixed horizontal datum, g is the acceleration due to gravity, and S fx and S fy are the bed shear stress components, which are defined as where ρ is the water density and C b is the bed friction coefficient, which may be estimated from C b � (gn 2 M /h 1/3 ), where n M is the Manning coefficient.
where W is the vector of dependent variables, F x and F y are the inviscid flux vectors, and S is the vector of source terms.
In full, the vectors are

The Radial Basis Function Meshless Method
Let Ω be an open domain of R 2 .Suppose U � U(x, t) is a function to be approximated in a set of N pairwise distinct nodes x � x 1 , x 2 , . . ., x N  .In the RBF meshless scheme, the approximation of Uat the node x � (x i , y i ) can be written as a linear combination of N RBFs: 2 Modelling and Simulation in Engineering T is a RBF centered at x, ‖.‖ denotes the Euclidean norm between nodes x and x j , and Λ�[λ 1 ,λ 2 ,...,λ N ] T 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 as where ε 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 d m denotes the smallest nodal distance.e expansion coefficients Λ in ( 5) are obtained by solving the following linear system of N × N algebraic equations of U(x i , t) � U i : e expansion coefficients are then calculated by where 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 x can be calculated as where k � [1, 2] denotes the first-and second-order derivatives.In a compact matrix form, using ( 9), ( 11) can be written as where

1. Convective Flux and Bottom Topography Term Approximations.
Let us assume that F n ξ (W(x)) are known convective fluxes along the axes x and y at time t � nΔt, where ξ � (x, y).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 (zF n ξ (W)/zξ) can be written as en, the compact matrix form of the partial derivative of the flux vector is Applying the same procedure described above for the bottom topography function z b (x), it is easy to obtain its discrete form.e 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 S n z b along the axis ξ is given as

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: Modelling and Simulation in Engineering where U � (u, v).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 S n+1 b may be expressed using a Taylor series as Rearranging the above equation leads to the following formula for updating water discharge hu at the new time step: In the second step, the value hU is taken to be the initial condition when solving the second equation in (23).

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].e procedure to advance the solution from the time t n to the next time t n+1 is carried out as ΔtL W (1)   , where n represents the time level and Δt is the time step.To achieve stability, for this explicit scheme, the time step must meet the following criterion: where CFL is the Courant number, such that 0 < CFL < 1, and d min denotes the smallest nodal distance between collocation points.

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.

Implementation Algorithms.
For a given initial condition W 0 � W(X, t � 0), the time integration procedure is (1) ALGORITHM 2: Calculation of flux vector and right-hand side. 4 Modelling and Simulation in Engineering described in Algorithm 1. e right-hand side (RHS) used in the algorithm represents the convective ux and bottom topography that are calculated in Algorithm 2.

Artificial Viscosity for Shock Capturing
For particular hydraulic problems involving shock waves, such as a hydraulic jump or dam break ows, 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 arti cial di usion, it is possible to damp these oscillations [23,28].We therefore augment the right-hand side of the hyperbolic system (2) with the arti cial di usion term: where D h is a tunable viscosity coe cient.e Peclet number is de ned as follows: It often controls the stability of the numerical solutions [23,29].In the case of convection-dominated ow, the Peclet number is large but nite, and the e ect of the di usion term (29) becomes negligible.erefore, 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 As boundary conditions, a zero discharge and a free boundary are considered at the left and right ends of the channel.e 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 ow after the dam fails, at three di erent times in terms of water surface elevation and discharge.Here, the model domain consists of 800 points.An arti cial viscosity of 0.001 m 2 s −1 is applied corresponding to a Peclet number of 300. Figure 2 shows close-up views of the evolving free surface elevation and discharge pro les along the channel at di erent times.In general, satisfactory agreement is achieved between the numerical and analytical solutions, although a very small amount of numerical di usion occurs where the surface gradient is steepest.6.1.Two-Dimensional Partial Dam Break Simulation.For this second test case, we consider a hypothetical twodimensional dam break problem with nonsymmetric breach that is a typical validation made in many presented papers, for example, [30][31][32][33].An illustration of this problem is shown in Figure 3, in which the domain is de ned 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.e initial water discharge and depth are given by hu(t 0, x, y) hv(t 0, x, y) 0, 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 ow after the dam fails, at three di erent 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 arti cial viscosity of 0.001 m 2 s −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.

Circular Dam Break.
is second test case consists of the instantaneous breaking of a cylindrical tank of diameter 20 m (Figure 5), initially lled with 2 m of water at rest.e wave generated by the breaking of the tank propagates into still water with an initial depth of 0.5 m. is 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 arti cial viscosity of 0.04 m 2 s −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 6 Modelling and Simulation in Engineering 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 di erent methods, that is, the RBF meshless method and the triangular nite volume method, as a based mesh method and a 1D reference solution are in an agreement (Figure 7).No singular corner e ects on smoothness of the solution can be noticed.

Dam Break over ree
Humps. e nonlinear shallow water equations are not too many analytical solutions for the problem of the ow over wetting and drying topography.In this context, Antuono et al. [36] studied an analytical solution of the shallow water ow over uneven topography, and this is for analysis of wave run-up on a beach.e current shallow water model is applied to simulate a dam break over an initially dry oodplain with three humps, recommended by Kawahara and Umetsu [37] for this challenging problem that involves complex ow hydrodynamics, bottom topography, and wetting and drying.e 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.e bed topography of the domain is de ned as e oodplain is initially dry, and a constant Manning coe cient of 0.018 s/m 1/3 is used throughout the domain.
e domain walls are assumed to be solid.e 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 oodplain.In Figure 9, we present the pro le of water depth and contours at di erent  Modelling and Simulation in Engineering 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 ow becomes steady due to the dissipation caused by bed friction, as shown at t 300 s when the ow is nearly motionless and the peaks of the humps are no longer submerged.Modelling and Simulation in Engineering e overall ow 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.e RBF meshless method performs well for this test problem since it does not di use the moving fronts, and no spurious oscillations have been observed when the water ows over the humps.

Conclusions
is paper has investigated a simple RBF meshless method for numerical simulation of dam break ows by solving the shallow water equations including bed friction and nonuniform bed elevation terms.e 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.e RBF meshless method has been found to be exible and straightforward to implement.e 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.
e 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.e numerical model gives promising predictions in comparison with previously references and published results.Modelling and Simulation in Engineering

Figure 1 :
Figure 1: Numerical and exact solutions of a dam break on a wet bed, at di erent times, using a grid of 800 points: (a) free surface elevation pro les; (b) pro les of water discharge per unit breadth.

Figure 3 :Figure 2 :
Figure 3: e computational domain and RBF collocation points of the 2D dam break problem.

Figure 5 :
Figure 5: e computational domain, RBF collocation points, and location of the circular dam break problem.

Figure 4 :
Figure 4: e pro le of water depth and velocity vectors of the 2D dam break at di erent output times: t 3 s, t 5 s, and t 7.2 s.

Figure 6 :
Figure 6: e water depth contours and pro le of the circular dam break at di erent output times: t 1 s and t 2.5 s.

Figure 7 :
Figure7: Cross section comparison between the numerical solution using the RBF meshless method and 2D FV-Roe scheme and the 1D reference solution.

Figure 8 :
Figure 8: e computational domain, RBF collocation points, and bed contours of the dam break problem over humps.

Figure 9 :
Figure9: e pro le of water depth and depth contours of the dam break over three humps at di erent output times: t 1 s, t 6 s, t 12 s, t 30 s, and t 300 s.