This paper presents a new meshless numerical scheme to overcome the problem of shock waves and to apply boundary conditions in cases of dam-break flows in channels with constant and variable widths. The numerical program solves shallow water equations based on the discrete mixed subdomain least squares (DMSLS) meshless method with collocation points. The DMSLS meshless method is based on the minimization of a least squares functional defined as the weighted summation of the squared residuals of the governing equations over the entire domain and requiring the summation of residual function to be zero at collocation points in boundary subdomains. The collocated discrete subdomain meshless method is applied on the boundary, whereas the collocated discrete least squares meshless technique is applied to the interior domain. The meshless scheme extends for dam-break formulation of shallow water equations. The model is verified by comparing computed results with analytical and experimental data for constant and varying width channels. The developed model is also used to study one-dimensional dam-break problems involving different flow situations by considering changes to the channel width, a bumpy channel with various downstream boundary conditions, and the effects of bed friction and bed slope as source terms on wave propagation. The accuracy of the results is acceptable.

INTRODUCTION

Dam-break flows propagate along rivers and can lead to devastating floods, damage to property, and loss of human life. Dam-break problems have been the subject of scientific research by many hydraulic scientists and engineers. Many efforts have been made to understand the mechanism of dam-break flows by conducting analytical, experimental, and numerical studies. Whereas dam break is associated with large free-surface deformation and fragmentation with the interfacial flow and mixing process, it makes modeling of the flow structure very complicated (Biscarini et al. 2010; Bellos et al. 2012; Ozmen-Cagatay & Kocaman 2014).

The most important feature of dam-break equations is that they allow discontinuities and smooth solutions. These discontinuities, often called ‘shock’ or ‘shock waves’, initiate the failure of a number of classical numerical methods. The results of dam-break modeling show that shallow water equations are logically appropriate for capturing discontinuous free surface, and methods for shallow water problems fit sufficiently with experimental and analytical results (Alcrudo & Garcia-Navarro 1993). Shallow water equations are typically used to represent the hydrodynamics and behavior of dam-break problems.

Earlier studies have been primarily based on analytical solutions of shallow water equations for idealized conditions. For example, Stoker (1957) developed an analytical solution to predict dam-break flows in an idealized channel in which the bed slope was assumed to be zero and the friction term was ignored. Chanson (2006) proposed an analytical solution of dam-break waves with flow resistance, and then it was applied to simulate tsunami surges on dry coastal plains. However, analytical solutions for shallow water equations which are nonlinear hyperbolic partial differential equations (PDEs) cannot be obtained except in special and simple cases. With the development of computer technology, hydrodynamic models based on one-dimensional (1-D) and two-dimensional (2-D) shallow water equations are increasingly being used in numerical solutions for dam-break flows. Numerous studies in the literature solved these equations using different numerical techniques, such as finite difference, finite volume, and finite element methods based on the mesh system (Bellos et al. 1991; Stansby et al. 1998; Mohamed et al. 2003; Janosi et al. 2004; Bukreev & Gusev 2005; Liang et al. 2006; Ghadimi & Reisinezhad 2012; Duricic et al. 2013; Di Francesco et al. 2015). However, mesh-based numerical methods have difficulty dealing with complicated problems such as in dam-break flows with violent free surface and interfaces. The main difficulty in using such methods is the capturing of sharp interfaces or shock waves. Furthermore, the complicated procedure needs to specify those cells containing the free surface. Also, in mesh-based methods, the problem of numerical diffusion is common due to advection terms (Zhihua 1989).

Recently, the hottest interest in the field of computational fluid mechanics has been concentrated on development of the next generation of numerical meshless methods, which are advantageous in applications with discontinuous free surface or large interfacial deformations compared to mesh-based methods. As their name implies, one common characteristic of all these methods is that they do not require the traditional mesh to construct the numerical formulation; they require node generation instead of mesh generation. In other words, there is no pre-specified connectivity or relationship among the nodes; thus, the computational costs associated with mesh generation are greatly reduced. The computational advantages of a meshless method suggest that they have the potential to solve a broad class of scientific and engineering problems (Liu 2002; Liu & Gu 2005). Over the last decade, meshless methods have increasingly been used for the solution of PDEs to solve both shallow water equations and Navier–Stokes equations. Ferrari et al. (2010), Kao & Chang (2012), and Ataie-Ashtiani et al. (2008) used the smooth particle hydrodynamic (SPH) concept to solve both shallow water equations and Navier–Stokes equations. The moving particle semi-implicit (MPS) method has been applied for simulation of fixed-bed dam-break flow on dry-bed and wet-bed downstream channels (Shakibaeinia & Jin 2010; Khayyer & Gotoh 2010). Darbani et al. (2011) used the natural element method to simulate 2-D shallow water equations in the presence of strong gradients. Despite their advantages, meshless methods have the difficulties of imposing boundary and treatment of boundary values, whereas the unreliable selection of weight functions and the complexity of algorithms for computing interpolation functions are all serious technical problems in such methods. This is in contrast to mesh-based methods which often have this property. Consequently, the imposition of boundary conditions requires certain attention in mesh-free methods and may degrade the convergence of the method (Firoozjaee & Afshar 2010a). Special techniques, such as the penalty method by Zhu & Atluri (1998), the transformation of approximate nodal values to actual nodal values by Cai & Zhu (2004), the nodal integration method by Beissel & Belytschko (1996), and the efficient computation of shape functions by Beitkopf et al. (2000) have been proposed to overcome these problems.

