Discretization of Multidimensional Mathematical Equations of Dam Break Phenomena Using a Novel Approach of Finite Volume Method

This paper was concerned to simulate both wet and dry bed dam break problems. A high-resolution finite volume method (FVM) was employed to solve the one-dimensional (1D) and two-dimensional (2D) shallowwater equations (SWEs) using an unstructured Voronoi mesh grid. In this attempt, the robust local Lax-Friedrichs (LLxF) scheme was used for the calculating of the numerical flux at cells interfaces. The model named V-Break was run under the asymmetry partial and circular dam break conditions and then verified by comparing the model outputs with the documented results. Due to a precise agreement between those output and documented results, the V-Break could be considered as a reliablemethod for dealing with shallowwater (SW) and shock problems, especially those having discontinuities. In addition, statistical observations indicated a good conformity between the V-Break and analytical results clearly.


Introduction
Floods induced by dam failures can cause significant loss of human life and property damages, especially when located in highly populated regions.These entail numerical and laboratory investigations of dam break flows and their potential damage.The shallow water equations (SWEs) are conventionally used to describe the unsteady open channel flow such as dam break.These equations are named as Saint Venant equations for one-dimensional (1D) problem and also include the continuity and momentum equations for twodimensional (2D) studies.
Recently, Shamsai and Mousavi evaluated the dam break parameters such as breach width, side slop, time of failure and peak outflow, using 142 case studies of previous researches.Then, using those evaluations and also the BREACH and FLDWAV application software, the studies of Aidoghmosh earth dam breach were examined [17].Xia et al. developed a 2D morphodynamic model for predicting dam break flows over a mobile bed [18].Shamsi et al. developed a 2D flow model based on SWEs.They used finite volume method (FVM), total variation diminishing (TVD), and weighted average flux (WAF) schemes as well as Harten-Lax-van Leer-Contact (HLLC) Riemann solver [19].Erpicum et al. presented a 2D finite volume (FV) multiblock flow solver, which was able to deal with the natural topography variation [20].Baghlani utilized a combination of the robust and effective flux-difference splitting (FDS) and flux-vector splitting (FVS) methods to simulate dam break problems based on FVM on a Cartesian grid.The method combined the effectiveness and robustness of the FDS and FVS methods to precisely estimate the numerical flux at each cell interface [21].Zhang and Wu developed a hydrodynamic and sediment transport model for dam break flows.The 2D SWEs were solved based on the FVM with an unstructured quadtree mesh grid [22].Singh et al. developed a 2D numerical model to solve the SWEs for the simulation of dam break problems [23].Chang et al. proposed a meshless numerical model to investigate the shallow water (SW) dam break in 1D open channel.A numerical model was used to solve the SWEs based on smoothed particle hydrodynamics (SPH).The concept of slice water particles was adapted in the SPH-SWE formulation [24].Shakibaeinia and Jin developed a new mesh-free particle model based on the weakly compressible MPS (WC-MPS) formulation for modeling the dam break problem over a mobile bed [25].Sarveram and Shamsai investigated the dam break problem in converge and diverge rectangular channels in the unsteady stance using Saint Venant equations and a quasi-Lagrangian method [26].
This paper attempts to present a novel development for 1D and 2D dam break problems in both wet and dry beds.A high-resolution FVM is employed to solve the SWEs on unstructured Voronoi mesh.The local Lax-Friedrichs (LLxF) scheme is used for the estimation of fluxes at cells and the numerical approximation of hyperbolic conservation laws.

Research Methodology
2.1.Governing Equations.The continuity and momentum equations of the SW can be written in different forms depending upon the requirements of the numerical solution of governing equations.The 2D SWEs with source terms are given in the vector form considering a rigid bed channel [1] as follows: In this study, the friction slopes were estimated using the Manning's formulas.In the case of the dam break flow, the influence of bottom roughness prevailed over the turbulent shear stress between cells.Therefore the effective stress terms were neglected in the computation.

