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

Abstract

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) shallow water 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 reliable method for dealing with shallow water (SW) and shock problems, especially those having discontinuities. In addition, statistical observations indicated a good conformity between the V-Break and analytical results clearly.

1. 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 two-dimensional (2D) studies.

Many researchers studied the dam break problem, such as Toro [1], Wu et al. [2], Wang and Shen [3], Mohapatra et al. [4], Zoppou and Roberts [5], Wang et al. [6], Wang and Liu [7], Venutelli [8], Ponce et al. [9], Zhou et al. [10], Liang et al. [11], Quecedo et al. [12], Begnudelli and Sanders [13], Loukili and Soulaïmani [14], Díaz et al. [15] and Aliparast [16], especially using computational fluid dynamic (CFD) methods.

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.

2. 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: 𝜕𝑈+𝜕𝑡𝜕𝐹+𝜕𝑥𝜕𝐺ℎ𝑢𝜕𝑦+𝑆=0,𝑈=𝑢ℎ𝑣ℎ,𝐹=𝑢ℎ₂1ℎ+2𝑔ℎ²𝑣𝑢𝑣ℎ,𝐺=𝑣ℎ𝑢ℎ₂1ℎ+2𝑔ℎ²0𝑆,𝑆=−𝑔ℎ_0𝑥−𝑛²ℎ^4/3𝑢√𝑢²+𝑣²−𝑔ℎ𝑆_0𝑦−𝑛²ℎ^4/3𝑢√𝑢²+𝑣². (1)(2)(3)(4)(5)

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.

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

2.3. 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 𝑝 (see Figure 1) [1] as follows: _𝐴𝜕𝑈𝜕𝑡𝑑𝐴+_𝐴∇⋅𝐻𝑑𝐴=_𝐴𝑆𝑑𝐴. (6) By implementing divergence theorem, (6) is yielded to (7) as follows: _𝐴𝜕𝑈𝜕𝑡𝑑𝐴+_𝑙𝐻⋅𝑛𝑑𝑙=_𝐴𝑆𝑑𝐴, (7) where 𝐻=𝐹𝑖+𝐺𝑗. Equation (7) can be written as (8) by approximating the line integral for all control volumes and nodes, generally, 𝑑𝑈_𝑖1𝑑𝑡=−𝐴_𝑖𝑗𝐻⋅𝑛_𝑖𝑗Δ𝑙_𝑖𝑗+𝑆_𝑖. (8)

Figure 1 

The 2D schematic Voronoi mesh cell used for describing the discretization of the governing equations.

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: 𝜕𝑈1𝜕𝑡=−𝐴_𝑓𝐻⋅𝐴𝑒_𝜀𝑓+𝑆,𝐻=𝐻_𝜀⋅𝑒_𝜀+𝐻_𝜂⋅𝑒_𝜂. (9)(10) Then, 𝜕𝑈1𝜕𝑡=−𝐴_𝑓𝐻_𝜀𝐴_𝑓+𝑆,_{𝑡𝑡+Δ𝑡}𝜕𝑈𝜕𝑡𝑑𝑡=_𝑡^𝑡+Δ𝑡−1𝐴_𝑓𝐻_𝜀𝐴_𝑓𝑑𝑡+_𝑡^𝑡+Δ𝑡𝑈𝑆𝑑𝑡,_𝑝𝑛+1−𝑈^𝑛_𝑝1=−𝐴_𝑓𝐻_𝜀𝐴_𝑓𝐻Δ𝑡+𝑆Δ𝑡,_𝜀=𝑛₁𝐹_𝜀+𝑛₂𝐺_𝜀,𝑛₁=𝑥_nb−𝑥_𝑝𝑥_nb−𝑥_𝑝²+𝑦_nb−𝑦_𝑝²,𝑛₂=𝑦_nb−𝑦_𝑝𝑥_nb−𝑥_𝑝²+𝑦_nb−𝑦_𝑝². (11) The discrete equation can be written as followS: 𝑈_𝑝^𝑛+1=𝑈^𝑛_𝑝−Δ𝑡𝐴_𝑓𝑛₁𝐹_𝜀𝐴_𝑓+𝑛₂𝐺_𝜀𝐴_𝑓^𝑛+𝑆^𝑛Δ𝑡. (12)

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 vanishing-viscosity 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: 𝑈_𝑝^𝑛+1=𝑈^𝑛_𝑝−Δ𝑡𝐴_𝑓𝐹_𝑓nb⋅out𝐴_𝑓𝑛_1𝑓+_𝑓𝐹_𝑓nb⋅in𝐴_𝑓𝑛_1𝑓+_𝑓𝐺_𝑓nb⋅out⋅𝐴_𝑓⋅𝑛_2𝑓−_𝑓𝐺_𝑓nb⋅in𝐴_𝑓𝑛_2𝑓+𝑆Δ𝑡. (13) Intercell fluxes can be estimated by implementing the following equations: 𝐹_𝑓nb⋅out=𝐹𝑢^𝑛_nb+𝐹𝑢^𝑛_𝑝2−12𝑢+𝑔ℎ𝑢^𝑛_nb+𝑢^𝑛_𝑝2𝑢^𝑛_nb−𝑢^𝑛_𝑝,𝐹_𝑓nb⋅in=𝐹𝑢^𝑛_nb+𝐹𝑢^𝑛_𝑝2−12𝑢+𝑔ℎ𝑢^𝑛_nb+𝑢^𝑛_𝑝2𝑢^𝑛_𝑝−𝑢^𝑛_nb,𝐺_𝑓nb⋅out=𝐺𝑣^𝑛_nb+𝐺𝑣^𝑛_𝑝2−12𝑣+𝑔ℎ𝑣^𝑛_nb+𝑣^𝑛_𝑝2𝑣^𝑛_nb−𝑣^𝑛_𝑝,𝐺_𝑓nb⋅in=𝐺𝑣^𝑛_nb+𝐺𝑣^𝑛_𝑝2−12𝑣+𝑔ℎ𝑣^𝑛_nb+𝑣^𝑛_𝑝2𝑣^𝑛_𝑝−𝑣^𝑛_nb. (14) 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: Δ𝑡=CFL×min𝑥_nb−𝑥_𝑝²+𝑦_nb−𝑦_𝑝²×1max𝑢+𝑔ℎ+1max𝑣+𝑔ℎ. (15) The CFL should range over [0,1] for achieving to the stability (0<CFL<1).

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