Recently, a family of collocation-based meshless methods has been emerging in the literature. The advantage of collocation methods is their efficiency in constructing the final system of equations because integration is not required, and shape functions are evaluated at nodal positions only. However, the accuracy and robustness of collocation approaches are weak points, especially if the approximation is based on a set of randomly scattered points. Collocation methods may lead to large numerical errors in these cases and involve numerical stability problems. An interesting discussion of these aspects may be found in Fries & Matthies (2004), where certain conditions are proposed which are often not met by standard collocation mesh-less methods. Furthermore, boundary conditions are an issue for collocation methods. The definition and location of the boundary surface may not be easy, e.g., for free surface flow, and methods of applying boundary conditions are not always straightforward. Firoozjaee & Afshar (2010b) proposed the collocated discrete least squares meshless (CDLSM) method to solve steady-state incompressible Navier–Stokes equations. Shobeyri & Afshar (2009) discussed simulating free surface problems using the discrete least squares meshless method. Afshar et al. (2009) used the CDLSM method to solve transient and steady-state problems.

In CDLSM, the penalty method is used and can only approximately satisfy the boundary conditions. Accuracy is affected by selection of the penalty coefficient; these boundary conditions are only approximately imposed, depending on the magnitude of the penalty coefficients. Theoretically, the larger the penalty coefficients, the more accurate the enforcement of boundary conditions will be. It is difficult to choose a set of penalty factors that are universally applicable for all kinds of problems. One hopes to use the largest possible penalty factors, but penalty factors that are too large often give numerical problems, and trials may be needed to choose a proper penalty factor (Liu 2002; Liu & Gu 2005).

In this study, the authors propose a new method named collocated discrete subdomain meshes (CDSM) method on the boundary for more accuracy as explained in the section ‘Implementation of boundary conditions’. In this approach, the residual functional can be written in summation form and minimized at collocation points; the sum of the residuals must be zero for each boundary subdomain. It is noteworthy that the boundary subdomains are usually centered by nodal points on the boundary.

This study investigated the discrete mixed subdomain least squares (DMSLS) meshless method using collocation points. To generalize this method, an attempt was made to check the balance of this approach for some dam-break flow simulations (where the definition and location of the boundary surface is not easy) for different bed and geometry conditions with maximum variables in boundary conditions (where applying boundary conditions are not straightforward). To verify the prediction results, this technique was tested for different experimental dam-break flows considering bed friction, bed slope, bumpy channels with various downstream conditions and varying width channels, and well-known analytical solutions for rectangular channels of constant width. Finally, the results were summarized and conclusions were drawn. The results of this study illustrate the potential ability of the DMSLS method to reproduce the detailed features of flow structure for highly transient flow with the large geometrical changes of the computational domain.

GOVERNING EQUATIONS

The shallow water theory is based on the hypothesis that a layer of water flows over a horizontal, flat surface with elevation Z (x, y, t). This means that the pressure distribution along each vertical is hydrostatic. If it is assumed that the horizontal scale of flow features is large compared to the depth of water, the flow velocity is independent of depth and the water within the layer is in hydrostatic balance. Thus, 1-D continuity and momentum equations, which were modified to provide a conservation formulation applicable to rectangular channels of varying width, which can be written in the following vector form as: 
formula
1
where 
formula
2
is the vector of conserved variables or the solution vector, 
formula
3
is the flux vector, and 
formula
4
represents the source term vector.
The non-conservation formulation of the 1-D Saint-Venant equations can be written in the following vector form as: 
formula
5
in which 
formula
6
is known as the convection matrix where g is the acceleration due to gravity, A is the cross-sectional area perpendicular to the flow, Q is the discharge, H is the water depth, U is the averaged longitudinal flow velocity, B is the width of the rectangular channel, is the longitudinal channel bed slope, t is the temporal coordinate, and is the longitudinal coordinate. The longitudinal boundary friction slope was estimated using the Manning resistance law: 
formula
7
in which n is the Manning resistance coefficient and Rh is the hydraulic radius.

NUMERICAL SCHEME

