Journal of Fluid Flow, Heat and Mass Transfer (JFFHMT)

ISSN: 2368-6111

1. Introduction

Dam-break flows over irregular bed often experience the cases as transcritical flows, steep bed slope, very small water depth, cells' wetting and drying, wave propagation. Inappropriate hydraulic simulation may leads to inaccurate and unstable prediction. Hence, it's necessary to apply a robust, accurate and effective hydraulic model to predict Hydraulic parameters of dam-break flow.

By using Godunov-type scheme, complicated shallow water flow phenomena such as transcritical flows, shock-type flows and moving wet-dry interface of a water wave front can be appropriately simulated. To solve a Godunov-type scheme, approximate Riemann solver is normally adopted to estimate the numerical fluxes of hyperbolic system and various numerical schemes have been proposed to solve Riemann problems. As explained in the next section, the Central-upwind method which is applied in current study, does not require computationally expensive decomposition of numerical flux on the basis of eigenvalues and furthermore, only an estimation of largest and smallest eigenvalues of the Jacobian matrix, which leads to a significant reduction of computational cost.

Flow over initially dry bed is very common in dam break flows, which involves complicated boundary conditions. For irregular topography, both positive and negative bed slopes generally exist, which may leads to cell drying and wetting with moving fronts which could not be easily solved by horizontal boundary condition. Some techniques have been developed in using finite volume method and shallow water equations. Zhao et al. [19] and Sleigh et al.[14] introduced two similar schemes to track the wetting and drying fronts, in which cells are divided into wet, dry and partially dry types according to two tolerances. Brufau et al. [2, 3] proposed a technique using unsteady wetting and drying conditions in flows and claimed that the method gave zero mass error, which is valid for an FVM with only first-order accuracy. Falconer et al. [5] and Falconer et al. [6] developed a wetting and drying method for regular grid finite difference model, which is recently refined for triangular grids by Xia et al. [17]. In present study, a technology for tracking wet-dry front is developed combing with method of Brufau et al. [3] to achieve zero mass error.

It is well-known that accuracy is the most important aspect for flow solver since it has a direct influence on the number of computational cells required. This means that a higher-order implementation involving a piecewise linear reconstruction is necessary. Higher order schemes often produce nonphysical oscillations which can be effectively suppressed by using limiters. In the present study, the multidimensional limiter proposed by Jawahar et al. [9] is adopted to calculate the limited gradient for reconstruction of variables. Moreover, a second order accuracy in time could be obtained in current model by apply Runge-Kutta scheme.

Based on the efficient divergence form of the slope source term proposed by Valiani et al. [16], Hou et al. [8]developed a novel slope source term treatment which is devised to transform the slope source of a cell into a flux form. By splitting the integral of the bed slope source term over a cell into those of the sub-cells, higher accuracy can be achieved by the novel treatment than that proposed by Valiani et al. [16]. This method can strictly preserve the well-balanced property and can be conveniently employed with second order or even higher order schemes. In addition, this treatment is able to handle the occurrence of wet-dry fronts, in conjunction with the non-negative water depth reconstruction.

2. Numerical Scheme

The 2-D shallow water is adopted in the current study as governing equations which can be described as:

The 2-D shallow water equations constitute a hyperbolic system which can be presented in the following vector form:

in which t = time, h = the water depth, x and y are horizontal coordinates, u and v are the depth-averaged velocity components in x and y directions respectively, Z is the bed elevation, g is the gravitational acceleration, v_tis eddy viscosity coefficient given by , n_b is Manning's roughness coefficient. and are friction slope terms in the x, y directions, respectively, which can be determined by conventional formulas involving Manning roughness coefficient n_b, , . Besides, in this work, water surface level η is used in the second order spatial reconstruction and the non-negative water depth reconstruction. η can be calculated as .

2.1. Discretization of Flow and Sediment Governing Equations

The coupled system is discretized on an unstructured triangular grid by finite volume method.

in which n_x and n_y are components of unit normal vector in the x and y directions, l_ik is the length of the kth edge of control volume i, F_ik is numerical convective flux across l_ik, is the diffusive flux across l_ik.

2.2. Central-upwind scheme for the interface flux