Numerical Modeling Algorithm.
The main advantage of the FVM is that volume integrals in a partial differential Equation (PDE) containing a divergence term are converted to surface integrals using the divergence theorem.These terms are then evaluated as fluxes at the surfaces of each FV.Because the flux entering a given volume is identical to that leaving the adjacent volume, these methods are conservative.Another advantage of the FVM is that it is easily formulated to allow for unstructured meshes.Unstructured grid methods utilize an arbitrary collection of elements to fill the domain.These types of grids typically utilize triangles in 2D and tetrahedrals in 3D, although quadrilateral, hexahedral, Voronoi, and Delaunay meshes can also be unstructured.In the Voronoi mesh, the chosen point has lower distance in the devoted domain rather than other points.If one point has the same distance from several domains, it will be divided between domains.Indeed, these points create Voronoi cell boundaries.Consequently, internal sections of the Voronoi mesh consist of nodes belonging to one domain and boundaries include nodes that belong to several domains [27].
In this paper, the studied domain was discretized using unstructured Voronoi meshes.Delaunay triangulation was created and then the Voronoi mesh was established using the Qhull program in MATLAB software.The governing equation was discretized applying the FVM.In this approach, the studied domain was divided into several separated control volumes without any overlapping.By integrating the governing differential equation over every control volume, the system of algebraic equations was created so that each of its formulations belonged to one control volume and each equation linked a parameter in the control volume node to different numbers of the parameter in adjacent nodes.This consequently led to the computation of the parameter in each node [28].
In order to solve discrete equations, the parameter in each node was computed considering its discrete equation and newest adjacent nodes' values.Solution procedure can be expressed as follows.
(1) Assuming an initial value in each node as an initial condition.
(2) Calculating the value in a node considering its discrete equation.(3) Performing pervious step for all nodes over the studied domain, one cycle is performed by repetition this step.(4) Verifying the convergence clause.If this clause is satisfied, the computing will end otherwise the computations will be repeated from the second step.
As exact values of boundary conditions were not distinct, the Riemann boundary condition was utilized for computing the investigating parameter.Therefore, by assuming a layer which is close to the boundary layer, / = 0 and / = 0 were defined for the investigating parameter and then calculated values for boundary adjacent nodes transform to related boundary nodes.This procedure will continue until the results difference is converged to zero.

Discretization of Governing Equations
2.3.1.FV Discretization.Various methods can be used to discretize the governing equations, among which the FVM due to its ability to satisfy mass and momentum conservation is frequently adopted.In this research, the discretization of (1) was performed using the FVM with unstructured Voronoi mesh, as shown in Figure 1.It should be mentioned that  face vertexes nomination direction is counterclockwise from  to  centering p (see Figure 1) [1] as follows: By implementing divergence theorem, ( 6) is yielded to (7) as follows: where ⃗  = ⃗  + ⃗ .Equation ( 7) can be written as ( 8) by approximating the line integral for all control volumes and nodes, generally, 2.3.2.Voronoi Mesh.By applying the Voronoi mesh, ( 8) is yielded to (9) for an investigated control volume; however, Figure 1 illustrates the 2D Voronoi mesh grid used for describing these equations as follows: Then, The discrete equation can be written as followS: 2.4.The Local Lax-Friedrichs (LLxF) High-Order Scheme.In shock capturing schemes, the location of discontinuity is captured automatically by the scheme as a part of the solution procedure.These slope-limiter or flux-limiter methods can be extended to systems of equations.In this paper, the algorithm is based upon central differences with comparable performance to Riemann type solvers when used to obtain a solution for PDE's describing systems.FV and Finite Difference (FD) methods are closely related to central schemes like the most shock capturing schemes [29].Rusanov scheme is often called the LLxF method, because it has the same form as the Lax-Friedrichs (LxF) method but the viscosity coefficient is chosen locally.It can be shown that this is a sufficient viscosity to make the method converge to the vanishingviscosity solution.It means that it is less diffusive than normal LxF, since it locally limits the numerical viscosity instead of having a uniform viscosity on the entire domain.LLxF formulates the FVM in space and the explicit Euler method in time.Many researchers (e.g., Lin et al. [30], van Dam and Zegeling [31], and Lu et al. [32]) utilized LLxF splitting scheme in different problems such as 2D SWEs, 1D adaptive moving mesh method and its application to hyperbolic conservation laws from magnetohydrodynamics (MHD), and the performance of the weighted essential nonoscillatory (WENO) method.
In this research, LLxF is used as a flux calculator.By expanding ( 12), ( 13) can be written as follows: Intercell fluxes can be estimated by implementing the following equations: After computing intercell fluxes by utilizing the LLxF scheme in Voronoi mesh, equations can be solved and the final result can be calculated for each time step.The Δ can be computed using the Courant Friedrichs Lewy (CFL) for each time step as follows: × The CFL should range over [0, 1] for achieving to the stability (0 < CFL < 1).