There are a variety of numerical techniques for approximating Equation (1) or shallow water equations, e.g., finite difference methods, finite element methods, finite volume methods, etc. As mentioned in the Introduction, the DMSLS meshless method with collocation points is used to solve Equation (5). The DMSLS meshless method was extended for solving 1-D dam-break problems, and its performance for solving transient problems on regular distribution of nodes was investigated. In this method, the least squares functional and residual functions are formed at the collocation points. The aim of the mixed subdomain-least squares meshless method is to find an approximate solution which minimizes the least squares function over the entire domain and causes the residual function to be zero in their boundary conditions at nodal points. As shown in Figure 1, the physical domain and boundaries are discretized using two sets of nodes called nodal points and collocation points. Collocation points are used in the physical domain and on its boundaries to make the least squares functional and subdomain residual function.
Figure 1

Local support domains in the meshless method.

Figure 1

Local support domains in the meshless method.

Moving least squares method

Several techniques have been developed to construct meshless shape functions. The moving least squares (MLS) approximation, the radial point interpolation method, and kriging interpolation are only some examples of current methods. Among these methods, the MLS method has gained more popularity in the meshless community. If it is assumed that u is a dependent variable function of , then the MLS approximation of is as shown in Figure 2 and can be defined at as: 
formula
8
where is the vector of polynomials (basis functions) of the spatial coordinates, is the vector of the spatial coordinates of nodes, and m is the number of the polynomial. The basis function is often built using monomials from the Pascal triangle to ensure minimum completeness. In some special problems, enhancement functions can, however, be added to the basis to improve the performance of the MLS approximation. In this paper, only pure polynomial basis was used. In Equation (8), is the coefficients vector given by: 
formula
9
Figure 2

Approximate function and nodal parameters in the MLS approximation.

Figure 2

Approximate function and nodal parameters in the MLS approximation.

As shown in Figure 2, the coefficients vector a(x) in Equation (9) is a function of . The coefficients can be obtained by minimizing the following weighted discrete norm: 
formula
10
where n is the number of nodes in the support domain of x, for which the weight function , and is the nodal parameter of u at . Equation (10) is a functional, weighted residual that is constructed using the approximated values and the nodal parameters of the unknown field function. Because the number of nodes used in the MLS approximation is usually much larger than the number of unknown coefficients , the approximated function does not pass through nodal values. In this research, the cubic spline weight function explained in more details by Liu (2002) and Liu & Gu (2005) was considered.
The approximation of the unknown function can be written in the familiar form of: 
formula
11
where denotes the shape function of node defined as: 
formula
12

MLS shape functions generally do not satisfy the Kronecker delta condition. Hence, cannot be treated like nodal values of the unknown function

The accuracy of interpolation for the point of interest depends on the nodes in the support domain, as shown in Figure 1. Therefore, a suitable support domain should be chosen to ensure an efficient and accurate approximation. For a point of interest at , the dimension of the support domain is determined by 
formula
13
where is the dimensionless size of the support domain, and is the nodal spacing near the point at . Generally, leads to good results for many of the problems that have been studied (Liu 2002; Liu & Gu 2005). It should be noted that the support domain is usually centered by a point of interest. It should also be noted that the polynomial basis of the second order was used in this study for 1-D problems to construct MLS shape functions, such as 1-D dam-break cases. Therefore, the accuracy of this approach is the second order in space.

Mathematical formulations

Using a semi-discretization approach, Equation (5) was discretized by the method in time as follows: 
formula
14
 
formula
15
 
formula
Assuming , the method is first-order accurate in time, and Equation (14) becomes: 
formula
16
The approximate value of A at a collocation point can be obtained through interpolation: 
formula
17
where is the number of nodal points having spatial coordinates in their support domain.
The linearized residuals in the problem domain and its boundaries can now be defined as follows: 
formula
18
Substituting Equation (17) into Equation (18) leads to the differential equation residual defined as: 
formula
19
 
formula
20
 
formula
21
 
formula
22
The total number of collocation points is N comprising internal collocation points, collocation points on the Neuman boundary, and collocation points on the Dirichlet boundary, i.e., 
formula
23
when there is no Neuman boundary condition, then .
Now, the least squares functional of all residuals at all collocation points in inner (field) domain can be constructed as: 
formula
24
Minimization of Equation (24) with respect to the nodal parameters leads to: 
formula
25
Substituting Equation (22) into Equation (25) yields the domain system of equations: 
formula
26
The typical components of the matrix and right-hand side vector are defined as: 
formula
27
 
formula
28
In the same way, the momentum part of Equation (5) becomes: 
formula
29
The approximate value of function Q at collocation point can be obtained through interpolation: 
formula
30
The typical components of the matrix and right-hand side vector in Equation (29) are defined as: 
formula
31
 
formula
32
 
formula
33
 
formula
34

Implementation of boundary conditions

This paper introduces the collocated subdomain meshless method to implement boundary conditions accurately. In this method, the summation of residuals must be zero for each boundary subdomain. It is noteworthy that the boundary subdomains are usually centered by nodal points on the boundary. In this paper, the boundary conditions are Dirichlet boundary conditions; therefore: 
formula
35
where is residual at collocation points inside the Dirichlet subdomain. Each Dirichlet subdomain is usually centered by a nodal point on the Dirichlet boundary, and is the number of collocation points inside each Dirichlet subdomain.
Substituting Equation (17) into Dirichlet boundary conditions leads to the Dirichlet boundary condition residual defined as: 
formula
36
subject to Dirichlet boundary .
The substitution of Equation (36) into Equation (35) for each Dirichlet subdomain leads to the boundary system of equations: 
formula
37
 