Central-upwind schemes on general triangular grids for solving two-dimensional systems of conservation laws are developed by Kurganov et al. [11], which enjoy the main advantages of the Godunov-type central schemes, i.e. simplicity, universality and robustness and can be applied to problems with complicated geometries. The triangular central-upwind schemes are based on the use of the directional local speeds of propagation and are a generalization of the central-upwind schemes on rectangular grids, introduced in [10].

Applying the Central-upwind scheme, the convective fluxes in Eqs. (7) could be estimated by:

in which and are normal fluxes on the left and right sides of the kth edge, respectively. S_Rand S_L are one-sided local speeds of propagation on right and left sides of kth edge, respectively, and can be determined by

with being the N eigenvalues of the Jacobian matrix using reconstructed variables or .

It can be deduced that, when a cell is dry or nearly dry, exists in denominator of equation Eqs. (8) which may cause numerical instability. In order to handle this situation that both S_Rand S_L are zero (or very close to zero), following the suggestion in Bryson et al. [4], the scheme Eqs. (8) reduces to:

2.3. Spatial linear Reconstruction

In this paper, the following 2-D linear reconstruction is employed:

in which is the reconstructed value of variables at point inside of cell i, and are component-wise approximation of numerical derivatives which are computed via a nonlinear limiter, used to minimize the oscillations of the reconstructions. In current study, the multidimensional limiter proposed by Jawahar et al. [9] is adopted to calculated limited gradient within a cell by taking the weighted average of three representative unlimited gradients.

In order to preserve the well-balanced property for second order schemes, as suggested by Audusse et al. [1], surface levels , , water depths , , flow discharges , , and are reconstructed at the midpoint of considered edges. The reconstructed bed levels at M are given by

Besides, flow velocities required at M are computed as

A threshold ε will be introduced for defining wet and dry cells. The second order reconstructions are only applicable to the wet cells.

2.4. Non-negative water depth reconstruction

A robust and efficient reconstruction approach to preserve non-negative water depth which is suggested by Liang et al. [12] is adopted in present study. The Z^M bed elevation at the midpoint of the considered edge is calculated as:

Then the non-negative water depth values on both sides are reconstructed as:

The discharges at M are in turn reconstructed as:

2.5. Treatment of the source term

2.5.1. Well-balanced treatment of the slope source term

Hou et al. [8] developed a novel slope source term treatment which is devised to transform the slope source of a cell into a flux form, which can strictly preserve the well-balanced property and can be conveniently employed with second order or even higher order schemes. In addition, this treatment is able to handle the occurrence of wet-dry fronts, in conjunction with the non-negative water depth reconstruction. The vector of source term at the considered faces k in cell i becomes

where Z^M and are obtained by spatial linear reconstruction.

2.5.2. Calculation of derivatives in source terms

The approach adopted by Mohammadian et al., [13] is used to calculate the unlimited gradient of variables. A similar approach may be also applied to calculate the diffusive terms , , and .

2.5.3. Friction source term treatment

In current study, a simple semi-implicit treatment suggested by Yoon et al. [18] is adopted to deal with very shallow water depth.

in which the superscript n denotes the time step.

2.6. Time integration and CFL condition

In order to obtain the second-order accuracy in time and retain the stability of the model, the two-stage explicit TVD Runge-Kutta method is applied in current study. The time step is limited by the Courant-Friedrichs-Lewy (CFL) condition.

2.7. Treatment of wetting and drying fronts

Wetting and drying fronts need to be specially treated to retain the numerical stability when very tiny water depth is introduced. In current model, a scheme of wetting and drying treatment is proposed which can be summarized as:

1. A tolerance water depth ε is introduced to classify the wet and dry cell. In present model, is used. In the dry cell, first order reconstruction will be applied to maintain the numerical stability.

2. For dry cell with , only continuity equation will be solved to deal with potentially wetting.

3. To deal with the cell-drying, if a cell with , it's treated as completely dry cell in next time-step and all hydraulic parameters are set to zero in this cell.

3. Numerical Tests

3.1. Quiescent water around a hump with sediment deposition

The first test is applied to verify the well-balanced property preserving of present model involving wet-dry interfaces. To implement this test, the elevation of a hump on flat bed is defined as:

