A well-balanced unsplit finite volume model with geometric flexibility

3.1. Finite volume method

The governing equations are discretized based on the cell centered finite volume method. By integrating Eq. (1) over an arbitrary computational cell, the basic equation of finite volume method can be obtained as:

where $\mathrm{\Omega }$ is the boundary of the cell, $\mathbf{n}$ is the outward unit vector that is normal to the boundary $\mathrm{\Omega }$, and $\mathbf{E}·\mathbf{n}$ is the normal flux vector and can be determined by using the rotational invariance property [11]:

where $\theta$ is the angle between the outward unit vector $\mathbf{n}$ and the $x$-axis, $\widehat{\mathbf{U}}=\mathbf{T}\left(\mathbf{U}\right)$ represent the transformed conserved variables, $\mathbf{E}\left(\widehat{\mathbf{U}}\right)=\mathbf{E}\left(\mathbf{T}\left(\mathbf{U}\right)\right)$ is the transformed normal flux, and $\mathbf{T}\left(\theta \right)$ and ${\mathbf{T}\left(\theta \right)}^{-1}$ are the transformation matrix and its reverse, respectively:

therefore, Eq. (3) can be discretized as [12]:

where $A$ is the area of an integration cell, $M$ is the number of cell interfaces, $m$ is the index of the cell interface, and $L^m$ is the width of each interface.

3.2. HLLC approximate Riemann solver

The HLL approximate Riemann solver [13] cannot take into account intermediate waves, such as shear waves or contact discontinuities [11]. To correct the defects of the HLL approach, Toro et al. [14] proposed a modified method called the HLLC scheme, in which C stands for Contact. The HLLC scheme provides results that satisfy the entropy condition [15]. If wave speeds are properly calculated and applied, the HLLC scheme can easily handle dry beds as well [16]. It also provides successful approximate solutions to practical applications [17].

Fig. 1Wave structure of the HLLC approximate Riemann problem

Fig. 1 shows the structure of waves assumed in the HLLC approximate Riemann solver. The conserved variable $\mathbf{U}$ is divided into four parts according to wave speeds$\mathrm{}{\mathrm{S}}_L$, ${\mathrm{S}}_R$ and ${\mathrm{S}}_{\mathrm{*}}$. If the conserved variables in each region have been determined, the HLLC numerical fluxes are evaluated from [11]:

where ${\mathbf{F}}_L=\mathbf{F}\left({\mathbf{U}}_L\right)$, ${\mathbf{F}}_R=\mathbf{F}\left({\mathbf{U}}_R\right)$, and ${\mathbf{U}}_L$ and ${\mathbf{U}}_R$ are the left and right Riemann states at each interface, respectively. ${\mathbf{F}}_{\mathrm{*}L}$ and ${\mathbf{F}}_{\mathrm{*}R}$ are the numerical fluxes in the left and right parts of the star region of the Riemann solution separated by a contact wave, respectively, and $S_L$, $S_{\mathrm{*}}$ and $S_R$ are the speeds of the left, contact and right waves, respectively.

To properly apply the HLLC approach to shallow water equations, it is important to calculate correct wave speeds $S_L$, $S_{\mathrm{*}}$ and $S_R$ [18]. In this study, a computational cell and a neighboring cell are set as the left and right hand sides of the Riemann interface for correctly calculating wave speeds using the unsplit method (Fig. 2(b)). At the same time, the wave speeds $S_L$, $S_{\mathrm{*}}$ and $S_R$ are estimated by taking into account the wet and dry bed conditions. The combinations of existing relations among wave speeds $S_L$ and $S_R$ are shown in Table 1.

In Table 1, $u_L$, $u_R$, $h_L$ and $h_R$ are the left and right initial values for a local Riemann problem. The contact wave speed $S_{*}$, water depth $h_{*}$ and $q_k$ ($k=L,\mathrm{}R$) are calculated from [11]:

3.3. Linear reconstruction of unsplit finite volume scheme

The MUSCL-Hancock scheme is a second-order extension of the Godunov upwind method. This scheme reconstructs conserved variables at the interface between cells prior to the calculation of fluxes for second-order accuracy in space and predicts fluxes at $\mathrm{}∆t/2$ for second-order accuracy in time. To achieve second-order accuracy in space and time, the two-step MUSCL-Hancock method consisting of prediction and correction procedures is employed.

