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

Authors: Elmiloud Chaabelasri
Journal: Modelling and Simulation in Engineering
 (2018)

---

## 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.


---

## Body

## 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, [1–3]. 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 1+εx−xi2are 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: Section2 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 are(1)∂h∂t+∂hu∂x+∂hv∂y=0,∂hu∂t+∂hu2+gh2/2∂x+∂huv∂y=−gh∂zb∂x−Sfxρ,∂hv∂t+∂huv∂x+∂hv2+gh2/2∂y=−gh∂zb∂y−Sfyρ,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, zb is the bed elevation above a fixed horizontal datum, g is the acceleration due to gravity, and Sfx and Sfy are the bed shear stress components, which are defined as(2)Sfx=ρCbuu2+v2,Sfy=ρCbvu2+v2,where ρ is the water density and Cb is the bed friction coefficient, which may be estimated from Cb=gnM2/h1/3, where nM is the Manning coefficient.(3)∂W∂t+∂FxW∂x+∂FyW∂y=SW,where W is the vector of dependent variables, Fx and Fy are the inviscid flux vectors, and S is the vector of source terms. In full, the vectors are(4)W=hhuhv,Fx=huhu2+gh22huv,Fy=hvhuvhv2+gh22,S=0−ghSfx−gh∂zb∂x−ghSfy−gh∂zb∂y.


## 3. The Radial Basis Function Meshless Method

LetΩ be an open domain of R2. Suppose U=Ux,t is a function to be approximated in a set of N pairwise distinct nodes x=x1,x2,…,xN. In the RBF meshless scheme, the approximation of Uat the node x=xi,yi can be written as a linear combination of N RBFs:(5)Ux,t=∑j=1Nλjtφx−xj,ε,where Ux,t=Ux1,t,Ux2,t,…,UxN,tT are the function values at node x, Φx=φx−x1,ε,φx−x2,ε,…,φx−xN,εT is a RBF centered at x, . denotes the Euclidean norm between nodes x and xj, and Λ=λ1,λ2,…,λNT 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(6)φx−xj,ε=1+ε2x−xj2,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]:(7)ε=0.8Ndm,where dm denotes the smallest nodal distance. The expansion coefficients Λ in (5) are obtained by solving the following linear system of N×N algebraic equations of Uxi,t=Ui:(8)ΦΛ=U.The expansion coefficients are then calculated by(9)Λ=Φ−1U,where U=U1,U2,…,UNT is the vector of approximate solutions and Φ=φxi−xj,ε,1≤i,j≤N is an N×N matrix of RBFs given as(10)φx1−x1,ε…φx1−xN,ε⋮⋮φxN−x1,ε…φxN−xN,εN×N.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(11)∂kU∂xk=∑j=1Nλjt∂kφ∂xkxi−xj,ε,where k = [1, 2] denotes the first- and second-order derivatives. In a compact matrix form, using (9), (11) can be written as(12)∂kU∂xk=ΦxkΛ=ΦxkΦU,where(13)∂kU∂xk=∂kUx1,t∂xk,∂kUx2,t∂xk,…,∂kUxN,t∂xkT,Φxk=∂kΦx1∂xk,∂kΦx2∂xk,…,∂kΦxN∂xkT.


## 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 thatFξnWx 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(14)FξnW=∑j=1Nλξ,jtφx−xj,ε.In the matrix form, this equation becomes(15)Fξn=ΛξnΦ.However, the expansion coefficients are then calculated by(16)Λξn=Φ−1Fξn.From (14), the partial flux derivative ∂FξnW/∂ξ can be written as(17)∂FξnW∂ξ=∑j=1Nλξ,jnt∂φx−xj,ε∂ξ.Then, the compact matrix form of the partial derivative of the flux vector is(18)∂FξnW∂ξ=ΦξΛξn=ΦξΦ−1FξnW.Applying the same procedure described above for the bottom topography functionzbx, it is easy to obtain its discrete form. The topography function can be approximated by(19)zbnx=∑j=1Nαjtφx−xj,ε.Its partial derivative is then expressed as(20)∂zbnx∂ξ=∑j=1Nαjt∂φx−xj,ε∂ξ.In the same way, the compact matrix form of the bottom topography termSzbn along the axis ξ is given as(21)Szbnξ=−ghn∂zbnx,t∂ξ=−ghnΦξΦ−1zbnx,t.