The hump is located at the center of a 8m×8m computational domain, the height of the bump is 2m. A quiescent lake around the dump is defined with initial water surface elevation of 1m, so the test could involve the wet-dry interface. Fig. 1 shows the 3d views of simulated bed profile and still water surface at t=50s. The undisturbed water surfaces are observed through the whole simulating process.

3.2. 2-D shorelines tracking in a parabolic bowl

This test is adopted here to investigate the hydraulic model and the accuracy of tracking the wetting and drying fronts. The bottom topography with the center (x₀,y₀) is defined as:

in which h₀ is the water depth at the center of the domain, a is the distance from the center to the edge of the shoreline. Since the topography is set to be frictionless in this test, as described in [15], the periodic analytical solution of the evolutions of surface elevation, water depth and velocities can be computed using following equations:

in which, σ is a constant that determines the amplitude of the motion; is the frequency of the pool's circulation around the center of the bowl.

In this test for current model, a 4m×4m domain is chosen, which is discretized to 14836 triangular cells, as shown in Fig. 2. The four boundaries are set to wall-boundary. According to the dimension, the parameters h₀=0.1 m, a=1.0m and σ=0.5m. The initial flow states are defined by Eqs. (20) and Eqs. (21) at t=0s.

Fig. 3 shows the comparison between simulated water surface profiles and analytical solutions at t=3T and 3.5T, respectively, where T represents the circulation period of the pool. It can be observed that the calculated free surface from current model agrees well with the analytical solution, no obvious distortion is observed near the shorelines, the treatment of wetting and drying boundaries successfully handle the task of tracking moving wet-dry fronts. Fig. 4 shows the velocity of circulating pool at t=3.5T, no unphysical high velocity is observed near the shoreline where very shallow water exists. The around area preserve zero velocity when it becomes dry.

3.3. 2-D partially dam-break flow on initially dry bed

For examining the numerical performance of the present scheme, a 2-d dam break problem with rapid varying unsteady flow is chosen as test cases. This test case is firstly introduced by Fennema et al. [7] in a numerical method study, which has been widely used by many researchers. The original case is mostly applied on an initially wet and fixed bed. It is adopted here for the aim of testing the capacity of current model to simulate the front wave propagation over an initial dry bed with wet-dry interface tracking, with a particular attention to the 2D aspects of the flow motion; and the model performance in prediction of fast erosion under high energetic flow in large scale. The modeling domain area is a 200m×200m basin over a flat and dry bed. A 10m-thick dam splits the basin into two equal-sized regions. The water depths are 10m and 0m on the left and right sides of the dam, respectively. At t = 0s, a 75m wide breach centered at y = 125m is assumed to form instantaneously. The duration of the simulation is 12s. The initial velocity of the whole modeling area is 0m/s and the outlet boundary at x = 200m is specified with a free out flow boundary condition, meanwhile all other boundaries are set to be standard wall condition. The simulating domain is discretized into 6806 triangular cells, as shown in Fig. 5.

Fig.6 shows 3D views of the water front-wave propagations and predicted bed erosion at varying times (first row). It can be seen that a shock wave forms and propagates downstream and a depression wave spreads upstream during the simulating process, no unphysical high velocities and oscillations are generated at the front-wave locations where wetting and drying interfaces are existing. The velocity field with contour of water depth is present in second row.

4. Conclusion

A robust two-dimensional numerical model for dam-break flow with wetting and drying is proposed in current study, based on finite volume method using unstructured triangular grids. The Central-upwind scheme can accurately estimate the convective fluxes. Current model can strictly preserve the well-balanced property with vector-form discretization of slope source term in conjunction with nonnegative water depth reconstruction. The semi-implicit method on friction slope term maintains the stability of present model. The adopted wetting and drying scheme could efficiently track the wet-dry fronts. The testing results confirm the capacity and accuracy of current model in dealing with various cases of dam-break flow over irregular bed in practical conditions.

Acknowledgements

This publication was made possible by NPRP grant #4-935-2-354 from the Qatar National Research Fund (a member of Qatar Foundation). The statements made herein are solely the responsibility of the authors.

References