https://orcid.org/0000-0001-6365-6849
Abstract
The Saint-Venant equations are numerically solved to simulate free surface flows in one dimension. A Riemann solver is needed to compute the numerical flux for capturing shocks and flow discontinuities occurring in flow situations such as hydraulic jump, dam-break wave propagation, or bore wave propagation. A Riemann solver that captures shocks and flow discontinuities is not yet reported to be implemented within the framework of a meshless method for solving the Saint-Venant equations. Therefore, a wide range of free surface flow problems cannot be simulated by the available meshless methods. In this study, a shock-capturing meshless method is proposed for simulating one-dimensional (1D) flows on a highly variable topography. The Harten–Lax–van Leer Riemann solver is used for computing the convective flux in the proposed meshless method. Spatial derivatives in the Saint-Venant equations and the reconstruction of conservative variables for flux terms are computed using a weighted least square approximation. The proposed method is tested for various numerically challenging problems and laboratory experiments on different flow regimes. The proposed highly accurate shock-capturing meshless method has the potential to be extended to solve the two-dimensional (2D) shallow water equations without any mesh requirements.
HIGHLIGHTS
A shock-capturing meshless method is presented for solving the 1D Saint-Venant equations on a highly variable topography.
The points are irregularly distributed along the channel to define minute topographical features.
The meshless method can accurately capture shocks and flow discontinuity in open channel flows.
INTRODUCTION
The free surface flow resulting from dam-break, levee breaching, or tidal bore propagation falls in the category of transcritical flow, wherein the flow regime changes from subcritical to supercritical and vice versa. The flow transition can be accurately described by a shock-capturing numerical scheme. Such flows may be treated as one-dimensional (1D) or two-dimensional (2D), depending on the flow dynamics and practical requirements. Although 2D flow models have gained tremendous popularity in the last two decades due to advancements in computing power and robust numerical schemes, 1D models are still in use for simulating flows dominant in one direction, especially for large computational domains. Moreover, the coupled 1D–2D models (Ghostine et al. 2015) are nowadays suitably used to reduce overall computation time and ensure enough lead time. This assists emergency action planning in cases of disasters such as dam-break flow, bore and flood wave propagation, and discontinuous flow (i.e., flows with shocks) that travels downstream of hydraulic structures. The Saint-Venant equations (i.e., shallow water equations in unidirectional form) are numerically solved to simulate 1D flows with discontinuity and shocks. Several numerical schemes have been developed for solving the Saint-Venant equations, such as the finite difference method (FDM) (Fennema & Chaudhry 1990; Molls & Chaudhry 1995), the finite element method (FEM) (Liang et al. 2008) and the finite volume method (FVM) (Eymard et al. 2000; LeVeque 2002; Yoon & Kang 2004; George 2008; Kuiry et al. 2008; Haleem et al. 2015). However, meshless methods, alias gridless or meshfree, or particle methods, are also becoming increasingly popular in recent years for simulating fluid flow problems. Several meshless methods have been developed, such as the smoothed particle hydrodynamics (SPH) (Gingold & Monaghan 1977), the finite pointset method (FPM) (Tiwari & Kuhnert 2003a, 2003b), the element-free Galerkin (EFG) method (Belytschko et al. 1995), the generalized finite difference method (GFDM) (Li & Fan 2017), and the discrete mixed subdomain least squares method (Fazli Malidareh et al. 2016). However, meshless methods are relatively new to free surface flow simulations. In particular, shock-capturing meshless methods for simulating flows with shocks and discontinuity are scarce in the literature.
The oldest and the most popularly used mesh-free method is the SPH. Although SPH was initially developed for astrophysical fluid dynamics (Gingold & Monaghan 1977), it was later extended to hydrodynamic modeling by Monaghan (1992) for solving the Navier–Stokes equations for incompressible free surface flows. The numerical solution of 1D shallow water equations for open channel flow by SPH was described by Chang et al. (2011). The main difficulty of the SPH is the incorporation of boundary conditions, which are necessary for solving the Saint-Venant equations that describe real-world flow dynamics. However, the latest modification of the SPH method for shallow water equations by Vacondio et al. (2012a, 2012b, 2013) can deal with capturing shocks, open boundaries, and balancing discontinuous bed slopes. Later, the EFG method was proposed by Belytschko et al. (1995). This method came into existence as an extension of the diffuse element method proposed by Nayroles et al. (1992). In the EFG method, the generalized moving least square interpolation was used to define the local approximation of spatial derivatives. The FPM developed by Tiwari & Kuhnert (2003a) is extensively applied to various problems of fluid dynamics (Tiwari et al. 2007; Tiwari & Kuhnert 2002, 2007; Kuhnert 2003; Kuhnert & Ostermann 2014; Jefferies et al. 2015; Suchde et al. 2017). The FPM uses a weighted least square (WLS) method to approximate spatial derivatives at each point of the domain. The FPM is a Lagrangian model for fluid problems involving rapidly changing flow domains with respect to time (Tiwari & Kuhnert 2003a; Jefferies et al. 2015). In addition, it is relatively simple to implement and set the boundary conditions (Tiwari & Kuhnert 2003a).
Wang & Shen (1999) were the first to develop a meshless 1D shallow water model for the dam-break problem over a wet bed but without the source terms such as bed and friction slopes. As a result, the model application is limited to ideal shallow water flow problems. The presence of source terms creates a numerical imbalance between the driving forces (e.g., source terms) and convective fluxes and artificially accelerates the flow as the solution advances with time. The study by Rogers et al. (2003) illustrated that such imbalance is created when different terms are discretized using different methods. Therefore, a consistent discretization method is required to avoid artificial acceleration of the flow. Nevertheless, some notable meshless shallow water solvers such as GFDM by Li & Fan (2017), finite point method by Buachart et al. (2014) and radial basis function (RBF) by Chaabelasri (2018) are developed over time. However, the GFDM for shallow water equations results in wiggles near the discontinuities, especially for the 1D dam-break problems, as reported by Li & Fan (2017). The finite point method for shallow water equations by Buachart et al. (2014) has difficulties in handling complex source terms. The RBF-based meshless method for solving the shallow water equations uses artificial viscosity to capture shocks (Chaabelasri 2018). On the other hand, the FVM-based discretization for the Saint-Venant equations uses the Riemann solvers to accurately capture shocks and flow discontinuity without the need for artificial viscosity. Therefore, Riemann solvers may be incorporated within the framework of meshless methods to enhance their capability to capture physical features of shallow water flow dynamics such as discontinuities, shocks, abrupt variations in the bed topography, and flow variables (e.g., depth and velocity). The extensive literature review by the authors suggests that enough studies have not been carried out in the direction of using Riemann solvers in meshless methods. This is especially missing for simulating flows on highly variable topography using the meshless method. Hence, implementation of, for example, the Harten–Lax–van Leer (HLL) Riemann solver proposed by Harten et al. (1983) with a two-wave model resolving three constant states can enhance the capability of meshless methods to solve practical shallow water flow problems.
The numerical solution of the shallow water equations for domains like rivers and prismatic channels does not need to change the domain geometry rapidly. Due to the nature of its relatively stable geometry, a meshless method that can represent the geometry by a fixed number of irregularly distributed points at the beginning, but need not be maintained later, may be adopted for solving the shallow water or the Saint-Venant equations. In this regard, the FPM is proven to be a robust method for complex fluid flow problems (Uhlmann et al. 2013; Tiwari et al. 2007; Jefferies et al. 2015). The FPM already serves as a numerical basis for two commercially used meshfree simulation tools, NOGRID by Moller & Kuhnert (2007) and VPS-PAMCRASH by Tramecon & Kuhnert (2013). Due to the Lagrangian nature, the FPM is suitable for handling flow problems with complicated and even rapidly changing geometry with time. However, in FPM, the points move with physical properties such as mass, momentum, and velocity. Hence, in every few time steps, the points need to be managed through the neighbor searching algorithm. Also, a few points need to be removed if the points get clustered in one place. On the other hand, new points need to be introduced if holes are found in the computational domain. Therefore, FPM needs to be implemented in the Eulerian framework to represent the geometry by a fixed number of irregularly distributed points at the beginning that need not be maintained later during the simulation. Such modification is required, especially for simulating flows on rigid bed topography. However, there is a lack of research on modeling the Saint-Venant equations using the Riemann solver-based meshfree method.
To fill this gap, in this study, we have developed a shock-capturing meshless method for solving the 1D Saint-Venant equations for simulating flows in nonprismatic open channels with abrupt variation in topography. The HLL approximate Riemann solver is implemented within the framework of a meshless method for capturing shocks and discontinuity. The proposed numerical scheme is well balanced due to the adoption of the Saint-Venant equations in the form given by Ying & Wang (2008), which automatically incorporates the surface gradient method (SGM) of Zhou et al. (2001). The physical domain is represented by a set of irregularly distributed points. A local cloud of neighboring points, called satellites, is identified for each of these points. Local approximations for the spatial derivatives are performed by using the WLS method. In the WLS method, the weight function can be quite arbitrary, but in our computations, a Gaussian weight function, as described in the FPM method by Tiwari & Kuhnert (2003b), is considered. The water surface gradient in the source term is treated as the average water surface level at the midpoints of the center and satellite points of the local cloud by using the piece-wise linearly reconstructed variables. The driving forces and fluxes are also computed using the piece-wise linearly reconstructed variables. This ensures that the solution scheme is consistent and well balanced. The HLL Riemann solver and consistent discretization enable the proposed meshless method to effectively handle discontinuous and shock flows on a variable bed topography. The performance of the proposed meshfree method is verified by solving various analytical test cases, including water at rest condition over an irregular bed and variable width, dam-break flow on a wet and dry horizontal bed, steady flow over a hump with a hydraulic jump, dam-break flow over a hump, and tidal flow over a vertical step. The proposed meshless method is then validated by simulating laboratory experiments on hydraulic jump in a diverging channel and dam-break wave propagation over a triangular hump. The results show that the proposed method is highly accurate and does not diffuse at the wavefront. The developed meshless method with the HLL approximate Riemann solver eliminates the formation of wiggles at the wavefront and accurately handles the complex source terms, unlike the existing meshless methods (e.g., GFDM and FPM) (Buachart et al. 2014; Li & Fan 2017). Moreover, the proposed meshless method does not require any artificial viscosity to dampen numerical oscillations observed in earlier implementations (e.g., Chaabelasri 2018). Thus, the novelty of the study is to successfully implement the shock-capturing scheme within the meshless method and analyze its accuracy. It should be noted that the true advantages of a meshless method can be realized in 2D applications. In that context, the present study can be considered as a foundation for advancing research on the implementation of the proposed meshless shock-capturing method for solving the 2D shallow water equations for real-world problems and thereby realizing its full advantages.
GOVERNING EQUATIONS
is the water surface level from a datum as shown in Figure 1; in which z is the bottom topography, n denotes the Manning coefficient, 
is the hydraulic radius; in which P is the wetted perimeter of the channel. In Figure 1, B is the top width of the flow through the cross-section. 
Definition sketch of the cross-section.
Ying et al. (2004) argued that the above form of the Saint-Venant equations (Equation (1)) is common in engineering practice for finding out the required flow variables h and Q. In the other form of the equations, the surface gradient 
is split into hydrostatic pressure and bed slope source terms. Rogers et al. (2003) observed that the simulated flow artificially accelerates in the case of still-water conditions due to different discretization methods for the split terms. Zhou et al. (2001) proposed the SGM and pointed out that the main source of error is caused by inaccurate reconstruction of water depth. Therefore, in this study, the driving forces are taken into a single term (i.e., water surface gradient) in Equation (1), so it avoids the numerical imbalance and results in a well-balanced scheme, as discussed in the following section.
NUMERICAL METHOD
The meshless method based on FPM is proposed in this study to solve the 1D Saint-Venant equations. The physical domain is represented by a set of an irregularly distributed finite number of points, which possess fluid properties such as depth, velocity, and momentum. For each of these points, a local cloud of neighboring points is identified. The cloud points around the considered point are called satellites. Local approximations of the spatial derivatives are performed using the WLS method. The standard FPM (Tiwari & Kuhnert 2003b) is Lagrangian in nature, and points move with fluid properties along the flow domain. Due to its Lagrangian nature, the FPM is suitable for handling flow problems with complicated and even rapidly changing geometry with time. However, for every few time steps, the points need to be managed through a neighbor searching algorithm so that points can be removed if they are clustered in one place or filled with new points if holes appear during the simulation.
In this study, due to the nature of rigid bed geometry in open channel flows, we have adopted the FPM but in the Eulerian framework for solving the 1D Saint-Venant equations (Equation (1)). In the proposed method, the points distributed at the beginning need not be maintained later during the simulation. The convective flux term is computed using the HLL Riemann solver to simulate the discontinuous flow. Therefore, using the HLL Riemann solver within the meshless method is a major contribution of this study. The source terms are also considered for solving real-world problems. The numerical imbalance between bed slope and convective flux terms is avoided following the SGM of Zhou et al. (2001). The well-balanced and shock-capturing meshless method is presented below.
Least square approximation of derivatives
In many practical applications, the mesh-based numerical schemes, such as FDM, FEM, and FVM, play a very important role in determining the accuracy of a numerical solution. However, they may lose accuracy if the grid points are poorly distributed. The method presented herein does not require regular grids to approximate spatial derivatives of a function because the proposed meshless method uses the cloud of grid points instead of neighboring grid points alone. In this section, we present the basic theory of the WLS method to approximate spatial derivatives of a function.
be a domain containing N number of scattered points and 
be a scalar function and 
be its discrete values at the points 
for 
Figure 2 shows an example of the typical distribution of points in the domain. These points are numerical grid points 
. The neighboring cloud of points 
for the ith grid point are considered as shown in Figure 3. In the cloud, the considered point 
is the center, and the other points are called the satellites. Let 
be the total number of satellites in the cloud 
. 
A typical distribution of points in the domain 
.
A typical meshless cloud of points.
and the center 
can be expressed as
start from 
to 
, and its length is given by
based on the discrete values of its neighboring points, a WLS method (Tiwari & Kuhnert 2003b) is adopted in the present work like the earlier studies of meshfree methods (Tiwari & Kuhnert 2003a; Ma et al. 2008). The weight function can be quite arbitrary; however, we consider a Gaussian weight function of the form:
is a positive constant and is considered as 6.25 in our numerical computations as it gives the best approximation of derivatives. The derivatives of the function f at the point 
are determined using the Taylor series about the point 
and can be expressed in the following form:
.
in Equation (8) can be computed by minimizing the following norm over the neighboring cloud 
is given by 
. This yields
can also be estimated using the midpoint between i and j as
and 
is estimated at the midpoint between i and j (Chen & Shu 2005; Ma et al. 2008).Generating grid points
The grid points may cluster or scatter if they are generated randomly. To avoid this, we consider two parameters 
and 
such that 
The algorithm starts with initializing the boundary points. Then, the first interior point is filled having a minimal distance 
and a maximal distance 
from the current point, where 
is the average distance between the points. The same procedure is followed until the entire domain is filled with grid points. In our numerical computations, we set the values as 
and 
.
Spatial discretization
through the midpoint of 
is evaluated by the two states on the left 
and right 
of the midpoint
and 
are chosen for computing the flux through the midpoint using Equation (15) instead of just at the left (i.e., 
) and right (i.e., 
) of the midpoint (i.e., where the discontinuity is assumed to exist), the numerical scheme results in a first-order numerical scheme, which is often diffusive in nature. For higher order spatial accuracy, the variables are reconstructed at the midpoint. Since Z and B are also required to be defined at the midpoint, like the variables A and Q, they all are reconstructed in the same way. For this purpose, we define the vector, 
and the left and right states are reconstructed as
using the WLS method as follows:
and 
are left- and right-side points of the center of the cloud 
.The HLL approximate Riemann solver
, because of its robustness and ease of implementation. The approximate solution may be expressed by the three averaged states as
and 
define the left and right wave speeds of the two states and are given by
and 
are velocities of the left and right states, respectively, and the variables at the intermediate state are defined as
and 
are the velocity and flow depth of the intermediate state, and 
and 
are the averaged water depths (
) of the left and right states.
at the intermediate state is computed as described by Ying & Wang (2008), given by
and 
are two components of the flux 
, that is 
.
and 
are estimated according to the following expressions:
Source terms
is the average water surface level at the midpoint.It is important to note that Equation (27) relates the water surface gradient to the piece-wise linearly reconstructed water surface level. This makes the numerical scheme consistent and well balanced, and artificial acceleration of flow is eliminated (Ying & Wang 2008).
The friction source term is explicitly evaluated in each cloud 
using the averaged values of flow variables Q and A over the cloud 
Stability criteria
at the start of each time step:
is the Courant number (
). In this study, 
is used for all the test cases for obtaining convergent solutions.Boundary conditions
Boundary at inflow and outflow; new points introduced at both boundaries.
Temporal discretization
. The forward discretization of Equation (29) is written as
denote the physical quantities at the nth and (n + 1)th time steps. Therefore, the vector of conservative variables at (n + 1)th time step can be obtained by
A second-order Runge–Kutta time discretization was also implemented. However, the second-order time discretization did not show any noticeable difference in the results for the cases presented in the following section. On the contrary, the computation time was found to be significantly longer than the Euler time discretization. Therefore, the simple explicit time discretization described above is implemented in this study.
RESULTS AND DISCUSSION
The proposed meshless method for solving the 1D Saint-Venant equations is tested rigorously to examine its accuracy for various flow conditions. The problems are selected on subcritical, supercritical, transcritical flows, flows with shocks and discontinuity, variable width channels, and wave propagation on a dry bed, as discussed below.
Water at rest over an irregular bed and variable width
Bed elevation (in meters)
 | 0 | 50 | 100 | 150 | 250 | 300 | 350 | 400 | 425 | 435 | 450 | 475 | 500 | 505 |
 | 0 | 0 | 2.5 | 5 | 5 | 3 | 5 | 5 | 7.5 | 8 | 9 | 9 | 9.1 | 9 |
 | 530 | 550 | 565 | 575 | 600 | 650 | 700 | 750 | 800 | 820 | 900 | 950 | 1,000 | 1,500 |
 | 9 | 6 | 5.5 | 5.5 | 5 | 4 | 3 | 3 | 2.3 | 2 | 1.2 | 0.4 | 0 | 0 |
