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 ﬁxed-bed dam-break ﬂow 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

*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: in which

*n*is the Manning resistance coefficient and

*R*is the hydraulic radius.

_{h}## NUMERICAL SCHEME

### Moving least squares method

*u*is a dependent variable function of , then the MLS approximation of is as shown in Figure 2 and can be defined at as: 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:

*a*(

*x*) in Equation (9) is a function of . The coefficients can be obtained by minimizing the following weighted discrete norm: 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.

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

### Mathematical formulations

*A*at a collocation point can be obtained through interpolation: where is the number of nodal points having spatial coordinates in their support domain.

*N*comprising internal collocation points, collocation points on the Neuman boundary, and collocation points on the Dirichlet boundary, i.e., when there is no Neuman boundary condition, then .

### Implementation of boundary conditions

*n*is the number of total nodal points.

*K*and vector

*F*are formed with the assembly process using . 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: where and are discharges in left and right boundaries, respectively.

*K*in Equations (41) and (46) can be seen to be symmetric and positive-definite.

## NUMERICAL CASES

*i*point, respectively.

^{th}### Idealized dam-break problem in rectangular channel

*H*and

_{up}*H*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): where is the initial water depth of the dam break,

_{down}*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).

*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 difﬁculty 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.

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

Reference | DMSLS | MPS |
---|---|---|

norm | 0.0488 | 0.0584 |

Reference | DMSLS | MPS |
---|---|---|

norm | 0.0488 | 0.0584 |

Water depth ratio | ||||||
---|---|---|---|---|---|---|

Case | No. of nodal points | No. of collocation points | CPU time (s) | L _{2} | CPU time (s) | L _{2} |

A | 101 | 201 | 6.62 | 0.083 | 7.26 | 0.306 |

B | 401 | 8.36 | 0.052 | 9.06 | 0.053 | |

C | 1,001 | 13.59 | 0.033 | 15.05 | 0.025 | |

D | 201 | 401 | 59.42 | 0.049 | 64.71 | 0.051 |

E | 1,001 | 78.09 | 0.031 | 84.77 | 0.032 | |

F | 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) | L _{2} | CPU time (s) | L _{2} |

A | 101 | 201 | 6.62 | 0.083 | 7.26 | 0.306 |

B | 401 | 8.36 | 0.052 | 9.06 | 0.053 | |

C | 1,001 | 13.59 | 0.033 | 15.05 | 0.025 | |

D | 201 | 401 | 59.42 | 0.049 | 64.71 | 0.051 |

E | 1,001 | 78.09 | 0.031 | 84.77 | 0.032 | |

F | 2,001 | 114.06 | 0.022 | 122.59 | 0.024 |

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

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

### Straight channel with a triangular bump in bed

^{⅓}. In the experiments, three cases were considered. The ﬁrst 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 (reﬂective) boundary was imposed. In all numerical computations, 153 nodal and 603 collocation points were used.

### Varying width dam-break problems

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

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