formula
38
where and are cross-sectional areas of the left and right boundaries, respectively.
The typical components of the matrix and the right-hand side vector are defined as follows: 
formula
39
 
formula
40
where is number of nodal points on the Dirichlet boundary, and n is the number of total nodal points.
Matrix K and vector F are formed with the assembly process using . 
formula
41
The assembly matrix of this method is similar to the assembly process in the finite element method. In the same way, the system of equations for discharge in the Dirichlet boundary becomes: 
formula
42
 
formula
43
where and are discharges in left and right boundaries, respectively.
The components of matrix and right-hand side vector in Equation (42) are defined as follows: 
formula
44
 
formula
45
subject to Dirichlet boundary .
Matrix K and vector F are formed with the assembly process using : 
formula
46
Now, the system of equations can be solved at each time step in a time-marching manner. The stiffness matrix K in Equations (41) and (46) can be seen to be symmetric and positive-definite.

NUMERICAL CASES

In this section, the numerical accuracy and applicability of the proposed DMSLS method are further examined by analytical and laboratory results of dam-break flows in open channels, including dam-break flows through an idealized and rough, flat channel, and a rough bumpy channel with various downstream boundary conditions and varying width channel. These validation cases encompass a number of flow phenomena in open channels, such as shock discontinuities, shock front motion, hydraulic jumps, dry/wet bed flow, contraction flow, overtopping flow, and partial reflections. For the numerical stability of the proposed scheme, the Courant–Friedriches–Lewy (CFL) condition should be satisfied: 
formula
47
where is the celerity is the Courant number that should be less than 1.
To evaluate overall accuracy, the following norm is used: 
formula
48
where and are the simulated water depth and the exact solution at the ith point, respectively.

Idealized dam-break problem in rectangular channel

The first numerical test is the 1-D idealized dam-break problem. The idealized channel is rectangular, frictionless, and horizontal, where Hup and Hdown are initial water depths upstream and downstream of the dam, respectively. The channel is 1,000 m in length, and a dam is located in its middle. At the beginning , the dam is broken instantaneously. A positive wave propagating downstream and a depression wave traveling upstream can be observed. The initial upstream water depth of the dam is 10 m, and the corresponding initial downstream water depths are given as 0.0 m (dry bed) and 2.0 m (wet bed) for each case. The analytical solution for dry bed on slope channel is given as (Chanson 2006): 
formula
49
where is the initial water depth of the dam break, H is the water depth, is the distance from the dam break, and t is the time. More details about the analytical solution on a wet, friction bed can be found in the literature (Stoker 1957; Chanson 2006).
In all the numerical computations, a uniform mesh of 201 nodal points with 801 distributed collocation points was used. The depths and velocities at t = 23 s after the dam failure for dry and wet bed were modeled with proposed methods and compared with the analytical solutions in Figures 3 and 4, respectively. Figure 3(b) reveals that solutions of velocity are evidently underestimated near the wave front edge, but the velocity distribution by the DMSLS method is more accurate than CDLSM. In fact, most of the numerical schemes have difficulty in accurately predicting velocity near the wave front edge for dry bed dam-break problem because of the low water depth in the area (Ying et al. 2004). There is little discrepancy in the water surface profiles for these two approaches, and they are approximately compatible together.
Figure 3

Numerical solutions of (a) water depths and (b) velocities for frictionless, dry bed.

Figure 3

Numerical solutions of (a) water depths and (b) velocities for frictionless, dry bed.

Figure 4

Numerical solutions of (a) water depths and (b) velocities for frictionless, non-dry bed.

Figure 4

Numerical solutions of (a) water depths and (b) velocities for frictionless, non-dry bed.

To study the effect of slope on the accuracy of DMSLS, the dam break was modeled on a sloping, frictionless bed as shown in Figure 5. The comparison between numerical and analytical solutions shows an acceptable agreement.
Figure 5

Numerical solutions of water depths for dam break on sloping bed.

Figure 5

Numerical solutions of water depths for dam break on sloping bed.

To examine the accuracy of the model compared with other meshless methods, an idealized dam break on a dry bed was modeled by these methods and compared with MPS and experimental observations (shown in Figure 6) (Fu & Jin 2014). Experiments were conducted by LaRocque et al. (2013) in a smooth wooden flume with a length of , height of , width of , a bottom slope of 0.93%, initial length of water column , and initial heights of water column upstream , respectively. Table 1 summarizes the norms between DMSLS and MPS with respect to experimental data. As can be seen, the results of this new approach (DMSLS) are in good agreement with the experimental data and have greater accuracy compared with the MPS method.
Table 1

