Discrete mixed subdomain least squares (DMSLS) meshless method with collocation points for modeling dam-break induced flows

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. doi: 10.2166/hydro.2016.116 om https://iwaponline.com/jh/article-pdf/18/4/702/390225/jh0180702.pdf 2020 Babak Fazli Malidareh (corresponding author) Seyed Abbas Hosseini Department of Civil Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran E-mail: fazli.babak@baboliau.ac.ir Ebrahim Jabbari School of Civil Engineering, Iran University of Science and Technology, Tehran, Iran

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 ). 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 () 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 fric-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 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 ).
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 ; Liu & Gu ). Over the last decade, meshless methods have increasingly been used for the solution of

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: is the vector of conserved variables or the solution vector, is the flux vector, and 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: 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, S 0 is the longitudinal channel bed slope, t is the temporal coordinate, and x is the longitudinal coordinate. The longitudinal boundary friction slope S f was estimated using the Manning resistance law: in which n is the Manning resistance coefficient and R h is the hydraulic radius.

NUMERICAL SCHEME
There are a variety of numerical techniques for approximating Equation (1)  assumed that u is a dependent variable function of x, then the MLS approximation of u x ð Þ is u h x ð Þ as shown in Figure 2 and can be defined at x as: where P x ð Þ is the vector of polynomials (basis functions) of the spatial coordinates, x T ¼ x ½ is the vector of the spatial coordinates of nodes, and m is the number of the polynomial. The basis function P x ð Þ 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), a x ð Þ is the coefficients vector given by: As shown in Figure 2, the coefficients vector a(x) in Equation (9) is a function of x i . The coefficients a x ð Þ can be obtained by minimizing the following weighted discrete L 2 norm: where n is the number of nodes in the support domain of x, for which the weight functionŴ x À x i ð Þ≠ 0, and u i is the nodal parameter of u at x ¼ x i . Equation (10)  The approximation of the unknown function can be written in the familiar form of: where N i x ð Þ denotes the shape function of node i defined as: MLS shape functions generally do not satisfy the Kronecker delta condition. Hence, u i cannot be treated like nodal values of the unknown function u h x i ð Þ ≠ u i À Á : 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 X k , the dimension of the support domain d s is determined by where α s is the dimensionless size of the support domain, and d c is the nodal spacing near the point at X Q . Generally, α s ¼ 2 ∼ 3 leads to good results for many of the problems that have been studied (Liu ; Liu & Gu ). 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 P ¼ 1; x x 2 Â Ã À Á 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: Assuming θ ¼ 0, the method is first-order accurate in time, and Equation (14) becomes: The approximate value of A at a collocation point x k can be obtained through interpolation: where n k is the number of nodal points having spatial coordinates (x k ) in their support domain.
The linearized residuals in the problem domain and its boundaries can now be defined as follows: Substituting Equation (17) into Equation (18) leads to the differential equation residual R Ω defined as: The total number of collocation points is N comprising N d internal collocation points, N t collocation points on the Neuman boundary, and N u collocation points on the Dirichlet boundary, i.e., when there is no Neuman boundary condition, then N t ¼ 0.
Now, the least squares functional of all residuals at all collocation points in inner (field) domain can be constructed as: Minimization of Equation (24) with respect to the nodal parameters u i leads to: Substituting Equation (22) into Equation (25) yields the domain system of equations: The typical components of the matrix K Ω and right-hand side vector F Ω are defined as: In the same way, the momentum part of Equation (5) becomes: The approximate value of function Q at collocation point x k can be obtained through interpolation: The typical components of the matrix K Ω and right-hand side vector F Ω in Equation (29) are defined as:

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: Substituting Equation (17) into Dirichlet boundary conditions leads to the Dirichlet boundary condition residual R u defined as: subject to Dirichlet boundary A ¼ A.

The substitution of Equation (36) into Equation (35) for
each Dirichlet subdomain leads to the boundary system of equations: where A L and A R are cross-sectional areas of the left and right boundaries, respectively.
The typical components of the matrix K u and the righthand side vector F u are defined as follows: where n u 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 K Ω , K u , F Ω , and F u .
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 Q L and Q R are discharges in left and right boundaries, respectively.
The components of matrix K u and right-hand side vector F u in Equation (42) are defined as follows: subject to Dirichlet boundary Q ¼ Q.
Matrix K and vector F are formed with the assembly process using K Ω , K u , F Ω , and F u : 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 (CFL) condition should be satisfied: where C is the celerity C ¼ ffiffiffiffiffiffi ffi gH p À Á and C t is the Courant number that should be less than 1.
To evaluate overall accuracy, the following L 2 norm is used: where H i sim and H i exact are the simulated water depth and the exact solution at the i th point, respectively. The analytical solution for dry bed on slope channel is   given as (Chanson ): where H 0 is the initial water depth of the dam break, H is the water depth, x is the distance from the dam break,  Table 1 summarizes the L 2 norms between DMSLS and MPS with respect to      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 L 2 norms between the numerical and exact solutions of different NNP and NCP, and results demonstrate that the L 2 norm decreases as NNP increases. As shown in Table 2, collocation points instead of nodal points can be increased to obtain better results.
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 m À 1 3 :s was considered. Other computational conditions were the same as in the dry bed test case discussed in the previous section. Figure    For these three cases, water depths were compared at four gauges as shown in Figures 13-15. 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 t ¼ 13 s occurs as the incident bore passes, and a second abrupt rise happens at t ¼ 35 s, because of the partial reflection of the overtopping bore from the upstream face of the bump.  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 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 con-  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