In the prediction step, the reconstruction of conserved variables is implemented at the interface between cells by using the data of both computational cells and neighboring cells. The simultaneous reconstruction of conserved variables at each interface is required to predict fluxes using the unsplit method. Thus, information about a computational cell and its neighboring cells, including cell numbers and primitive variables, must be known, and the order of the neighboring cells needs to be determined for two-dimensional calculations. To obtain the information necessary for calculation, the following algorithm is introduced in this study:

1) The angles produced by the centroid coordinate of a computational cell and the coordinates of nodes consisting of a computational cell are calculated using Eq. (11) for all cells in the computational domain:

2) As shown in Fig. 2(a), the index of nodes is determined in ascending order based on the angles calculated in the previous step.

3) Two nodes from the lowest number are read counter-clock wise for a computational cell. At the same time, the node numbers for the other cells, except the computational cell, are read in the same manner. The process is repeated to find the nodes matched with the two previously read nodes.

4) The cells with two shared nodes become the computational cell and one of its neighboring cells. The order of neighboring cells is determined in the order of angles produced between the centroid and the nodes of the computational cell.

The proposed algorithm can be easily applied without any limits to the configuration of cells, even in complex flow domains such as natural rivers.

Fig. 2Schematic diagram of topological relationships between cells a) The order of neighboring cells and b) Linear reconstruction at interface for unsplit scheme

a)

b)

The SGM [3] is coupled with the MUSCL-Hancock scheme so as to keep the balance between flux gradients and source terms. The SGM uses the water level $\eta$ instead of the water depth $h$ of depth gradient method (DGM) for reconstruction of conserved variables at interfaces in the prediction procedure. Accordingly, $\delta {\mathbf{U}}^{\mathrm{\text{'}}}$, which is the variation of ${{\mathbf{U}}^{\mathrm{\text{'}}}=[\mathrm{}\eta ,\mathrm{}\mathrm{}hu,\mathrm{}\mathrm{}hv]}^T$, instead of the variation of ${\mathbf{U}=[\mathrm{}h,\mathrm{}\mathrm{}hu,\mathrm{}\mathrm{}hv]}^T,$ is estimated as:

where $l$ is a computational cell and $k=$1, 2, 3 and 4 represent the index of neighboring cells composing cell $l$ (Fig. 2(a)). To prevent numerical oscillations that may occur near discontinuities with the use of the second-order-accurate MUSCL-Hancock method, the TVD scheme, which limits the slopes of conserved variables to be reconstructed at the interface, is applied in this study. The Superbee limiter [11] is extended to the $x$ and $y$ directions for unsplit scheme application as Eq. (13):