| | 0 | 50 | 100 | 150 | 250 | 300 | 350 | 400 | 425 | 435 | 450 | 475 | 500 | 505 |
| | 0 | 0 | 2.5 | 5 | 5 | 3 | 5 | 5 | 7.5 | 8 | 9 | 9 | 9.1 | 9 |
| | 530 | 550 | 565 | 575 | 600 | 650 | 700 | 750 | 800 | 820 | 900 | 950 | 1,000 | 1,500 |
| | 9 | 6 | 5.5 | 5.5 | 5 | 4 | 3 | 3 | 2.3 | 2 | 1.2 | 0.4 | 0 | 0 |
Width of the channel (in meters)
 | 0 | 50 | 100 | 150 | 200 | 250 | 300 | 350 | 400 | 425 | 435 | 450 | 470 | 475 | 500 |
 | 40 | 40 | 30 | 30 | 30 | 30 | 30 | 25 | 25 | 30 | 35 | 35 | 40 | 40 | 40 |
 | 505 | 530 | 550 | 565 | 575 | 600 | 650 | 700 | 750 | 800 | 820 | 900 | 950 | 1,000 | 1,500 |
 | 45 | 45 | 50 | 45 | 40 | 40 | 30 | 40 | 40 | 5 | 40 | 35 | 25 | 40 | 40 |