The norms between the simulated water depths and exact solutions

Reference DMSLS MPS 
norm 0.0488 0.0584 
Reference DMSLS MPS 
norm 0.0488 0.0584 
Figure 6

Comparison of DMSLS meshless method with MPS and experimental data (Ying et al. 2004).

Figure 6

Comparison of DMSLS meshless method with MPS and experimental data (Ying et al. 2004).

To study the effect of number of nodal and collocation points on the accuracy and point dependency of the proposed method, a sensitivity analysis was performed for the DMSLS meshless method. For this purpose, two different numbers of nodal points (NNP), NNP = 101, 201 with three different numbers of collocation points (NCP) were considered. The non-dry bed dam-break problem with NNP = 201 and a water depth ratio equal to 0.2 is evaluated in Figure 7. Good agreement of the wave-front tracking process between the numerical and analytical solutions was obtained using NCP = 1,001, 2,001. Table 2 summarizes the norms between the numerical and exact solutions of different NNP and NCP, and results demonstrate that the norm decreases as NNP increases. As shown in Table 2, collocation points instead of nodal points can be increased to obtain better results.
Table 2

The norms between the simulated and exact water depths for different numbers of nodal points (NNP) and numbers of collocation points (NCP) and water depth ratios

      Water depth ratio
 

 

 
Case No. of nodal points No. of collocation points CPU time (s) L2 CPU time (s) L2 
101 201 6.62 0.083 7.26 0.306 
401 8.36 0.052 9.06 0.053 
1,001 13.59 0.033 15.05 0.025 
201 401 59.42 0.049 64.71 0.051 
1,001 78.09 0.031 84.77 0.032 
2,001 114.06 0.022 122.59 0.024 
      Water depth ratio
 

 

 
Case No. of nodal points No. of collocation points CPU time (s) L2 CPU time (s) L2 
101 201 6.62 0.083 7.26 0.306 
401 8.36 0.052 9.06 0.053 
1,001 13.59 0.033 15.05 0.025 
201 401 59.42 0.049 64.71 0.051 
1,001 78.09 0.031 84.77 0.032 
2,001 114.06 0.022 122.59 0.024 
Figure 7

Numerical solutions of water depths for various numbers of collocation points under water depth ratio of 0.2.

Figure 7

Numerical solutions of water depths for various numbers of collocation points under water depth ratio of 0.2.

It is notable that increasing nodal points to obtain better results yields an increase in the dimension of the coefficient matrix, but increasing the number of collocation points yields better results without increasing the dimension of the coefficient matrix. Thus, this method can be easily used in practical cases with high efficiency and less computational effort, and it can simulate dam-break problems in real size (in the field) rather than on a small scale with good accuracy and less simulation time.

Real dam-break problem in a rectangular channel

In this test case, a real dam break with a Manning roughness coefficient n = 0.03 was considered. Other computational conditions were the same as in the dry bed test case discussed in the previous section. Figure 8 represents the numerical and analytical solution of water depth H at 23 s after the dam failure. Good agreement between numerical and analytical solutions can be observed. To study the effect of bed friction on wave propagation, downstream channels were considered to be wet, and the initial downstream water depth was 2.0 m. The results of the wave propagated downstream are shown in Figure 9. It can be seen that the effect of the bed friction on dam break with a wet bed is negligible, and it can be neglected in numerical modeling. This result is in accordance with previous results published by other researchers (Chanson 2006) and can be a suitable case for DMSLS verification.
Figure 8

Numerical solutions of water depths for frictional, dry bed.

Figure 8

Numerical solutions of water depths for frictional, dry bed.

Figure 9

Comparison of idealized and real numerical solutions of water depths for non-dry bed.

Figure 9

Comparison of idealized and real numerical solutions of water depths for non-dry bed.

Straight channel with a triangular bump in bed

The third test case selected from the CADAM studies was a dam-break flow into a channel with a triangular bump in the bed (Morris 2000). The channel was 38 m long and 0.75 m wide. A gate was located 15.5 m from the upstream end of the channel at the exit of a reservoir, and a triangular bump 0.4 m high and 6 m long was located 13 m downstream of the gate (Figure 10). The initial water depth at the reservoir was 0.75 m. Manning's friction coefficient for the channel bed was 0.0125 s/m. In the experiments, three cases were considered. The first case (case I) was a dry bed with open conditions at the downstream end, in which the water depth in the downstream channel was zero; two wet bed cases were then considered. In both of them, the depth of water in the channel after the bump was set to 0.15 m, but two different boundary conditions were considered at the end of the channel: in case II an open (transmissive) condition was applied, and in case III a closed (reflective) boundary was imposed. In all numerical computations, 153 nodal and 603 collocation points were used.
Figure 10

Straight channel with a triangular bump in the bed. Side elevation showing gauge positions indicated with a suffix number after G which means distance from the gate, i.e., G4 indicates that it is 4 m away from the gate (Morris 2000).