For linear reconstruction of flow variables at interface, the method which uses the distance between the center of a computational cell and midpoint of edges is employed (Fig. 2(b)). This technique is efficient to treat uniform or nor-uniform mesh since it can capture the correct behavior of linear function within a cell [13]. Therefore, ${{\mathbf{U}}^{\mathrm{\text{'}}}=[\mathrm{}\eta ,\mathrm{}\mathrm{}hu,\mathrm{}\mathrm{}hv]}^T$ for each interface between cells is extrapolated as:

where $L_{ck}$ ($k=$1, 2, 3 and 4) is the distance from the center of cell $l$ to each midpoint of interface (Fig. 2(b)). To estimate the conserved variable ${\mathbf{U}=[\mathrm{}h,\mathrm{}\mathrm{}hu,\mathrm{}\mathrm{}hv]}^T$ needed to calculate the fluxes at the interfaces, $hu$ and $hv$, both of which are conserved variables of momentum equation calculated in Eq. (14), and $h$, which is the conserved variable of the continuity equation calculated in Eq. (15), are used [13]:

where $h_{l,k}$, ${\eta }_{l,k}$ and $z_{l,k}$ refer to the water depth, water level and average bed elevation at interfaces between a computational cell ($l$) and neighboring cells ($k$). Audusse et al. [20] suggested $h_{l,k}=\text{max}\left(0,{\eta }_{l,k}-z_{l,k}\right)$ instead of Eq. (15) to preserve nonnegative water depth. However, this suggestion can cause imbalance between fluxes and sources terms in dry bed [21]. In this study, a simple technique propose to preserve nonnegative water depth and the well-balanced solutions in a dry bed application as Eq. (16):

where $k\text{'}$ is the counterpart number of $k$ (e.g., $k=$1, $k^{\text{'}}=$3).

The prediction fluxes ${\mathbf{F}}^p$ and ${\mathbf{G}}^p$ are obtained by substituting the conserved variables at interfaces into Eq. (2). For second-order accuracy in time, a conserved variable in $\mathrm{\Delta }t/2$ is computed as:

where $A_l$ is the area of the cell $l$, $\mathrm{\Delta }t$ is the computational time step, $N_l$ is the number of the neighboring cells of cell $\mathrm{}l$, ${\mathbf{n}}_{l,k}^x$ and ${\mathbf{n}}_{l,k}^y$ are the $x$ and $y$ components of the normal outward unit vector of the interface ($l,k)$ and $w_{l,k}$ is the width of the interface.

In the correction step, a conserved value $\mathbf{U}$ at the time level $n+1$ is achieved by solving a series of Riemann problems at the interface based on data from the prediction step [3, 19]:

where the correction fluxes ${\mathbf{F}}^c$ and ${\mathbf{G}}^c$ in the $x$ and $y$ directions are computed using the HLLC approximate Riemann solver.

3.4. Treatment of source terms

In developing a numerical model to compute shallow water flows with an irregular bed topography, the numerical treatment of source terms largely affects the accuracy of the calculated results [17, 22-23]. The SGM is accurate and effective in treating source terms and can be easily coupled with the MUSCL-Hancock scheme [3]. Accordingly, this study improved the accuracy and stability of the developed model by combining the two schemes. The source terms consist of two parts, including the bed slope and the friction slope.

For the discretization of the bed slope terms, the application of the divergence theorem yields the rate of change of bottom elevation $z$ in the $x$ and $y$ directions of cell $l$ as follows [22]:

where subscripts 1, 2, 3 and 4 represent the index of nodes composing the computational cell, as shown in Fig. 2(b). The friction slopes in the $x$ and $y$ directions are estimated using Manning’s formula with the wide channel assumption [24]:

where $n$ is Manning coefficient, and the primitive variables $h$, $u$ and $v$ are estimated using the average values over the computational cell $l$.

3.5. Boundary conditions

Fig. 3Assignment of ghost cell numbers for boundary condition treatment

The finite volume method requires ghost cells for the treatment of the boundary condition. Fig. 3 depicts the computational domain and ghost cells. In this study, the numbers to be assigned to each ghost cell are determined using boundary cell numbers as shown in Fig. 3. This technique delivers an efficient treatment of complex topographies at the boundary. For general shallow water problems, two types of boundary conditions including the transmissive boundary condition and the reflective boundary condition are taken into account without considering the flow regime [11]. The Eq. (21) and (22) are used for a transmissive, and a reflective boundary condition, respectively:

where the subscripts $L$ and $\mathrm{*}$ refer to the left hand side and right hand side of the Riemann states at the interface between boundary cells and ghost cells, respectively.

4.1. 1D Steady flow over a bump

If boundary conditions remain constant in a computational domain in which stationary flow is taken as the initial state, the flow will become steady-state despite the presence of a bed slope. However, nonphysical oscillations occur around the discontinuity of the bed slope if the source terms are not treated appropriately. Simulations investigating the convergence to a steady state are generally used to verify the balance between the flux gradients and source terms of a numerical model. In order to demonstrate convergence to a steady state, this study takes into account a channel proposed by Goutal and Manual [25]. The channel is frictionless and rectangular, 25 m in length and 1 m in width. The bed elevation is defined as:

According to the initial and boundary conditions, the steady flow over a bump can be subcritical or transcritical flow with or without a shock [25]. Of the three cases, subcritical flow and steady transcritical flow with a shock are selected for model verification, because the two cases are considered most challenging for a numerical model. We refer to Goutal and Maurel [25] for analytical solutions.

4.1.2. Transcritical flow with a shock

For the transcritical flow with a shock, a discharge per unit width of $q=$0.18 m²/s is imposed at the upstream boundary and a water level of $\eta =$0.33 m is fixed at the downstream boundary. The initial water level is $\eta =$0.33 m throughout the entire computational domain, in which case the hydraulic jump occurs at the end of the bump.

Fig. 4Steady subcritical flow over a bump a) Water level and b) Discharge per unit width