### 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:(22)∂hU∂t=Sbn+1W=−gn2UUh−1/3n+1,∂hU∂t+LW,Zb=0,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]:(23)hU¯−hUnΔt=Sbn+1W,where the friction term Sbn+1 may be expressed using a Taylor series as(24)Sbn+1W=SbnW+∂Sb∂hUnhU¯−hUn+0ΔhU2.Rearranging the above equation leads to the following formula for updating water dischargehu¯ at the new time step:(25)hU¯=hUn+ΔtSbnW1+2ΔtCfUh−1.In the second step, the valuehU¯ 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 tn to the next time tn+1 is carried out as(26)W1=Wn+ΔtLWn,W2=34Wn+14W1+14ΔtLW1,Wn+1=13Wn+23W2+23ΔtLW2,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:(27)Δt=CFLdminmaxui+ghi,where CFL is the Courant number, such that 0 < CFL < 1, and dmin 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 conditionW0=WX,t=0, 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.Algorithm 1: Time integration procedure.
(1)
For
i
=
1
,
2
,
3
do
(2)
Wi=W0i
(3)
For
t
=
Δ
t
,
2
Δ
t
,
3
Δ
t
,
…
do
(4)
Fori=1,2,3do
(5)
K1=Wi+ΔtRHSWi
(6)
K2=3/4Wi+1/4K1+1/4ΔtRHSK1
(7)
Wi=1/3Wi+2/3K2+2/3ΔtRHSK2
(8)
Fork=1,2,…,NBoundarydoForallboundarypoints
(9)
Ifxk∈ΓopenThenBoundaryconditionsforinletandoutletdepending  on  the  case
(10)
Wk=Wknowed∂Wk/∂n=0
(11)
ElseBoundaryconditionsforsolid walls
(12)
huk=hvk=0&∂hk/∂n=0Algorithm 2: Calculation of flux vector and right-hand side.
(1)
F
x
1
=
W
2
,
F
x
2
=
W
2
2
/
W
1
+
0.5
g
W
1
2
,
F
x
3
=
W
2
×
W
3
/
W
1
(2)
F
y
1
=
W
3
,
F
y
2
=
W
2
×
W
3
/
W
1
,
F
y
3
=
W
3
2
/
W
1
+
0.5
g
W
1
2
(3)
For
i
=
1
,
2
,
3
do
(4)
RSHi=MxM−1Fxi+MyM−1Fyi
(5)
For
i
=
1
,
2
,
3
do
(6)
RSHi=RSHi−maxU(i)+ghi/PedminMxxM−1Wi


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

Let us assume thatFξnWx 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(14)FξnW=∑j=1Nλξ,jtφx−xj,ε.In the matrix form, this equation becomes(15)Fξn=ΛξnΦ.However, the expansion coefficients are then calculated by(16)Λξn=Φ−1Fξn.From (14), the partial flux derivative ∂FξnW/∂ξ can be written as(17)∂FξnW∂ξ=∑j=1Nλξ,jnt∂φx−xj,ε∂ξ.Then, the compact matrix form of the partial derivative of the flux vector is(18)∂FξnW∂ξ=ΦξΛξn=ΦξΦ−1FξnW.Applying the same procedure described above for the bottom topography functionzbx, it is easy to obtain its discrete form. The topography function can be approximated by(19)zbnx=∑j=1Nαjtφx−xj,ε.Its partial derivative is then expressed as(20)∂zbnx∂ξ=∑j=1Nαjt∂φx−xj,ε∂ξ.In the same way, the compact matrix form of the bottom topography termSzbn along the axis ξ is given as(21)Szbnξ=−ghn∂zbnx,t∂ξ=−ghnΦξΦ−1zbnx,t.


## 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:(22)∂hU∂t=Sbn+1W=−gn2UUh−1/3n+1,∂hU∂t+LW,Zb=0,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]:(23)hU¯−hUnΔt=Sbn+1W,where the friction term Sbn+1 may be expressed using a Taylor series as(24)Sbn+1W=SbnW+∂Sb∂hUnhU¯−hUn+0ΔhU2.Rearranging the above equation leads to the following formula for updating water dischargehu¯ at the new time step:(25)hU¯=hUn+ΔtSbnW1+2ΔtCfUh−1.In the second step, the valuehU¯ 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 tn to the next time tn+1 is carried out as(26)W1=Wn+ΔtLWn,W2=34Wn+14W1+14ΔtLW1,Wn+1=13Wn+23W2+23ΔtLW2,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:(27)Δt=CFLdminmaxui+ghi,where CFL is the Courant number, such that 0 < CFL < 1, and dmin 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 conditionW0=WX,t=0, 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.Algorithm 1: Time integration procedure.
(1)
For
i
=
1
,
2
,
3
do
(2)
Wi=W0i
(3)
For
t
=
Δ
t
,
2
Δ
t
,
3
Δ
t
,
…
do
(4)
Fori=1,2,3do
(5)
K1=Wi+ΔtRHSWi
(6)
K2=3/4Wi+1/4K1+1/4ΔtRHSK1
(7)
Wi=1/3Wi+2/3K2+2/3ΔtRHSK2
(8)
Fork=1,2,…,NBoundarydoForallboundarypoints
(9)
Ifxk∈ΓopenThenBoundaryconditionsforinletandoutletdepending  on  the  case
(10)
Wk=Wknowed∂Wk/∂n=0
(11)
ElseBoundaryconditionsforsolid walls
(12)
huk=hvk=0&∂hk/∂n=0Algorithm 2: Calculation of flux vector and right-hand side.
(1)
F
x
1
=
W
2
,
F
x
2
=
W
2
2
/
W
1
+
0.5
g
W
1
2
,
F
x
3
=
W
2
×
W
3
/
W
1
(2)
F
y
1
=
W
3
,
F
y
2
=
W
2
×
W
3
/
W
1
,
F
y
3
=
W
3
2
/
W
1
+
0.5
g
W
1
2
(3)
For
i
=
1
,
2
,
3
do
(4)
RSHi=MxM−1Fxi+MyM−1Fyi
(5)
For
i
=
1
,
2
,
3
do
(6)
RSHi=RSHi−maxU(i)+ghi/PedminMxxM−1Wi