Figure 10

Straight channel with a triangular bump in the bed. Side elevation showing gauge positions indicated with a suffix number after G which means distance from the gate, i.e., G4 indicates that it is 4 m away from the gate (Morris 2000).

In all three cases, the flow initially generated a bore wave which propagated along the channel. Upon reaching the bump, the bore ran up, overtopped the bump, and then accelerated downstream. During the run up, the bore was partially reflected and generated a small wave which propagated upstream. After the initial wave overtopped the bump, it began to act like a weir with supercritical flow over the downstream face, causing a hydraulic break. In the wet-bed cases, this supercritical flow caused a hydraulic jump to form (Figure 11). When the end of the channel was closed, this bore wave was reflected from the end wall and propagated upstream (Figure 12). The results indicate that the flow is potentially very complex, involving multiple wave interactions, thus providing a demanding test for verification and validation of a numerical method.
Figure 11

Straight channel with a triangular bump in the bed, wet bed with open end. Water surface profile at s.

Figure 11

Straight channel with a triangular bump in the bed, wet bed with open end. Water surface profile at s.

Figure 12

Straight channel with a triangular bump in the bed, wet bed with closed end. Water surface profile at s.

Figure 12

Straight channel with a triangular bump in the bed, wet bed with closed end. Water surface profile at s.

For these three cases, water depths were compared at four gauges as shown in Figures 1315. In almost all cases, there was a good agreement between numerical and experimental results. Figure 13 shows the comparison of CDLSM, DMSLS, and experimental data for case I (dry bed). As seen in Figure 13, the DMSLS method provides more accurate results compared with the CDLSM approach. As can also be observed in Figures 14 and 15, near the reservoir (at gauge 4), the sudden rise in depth at s occurs as the incident bore passes, and a second abrupt rise happens at s, because of the partial reflection of the overtopping bore from the upstream face of the bump. Results at gauge 4 are virtually identical in all three cases over the plotted time since the effects of different downstream conditions are separated from this portion of the flow domain by a hydraulic break. The results at gauge 10, which is located at the upstream toe of the bump, show similar characteristics to those at gauge 4. Numerical results for gauge 13 at the apex of the bump follow those of experiments, but they consistently underestimate the water depth. This may be due to the lack of vertical acceleration terms in the shallow water model and the local discontinuity in the bed slope.
Figure 13

Straight channel with a triangular bump in the bed, dry bed with an open end. Comparisons of experimental and numerical water depths at selected gauges.

Figure 13

Straight channel with a triangular bump in the bed, dry bed with an open end. Comparisons of experimental and numerical water depths at selected gauges.

Figure 14

Straight channel with a triangular bump in the bed, wet bed with an open end. Comparisons of experimental and numerical water depths at selected gauges.

Figure 14

Straight channel with a triangular bump in the bed, wet bed with an open end. Comparisons of experimental and numerical water depths at selected gauges.

Figure 15

Straight channel with a triangular bump in the bed, wet bed with a closed end. Comparisons of experimental and numerical water depths at selected gauges.

Figure 15

Straight channel with a triangular bump in the bed, wet bed with a closed end. Comparisons of experimental and numerical water depths at selected gauges.

Varying width dam-break problems

To further verify the proposed meshless scheme, dam-break surge propagation in a rectangular channel with varying width was considered. This case is important, because different flow regimes can occur. Analytical solutions for this case are not readily obtainable. However, the experimental results of Bellos et al. (2012) involving dam-break flows in a non-prismatic channel provide suitable data for the verification of the proposed meshless scheme. Figure 16 shows the geometry of the flume used by Bellos et al. (Hicks et al. 1997). The rectangular flume was made of steel and glass with a gate fixed at the section of minimum width. During each dam-break simulation, water levels were monitored at eight different locations along the center line of the flume as shown in Figure 16. The data, which were collected from five of these observational stations (x = 0.0 m, 4.5 m, 8.5 m, 13.5 m, and 18.5 m), are used for the purpose of this comparison. Figures 17 and 18 show the comparison of the DMSLS meshless method results with the experimental data for (upstream depth and downstream depth ) and ( and ). The results are presented in five locations along the channel where data were provided by Bellos et al. As Figures 17 and 18 illustrate, this numerical scheme provided accurate and comparable water level predictions. Only numerical results differ from the experimental measurements at x = 18.5 m after 12 sec. As data describing the weir dimensions at the downstream end of the experimental channel were not available, the downstream boundary condition was assumed as the constant water depth. As a result, the numerical modeling was unable to predict waves reflected off the downstream boundary of the channel. It seems the numerical solution underestimated the stage at x = 8.5 m (the point of maximum plan flow curvature) and may also reflect 2-D effects that are not considered in 1-D modeling.
Figure 16

Top view of experimental set-up of Bellos et al. (2012).

Figure 16

Top view of experimental set-up of Bellos et al. (2012).