a)

b)

To investigate accuracy dependence on cell size, simulations for 100 and 200 cells were carried out, respectively. Fig. 5 shows a comparison between the numerical predictions and analytical solutions for the water level and discharge per unit width estimated after 360 seconds. In both cases, numerical errors occur at the end of the bump. When the number of cells is 100, discharge near the point where the hydraulic jump occurs cannot meet steady-state conditions. On the other hand, when the number of cells is increased to 200, the numerical prediction is shown to have good agreement with analytical solutions. The maximum deviation between the discharge predictions and the analytical solutions is 15 % and 2 % for 100 and 200 cells, respectively. It is clearly shown that the error between the numerical solution and analytical solutions decreases as the number of cells increases.

Fig. 5Steady transcritical flow over a bump with a shock a) Water level and b) Discharge per unit width

a)

b)

The simulated results of the steady flow over a bump show close agreement with analytical solutions without spurious oscillations. This suggests that the unsplit method applied to the proposed algorithm in this study and the source-term treatment technique combined with the SGM maintains a numerical balance between fluxes and source terms.

In order to assess the efficiency and accuracy of the source term treatment with the SGM in the developed model, a case with the DGM is employed for comparison. The SGM is the same as the DGM in procedure, except that it uses the water level instead of the water depth for the reconstruction of conserved variables at the interfaces between cells. The SGM and DGM are applied to the previous two test cases of steady flow over a bump, respectively. In those applications, 200 cells are used. Fig. 6 and 7 show a comparison of the unit discharge predictions of the SGM and DGM to subcritical flow and transcritical flow with a shock, respectively. In the comparison results, the numerical treatment of source terms using the DGM resulted in far larger numerical errors than that of source terms using SGM. This occurs because the DGM does not reflect the real variation of water depth [3]. As a result, the DGM cannot accurately reconstruct water depth at the interface. Thus, the incorrect reconstruction of water depth affected the entire calculation by causing an inaccurate estimation of fluxes at the interface. This demonstrates that the application of the SGM for the treatment of source terms can improve the accuracy and efficiency of the developed model.

Fig. 6Comparison of the DGM and SGM using the proposed model for steady subcritical flow over a bump a) DGM and b) SGM

a)

b)

Fig. 7Comparison of the DGM and SGM using the proposed mode for steady transcritical flow over a bump with a shock a) DGM and b) SGM

a)

b)

4.2. Dam-break flow in an adverse slope channel

In this section, the accuracy of the proposed algorithm is verified through simulating dam-break experiments designed by Aureli et al. [28]. Laboratory experiments set up with different upstream depths and adverse slopes were intended to induce shock origination and propagation, a reverse flow and wetting and drying conditions. As shown in Fig. 8, the experimental domain was composed of a rectangular reservoir and a channel with an adverse slope connected to it. The length and width of this channel were set at 7 m and 1 m, respectively. To apply the numerical model, a computational domain consisting of 140 cells with $∆x=$0.05 m and $∆y=$1.0 m was composed. The Manning roughness coefficient is set to 0.010 m^-1/3s, and the reservoir has an initial water depth of 0.25 m. The initial water depth downstream of the dam is zero for the dry bed test, and the total simulation time is set to 40 seconds. The side walls and upstream region of the channel have closed boundary conditions, and the downstream end has transmissive boundary conditions. Under these same conditions, Aureli et al. [28] and Tseng [29] validated a numerical model using the TVD-Maccormack scheme and the flux-limited TVD scheme, respectively. Fig. 8 shows the experimental configurations and simulation conditions. The bed topography is defined as:

Fig. 8Schematic view of the dam-break experiment in an adverse slope channel