## 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:(28)DW=Dh∂2W∂x2,where Dh is a tunable viscosity coefficient. The Peclet number is defined as follows:(29)Pe=maxUi+ghiDhdmin.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(30)DW=maxUi+ghiPedminΦxxΦW.


## 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(31)hux,0=0andhx,0=10,ifx≤10,2,otherwise.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. Figure1 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.Figure 1
Numerical and exact solutions of a dam break on a wet bed, at different times, using a grid of 800 points: (a) free surface elevation profiles; (b) profiles of water discharge per unit breadth.
(a)
(b)Figure 2
Zoom-in view of Figure 1, focusing on discontinuities in the water level.
(a)
(b)


### 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, [30–33]. 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(32)hut=0,x,y=hvt=0,x,y=0,ht=0,x,y=10,ifx≤95,5,otherwise.Figure 3
The computational domain and RBF collocation points of the 2D dam break problem.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. Figure4 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.Figure 4
The profile of water depth and velocity vectors of the 2D dam break at different output times:t = 3 s, t = 5 s, and t = 7.2 s.


### 6.2. Circular Dam Break

This second test case consists of the instantaneous breaking of a cylindrical tank of diameter 20 m (Figure5), 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 5
The computational domain, RBF collocation points, and location of the circular dam break problem.Figure6 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.Figure 6
The water depth contours and profile of the circular dam break at different output times:t = 1 s and t = 2.5 s.Figure 7
Cross section comparison between the numerical solution using the RBF meshless method and 2D FV-Roe scheme and the 1D reference solution.
(a)
(b)


### 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(33)zx,y=max0,1−18x−302+y−62,1−18x−302+y−242,1−310x−47.52+y−152.Figure 8
The computational domain, RBF collocation points, and bed contours of the dam break problem over humps.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 Figure9, 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.Figure 9
The profile of water depth and depth contours of the dam break over three humps at different output times:t = 1 s, t = 6 s, t = 12 s, t = 30 s, and t = 300 s.
(a)
(b)
(c)
(d)
(e)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.


## 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(31)hux,0=0andhx,0=10,ifx≤10,2,otherwise.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. Figure1 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.Figure 1
Numerical and exact solutions of a dam break on a wet bed, at different times, using a grid of 800 points: (a) free surface elevation profiles; (b) profiles of water discharge per unit breadth.
(a)
(b)Figure 2
Zoom-in view of Figure 1, focusing on discontinuities in the water level.
(a)
(b)


## 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, [30–33]. 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(32)hut=0,x,y=hvt=0,x,y=0,ht=0,x,y=10,ifx≤95,5,otherwise.Figure 3
The computational domain and RBF collocation points of the 2D dam break problem.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. Figure4 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.Figure 4
The profile of water depth and velocity vectors of the 2D dam break at different output times:t = 3 s, t = 5 s, and t = 7.2 s.


## 6.2. Circular Dam Break

This second test case consists of the instantaneous breaking of a cylindrical tank of diameter 20 m (Figure5), 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 5
The computational domain, RBF collocation points, and location of the circular dam break problem.Figure6 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.Figure 6
The water depth contours and profile of the circular dam break at different output times:t = 1 s and t = 2.5 s.Figure 7
Cross section comparison between the numerical solution using the RBF meshless method and 2D FV-Roe scheme and the 1D reference solution.
(a)
(b)


## 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(33)zx,y=max0,1−18x−302+y−62,1−18x−302+y−242,1−310x−47.52+y−152.Figure 8
The computational domain, RBF collocation points, and bed contours of the dam break problem over humps.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 Figure9, 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.Figure 9
The profile of water depth and depth contours of the dam break over three humps at different output times:t = 1 s, t = 6 s, t = 12 s, t = 30 s, and t = 300 s.
(a)
(b)
(c)
(d)
(e)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.



---

*Source: 4245658-2018-05-15.xml*