Preparation and Validation of the Numerical Algorithm.
The CFD code named V-Break was prepared on the unstructured Voronoi mesh grid using MATLAB programming.Figure 2 illustrates the running process of the V-Break as a flowchart.

1D Dam Break Test.
At the instant of the dam break, water is released through the breach, forming a positive wave that propagates downstream and a negative wave that moves upstream.Here, Stoker's analytical solution of dam break problem can be used to illustrate the accuracy of the numerical schemes.Stoker derived this theory just for the wet bed dam break problem, but it is possible to develop it for dry bed considering a downstream water depth very close to zero [33].A 1D horizontal rectangular channel having 1 m in length with wall at either ends and no roughness was considered, however, the source term equaled to zero.The initial velocity is zero and a barrier is present at  = 0.5 m, which is removed at  = 0 s.The water depths of the upstream and downstream are 1 and 0.5 m, respectively.Also, the CFL parameter is equal to 0.9. Figure 3 illustrates the 1D dam break domain [34].Figures 4 and 5 show the water depth and velocity of the flow obtained using V-Break and Stoker's analytical solution at  = 0.02 and 0.1 s.Furthermore, Table 2 presents the results of the Mann-Whitney test for a statistical comparison between the obtained results using V-Break and Stoker's analytical solution.
In terms of comparing two groups of data statistically, the null hypothesis ( 0 ) is usually a hypothesis of "nondifference." It means that there is no difference between the ranks of the two comparing groups.These two groups can be defined as V-Break and Stoker's analytical results.When computing the value of Mann-Whitney  test, the number of comparisons equals to the product of the number of data in the first group (  ) times the number of data in the second group (  ), which is equal to   ×   .If the null hypothesis is true, then the value of Mann-Whitney  test should be about half thit value.The smallest possible value of Mann-Whitney  test is zero.The largest possible value is equal to (  ×   )/2.If this value is much smaller than that, the -value will be small.The  value or calculated probability is the estimated probability of rejecting the null hypothesis of a study question when that hypothesis is true, while the significance level is usually considered as  = 0.05.So, if the  value is more than 0.05, the null hypothesis will be true, especially when  value is closer to 1. Also, the null hypothesis is true when RS1 ≤ Wilcoxon W ≤ RS2.Here, RS1 is the sum of ranks for the first group and RS2 is the sum of ranks for  the second group.Similarly, the null hypothesis will be true if || ≤  0 , while  0 is the critical value extracted from the -distribution graph statistically.This value for the present study equal to 1.96.
In this research, all the above conditions were satisfied in Table 2. So, the null hypothesis is true and there is no difference between V-Break and Stoker's analytical solution results.Consequently, V-Break results are acceptable having a good agreement with Stoker's analytical solutions.

V-Break.
For this case study, a 2D partial dam break problem with asymmetrical breach was considered.The computational domain was defined in a channel with 200 m in length and 200 m in width.The breach is 75 m in length and the dam is 15 m in height.The initial upstream water depth is 10 m.The downstream water depth is 5 m in a wet bed and 0.1 m in a dry bed.The roughness coefficient was assumed zero implying a frictionless surface.In this example, a domain with 40 × 40 node points, it was proposed and a time of 7.2 s was considered as the total time for the calculation procedure.Also, the CFL parameter is equal to 0.9. Figure 6 illustrates the 2D dam break domain.Figure 7 shows the results of previous studies.
Different conditions of the domain were simulated using V-Break.Figure 8 shows the output results of V-Break.

Mesh Configuration Comparison. The output results
were then compared with the previous studies.A section was considered at  = 130 m and the leading points of water depth contours were exported as a function of distance ().Figure 9 shows the comparison among results obtained using V-Break on Voronoi, results obtained by Liang et al. [11] on rectangular, and results obtained by Loukili and Soulaïmani [14] on triangular mesh grids in both wet and dry beds.