Fig. 9 shows a comparison of the calculated water depth and measured water depth at four different points. The simulated results for predicting the propagation time of the forward and reverse waves, wet/dry fronts and water depth are in close agreement with the observed values. From these simulated results, it can be demonstrated that a new method for a linear reconstruction of conserved variables in effectively applying an unsplit scheme and a simple technique to preserve balance between fluxed and source terms in a dry bed can handle reverse flow and wetting and drying beds very well.

Fig. 9Dam-break flow in an adverse slope channel: comparison of water depths at a) G1, b) G2, c) G3 and d) G4

a)

b)

c)

d)

4.3. Malpasset dam-break flow

To verify the applicability of the developed model to a natural river with an irregular bed elevation and dry terrain, it is applied to the Malpasset dam-break simulation. The dam was located in the valley of the Reyran River about 12 km upstream of the Frejus city in France and was constructed to provide irrigation and drinking water for that city in 1954. It was an arch dam with a maximum height of 66.5 m, a crest length of 223 m and a storage capacity of 55×10⁶ m³. It collapsed due to a sudden rise of water level in the reservoir resulting from exceptional rainfall on December 2, 1959, and the resulting flood waves were propagated along the Reyran valley to Frejus, leaving 421 casualties [30]. This dam-break case was simulated by other researchers [30-33] to verify their numerical model.

The high water marks of flood waves resulting from the dam-break were surveyed by the police along the left and right banks of the Reyran River. Fig. 10 shows the surveyed points, including P1-P17. In 1964, a laboratory experiment on the Malpasset dam-break was conducted using a physical model scaled at 1:400, made by LNH-EDF. In the experiment, the arrival time of the flood wave and the maximum free surface elevations were measured in S6-S14 gauges as shown in Fig. 10 [32]. These data are used to validate the numerical model in this study. The boundary condition for the upstream region is an imposed reflective boundary condition without any inflow into the reservoir because it was very small in comparison to the discharge resulting from the dam-break [30]. The downstream boundary conditions are set as a free flow which flood waves run into the sea. For the initial conditions, the water level of the reservoir is imposed at 100 m, and the downstream region of the dam is considered as a dry bed. The 13,550 quadrilateral meshes composing the computational domain are used and Manning coefficient is specified at 0.033 m^-1/3s for the entire area.

In Fig. 11, the predicted arrival time of the flood front and the maximum water surface elevation are compared with the measurements of the physical model. Fig. 12 shows a comparison of the calculated maximum water surface elevation with field data surveyed in the left and right banks. The predictions in the proposed model demonstrate relatively good agreement with the field and experimental data. However, there are small discrepancies between the simulation results and high water marks at surveyed points, P4, P5, P7, P8, P13 and P16. It should be noted that, in the simulation, bottom elevation data digitized from maps dating back to 1931 are used [30]. However, the field data were collected after the dam-break, which was accompanied by abrupt changes in topography. This may be the reason for differences between the numerical results and the high water marks.

Fig. 10Topography and locations of field surveyed points (P1-P17) and gauges at physical model (S6-S14)

Fig. 11Comparison between predicted results and experimental data a) Propagation time at gauges, b) Maximum water level at gauges

a)

b)

Fig. 12Comparison between predicted results and field data a) Maximum water level at right banks, b) Maximum water level at left banks

a)

b)

Fig. 13 illustrates the three-dimensional flood inundation extent and water depth distributions at 0, 3, 12, 21, 30 and 35 minutes after the dam-break. As shown in Fig. 13, the main channel in the center of the river valley is included in the topography and flood waves due to dam-failure are rapidly propagated along the steep and narrow valley, starting at the dam.

Fig. 133D inundation water depth and extent at a) 0, b) 3, c) 12, d) 21, e) 30 and f) 35 min

a)

b)

c)

d)

e)

f)

However, as the flood waves travel downstream, where the floodplain becomes wider, they are propagated more slowly and the water depth decrease. Additionally, the wave front is also propagated over a wider area. The water level, which is initially set to 100 m in the reservoir, dramatically drops for 35 minutes after the dam-break. Fig. 13 shows that the main characteristics of the flood resulting from the Malpasset dam-break can be well reproduced by the proposed model. It also reveals that the techniques presented in this study can be simply and efficiently applied to irregular topographies such as natural channels.