Figure 17

Comparison of the DMSLS solution with experimental data for : (a) x = 0.0 m, (b) x = 4.5 m, (c) x = 8.5 m, (d) x = 13.5 m, (e) x = 18.5 m.

Figure 17

Comparison of the DMSLS solution with experimental data for : (a) x = 0.0 m, (b) x = 4.5 m, (c) x = 8.5 m, (d) x = 13.5 m, (e) x = 18.5 m.

Figure 18

Comparison of the DMSLS solution with experimental data for : (a) x = 0.0 m, (b) x = 4.5 m, (c) x = 8.5 m, (d) x = 13.5 m, (e) x = 18.5 m.

Figure 18

Comparison of the DMSLS solution with experimental data for : (a) x = 0.0 m, (b) x = 4.5 m, (c) x = 8.5 m, (d) x = 13.5 m, (e) x = 18.5 m.

CONCLUSIONS

The moving least squares (meshless) scheme has been developed to solve the 1-D shallow water equations for the computation of dam-break induced flows. A key feature of the model is the use of collocated points to minimize the sum of the squared residual functional. The problem domain is discretized using some nodes. Then, using the MLS shape function with a least squares technique, an asymmetric stiffness matrix is constructed. Sampling (collocated) points are defined to minimize the sum of the squared residual functional over interior domain and to be zero residual function at their boundaries, and field nodes (nodal points) are defined to construct MLS shape functions. Considering the number of sampling (collocated) points to be more than the number of field nodes will improve the results and guarantee stability.

Nonlinear 1-D dam-break problems have been solved using the newly proposed CDSM method for boundary conditions and the CDLSM method for interior domain, and the results have been presented and compared to the analytical and experimental results. The variety of dam-break problems involving different flow regimes by changing the width of the channel, considering the effects of bed friction and bed slope as source terms on wave propagation, and bumpy channel with reflective boundary were modeled. The results clearly indicate that the proposed DMSLS method with collocation points for boundary conditions is capable of obtaining stable and more accurate results in dam-break problems, especially for wet/dry interfaces, different bed geometry conditions, and shock and hydraulic break waves.

The formulation presented in this paper offers a number of advantages over many existing meshless methods used to solve PDEs. In CDLSM, the penalty method has been used for boundary conditions; these boundary conditions were imposed only approximately, depending on the magnitude of the penalty coefficients. It is difficult to choose a set of penalty factors that are universally applicable for all kinds of problems. One hopes to use the largest possible penalty factors, but penalty factors that are too large often give numerical problems, and trials may be needed to choose a proper penalty factor. Thus, this paper introduces the CDSM method to enforce boundary conditions in moving least square mesh-free methods. Another advantage of this method is the use of collocated points. It is notable that increasing nodal points for obtaining better results yields to an increase in the dimension of the stiffness matrix, but increasing the number of collocation points yields better results without increasing the dimension of the stiffness matrix. This issue prevents a large stiffness matrix, increases the efficiency of the scheme, and consumes less time to solve the system of equations. These advantages can help researchers simulate flood events in the field induced from prototype dam break. This method can also be used for a moving bed such as dam-break modeling over an erodible bed, an issue currently under research.

REFERENCES

