Abstract
The last two decades have seen great progress in mathematical modeling of fluvial processes and flooding in terms of either approximation of the physical processes or dealing with the numerical difficulties. Yet attention to simultaneously taking advancements of both aspects is rarely paid. Here a well-balanced and fully coupled noncapacity model is presented of dam-break flooding over erodible beds. The governing equations are based on the complete mass and momentum conservation laws, implying fully coupled interactions between the dam-break flow and sediment transport. A well-balanced Godunov-type finite volume method is used to solve the governing equations, facilitating satisfactory representation of the complex flow phenomena. The well-balanced property is attained by using the divergence form of matrix related to the static force for the bottom slope source term. Existing classical tests, including idealized dam-break flooding over irregular topography and experimental dam-break flooding with/without sediment transport, are numerically simulated, showing a satisfactory quantitative performance of this model.
1. Introduction
Mathematical modeling of fluvial flows and flooding has been a routine work in support of flood risk management and river management. Nevertheless, efforts to improve the quality of mathematical modeling for fluvial flows and flooding have never stopped due to their complexities: hydraulic jump/drop, wet/drying, sediment transport, bed deformation and heterogeneous sediment sizes, and so forth. This work is motivated by the general recognition that physically meaningful results in agreement with observations depend on not only accurate numerical algorithms, but also physically well-grounded model formulations [1].
In terms of numerical algorithm, a mathematical model for dam-break flooding should be capable of capturing the transitions between subcritical and supercritical flow regimes such as hydraulic jumps/drops (i.e., shock) and the moving wet-dry fronts. One of the widely used methods dealing with these is the Riemann solver-based technique: Godunov-type finite volume method [2]. This, at the same time, gives rise to an important issue: the well-balanced property of a model, which refers to the ability of the model to reproduce a body of static water over irregular topographies [3, 4]. It is because the Riemann solver-based technique estimates the hydrostatic pressure term (i.e., flux gradients) and the bed slope source term differently. This subsequently may make the two otherwise equal terms in the case of a body of static water become unequal/unbalanced over irregular topography. Fortunately, the last several decades have seen numerous novel techniques resolving this issue and wide applications can be found. Those include the surface gradient method [4], the flux correction method [5], upwind discretization of bed slope [3], divergence form for bed slope source term [6], wave propagation algorithm/augmented Riemann solver [7], prebalanced shallow water equations [8], and hydrostatic reconstruction technique [9].
Indeed, the development of these well-balanced numerical techniques for Godunov-type finite volume method greatly improved confidence in mathematical representation of the complex phenomena of the fluvial flows and the dam-break flooding. However, most of those well-balanced models focus on clear water flow [3–8, 10, 11], in which the potential sediment transport and high erodibility of the flow on the bed are ignored. For these models that have taken into account sediment transport, a capacity description is often adopted, in which the sediment transport rate/concentration is directly estimated by empirical relations [9, 12, 13]. In fact, however, sediment transport is governed also by physical laws, that is, the mass conservation law. In the noncapacity model, sediment concentration at a specific control volume is computed as a result of the net flux across the control volume and the sediment exchange between the flow and the bed. The advantages of the noncapacity modelling have been recently well recognized [14].
This paper presents a well-balanced and fully coupled noncapacity model for dam-break flooding over erodible beds. The governing equations are numerically solved using a second-order Godunov-type finite volume method: a predictor-corrector time stepping along with the HLLC approximate Riemann solver for flux estimation. The bed slope source term is rewritten in a divergence form of matrix related to the static force due to bottom slope, facilitating straightforward satisfaction of the (conservation)-property. Unstructured triangular grid system is used to represent the computational domain. The quantitative performance of the model was tested against classical idealized and experimental dam-break flooding benchmark tests.
2. Numerical Model
2.1. Governing Equations
The governing equations of the model comprise the mass and momentum conservation equations for the water-sediment mixture flow and the mass conservation equations, respectively, for sediment and bed material. Two-dimensional governing equations are written in a matrix form as follows: where = vector of conserved variables; , = vectors of flux variables; , = the diffusive vectors; = vectors of bed slope source term; = vectors of friction source term; = vectors of source terms representing the feedback impacts of sediment transport on the flow; = time; = streamwise coordinate; = cross channel coordinate; = flow depth; , = depth-averaged velocities in the - and -directions; = depth-averaged volumetric sediment concentration; = gravitational acceleration; is the turbulent viscosity coefficient, = bed shear velocity, = bed elevation; , = friction slopes in - and -directions; , = sediment entrainment and deposition fluxes across the bottom boundary of the flow; = bed sediment porosity; , = densities of water and sediment; = density of water-sediment mixture; and = density of saturated bed.
Based on (1), (2a), (2b), (2c), (2d), (2e), (2f), (2g), (2h), and (3), the model is categorized as a coupled model. By “coupled,” the meaning is twofold. The first feature is the nonequilibrium description of sediment transport by the mass conservation of sediment carried by the flow (the fourth component of (1) and (2a), (2b), (2c), (2d), (2e), (2f), (2g), and (2h)), which determines the sediment concentration. This is in contrast with the commonly used equilibrium description in the existing well-balanced morphological models, which assume sediment concentration equal to the sediment transport capacity of the flow. The second feature refers to simplifications in the mass and momentum conservation of the flowing water-sediment mixture: the first, second, and third components of the source term were set equal to zero in previous well-balanced morphological models.
2.2. Numerical Algorithm
Using the unstructured triangular mesh system, the governing equations are discretized by the finite volume method, and the interface numerical fluxes are estimated by the HLLC approximate Riemann solver. In order to achieve second-order accuracy in both space and time, the monotone upstream schemes for conservation laws (MUSCL) reconstruction are implemented before the numerical fluxes are estimated. In the MUSCL, the water surface level is used following the surface gradient method [4]: the flow depths are calculated by applying the MUSCL to the water surface level and then subtracting the bed elevation. The well-balanced property is achieved by using the DFB form for the bed slope source term.
2.2.1. Well-Balanced Property
Essentially, the issue of the well-balanced property arises from two distinct terms of the momentum equation: the flux gradient/hydrostatic pressure term and the bed slope source term. Take the case of a body of static water as an example: the two terms in the momentum equation of the -direction have the relation . While the left-hand-side (LHS) term (the flux gradient/hydrostatic pressure term) is computed by an approximate Riemann solver, the right-hand-side (RHS) term [the bed slope source term] is often estimated by centered discretization. The different estimation ways would easily give rise to inconsistency: the computed is not equal to in the case of a body of static water especially when the topography is irregular.
The present paper makes use of the method of DFB (the divergence form for bed slope source term) to resolve this issue. The DFB method was proposed by Valiani and Begnudelli [6] following the recognition that “the simplest way towards the numerical closure of physical balances is the divergence form of physical laws, which leads to conservation of physical properties such as mass, momentum, energy, and so on.” By DFB, the bed slope source term is written as the divergence form of a proper matrix related to the static force (the DFB method). Following Valiani and Begnudelli [6], one haswhere is the free water level and is a constant value for the free water level. It has been proved that reformulating the bed slope source term in the form of (4), (5a), (5b), and (5c) facilitates straightforward satisfaction of the well-balanced property irrespective of the specific numerical algorithm. It is noted that the governing equations of the present model involve sediment transport. However, the extension of this DFB method from clear water flow to the present sediment-laden flow is straightforward and is justified, because in the simplest yet challenging case of a body of static water over an irregular topography, governing equations for sediment-laden flow become the same as those for clear water flow.
Make use of (4), (5a), (5b), and (5c) and let and ; the governing (1) can be rewritten as
2.2.2. Godunov-Type Finite Volume Discretization
Figure 1 shows a sketch for the unstructured triangular mesh system. Bed elevation is defined at nodes of a cell. Flow variables are stored at the centre of each cell. Integrating (6) over the area of an arbitrary cell giveswhere = average free water level of the cell which is set equal to the flow depth plus the average bed elevation of the cell : , where = bed elevation of the th node of the cell , the subscript denotes the cell number, and the subscript denotes the edge number or node number of a cell. For a neighborhood description between cells, nodes of a cell, and neighboring cells of a cell, the model follows the integer mapping arrays used in Begnudelli and Sanders [5].
Applying Green’s theorem to the two matrix terms, the second term in the LHS and the first term in the RHS, (7) becomeswhere = time step, the superscript is the time step index, are vector of the cell-averaged conservative variables, is the unit outward vector normal to the edges of the cell, = length of the th edge of the cell , and and = source term vectors integrated over the cell .
Estimation of the source vectors and is straightforward as they are simply algebraic function of the flow variables, except the spatial gradient of sediment concentration, which is discretized by a centered difference. A semi-implicit method is adopted for computing in the present model, following Liang and Borthwick [8]. The bed slope source term , that is, the first term in the RHS of (8), is estimated following Valiani and Begnudelli [6]:where is the water depth at the midpoints of the th edge of the cell and = average bed elevation of the th edge of the cell .
For estimation of the intercell fluxes, that is, the second term in the LHS of (8), the HLLC approximate Riemann Solver is used, of which the details can be found in many literature reviews [2]. Thus, details for the HLLC approximate Riemann solver are not introduced here. For numerical flux estimation, two sets of conservative variables, and at the left and right sides of the interface of two neighboring cells, respectively, need to be estimated. One can directly use the cell-averaged conservative variables, which, however, is first-order accuracy only. To achieve second-order accuracy in both time and space, two steps are implemented to obtain and . The first step advances the solutions from time to for second-order accuracy in time [5], and the second step achieves second-order accuracy in space by making slope-limited linear extrapolation of the half-time level solution from the cell center to the cell boundaries. For the details of the MUSCL method in unstructured triangular grid system, one can refer to Begnudelli and Sanders [5].
2.2.3. Boundary Conditions
The boundary conditions used in this model include two types, that is, closed and open boundaries. At a closed boundary, a slip condition is used, namely, = 0, , , and . If an open boundary is involved, the method of characteristics using the Riemann invariants will be implemented along with the boundary conditions to obtain flow variables at the boundaries [8, 15].
2.2.4. Bed Evolution
Cell-averaged bed evolution is computed from (3) with the state information due to (8): in which is the starting time level. The computed for each cell is distributed to each node of the cell so that bed elevations at each node can be updated [16].
2.2.5. Time Step
The time step is constrained by the Courant-Friedrichs-Levy condition where = cell numbers, = wave speed normal to the th face of the th cell, and = Courant number. The Courant number = 0.95 is used for the numerical cases here.
2.3. Empirical Relations
The friction slopes are determined using the Manning roughness : The sediment entrainment and deposition fluxes are estimated aswhere = settling velocity of a single sediment particle in tranquil clear water, which is calculated using Zhang’s formula [17]; the exponent , , is the kinematic viscosity of water, , = particle diameter, and = sediment concentration at capacity. The coefficient represents the difference between the near-bed concentration and the depth-averaged concentration, and usually larger than unity. Physically the near-bed concentration must not be larger than (). In the present model, it is specified as . The formulation of Zyserman and Fredsoe [18] is deployed for estimating sediment concentration at capacity aswhere is the Shields parameter, is the critical Shields parameter for incipient sediment motion, and = bed shear velocity.
3. Numerical Case Studies
In this section, numerical case studies are presented to demonstrate the performance of the model.
3.1. Steady Transcritical Flow over a Bump
This numerical case is used to test whether the present model can simulate steady flow over irregular topography, the main indication of satisfying the well-balanced property. The computational domain is a 25 m long one-dimensional channel. The bed elevation is defined as . Initially the flow is static in the channel with water level m. A steady unit-discharge of 0.18 m2/s is imposed at the upstream boundary (m) and the water depth is kept constant at the downstream boundary (m). Figure 2 presents a comparison of the computed water level by the present well-balanced model against the analytical solution. A good agreement of the numerical results by the well-balanced model with the analytical solutions is observed. Also included in Figure 2 is the computed water level by the model of Cao et al. [19], in which the well-balanced issue was not accounted for. Obviously, the present model is well-behaved in preserving the steady state, whereas the previous model [19] is not.
3.2. Idealized Dam-Break in a Closed Dry Channel with Irregular Rough Bed
The present model is now tested against a widely used idealized case of dam-break flow in a closed dry channel with irregular rough bed. Dam-break flow requires the model to be capable of capturing shock waves. The dry property of the channel requires the model to deal with wet/dry fronts satisfactorily. The irregular bed requires the model to be well-balanced. Dissipation of energy due to friction of the rough bed will finally lead to a stationary state of the dam-break flow, which combined with the irregular feature of the bed provides challenging testing for the well-balanced property.
The closed channel is 75 m long and 30 wide. A dam is arranged 16 m to the right of the left end ( = 16 m). Upstream of the dam is still clear water with surface elevation of 1.875 m, and downstream of the dam is initially dry bed. The bed topography is defined by By the above equation, there are two low humps and one high hump downstream of the dam. The manning roughness of the channel is set as = 0.018. Computational grids (in total 20480 cells) are symmetrical along centre lines. The computation starts with the instant break of the dam and ends with the establishment of the stationary state.
Figure 3 shows the 3D view of the evolution of the dam-break flow over the three humps at six selected times. Figure 4 shows the corresponding velocity distributions. The following are observed from Figures 3 and 4. Dam-break flow induces a bore propagating downstream. In response to the arrival of the bore, the two low humps are fully submerged and hydraulic jumps occur upstream of the humps ( = 5 s and 10 s). Along with further downstream propagation of the dam-break flow, the high hump becomes submerged and the high water around the two low humps recedes ( = 15 s and 20 s). When the whole channel is covered by water, the flow advances and recedes due to the wall reflections ( = 20 s and 30 s). During these processes, flow energy is dissipated and flow velocity decreases (Figure 4) until a stationary state is established. To test if the model simulation is stable, the computation is continued for sufficiently long time, showing that after about ten hours the stationary state is stable ( = 10 h). In contrast, as shown in Figures 5 and 6, the numerical solutions by the model of Cao et al. [19] that did not consider the well-balanced property do not converge to the final static state. These results suggest that the present model can capture shock waves and deal well with wet/dry fronts, satisfying the well-balanced property.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
3.3. Experimental Dam-Break in a Converging-Diverging Channel
Bellos et al. [20] presented a set of experiments for two-dimensional dam-break flood waves in a converging-diverging channel. Different flow conditions and bed slopes are considered in the experiments, and depth hydrographs at eight locations are measured. Figure 7 shows a sketch for the experimental flume. The length of the experimental flume is 21.2 m and maximum width is 1.4 m, and the Manning roughness is 0.012, the dam located at = 0 m. Initially, the depth upstream of the dam site is 0.15 m, and bed is dry downstream of the dam. The experiment simulated in the present paper is characterized by a bed slope of 0.002.
Figure 8 shows comparison of calculated and measured water depth. Upstream of the dam the water stage at P1 and P2 decreases after the dam failure, whereas those at P3 and P4 experience an increase before decreasing. The good agreements between the computation and the observation indicate that with sensibly specified parameters the quantitative accuracy of the present mathematical model in numerically reproducing dam-break flows is sufficiently high.
(a) P1
(b) P2
(c) P3
(d) P4
3.4. Experimental Dam-Break Flow in an Abruptly Widening Erodible Channel
At the Civil Engineering Laboratory of the University of Catholique de Louvain (UCL), Belgium, a series of experimental dam-break flows in an abruptly widening erodible channel are conducted [21, 22]. Here the experiment conducted in a 6 m long channel with a sudden asymmetric enlargement of the channel width from 0.25 m to 0.5 m is numerically simulated. A plane view of the experimental configuration is shown in Figure 9. A layer of saturated sand was put in the flume bottom to form a mobile channel: 0.1 m thick, particle diameter = 1.82 mm, sediment density = 2680 kg/m3, and sediment porosity = 0.47. A dam is constructed at = 3 m (see Figure 9). Upstream of the dam is initially still water with a water stage of 0.25 m; downstream of the dam is initially dry bed. Experiment is initiated by rapidly moving down the dam gate so that the dammed water is poured into the downstream dry bed. At the downstream end of the flume, free outflow condition is achieved by installing a weir there. Measurements of the water stage are made at six selected positions, indicated as P1 ( = 3.75 m, = 0.125 m), P2 ( = 4.2 m, = 0.125 m), P3 ( = 4.45 m, = 0.125 m), P4 ( = 4.95 m, = 0.125 m), P5 ( = 4.2 m, = 0.375 m), and P6 ( = 4.95 m, = 0.375 m). Measurements of the final bed profile are made at two cross sections, indicated as CS1 ( = 4.1 m) and CS2 ( = 4.4 m), respectively. The manning roughness is set as = 0.024; the threshold Shields parameter for incipient sediment motion is set as = 0.04.
Figures 10 and 11 present the comparisons between the numerical solutions and the measured data for this case: time variations of the water stage at the six positions (Figure 10) and transverse profiles of bed deformation depth at the two cross sections (Figure 11). From Figures 10 and 11, the time histories of the water stage at the six positions are well reproduced as compared to the measurement, whereas the agreements between the computed and measured bed deformation are less satisfactory. Nevertheless the qualitative trend of bed deformation is well captured. The appreciable quantitative differences between the computed and measured bed deformation depth are understandable considering the complexity of sediment transport, which is similar to previous numerical simulations [23].
(a) CS1
(b) CS2
3.5. Experimental Partial Dam-Break in a Straight Erodible Channel
In this section the present model is further tested against an experimental partial dam-break flow in a straight erodible channel, which is also conducted at UCL, Belgium. This experiment differs from the one used in the previous section mainly in two points: the dam only collapses in the center part with a width of 1 m (see Figure 12 for the experimental configuration); the width of the channel is constant. As shown in Figure 12 of the plane view, the channel is 3.6 m wide and 36 m long. The coordinate (0, 0) is arranged at the center of the dam. Upstream of the dam is still water with a water stage of 0.47 m and downstream of the dam is initially dry bed (this paper only simulates the dry bed case as it is more difficult to be numerically reproduced as compared to the wet bed case; see Wu et al. [24]). Initially, a layer of sand sediment is put at the flume bottom extending from 1 m upstream of the dam to 9 m downstream of the dam so that an erodible bed is formed: 0.085 m thick, particle diameter = 1.61 mm, sediment density = 2630 kg/m3, and bed sediment porosity = 0.42. The roughness is set to be 0.0165 for the mobile part and 0.01 for the fixed bed portion [23, 24]. The experiment is initiated by rapidly lifting the center gate of the dam (1.0 m wide). The experiment is stopped 20 s after the dam-break to avoid reflection from the downstream end of the flume. Measurements of the water stage are conducted at 8 positions: P1 ( = 0.64 m, = −0.5 m), P2 ( = 0.64 m, = −0.165 m), P3 ( = 0.64 m, = 0.165 m), P4 ( = 0.64 m, = 0.5 m), P5 ( = 1.94 m, = −0.99 m), P6 ( = 1.94 m, = −0.33 m), P7 ( = 1.94 m, = 0.33 m), and P8 ( = 1.94 m, = 0.99 m). These positions distribute asymmetrically in the two sides of the dam. After the experiments, longitudinal bed profiles are measured at three positions ( = 0.2, 0.7, and 1.45 m). The experiment is repeated so that a second measurement is made to ensure the accuracy of the measurements.
Figure 13 presents the comparisons between the computed and measured time variations of the water stages at the eight positions. The measured data are averaged values of the two repeated experiments. Positions are arranged in pairs (P1 and P4; P2 and P3; P5 and P8; P6 and P7) in the comparison because the positions are symmetrical along the longitudinal center line of the flume. Not surprisingly, the computed water stages for each position of a pair are identical, indicating good performance of the model. While differences between the measurements and the numerical solutions are obvious, the computation reasonably captures the magnitude and trend of the time variations of the water stage. For positions closer to the dam (P1, P2, P3, and P4), the discrepancy between the computation and the measurement is more obvious than those positions further away from the dam (P5, P6, P7, and P8). This is due to the 3D effects of the abrupt change in the plane form of the channel, which may be not well accounted for in the model.
(a)
(b)
(c)
(d)
Figure 14 shows the longitudinal bed profiles ( = 0.2 m, 0.7 m, and 1.45 m) after the experiments including the numerical solutions and two repeated measurements. It is seen that the extent of the discrepancy between the computation and experiments is similar to the discrepancy between the two repeated measurements. This indicates good numerical reproduction of the bed deformation change.
(a)
(b)
(c)
4. Conclusions
A well-balanced and fully coupled noncapacity 2D depth-averaged model for sediment-laden dam-break flooding over erodible beds is developed. The satisfaction of the C-property (well-balanced property) of the model is demonstrated through a case of steady flow over a bump and a case of dam-break flow over three humps. The numerical accuracy of the model in terms of reproducing key flow variables (i.e., water level) and bed deformation is demonstrated against three sets of experimental dam-break flows over movable beds.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the Basic Research Program of Yangtze River Waterway Research Institute (2011-02-029), the Research Fund for Doctoral Program of Higher Education of China (20130101120152), Open Fund of the State Key Laboratory of Satellite Ocean Environment Dynamics (SOED1309), and the National Natural Science Foundation of China (11402231).