Riemann Solver Comparison.
As previously mentioned, Wang and Liu implemented four typical FVMs, including the Reo-MUSCL, Reo-upwind, HLL-MUSCL, and CFLF8 composite methods on unstructured triangular meshes to simulate a 2D dam break problem [7].In this research, the Riemann solver used in the current study was compared with solvers used by Wang and Liu considering the particular mesh grids used at each study [7].For a better visual comparison, a section was considered at  = 130 m and the leading points of water depth contours were exported, and then results are shown in Figure 10.

2D Circular Dam Break Test.
In this test, 200 computational cells were used where initial conditions consist of two states separated by a circular discontinuity.The radius of the circle is 50 m and it is centered at  = 100 m.Both components of the velocity, u and v, were set to zero everywhere and ℎ was set to 10 m within the circle and 1 m outside the circle.The computational grid consisted of 40 × 40 cells and the solution was sought after 2 s.The domain was defined for the V-Break, and then the output results were compared with the previous studies.Figures 11  and 12 illustrate output results of V-Break including the water surface formation and water depth contours, respectively.Also, Figure 13 shows results performed by Baghlani [21] in order to use them in comparison.Figures 11,12,and 13 show that the V-Break can simulate the 2D circular dam break well, as compared with previous study by Baghlani [21].

Conclusions
In the current research, a novel and friendly user code named V-Break was written showing that the LLxF scheme along with the FVM on the unstructured Voronoi grid is a suitable combination in order to simulate 1D/2D dam break problems.The advantages of this method are very promising, especially in reconstructing the conducted tests.For 1D dam break, Stoker's analytical solution was considered for validation.Illustrations and computed  values demonstrated that there are no significant differences between Stoker's analytical solution results and V-Break outputs.In addition to 1D partial dam break, the 2D partial dam break tests were done.Results were compared with some other previous studies for validation at  = 7.2 s.Each test can be done in the wet or dry bed conditions.In all scenarios considered here, no significant numerical dispersion problem or nonphysical alternation was observed in the results.The comparison showed a good agreement between V-Break results with previous studies.In terms of mesh grid comparison, it was seen that the Voronoi mesh grid results are closer to triangular mesh grid results compared with rectangular mesh grid results.In addition, the 2D circular dam break test was done and its results were compared with previous studies at  = 2s which seems accepted.Generally, the unstructured Voronoi mesh grid is able to model inlet and outlet fluxes in every direction of control volume faces.The results indicated a higher efficiency and precision of the discrete equations resulted from the Voronoi mesh.Thus, it could be recommended to utilize the Voronoi mesh in the numerical discrete equations.The Voronoi mesh grid is able to model complicated geometries, and also it could produce the final discrete equations leading

𝑛:
Denotes parameters belonging to time of   + 1: Denotes parameters belonging to time of  + Δ.

Figure 1 :
Figure 1: The 2D schematic Voronoi mesh cell used for describing the discretization of the governing equations.
: Verticalcoordinatecomponent ∑  : The sum over the all Voronoi cells : Central node of investigated Voronoi cell Δ: Time interval Δ: The distance between the central node of investigated Voronoi cell and the adjacent cells  0 : The null hypothesis : Thevalueof-Test  0 : The critical value of  extracted from the statistical -distribution graph RS1: Sum of ranks for the first comparing group RS2: Sum of ranks for the second comparing group : The significance level   : The number of data in the first comparing group   : The number of data in the second comparing group.Subscripts : Denotes the parameter at the Voronoi cells side's area  : Counts all central control volumes : Counts all nodes of the central control volumes nb: Denotes parameters at the central node of adjacent Voronoi cells : Denotes parameters at central node of investigated Voronoi cell in: Denotes the parameter outside the Voronoi cells side's area  out: Denotes the parameter outside the Voronoi cells side's area  : The parameter component in   direction : The parameter component in   direction.

Table 1 :
General summarization of the equipment used by some other researchers and present study.

Table 2 :
Mann-Whitney test parameters for comparison between V-Break and stoker's analytical solution results.