Figure 2 

The V-Break flowchart.

V-Break was then validated using Stoker’s analytical solution in 1D [33], and previous results were obtained by other researchers in 2D (e.g., Wang and Liu [7], Loukili and Soulaïmani [14], Liang et al. [11], and Baghlani [21]). Table 1 presents key solvers, algorithms, and methods used in these previous studies.

Table 1 

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

3. Results

3.1. 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].

Figure 3 

1D Studied domain for verification.

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.

Table 2 

Mann-Whitney test parameters for comparison between V-Break and stoker’s analytical solution results.

(a)

(b)

Figure 4 

1D Water depth values obtained using V-Break and Stoker’s analytical solution. (a) 1D water depth at 𝑡=0.02 s, (b) 1D water depth at 𝑡=0.1 s.

(a)

(b)

Figure 5 

1D Velocity values obtained using V-Break and stoker’s analytical solution. (a) Velocities at 𝑡=0.02 s, (b) velocities at 𝑡=0.1 s.

In terms of comparing two groups of data statistically, the null hypothesis (𝐻₀) 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 |𝑍|≤𝑍₀, while 𝑍₀ 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.

3.2. 2D Partial Dam Break Test

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

Figure 6 

2D studied domain for verification.

Figure 7 

2D asymmetrical partial dam break test in a frictionless, horizontal, and wet/dry bed domain at 𝑡=7.2 s by some other researchers. (A1) Liang et al.; Wet bed [11], (A2) Liang et al.; dry bed at 𝑡=5 s [11], (B1) Loukili and Soulaïmani-Wet Bed [14], (B2) Loukili and Soulaïmani-Dry Bed [14], and (C) Wang and Liu [7].

Different conditions of the domain were simulated using V-Break. Figure 8 shows the output results of V-Break.

(a)

(b)

Figure 8 

2D asymmetrical partial dam break test and in a frictionless horizontal domain at 𝑡=7.2 s obtained by V-Break. (a) Wet bed and (b) dry bed.

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

(a)

(b)

Figure 9 

The comparison among Voronoi, rectangular [11], and triangular [14] mesh grids. (a) Wet bed and (b) dry bed.

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

Figure 10 

The comparison among various Riemann solvers used in Wang and Liu’s code (Reo-MUSCL, Reo-Upwind, HLL-MUSCL and CLFL8) and V-Break (LLxF).

3.3. 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, 𝑢 and 𝑣, 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.

Figure 11 

2D Circular dam break problem at 𝑡=2 s using V-Break.

Figure 12 

2D Circular dam break water depth contours obtained using V-Break.

Figure 13 

2D Circular dam break problem at 𝑡=2 s by Baghlani [21].

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

4. 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 𝑡=2 s 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 to accurate results within a lower computational demand compared to other unstructured meshes. Furthermore, a high capability of LLxF scheme was demonstrated in the 2D dam break pneumonia evaluation compared with the other presented Riemann solutions.

Nomenclatures

𝐴:	Adjacent surface vector of the investigated Voronoi cell
𝐴:	Area of the adjacent surface vector of the investigated Voronoi cell (𝐴)
𝐴𝑥:	The component of 𝐴 in 𝑥 direction
𝐴𝑦:	The component of 𝐴 in 𝑦 direction
𝐹,𝐺:	Flux vector functions
𝐹𝜀,𝐺𝜀:	The Voronoi cell normal flux vectors
𝑈:	The vector of conserved variables
𝐻:	Input and output fluxes to a Voronoi cell
𝑆:	The vector of source terms
𝑆0𝑥:	Bed slope in the 𝑥 direction
𝑆0𝑦:	Bed slope in the 𝑦 direction
𝑎,𝑏:	Nodes of the both sides of the investigated Voronoi cell 𝑓
𝑒𝜉:	Unit outward normal vector in each Voronoi cell 𝑓
𝑒𝜂:	The unit tangent vector in each Voronoi cell 𝑓
𝑓:	Joint surface element between investigated cell and other adjacent cells
𝑔:	Gravity acceleration
ℎ:	The water depth
ℎ∶	The mean water depth
ℎ𝑢:	The upstream water depth at 𝑡=0
ℎ𝑑:	The downstream water depth at 𝑡=0
𝑙:	The boundary of the 𝑖th control volume
𝑛:	The Manning’s roughness coefficient
𝑛:	The outward unit vector normal to the boundary
nb:	The central node of adjacent cells
𝑡:	Time
𝑢:	Velocity vector component in 𝑥 direction
𝑢:	The mean velocity vector component in 𝑥 direction
𝑣:	Velocity vector component in 𝑦 direction
𝑣:	The mean velocity vector component in 𝑦 direction
𝑥:	Horizontal coordinate component
𝑦:	Vertical coordinate component
∑𝑓∶	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
𝑍:	The value of 𝑍-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.

Superscripts

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