Abstract
A second-order accurate, Godunov-type upwind finite volume method on dynamic refinement grids is developed in this paper for solving shallow-water equations. The advantage of this grid system is that no data structure is needed to store the neighbor information, since neighbors are directly specified by simple algebraic relationships. The key ingredient of the scheme is the use of the prebalanced shallow-water equations together with a simple but effective method to track the wet/dry fronts. In addition, a second-order spatial accuracy in space and time is achieved using a two-step unsplit MUSCL-Hancock method and a weighted surface-depth gradient method (WSDM) which considers the local Froude number is proposed for water depths reconstruction. The friction terms are solved by a semi-implicit scheme that can effectively prevent computational instability from small depths and does not invert the direction of velocity components. Several benchmark tests and a dam-break flooding simulation over real topography cases are used for model testing and validation. Results show that the proposed model is accurate and robust and has advantages when it is applied to simulate flow with local complex topographic features or flow conditions and thus has bright prospects of field-scale application.
1. Introduction
Solution to the two-dimensional (2D) hyperbolic conservation laws of the shallow-water equations (SWE) is relevant to many real-life hydrodynamic problems, such as free surface flows in shallow lakes, dam breaks, dyke breaches, wide rivers, flash floods, tidal bores, and coastal inundation, among others. Meanwhile, a variety of methods can be used to solve the nonlinear shallow-water equations, such as variational principle [1, 2], finite element method [3], and finite volume method (FVM) [4–8]. HE [1] deduced a semi-inverse method for establishment of a generalized variational principle which applies to any conservation equations of hydrodynamics, and this method is helpful to study the analytical solutions which are useful in shallow-water equations. In recent years, with the development of numerical methods, more and more attentions have been paid to the FVM which are based on the conservative form of the governing equations. In particular, the upwind Godunov-type schemes via the FVM are very popular because of their robustness and accuracy in capturing discontinuities.
When solving the shallow-water equations in the context of a Godunov-type framework, it is generally necessary and efficient to use the well-balanced methods, which employ a special treatment for the bed slope source terms that balanced the source terms and flux gradients. This well-balanced concept is also known by exact conservation property (C-property) [4]. Several authors have conducted this study; see, for example, [4, 6, 7]. Among these, the prebalanced SWEs derived by choosing the water level and discharge as dependent variables are derived by Liang and Borthwick [6]. However, near the wet/dry fronts, negative water depth would be produced which could cause the mass error. To deal with this problem, Liang and Marche [8] then proposed the nonnegative reconstruction of Riemann states and compatible discretization of the slope source term to maintain stable and well-balanced property for dry-bed simulations. Begnudelli and Sanders [9] adopted flux correction terms to preserve the well-balanced property in a fully wet cell; however, the correction terms were later proved to be a potential source of large unphysical velocity [10].
Moreover, it is important that the numerical shallow-flow models are able to simulate accurately the evolution of relatively small-scale features, such as fronts, within large complex domains. Therefore, various approaches to adaptive mesh resolution for improving computational efficiency have been suggested for the shallow-water equations in the published literatures. Hu et al. [11] briefly reviewed the different mesh patching methods and presented a new patching technique for solving the SWEs using a finite volume Godunov-type scheme. Rogers et al. [12] and Liang and Borthwick [6] presented typical adaptive quadtree solvers of the shallow-water equations to achieve local grid refinement. Sleigh et al. [13] developed a finite volume method model based on adaptive unstructured triangular grids, which are widely used in the domains with highly complicated boundary configurations due to their capability of robustness. However, it is noteworthy that most of the above adaptive grids generation techniques require data structure to store their neighbour information that increases the computational cost. All in all, it is reasonable to anticipate that the ideally grid generators should be robust and efficient and fit the flow boundaries correctly.
Various numerical methods developed for general systems of hyperbolic conservation laws have been applied to achieve second-order accuracy in space, on structured and unstructured meshes. These methods provide efficient approaches to accommodate sub-, super-, and transcritical flows as well as flows with discontinuities. The conventional depth gradient method (DGM) approach which evaluates intercell water depths starting from the extrapolation of the same conserved variable can obtain a more robust behavior in the front tracking. However, if a centered discretization is used for the bed slope source terms, the depth gradient method may not maintain the static condition. Zhou et al. [7] proposed surface gradient method (SGM), which is suitable for both steady and unsteady shallow-water flows and is easy to implement. A central-upwind scheme using the surface elevation instead of the water depth has been used by Kurganov and Levy [14]. The SGM is efficient, but near the wet/dry fronts on uneven topographies the SGM reconstruction can lead to very small or negative predicted water depths that may cause numerical instabilities. Therefore, there is a need to modify the numerical scheme in order to avoid unphysical results.
It is significant to design a numerical approach to handle wetting and drying, which is one of the central challenges that researchers are tackling. In Godunov-type finite volume schemes, the depth-averaged velocity components are normally determined by dividing the discharge per unit width by the local water depth which can inevitably lead to predictions of unacceptable negative water depth and numerical instability. Liang and Borthwick [6] presented a Godunov-type shallow-flow model to solve the well-balanced SWEs with wet/dry fronts by a local bed slope modification method. Brufau et al. [15] presented a numerical technique that temporarily modifies ground elevation to approximate wetting and drying for both steady and unsteady flows. Begnudelli and Sanders [9] tracked fluid volume and the free surface elevation in partially submerged cells for wetting and drying and developed a high-resolution, unstructured finite volume model for shallow-water flows over arbitrary topography. The fully implicit treatments of friction source terms are employed by Yoon and Kang [16] to alleviate the problems associated with numerical instabilities due to small water depths near the wet/dry fronts. Moreover, near the wet/dry interface, where the depth is vanishing, the friction terms may cause the stiff problem and significantly affect the stability of the numerical model [17]. Therefore, it is necessary to employ a robust scheme, such as the semi-implicit scheme [9, 17], for friction terms when modeling the real-world shallow flow over irregular terrain.
In this paper, we describe details of a computationally efficient numerical approach to simulating complex shallow flow. A structured, but nonuniform, Cartesian grid system which allows a simulation to be carried out in dynamic refinement grids is used for solving the SWEs. A weighted surface-depth gradient method (WSDGM) is adopted with the aim of combining the advantages of SGM and DGM approaches and avoiding their drawbacks. An efficient and robust technique is used to track the stationary or move the wet/dry fronts. The friction terms are treated by a semi-implicit scheme, which may relieve the problem of high velocity when the water depth is small and does not invert the direction of velocity components and thus would benefit model robustness. This paper is organized as follows. In Section 2, the prebalanced shallow-water equation is derived. Section 3 presents the nonuniform gridding technique. Section 4 describes the numerical solution to the equations on dynamic refinement grids using a second-order Godunov-type finite volume scheme. Finally, the results of model validation tests and a simulation of a flood wave are drawn in Section 5.
2. Shallow-Water Equations
The two-dimensional depth-integrated shallow-water equations are obtained by integrating the Navier-Stokes equations over the flow depth with the following assumptions: hydrostatic pressure distribution, small bottom slope, uniform velocity distribution in the vertical direction, incompressible fluid, and constant fluid density with free surface. In matrix form, the conservation laws of the two-dimensional nonlinear shallow-water equations are where is the time; and are Cartesian coordinates; , , , and are vectors of the conserved flow variables, fluxes in the - and -directions, and source terms, respectively. Traditionally, the vector terms are given by where is the water depth; and are the depth-averaged velocity components in the - and -directions, respectively; is the gravity acceleration; and are bed slopes in the - and -directions, respectively; and are friction slopes in the - and -directions, respectively. In this study, the friction slopes are estimated by using the Manning formulae in which is Manning’s roughness coefficient. In cases involving wetting and drying, any grid cell with is assumed to be dry.
Note that , in which is defined as the water surface elevation above datum. Then, can be rewritten as . Simple algebraic manipulation leads to A similar analysis can be applied in the -direction momentum equation. Based on this, (2) will be transformed into another formulation constituted by (1) and
Consider the -direction momentum equation at the initial steady state in which is a constant for a wet-bed application with . After eliminating terms with zero velocities, the momentum equation reduces to When using HLLC approximate Riemann solver based Godunov-type numerical scheme to calculate the interface fluxes, the term on the left-hand side of (6) is approximated by , where and are fluxes through the east and west cell interfaces, respectively. Since is a constant and , it is easy to prove that no matter which calculation formula in HLLC solver is used, the fluxes are Based on (7), the flux gradient is If the bed slope term is discretized as , (6) is balanced. A similar analysis can be applied in the -direction momentum equation. Therefore, (1) and (5) are in flux gradient and source terms balanced form without any need to use special numerical treatment for the source terms.
3. Computational Grid
3.1. Structured but Nonuniform Grids
Because of the complexity of computational domain, a simple structured rectangular mesh requires a large number of cells to resolve the detailed flow pattern near the dam breaking location, navigation channels, and in-stream structures. To optimize the use of computational resources, we use the multiple-level nonuniform rectangular mesh with local refinement proposed by Liang [18].
In this mesh system, various levels of fine cells are placed close to the dam-break locations and in-stream structures or inlet areas where the flow gradients are high, while coarse grids are used in the low-gradient regions. To simplify the mesh, a cell is refined by splitting into four equal child cells. Corresponding to this refining, the computational mesh will be simpler because any cell has one or two faces on each of its neighbors, as shown in Figure 1.
To obtain a structured but nonuniform grid is straightforward and only involves the following four simple steps.(1)Discretize the computational domain using a coarse uniform grid.(2)Divide the unit square into four equal quadrant cells according to specific subdivision criteria.(3)Check each cell in turn and subdivide it if necessary.(4)Regularize the grid to ensure that an arbitrary cell is at most two times bigger or smaller than any of its neighbors (2 : 1 rule).
The first step is simply to determine the number of cells in the - and -directions, respectively. In the second step of subdividing the unit uniform grid (root grid) into four quadrants (child cells), seeding points are adopted or certain criteria may be used in this work. Then, each cell (including child cells) is checked in turn and a prescribed highest subdivision level is allocated to those cells containing seeding points. The procedure essentially creates an irregular grid, where a cell may have a neighbor that is many times bigger or smaller than its own. The stability and accuracy of the solution may be affected if no action is taken. In order to reduce grid irregularity, the 2 : 1 rule is satisfied and provides a trade-off between grid flexibility and accuracy.
Figure 1 shows the different levels of cells , , and , where and represent the cell indexes in the - and -directions, respectively. Considering cell , as Figure 1 shows, the subgrid is indexed by and , where is the number of subgrid cells in the - or - direction. with lev being the subdivision level of cell . The advantage of this nonuniform grid system is that no data structure is needed to store the neighbor information, neighbors are directly specified by simple algebraic relationships, and details on implementing the neighbor interpolation are well documented in the literature [18].
3.2. Dynamic Grid Adaptation
In this work, the grid adaptation criterion is based on the averaged water surface gradient and the critical value of and should be prescribed for grid refinement and coarsening. The averaged water surface gradient at cell is defined as [6] The following procedures are implemented for grid refinement and coarsening.(1)Check every cell in turn; if any cell has and the subdivision level is less than the maximum, the cell will be marked for further subdivision and its subdivision level will be increased to lev +1. When a cell is marked for subdivision, its eight neighbors are also checked and subdivided if necessary in order to meet the 2 : 1 rule.(2)For grid coarsening, if the cell has and the subdivision level is greater than 0, the cell will be marked for coarsening and its subdivision level will be changed to lev −1. Again the 2 : 1 rule must remain satisfied after grid coarsening.
It is worth noting that the above adaptation procedure is performed at every time step during a simulation, and the flow information of a newly created subcell is obtained by the aforementioned interpolation method.
4. The Numerical Method
4.1. Finite Volume Method
Integrating (1) over an arbitrary cell, the integral form thereof is in which all the vector terms are given by (5). By applying Green’s theorem to the second term in (10), where is the boundary of and is the vector of flux functions through given by in which and are the Cartesian components of the unit normal vector . The FVM is based on (11), the discretization of which ensures that mass and momentum are conserved across a discontinuity. Discretizing (11), the equation becomes where superscript is the time level, subscript represents the cell index, is the time step, and are the fluxes through the east and west cell interfaces, respectively, and and are the fluxes through the south and north cell interfaces, respectively.
4.2. HLLC Approximate Riemann Solver
To capture shock waves, reduce the numerical damping effect and overcome the numerical oscillation problems; the Harten-Lax-van Leer-Contact (HLLC) Riemann solver is employed to solve the local Riemann problems and calculate the interface fluxes. An interface flux, such as , is computed [6]: where and are calculated from the left and right Riemann states and , respectively, for a local Riemann problem and and are the numerical fluxes in the left and right sides of the middle region of the Riemann solution, respectively, given by where and represent the left and right tangential velocity components of the Riemann states, respectively, and and are the first two components of the normal flux , which is calculated using the HLLC solver by where , , and represent the speeds of the left, contact, and right waves, respectively, which are determined by (17). The third flux component in the middle flux vector is calculated by or : where , , , and are the components of the left and right initial Riemann states for a local Riemann problem, with the middle states and calculated by Similar formulae are used to calculate , , and .
4.3. High-Order Reconstruction and WSDGM
To achieve higher-order accuracy, which is the most important aspect for any flow solver, better reconstruction techniques than a piecewise constant function should always be used. Overall second-order accuracy in space and time is achieved using a two-step unsplit MUSCL-Hancock method. The predictor and corrector steps are given by (19) and (20), respectively, No Riemann solutions are required at the predictor step and the interface fluxes are essentially evaluated within each cell using interpolated values at the extremities of the inner interfaces [6].
In the corrector step, the fluxes through the cell interfaces are evaluated by the HLLC approximate Riemann solver. Then, the Riemann states can be reconstructed using the WSDGM, taking the eastern interface as an example: where and are the predicted cell-centre values at cell and its eastern neighbor obtained in the predictor step; the slope limiter is minmod limiter function; r is the ratio of gradients given by Using WSDGM, the first component of (21) is reconstructed as where is a weighting parameter which is defined as and , , , and are the DGM and SGM reconstructions, respectively. They are labeled as where represents an estimation of the intercell bed elevation if a linear variation of the bed profile is assumed.
4.4. Friction Terms Treatments
When considering a case with complex topography and involving wetting and drying, the friction term actually dominates the stability of a numerical model. In particular, when the water depth becomes very small, the Manning formulation of friction may lead to an exaggerated force that can even reverse the flow, which is obviously physically incorrect. To deal with this problem, in this work, the friction terms are solved by the semi-implicit scheme which is given by where ; , and are the initial states for the friction source term treatment. This scheme does not invert the direction of velocity and prevent instability effectively because (26) means and . Moreover, when the water depth becomes very small, this method may relieve the problem of high velocity. For example, when , , , and , then the velocities computed by (26) are and . This indicates that the semi-implicit scheme (26) can well benefit the computational stability.
4.5. Time Stepping and Wet/Dry Treatments
The numerical model is explicit, and its stability is controlled by the Courant-Friedrichs-Lewy (CFL) condition. For a two-dimensional Cartesian grid, the CFL criterion is used for predicting an appropriate time step [6]: with and , where is the Courant number and = 0.8 is adopted in this study. It is worth mentioning that the local time stepping method [18] can be implemented to achieve higher computational efficiency.
When considering the application to wetting and drying processes on a complex terrain, there is a simple but effective method to capture wet/dry fronts. After the whole solution is updated by (20), the computed water depths lower than zero are set at zero; besides, the velocities of the dry cell and its neighbors are also set at zero. This simple method which prevents the fake fluxes through the dry cell interfaces could sharpen the wet/dry fronts of the new conserved flow variables .
5. Numerical Results and Discussion
In all the simulations undertaken in this section, and .
5.1. Still Water Test with Wet/Dry Front
This test is chosen to verify that the present model indeed preserves the stationary steady flow with wet/dry fronts over an uneven bottom. The numerical experiment takes place in a square pool . At the centre of the pool, a symmetric bump at the bottom is mathematically defined as the form
The still water as the initial condition covering partially the bump in a frictionless domain is fully closed by vertical walls; the water surface elevation is 0.1 m and 0.3 m. A simulation of s is performed on a grid of 50 × 50 cells. Figure 2 shows the computed 3D water surface in the different conditions. Similarly, Figure 3 presents the water level in cells versus the -coordinate of the cell centroid along the line m. In terms of [19, 20], the errors, where , , and denote the analytical solutions and is the number of cells with positive water depth, are listed in Table 1. The errors are sufficiently small to the level of round-off error of a double accuracy calculation. In addition, , represents the maximum relative error in global mass conservation over the simulation duration under the two conditions, in which and are the initial volume of water and the computed volume of water at the time , respectively, and represents the net volume transported out of the domain from the beginning of the simulation period to the time .
(a)
(b)
(a)
(b)
The results show that the stationary solutions are well maintained, and the current numerical scheme is able to preserve the lake at rest solution involving wet/dry interface and hence is well balanced.
5.2. Moving Shorelines in a Parabolic Bowl with Friction
Thacker [21] provides a good benchmark test case for the present shallow-flow model in dealing with wetting and drying and bed slope as well as the friction source term. This paper firstly deduced the analytical solution of the test case.
Consider the case where the motion of shallow water in a basin is governed by the equations [21] where is the height of the water surface above mean water level, is the bottom surface, is the averaged component of the water current to the east, is the depth-averaged velocity component of the water current to the north, is the bottom friction parameter which is considered to be a constant (related to the bed friction coefficient by ), is the acceleration due to gravity, and is the time.
Assuming that was function of only, , and , then (31) and (32) imply that
The bed profile of the domain is defined by where and are constants so that flow takes place above parabolic bottom topography.
Substituting (33), (34), (35), , and in (30) gives Following Thacker, equate the time-varying coefficients of the linearly independent terms, 1 and , respectively, Substituting (34) in (38), Substituting (34) in (37), The auxiliary equation for (39) is The roots of this equation are
In this test case, the computational domain has planned dimensions of ; the coefficients are , , and . So, considering the case with , a solution of (39) is where is a constant, obtained by using given values for and . Substituting (43) in (40) and integrating with respect to give Substituting (43) in (34) gives Substituting (44) and (45) in (33) gives Thus the analytical water level is given by
We compared the numerical solution of the model in this paper and the analytical solution; the simulation lasted 6000 s. Figure 4 shows a good agreement between the numerical predictions and analytical solutions of the free surface profile at different times; the moving shoreline is accurately simulated with no signs of amplified spurious oscillations and due to the friction effect, the sloshing amplitude of water level is damping with increasing simulation time. Figure 5 shows the comparison between numerical and analytical solutions for water depth at points −2750 m, 50 m, and 2750 m. The numerical predictions are found to match perfectly with those given by the analytical solutions and no obvious distortion is detected at the wet-dry interface. The results validate the wetting and drying method of the present model, as well as the source terms treatments.
(a)
(b)
(c)
(d)
(e)
(f)
5.3. 1D Steady Flows over a Steep Bump with WSDGM
The goal of this test case is to investigate the ability of the present model to converge to a nonstationary steady state. The test concerns 1D steady flows in a , which is discretized with 0.1 m grid spacing to make the test more severe, frictionless channel. The bottom profile, characterized by the presence of a steep bump, is described by the following equation [22]:
In order to demonstrate the capability of the WSDGM reconstruction, the results of supercritical flow, subcritical flow, and transcritical flow with a shock are displayed in Figures 6, 7, and 8. The three cases were simulated using a mesh of 400 () cells to steady states.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
For the case of supercritical flow over the bump, we impose a constant discharge per unit width at the upstream boundary and the downstream boundary set to be transmissive [22]; the range from 2.31 to 3.83. Figure 6, in the subdomain , shows the numerical results of the DGM (, (a)), SGM (, (b)), and WSDGM (, (c)) reconstruction, respectively. Figure 6(b) shows that the numerical scheme with converges toward the steady solution worse than . The results verified that pure SGM reconstruction is unsuitable when the flow is supercritical.
On the other hand, if the flow is subcritical (e.g., at the upstream and at the downstream), then the range from 0.14 to 0.41; Figure 7(a) shows a dip in the water elevation upstream and downstream of the bump, where the bottom profile is not smooth. Moreover, the pure DGM reconstruction does not satisfy the C-property.
In Figures 6(c), 7(c), and 8, excellent agreements are observed between the numerical predictions and analytical solutions under different conditions. All in all, the results prove that WSDGM performs better than the schemes based on pure DGM or SGM reconstructions. For this reason, the parameter was set at 1 also in the following tests.
5.4. 2D Circular Dam-Break Problem
This problem is based on the hypothetical test case studied by Alcrudo and Garcia-Navarro [23], which involves the breaking of a circular dam. In this problem, a circular discontinuity over a 200 m × 200 m is investigated. The circle’s radius is 50 m and it is centered at the middle of the square. The initial states are 10 m water depth inside the circle, and the water is at rest. The initial water depth outside the circle is considered for two cases, that is, 5 m and dry bed. The wall is then assumed to be removed completely and no slope or friction is considered. It is assumed that the circular dam is broken instantaneously and water surface profile, 2 s after dam failure, is studied.
In these examples, the computational grid consists of 100 × 100 cells with each cell 4 m2 as the background grid. During the simulation, the adaptation criteria and are used. Figure 9 shows the computed water surface profile 2 s after dam failure outside water depth equal to 5 m and dry, respectively. Obviously, the circular dam-break waves spread and propagate radically and symmetrically as the water drains from the deepest region. Figure 10 presents the final numerical outputs at 2 s after the dam failed in terms of adapted grid. From the 3D water surface and depth contours, it is conspicuous that the sharp shock front is well captured and the depth contours show that most of the depth contour lines appear to be circular. The case outside the water depth is equal to 5 m, and the flow regime changes from subcritical to supercritical. This case can be a good test for validation of a numerical scheme. In the case of dry-bed simulation, it can be seen that no bore forms, but instead a rarefaction wave extends into the dry region (Figure 11). The scheme is capable of handling the dry-bed problem. Similar results were obtained by Alcrudo and Garcia-Navarro [23] and Zoppou and Roberts [5].
(a)
(b)
(a)
(b)
In order to perform the advantage of the adaptive grid-based model, the simulations are also run on uniform grids with different resolutions. The results are shown in Figure 12 in terms of water surface profile. The structured adaptive grid-based prediction matches closely the one produced on the finest uniform grid with 200 × 200 cells, which indicates similar accuracy of the two simulations. On the adaptive grid system, the mesh is refined to those areas with sharp surface gradients. Therefore, the computational efficiency is greatly improved due to the fact that the total number of cells is reduced. As listed in Table 2, the new adaptive grid can save more than twice of the computational cost, compared with the simulation on a uniform grid with 200 × 200 cells.
5.5. Bore Wave Past a Dip
This test gives unsteady results of a right-going bore wave past a downward hump. The Froude number is , and the initial conditions to the left and the right of the bore are The initial discontinuity is located at . The bottom function is defined by The computational domain is a square . We use inflow boundary conditions at the west boundary and outflow boundary conditions for the other boundaries. Simulation is carried out for 2 s when the bore has passed the dip and travels towards the eastern boundary.
Instead of using a uniform fine mesh to decompose the entire computational domain, adaptive grid on the coarse background grid 160 × 200 with one level local refinement is used. and are prescribed for grid adaptation, as Figure 13(a) shows. The complex wave pattern is captured by the adapted grid with steep surface gradient; the dynamically adaptive grid ranges from 32000 to 44684. The results of the simulation are presented in Figures 13(b) and 13(c) in terms of 3D water surface and depth contours, which are compared very well with those presented by Liang [18] and Tang [24]. Comparing Figure 13(a) with Figure 13(c), it is obvious that the wave fronts can be partly predicted by the adapted grid.
(a)
(b)
(c)
5.6. Dam-Break Flooding over Three Humps in a Closed Channel
This test of a 2D dam-break wave traveling over an initially dry and rough floodplain with three hillocks has been widely accepted as a standard benchmark to assess the adequacy of a numerical model for realistic flood modeling applications. The dam break occurs in a 75 m long and 30 m wide closed container, with the bed topography defined by
The dam is initially placed at m to hold a tranquil water body with a surface elevation of 1.875 m. The rest of the domain is dry and a global Manning coefficient is set to be . In this test, grid enrichment and coarsening are undertaken according to and .
Figure 14(a) illustrates the bed topography within the computational domain and the upstream still water retained by the dam. Figures 14(b), 14(c), 14(d), 14(e), and 14(f) illustrate the 3D view of free surface elevation at s, 6 s, 12 s, 30 s, and 300 s, respectively. As illustrated in Figure 14(c), for s, while the violent dam-break flow continues on moving downstream and interacts with the three hills, a reflected shock is developed after the depression wave hit the western boundary wall and propagating downstream. At s, the flood water is passing either side of the large hill and starts to flood the lee of the hill. The upstream directed reflection shock waves from the small hills have coalesced into a nearly straight bore. The main flood wave has reflected against the large hill. By s, steady state is achieved, with the peaks of the small hills no longer submerged. Similar results were obtained by [6, 17]. Figure 15 shows the evolution of the number of cells in this adaptive grid, which initially increases from 2364 to 6126 and then decreases to a steady value of about 2700 in virtue of the flow becoming steadier, this essentially indicates the efficiency by using dynamically adaptive grids for flood inundation problems.
(a)
(b)
(c)
(d)
(e)
(f)
5.7. Flood Flows over Complex Topography
In order to test the applicability of the present model to real case studies characterized by an irregular topography, a practical application related to river flood over realistic topography is presented. The study area is a river in southern China. A free outgoing condition is imposed along all boundaries and the Manning coefficient is chosen to be 0.019 over the whole domain.
Figure 16(a) illustrates the plan view of the bed elevation together with the locations of the gauges; Figures 16(b) and 16(c) show that the peak-flood discharge is imposed at the inflow boundary and the field-measured water level is specified at the outflow boundary, respectively. As s, the velocity field is set to , whereas the initial water level is set to 17.366 m in the entire computational domain.
(a)
(b)
(c)
Simulations are executed for 300 hours after the flood as mentioned previously. Figure 17 presents the velocity field of regions A, B, and C. The close-up shows that the major stream flow is confined along the main river channel, which indicates that the simulated velocity field is quite reasonable. It is worth noting that this is a sandbar in region C where the bed elevation is much higher than the water surface elevation of its nearby area; thus an obviously detour flow phenomenon is shown in Figure 17(c), validating the wetting and drying method in the present model. Figure 18 shows the comparison of the measured and simulated water flow at three gauges, and the computed results match the measured water flow very well. Besides, the suddenly rising phenomena of water level are shown obviously at the start of the flood tide and then the water level began to drop. All velocities were monitored during the simulation, and we noticed that no unphysical velocities (bigger than 25 m/s) happened.
(a)
(b)
(c)
(a)
(b)
(c)
6. Conclusions
A Godunov-type numerical model based on the finite volume method has been developed and tested to simulate shallow-water flows over a wet or dry bed with complex domain topography. The physical domain has been discretized on uniform and nonuniform structured grids with dynamic refinement. The numerical model solves the two-dimensional shallow-water equations derived as conservation laws with the flux and source terms conveniently balanced so that there is no need for further numerical treatment. The weighted surface-depth gradient method (WSDGM) computes water depth at the cell boundaries through a weighted average, based on the local Froude number of the extrapolated values deriving from DGM and SGM reconstructions. The method is simple, accurate, robust, and easy to implement and can be used to solve both steady and unsteady shallow-water problems. The present results demonstrate the accuracy of the new finite volume method and its capability to simulate water flows in the hydrodynamic regimes considered.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is financially supported by the Key Program of the National Natural Science Foundation of China (no. 51239004), the National Science & Technology Pillar Program during the Twelfth Five-Year Plan Period (Grant no. 2012BAK10B04), the CRSRI Open Research Program (CKWV2014220/KY), and the Open Research Foundation of PRHRI (no. 2013KJ01).