| | 0 | 50 | 100 | 150 | 200 | 250 | 300 | 350 | 400 | 425 | 435 | 450 | 470 | 475 | 500 |
| | 40 | 40 | 30 | 30 | 30 | 30 | 30 | 25 | 25 | 30 | 35 | 35 | 40 | 40 | 40 |
| | 505 | 530 | 550 | 565 | 575 | 600 | 650 | 700 | 750 | 800 | 820 | 900 | 950 | 1,000 | 1,500 |
| | 45 | 45 | 50 | 45 | 40 | 40 | 30 | 40 | 40 | 5 | 40 | 35 | 25 | 40 | 40 |
Water at rest over an irregular bed and variable breadth: (a) bed topography and (b) channel breadth.
Water at rest over an irregular bed and variable breadth: (a) bed topography and (b) channel breadth.
Water at rest over an irregular bed and variable breadth: (a) water surface level and (b) discharge at time 
.
Water at rest over an irregular bed and variable breadth: (a) water surface level and (b) discharge at time .
error is defined as given below:
The errors for different number of grid points are summarized in Table 3.
The 
errors in the computed water surface elevation and velocity on three different points
| N | Maximum error for Z and u | |
|---|---|---|
 |  | |
| 100 | 0 | 3.99 × 10−16 |
| 200 | 0 | 2.01 × 10−16 |
| 500 | 0 | 6.62 × 10−18 |
| N | Maximum error for Z and u | |
|---|---|---|
| | |
| 100 | 0 | 3.99 × 10−16 |
| 200 | 0 | 2.01 × 10−16 |
| 500 | 0 | 6.62 × 10−18 |
Dam break on wet and dry beds
from the upstream end. The initial upstream water depth is 
, and downstream water depth is 
for wet and 
for dry beds. A total of 600 meshless points are irregularly distributed along the channel length. At time 
the dam is removed instantaneously. In the numerical computation, a dry bed tolerance of 
is used to distinguish between wet and dry bed, as discussed in Kuiry et al. (2012). The results are reported at 30 s after the removal of the dam. The computed water surface profiles and discharges compared with the analytical solutions are available in Toro (2001). Figures 7 and 8 show the wet bed and dry bed cases, respectively. It should be noted that the proposed meshless method can accurately capture the sharp changes in depth and velocity profiles, unlike the reported mesh-based finite volume solutions of Kuiry et al. (2012). Therefore, the proposed meshless model can accurately capture discontinuous flow with shocks over wet and dry beds. 
Dam break on the wet bed: (a) water depth and (b) discharge.
Dam break on the dry bed: (a) water surface elevation and (b) discharge.
error norm of h and u for wet bed case and root mean square error of h and Q for wet bed are investigated. The error functions are defined as
and 
are numerical and analytical values of the ith point; and 
and 
are the numerical and analytical values of the ith point. Tables 4 and 5 summarize the errors for different values of N. Norm of non-dimensional error 
and 
for the dam-break on the wet bed
| Points | Without limiter | With limiter | ||
|---|---|---|---|---|
| N |  |  |  |  |
| 80 | 3.89 × 10−2 | 3.97 × 10−2 | 2.26 × 10−2 | 2.45 × 10−2 |
| 150 | 3.55 × 10−2 | 3.26 × 10−2 | 1.59 × 10−2 | 1.52 × 10−2 |
| 300 | 2.14 × 10−2 | 2.11 × 10−2 | 9.41 × 10−3 | 9.43 × 10−3 |
| 600 | 1.41 × 10−2 | 1.38 × 10−2 | 6.05 × 10−3 | 5.82 × 10−3 |
| Points | Without limiter | With limiter | ||
|---|---|---|---|---|
| N | | | | |
| 80 | 3.89 × 10−2 | 3.97 × 10−2 | 2.26 × 10−2 | 2.45 × 10−2 |
| 150 | 3.55 × 10−2 | 3.26 × 10−2 | 1.59 × 10−2 | 1.52 × 10−2 |
| 300 | 2.14 × 10−2 | 2.11 × 10−2 | 9.41 × 10−3 | 9.43 × 10−3 |
| 600 | 1.41 × 10−2 | 1.38 × 10−2 | 6.05 × 10−3 | 5.82 × 10−3 |
Root mean square errors 
and 
for dam-break on the dry bed
| Points | Without limiter | With limiter | ||
|---|---|---|---|---|
| N |  |  |  |  |
| 80 | 2.409 × 10−1 | 2.467 × 10 | 1.181 × 10−1 | 1.217 × 10 |
| 150 | 1.804 × 10−1 | 1.561 × 10 | 7.636 × 10−2 | 5.199 × 10−1 |
| 300 | 1.092 × 10−1 | 1.113 × 10 | 3.113 × 10−2 | 3.017 × 10−1 |
| 600 | 6.842 × 10−2 | 7.171 × 10−1 | 1.451 × 10−2 | 1.279 × 10−1 |
| Points | Without limiter | With limiter | ||
|---|---|---|---|---|
| N | | | | |
| 80 | 2.409 × 10−1 | 2.467 × 10 | 1.181 × 10−1 | 1.217 × 10 |
| 150 | 1.804 × 10−1 | 1.561 × 10 | 7.636 × 10−2 | 5.199 × 10−1 |
| 300 | 1.092 × 10−1 | 1.113 × 10 | 3.113 × 10−2 | 3.017 × 10−1 |
| 600 | 6.842 × 10−2 | 7.171 × 10−1 | 1.451 × 10−2 | 1.279 × 10−1 |
Dam break flow on a variable bottom topography
and 
. The variable bed topography is defined as
, which separates the channel into two initial water columns. The initial conditions are defined as
. At the time 
, the dam is removed instantaneously, and results are noted at 
after the removal of the dam. The computed water surface profile and discharge along the channel are compared with the analytical solutions by Hubbard's approach in Figure 9. The analytical solution is discussed in Hudson (1999). It is very interesting to observe that the proposed meshless method can accurately capture the transition in flow depth and velocity. 
Dam break on a variable depth riverbed: (a) water surface elevation and (b) velocity.
Dam break on a variable depth riverbed: (a) water surface elevation and (b) velocity.
Steady flow over a hump
Depending on the initial and boundary conditions, the flow may be subcritical, transcritical with shock or without shock, or supercritical. For all the cases presented below, the domain is represented by 233 irregularly distributed points. The time step, 
, is used for a stable solution. Three different flow conditions are simulated and presented below.
Subcritical flow
and the constant water surface level of 
are set as upstream and downstream boundary conditions, respectively. The simulation is continued for a sufficiently long time to achieve the steady-state condition. The results are plotted in Figure 10 and show a good agreement between numerical and analytical solutions. The discharge at the point of bed transition downstream is deviated by 
However, such a deviation is insignificant compared with previous studies (e.g., Valiani & Begnudelli 2006). 
Subcritical flow over a hump: (a) water surface elevation and (b) discharge.
Transcritical flow without a shock
is imposed, and at the downstream, no boundary condition is needed since the flow is supercritical. Figure 11 compares the numerical and analytical solutions at the steady state. In this case, there is no deviation in discharge that can be observed. 
Transcritical flow without a shock: (a) water surface elevation and (b) discharge.
Transcritical flow without a shock: (a) water surface elevation and (b) discharge.
Transcritical flow with a shock
and a constant water surface elevation of 
are set as upstream and downstream boundary conditions, respectively. Figure 12 shows the comparison of numerical and analytical solutions of water surface elevation and unit discharge along the channel at the steady state. The discharge at the point of bed transition at downstream is deviated by 
Similar differences were reported as 
by Zhou et al. (2002), 
by Valiani & Begnudelli (2006), and 
by Brufau et al. (2002). It is worth stressing that an error on the discharge does not imply errors in mass conservation (Valiani & Begnudelli 2006), and hence the behavior of the system may still be considered very good. 
Transcritical flow with a shock: (a) water surface elevation and (b) discharge.
It can be concluded from the above set of tests that the proposed meshless method can successfully simulate different flow regimes in the presence of bed variation and shocks. It should also be noted that the results obtained from the meshless method are superior compared with many previous studies based on the FVM and approximate Riemann solvers.
Tidal wave flow over a step
with the analytical solutions in Figure 13. It can be observed from the comparison that the proposed meshless method can accurately capture the flow discontinuity due to the presence of vertical steps on the bed. 
Tidal wave over a step: (a) water surface elevation and (b) velocity.
Hydraulic jump in a divergent channel
with rectangular cross-sections. However, the channel width varies, as shown in Figure 14, following the below expression:
Hydraulic jump in a rectangular channel: diverging channel width variation.
and discharge of 
at the inlet, and a water depth of 
at the outlet are prescribed as boundary conditions. At steady state, a hydraulic jump is formed as shown in Figure 15. The proposed numerical scheme can accurately capture the location and profile of the hydraulic jump as compared with the experimental observation in Figure 15. The present numerical solution is also compared with the 2D model results of Younus & Chaudhry (1994) and no significant differences are observed. 
Hydraulic jump in a diverging channel.
Flow over a triangular obstacle
contained by a dam located at 
and a dry bed downstream within a rectangular channel. A symmetric triangular obstacle (
long and 
high) was placed on the channel, with its peak located 
downstream of the dam. To observe depth evolution, gauges were located at 
(G4), 
(G10), 
(G11), 
(G13), and 
(G20) from the dam as shown in Figure 16. The gauge point G13 was located at the vertex of the obstacle, and therefore it is a critical point. The fixed boundaries are walls except for the open outlet. The Manning's roughness coefficient was supplied as 
by the experimentalists (Brufau et al. 2002). 
Geometry and gauge locations in the experimental and model setup.
, and a dry bed tolerance, 
, are used. The simulation is carried out for 
. The depth variation with respect to time at the five gauges is compared with the experimental observations in Figure 17. It can be seen that the simulated flow depths closely follow the observed data with respect to depth variation and arrival time of the waterfront to the respective gauges. The transition from wet to dry and vice versa at the gauge point G13 is also well predicted. The flow depth comparison at the tip of the hump shows that the proposed meshless method has better agreement with the experimental measurement than the FVM-based solution reported by Kuiry et al. (2012). However, similar to earlier studies by Chaabelasri et al. (2019), a large deviation with respect to experimental data can be seen at G20, but the flow depth variation at this location during the entire simulation period is very small. 
Computed and measured water depth variations at different stations: (a) G4, (b) G10, (c) G11, (d) G13, and (e) G20.
Computed and measured water depth variations at different stations: (a) G4, (b) G10, (c) G11, (d) G13, and (e) G20.
CONCLUSIONS
A meshless shock-capturing method is presented for solving the 1D Saint-Venant equations describing open channel flow. The HLL approximate Riemann solver is implemented within the framework of the WLS method for capturing shocks and flow discontinuity. The first-order Euler time discretization results in similar level of accuracy to that of the second-order Runge–Kutta time discretization method but in significantly less computation time. A stability condition for the proposed explicit scheme is shown to produce convergent solutions. The paper describes the theoretical background, verification, and validation of the 1D shock-capturing meshless method for solving the Saint-Venant equations. The robustness of the proposed meshless method is verified by comparing the simulated results with several theoretical benchmark tests and laboratory experiments. The presented method is stable and does not artificially accelerate while solving still-water conditions in a closed domain. The predictions of the water surface elevation and discharge for an ideal dam break on the wet and dry bed are in very good agreement with the analytical solutions. The flow discontinuity over vertical steps is accurately captured by the developed method. The hydraulic jump formed in a laboratory scale divergent channel is reproduced with a similar accuracy level to a previous study. The water wave propagation due to a dam break on a dry bed over a triangular obstacle is well simulated, and the computed results closely follow the experimental observations. Also, the shifting of dry and wet conditions is well captured. The meshless method presented in this study can be extended to solve the 2D shallow water equations to capture sharp changes in topography and flow dynamics by distributing a greater number of localized points.
ACKNOWLEDGEMENTS
This research is financially supported by the Indian Institute of Technology Madras under the MHRD project (Project No. SB20210848MAMHRD008558).
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.
CONFLICT OF INTEREST
The authors declare there is no conflict.
%20cropped.png?versionId=7142)