REFERENCES
Afshar
M. H.
Lashckarbolok
M.
Shobeyri
G.
2009
Collocated discrete least squares meshless (CDLSM) method for the solution of transient and steady-state hyperbolic problems
.
International Journal for Numerical Methods in Fluids
60
,
1055
1078
.
Alcrudo
F.
Garcia-Navarro
P.
1993
A high-resolution Godunov type scheme in finite volumes for the two-dimensional shallow water equations
.
International Journal for Numerical Methods in Fluids
16
,
489
505
.
Ataie-Ashtiani
B.
Shobeiry
G.
Farhadi
L.
2008
Modified incompressible SPH method for simulating free surface problems
.
Fluid Dynamics Research
40
(
9
),
637
661
.
Beissel
S.
Belytschko
T.
1996
Nodal integration of the element-free Galerkin method
.
Computer Methods in Applied Mechanics & Engineering
139
,
49
74
.
Beitkopf
P.
Rassineux
A.
Gibert
T.
2000
Explicit form and efficient computation of MLS shape functions and their derivatives
.
International Journal for Numerical Methods in Engineering
48
,
451
466
.
Bellos
C. V.
Soulis
J. V.
Sakkas
J. G.
1991
Computation of two dimensional dam-break induced flows
.
Advances in Water Resources
14
(
1
),
31
41
.
Bellos
C. V.
Soulis
J. V.
Sakkas
J. G.
2012
Experimental investigation of two-dimensional dam-break induced flows
.
Journal of Hydraulic Research
30
(
1
),
47
63
.
Biscarini
C.
Di Francesco
S.
Manciola
P.
2010
CFD Modeling approach for dam break flow studies
.
Hydrology and Earth System Sciences
14
,
705
718
.
Bukreev
V. I.
Gusev
A. V.
2005
Initial stage of the generation of dam-break waves
.
Doklady Physics
50
(
4
),
200
203
.
Chanson
H.
2006
Analytical solutions of laminar and turbulent dam break wave
.
The University of Queensland
,
Brisbane
,
Australia
,
Taylor & Francis Group, London
.
Darbani
M.
Ouahsine
A.
Villon
P.
Naceur
H.
Smaoui
H.
2011
Meshless method for shallow water equations with free surface flow
.
Applied Mathematics and Computation
217
,
5113
5124
.
Di Francesco
S.
Biscarini
C.
Manciola
P.
2015
Numerical simulation of water free-surface flows through a front-tracking lattice Boltzmann approach
.
Journal of Hydroinformatics
17
(
1
),
1
6
.
Duricic
J.
Erdik
T.
Gelder
P. V.
2013
Predicting peak breach discharge due to embankment dam failure
.
Journal of Hydroinformatics
15
(
4
),
1361
1376
.
Ferrari
A.
Fraccarollo
L.
Dumbser
M.
Toro
E. F.
Armanini
A.
2010
Three-dimensional flow evolution after a dam break
.
Journal of Fluid Mechanics
663
,
456
477
.
Firoozjaee
A. R.
Afshar
M. H.
2010a
Solution of convection-dominated problems on irregular meshes by collocated discrete least squares mesh-less (CDLSM) method
.
Scientia Iranica
16
(
4
),
340
350
.
Firoozjaee
A. R.
Afshar
M. H.
2010b
Steady-state solution of incompressible Navier–Stokes equations using discrete least-squares meshless method
.
International Journal for Numerical Methods in Fluids
67
,
369
382
.
Fries
T. P.
Matthies
H. G.
2004
Classification and overview of meshfree methods
.
Technical Report Informatikbericht Nr 2003-3
,
Technical University Braunschweig
,
Germany
.
Fu
L.
Jin
Y. C.
2014
Simulating velocity distribution of dam breaks with the particle method
.
ASCE Journal of Hydraulic Engineering
140
(
10
),
04014048
.
Hicks
F. E.
Steffler
P. M.
Yasmin
N.
1997
One dimensional dam break solutions for variable width channel
.
Journal of Hydraulic Engineering
123
(
5
),
464
468
.
Janosi
I. M.
Jan
D.
Szabo
K. G.
Tel
T.
2004
Turbulent drag reduction in dam-break flows
.
Experiments in Fluids
37
(
2
),
219
229
.
Khayyer
A.
Gotoh
H.
2010
On particle-based simulation of a dam break over a wet bed
.
Journal of Hydraulic Research
48
(
2
),
238
249
.
LaRocque
L.
Imran
J.
Chaudhry
M.
2013
Experimental and numerical investigations of two-dimensional dam-break flows
.
Journal of Hydraulic Engineering
139
(
6
),
569
579
.
Liang
D.
Falconer
R. A.
Lin
B.
2006
Comparison between TVD-MacCormack and ADI-type solvers of the shallow water equations
.
Advances in Water Resources
29
(
12
),
1833
1845
.
Liu
G. R.
2002
Mesh Free Methods: Moving beyond the Finite Element Method
,
1st edn
.
CRC Press
,
Boca Raton, FL
,
USA
.
Liu
G. R.
Gu
Y. T.
2005
An Introduction to Meshless Methods and their Programming
,
1st edn
.
Springer
,
Berlin
,
Germany
.
Mohamed
S.
Ghidaoui
M. S.
Li
N.
2003
Generalized Boltzmann equation for shallow water flows
.
Journal of Hydroinformatics
5
(
1
),
1
10
.
Morris
M.
2000
CADAM: Concerted action on dam-break modeling–Final report
.
Rep. No. SR 571
,
HR Wallingford
,
Wallingford
,
UK
.
Ozmen-Cagatay
H.
Kocaman
S.
2014
Dam-break flow in the presence of obstacle: experiment and CFD simulation
.
Engineering Applications of Computational Fluid Mechanics
5
,
541
552
.
Shakibaeinia
A.
Jin
Y. C.
2010
A weakly compressible MPS method for simulation open-boundary free-surface flow
.
International Journal for Numerical Methods in Fluids
63
(
10
),
1208
1232
.
Stansby
P. K.
Chegini
A.
Barnes
T. C. D.
1998
The initial stages of dam-break flow
.
Journal of Fluid Mechanics
374
,
407
424
.
Stoker
J. J.
1957
Water Waves: The Mathematical Theory with Applications
.
Inter Science
,
New York
,
USA
.
Ying
X.
Khan
A.
Wang
S.
2004
Upwind conservative scheme for the Saint Venant equations
.
Journal of Hydraulic Engineering
130
(
10
),
977
987
.
Zhihua
O.
1989
An adaptive characteristics method for advective-diffusive transport
.
Applied Mathematical Modelling
13
(
12
),
682
692
.