Simulation of Dam Breaks on Dry Bed Using Finite Volume Roe-TVD Method

Abstract

Dams are one of the most important hydraulic structures. In view of unrecoverable damages occurring after a dam failure, analyzing a dams’ break is necessary. In this study, a dam located in Iran is considered. According to adjacent tourist and entertainment zones, the breaking of the dam could lead to severe problems for the area and bridges downstream of the river. To investigate the issue, a numerical FORTRAN code based on the 2D finite volume Roe-TVD method on a fixed bed is provided to assess the effects of the dam break. Turbulence terms and dry bed conditions were considered in the code. A numerical wave tank (NWT) with a triangular barrier in the bed was numerically modeled and compared with analytical models to verify the capability of the code. Comparing numerical, experimental and analytical results showed that estimated water level and mass conservation in the numerical model is in good agreement with the experimental data and analytical solutions. The 2D approach used has reduced the cost of computing compared to a 3D approach while obtaining accurate results. The code is finally applied to a full-scale dam-break flood. Six KM of the natural river downstream of the dam, including two bridges, B1 and B2, is considered. Flood flow hydrographs and water level variations at bridges B1 and B2 are presented. The results denoted that bridges B1 and B2 will be flooded after 12 and 21 min, respectively, and are at risk of the potential break. Thus, it is necessary to announce and possibly evacuate the resort area alongside the dam in order to decrease losses.

1. Introduction

In recent years, the statistics of the dam break phenomenon have been collected to improve the knowledge for the construction dams. The primary goals of the studies are considering and implementing the new findings in designs to have safe dams adjacent to urban areas and keeping the projects time and cost-efficient. On the other hand, the consequences of a dam break need to be considered in the designing stage in order to prevent probable damage to nearby infrastructure and resort areas.

Modeling and analyzing dam breaks have been carried out by several researchers. In some one-dimensional models, the non-conservative form of the shallow water equations has been solved [1,2]. Zoppou et al., in 2000 [3], described a two-dimensional numerical model for dam-break problems. The model was able to consider shocks, complex geometries, including steep bed slopes, and was capable of simulating the wetting and drying process. In this research, the performance of 20 numerical schemes used to solve the SWE for simulating the dam-break problem was examined, and some form of flux or slope limiter was used to eliminate oscillations.

Bradford et al., in 2002, developed a model based on the finite-volume method for 2D unsteady, shallow-water flow over arbitrary topography with moving lateral boundaries. They also introduced a new technique to prevent numerical truncation errors caused by the pressure and bed slope [4]. In the same year, Valiani et al. validated their code with the Malpasset dam-break [5]. In 2003 also, Ying et al. [6] developed a numerical model to simulate flood inundation due to a dam break. Sanders et al., in 2008, by applying Reynolds transport theorem to a finite control volume, derived shallow-water equations appropriate for urban flood modeling [7]. Ni et al., in 2018, presented an approximate solution to dam breaks [8]. The proposed method only worked on uniform slope channels and cannot be applied to real dam-break problems.

Fent et al. [9] investigated the vertical velocity distribution of a dam-break wave using experimental equipment. Using the OpenFOAM, a new code for modeling the two-phase flow in the dam break problem that was created by Park et al. [10], the accuracy of the new method was confirmed by comparing the results with laboratory data. The solution domain was wet in the proposed model. The effect of an obstacle in the river on the hydraulic characteristics of the dam-break flow was investigated by Issakhov et al. [11].

Different forms of obstacles were investigated numerically and experimentally. Different numerical methods for dam-break flow modeling were reviewed by Yang et al. in 2017 [12]. This study showed that Navier–Stokes equations with turbulence modeling have good accuracy. The finite volume method has been applied to a wide range of numerical methods [13]. In this method, volume integrals in a partial differential equation are converted into surface integrals, using the divergence theorem. In this regard, determining the passing flux through the surface is the primary goal of this method. However, estimation of the flux on the dry surface is more complicated. Haltas et al., for the failure problem, reviewed the calculation of the flux at the dry boundary [14]. An individual system is conducted to specify the dry, semi-dry and wet elements in any models. The wet cells remain in the simulation, and dry cells remove from the calculated boundary. The wet-dry method was applied in a three-dimensional finite-difference model by Medeiros and Hagen in 2013 [15]. In their model, different conditions have been used to ensure that the cell is wet and remains in the computational interval during each time step. The drought control was performed based on a length scale, defining the bed roughness. The points that must be removed from the calculated domain are specified after estimating the water depth in each cell side and cell center. In order to reconsider the cells into the calculated domain, the length scale is used, and it is compared with the water depth of the surrounding points. If it guarantees the conditions of the flow, then the cell returns to the next calculations time step. Casulli and Walters combined a three-dimensional difference-finite volume model with a dry-wet method [16]. In their method, the water depth was estimated, and the distance of vertical points was also updated based on each time step. By measuring the water depth as zero, the surface velocity and height were set to zero. Hou et al., in 2013 [17], experimentally investigated dam-break flow on a dry bed. The effect of topography was investigated in this study. Ji et al. used the Environmental Fluid Dynamics Code (EFDC) [18] and applied the dry-wet method presented by Hamrick and William to analyze flood and drought behaviors [19]. In this method, determining dry cells are the key, and if the cells are specified as dry, then the transmitted flux would be zero. In Brufau’s finite volume model, the dry side is similar to the boundary, and water flow is controlled by the depth [20,21]. In this light, to preserve the conservation of mass, bed elevation deference is redefined locally. The finite volume Roe-TVD scheme was applied to the 2D shallow water equation by Cea [22]. He considers three depth average turbulence models by considering the behavior of wet-dry fronts for river flow. Song et al. used the Godunov-type finite volume model in order to simulate two-dimensional dam-break floods over complex topography with wetting and drying [23]. Vichiantong et al., in 2019, used a well-balanced finite volume method for flood simulation [24]. In their work, a well-balanced scheme with bottom slope approximation was developed. The accurate simulation of dam-break problems was investigated numerically by Antunes et al. in 2019 [25]. They presented an efficient technique for switching from Serre [26] to Saint-Venant in the breaking zone. The majority of these models are only applicable to dam-break flows over fixed beds [27]. Spinewine and Zech [28] used two-layer models to simulate dam-break flows over mobile beds. These models were applied to the morphological changes caused predominately by the non-equilibrium transporting of bedload. The impact of a dam-break wave on an erodible embankment with a steep slope was studied by Di Cristo et al. [29]. The simulation was carried out using a two-phase depth-integrated model. Recently, smoothed particle hydrodynamics (SPH) modeling represents a valuable numerical method of CFD particularly suitable for the simulation of dam-break floods [30]. The internal boundary condition was used to apply the bridge piers boundary condition of 2D shallow water equations by Dazzi et al. [31]. They confirmed the capability of this boundary condition in different flow regimes. Different methods of considering the bridge piers on 2D shallow water equations were investigated by Ratia et al. [32]. They applied the headloss method (HL) and solid wall boundary (MD) to consider the effects of the bridge piers, and it was concluded that the MD approach achieves results close to reality even when regimes other than free surface flow are involved.

When a dam breaks, it can cause casualties. These casualties are dependent on the extent of the inundation area, the population in the flood zone and the amount of alert time available. An investigation of the history of dam failures in the world shows that a large number of dams have broken and caused great financial and human losses. For example, the failure of dams Vajont (1963) in Italy, Johnstown (1889) in the USA and Macho (1979) in India caused 2600, 2200 and 2000 deaths, respectively [33]. According to Coasta [34], the death toll when a dam breaks without a flood warning or weak warning is 19 times higher than if there is a proper flood warning system in the area. Due to the critical situation at the time of dam failure, field data of the flood flow characteristics, such as flow depth, flow velocity, etc., are very limited. Therefore, the use of numerical modeling can be useful to determine the flood zone or the duration of response to the flood wave due to dam failure [33]. In this regard, a dam located in East Asia is considered to investigate the consequences after breaching. Breaking of the dam can initiate a detrimental risk to downstream structures, highways, and infrastructure and flow could cause the failure of two bridges, B1 and B2, and also an important highway. In this paper, the 2D finite volume method (Roe-TVD) is used on a fixed bed in order to model the dam break and assess its potential risks. Therefore, a numerical code is developed to model the dry bed. The maximum time that is necessary to apply precautionary behaviors to prevent potential risks is also presented. Due to the limited availability of real-scale flood data caused by dam break, code validation has only been performed with laboratory-scale data [35].

2. Shallow Water Equations

Shallow water equations are used when the length of the flume is much longer than depth, i.e., rivers. These equations can be obtained by integrating the 3D Navier–Stokes equations over the flow depth, considering the incompressibility of fluid and the hydrostatic pressure distribution. The equations are applied to study a wide range of physical phenomena such as dam break, flow in an open channel, flood waves, forces on offshore and nearshore structures [36,37] and pollution transfer. The two-dimensional form is [22]:

$$\frac{\partial W}{\partial t}+\frac{\partial F_x}{\partial x}+\frac{\partial F_y}{\partial y}={\displaystyle \sum _{k=1}^3{G}_k},$$

(1)

where

$$W=\left(\begin{array}{c}h\\ q_x\\ q_y\end{array}\right);\hspace{1em}{F}_y=\left(\begin{array}{c}{q}_y\\ \frac{q_x{q}_y}{h}\\ \frac{q_y^2}{h}+\frac{g{h}^2}{2}\end{array}\right)F_x=\left(\begin{array}{c}{q}_x\\ \frac{q_x^2}{h}+\frac{g{h}^2}{2}\\ \frac{q_x{q}_y}{h}\end{array}\right);$$

(2)

where W is the vector of the conserved variables, including the water depth h, and the unit discharges in each direction are q_x and q_y as well as the vectors F_x and F_y account for the convective fluxes in the x and y directions, and g is the acceleration due to gravity, respectively. The vector G_k is a source term composed of the bed slope G₁, bed friction G₂, and turbulence terms G₃:

$$G_1=\left(\begin{array}{c}0\\ -gh\frac{\partial Z_b}{\partial x}\\ -gh\frac{\partial Z_b}{\partial y}\end{array}\right); G_2=\left(\begin{array}{c}0\\ -\frac{{\tau }_{b,x}}{\rho }\\ -\frac{{\tau }_{b,y}}{\rho }\end{array}\right); G_3=\left(\begin{array}{c}0\\ -\frac{\partial h \stackrel{}{\overline{u_i^{\prime }{u}_j^{\prime }}}}{\partial x}\\ -\frac{\partial h\stackrel{}{\overline{u_i^{\prime }{u}_j^{\prime }}}}{\partial y}\end{array}\right)\hspace{1em}\left(i,j=1,2\right)$$

(3)

where Z_b is the bed elevation, τ_b,_x and τ_b,y are the bed shear stresses due to the friction in the x and y directions, ρ is the fluid density and $\overline{u_i^{\prime }{u}_j^{\prime }}$ is Reynolds stress.

3. Numerical Solution for Equations by Finite Volume Method

Different methods have been used for generating a mesh in finite volume methods. In the unstructured method, triangular cells are used, and the center of the triangular cell is the basis of numerical computation. Bermudez et al., in 1998, introduced unstructured quadrilateral cells [38]. In this method, first, the domain is divided into some triangles (Figure 1a). Then the middle of each side of the triangle is considered as the center of a quadrilateral cell. The two corners of this latter cell are the start and the end of the triangular cell side, and the other two corners are the triangular cell centers around the line (Figure 1b).

A numerical code in the FORTRAN by using the finite volume Roe-TVD method is written to solve shallow water equations. A second-order method is applied by discretizing equations in time and using a semi-step. A multi-dimensional slope limiter is proposed to achieve second-order accuracy. Two depth-averaged turbulence models, including the $k-\epsilon$ model and algebraic stress model (ASM), are used to calculate turbulence terms (G₃). The flowchart of the code is presented in Figure 2.

By time discretization of the system (1) and simplification, the following equations are obtained with a second-order of accuracy in time, [39]:

$$\begin{gathered} W^{n+\frac{1}{2}}=W^n-\frac{\Delta t}{2}\left(\frac{\partial F_x}{\partial x}\left(W^n\right)+\frac{\partial F_y}{\partial y}\left(W^n\right)\right)+\frac{\Delta t}{2}{\displaystyle \sum _{k=1}^3{G}_k^n} \\ {W}^{n+1}=W^n-\Delta t\hspace{0.17em}\left(\frac{\partial F_x}{\partial x}\left(W^{n+\frac{1}{2}}\right)+\frac{\partial F_y}{\partial y}\left(W^{n+\frac{1}{2}}\right)\right)+\Delta t{\displaystyle \sum _{k=1}^3{G}_k^{n+\frac{1}{2}}} \end{gathered}$$

(4)

where $W^n$ is the vector of conserved variables at time t_n, and Δt is the time step. For spatial discretization, an upwind model may be implemented. In the upwind method, the numerical flux is defined as:

$$\begin{gathered} {\varphi }_{ij}=\frac{Z(W_i,n_{ij})+Z(W_j,n_{ij})}{2}-\frac{1}{2}\left|{\varphi }_{}\left(W_i,W_j,n_{ij}\right) \right|\hspace{0.17em}\left(W_j-W_i\right) \\ \left|\varphi \right|=X\left|D\right|X^{-1}\hspace{1em};\hspace{1em}Z=F_x{n}_x+F_y{n}_y \\ \left|D\right|=\left[\begin{array}{ccc}\left|{\overline{\lambda }}_1\right|& 0& 0\\ 0& \left|{\overline{\lambda }}_2\right|& 0\\ 0& 0& \left|{\overline{\lambda }}_3\right|\end{array}\right]\hspace{1em}; \\ X=\left[\begin{array}{ccc}0& 1& 1\\ -{\tilde{n}}_y& {\overline{U}}_x+\overline{C}\hspace{0.17em}{\tilde{n}}_x& {\overline{U}}_x-\overline{C}\hspace{0.17em}{\tilde{n}}_x\\ {\tilde{n}}_x& {\overline{U}}_y+\overline{C}\hspace{0.17em}{\tilde{n}}_y& {\overline{U}}_y-\overline{C}\hspace{0.17em}{\tilde{n}}_y\end{array}\right] \\ {\overline{\lambda }}_1=n_x\hspace{0.17em}{\overline{U}}_x+n_y\hspace{0.17em}{\overline{U}}_y\hspace{1em};\hspace{1em}{\overline{\lambda }}_2={\overline{\lambda }}_1+\overline{C}\hspace{0.17em}{L}_{ij}\hspace{1em}; \\ {\overline{\lambda }}_3={\overline{\lambda }}_1-\overline{C}{L}_{ij} \end{gathered}$$

(5)

The averaged values in the first-order scheme of Roe at each cell are defined as:

$${\overline{U}}_x=\frac{\sqrt{h_i}{U}_{x,}{}_i+\sqrt{h_j}{U}_{x,}{}_j}{\sqrt{h_i}+\sqrt{h_j}}{\overline{U}}_y=\frac{\sqrt{h_i}{U}_{y,}{}_i+\sqrt{h_j}{U}_{y,}{}_j}{\sqrt{h_i}+\sqrt{h_j}}\overline{C}=\sqrt{g\frac{h_i+h_j}{2}}$$

(6)

where U_x,i and U_y,i are the velocities of the flow in cell i in the x- and y- directions, respectively, ${\phi }_{ij}$ is the numerical flux at the cell face ij, hi is the water depth in cell i, L_i is the boundary of the cell i, and $\tilde{n}=\left({\tilde{n}}_{\mathrm{x}},\hspace{0.17em}{\tilde{n}}_{\mathrm{y}}\right)$ is the unit vector normal to the cell face, Figure 1a. To achieve a second order accuracy in the method of Roe, the conserved variables at the triangular cell faces are reconstructed using a spatial limiting technique [38], Figure 3b:

$$W_{I j}=W_i+r{\nabla }_i^{1}W$$

(7)

in which $W_{I j}$ is the value of W_i at the boundary with cell j, r is the distance vector between the cell area center i and the middle of $L_{ij}$, and ${\nabla }_i^1$ is the limited gradient of variables at cell i, defined by:

$${\nabla }_i^{1}W=w_a\nabla W_a+w_b\nabla W_b+w_c\nabla W_c$$

(8)

where $w_a$, $w_b$ and $w_c$ are weighting factors, and $\nabla W_a$, $\nabla W_b$ and $\nabla W_c$ are unlimited gradients of the three surrounding cells a, b and c. The unlimited gradient for cell i is computed using the area-weighted average gradients at the three faces:

$$\nabla W_i=\frac{A_{i1a2}{\left(\nabla W\right)}_1+A_{i2b3}{\left(\nabla W\right)}_2+A_{i3c1}{\left(\nabla W\right)}_3}{A_{i1a2}+A_{i2b3}+A_{i3c1}}$$

(9)

where $A_{i1a2}$ is the area of quadrilateral $i1a2$ (Figure 3a), and ${\left(\nabla W\right)}_m$ is the gradient of the variable W at the face m of cell i. This gradient may be computed from the divergence theorem and an area-weighted average of two triangles around each face. For example, for face 1 in Figure 1a, it can be written:

$$\begin{gathered} \left(\nabla W\right){ }_1=\frac{A_{1a2}\hspace{0.17em}\nabla W\left|{}_{1a2}\right.+A_{1i2}\hspace{0.17em}\nabla W\left|{}_{1i2}\right.}{A_{1a2}+A_{1i2}}; \\ \nabla W\left|{}_{1a2}=\right.\frac{1}{A_{1a2}}{\displaystyle {\oint }_{\hspace{0.17em}\mathsf{\Gamma }} Wnd\mathsf{\Gamma }} \end{gathered}$$

(10)

in which Γ is the integral path along the circumference of each sub-triangle (e.g., 1a2). The weighting functions $w_a$, $w_b$ and $w_c$ are defined as:

$$\begin{gathered} w_a=\frac{\left(g_b{g}_c+{\epsilon }^2\right)}{{\left(g_a^2+g_b^2+g_c^2+3{\epsilon }^2\right)}^{}};w_b=\frac{\left(g_a{g}_c+{\epsilon }^2\right)}{{\left(g_a^2+g_b^2+g_c^2+3{\epsilon }^2\right)}^{}} ; \\ {w}_c=\frac{\left(g_a{g}_b+{\epsilon }^2\right)}{{\left(g_a^2+g_b^2+g_c^2+3{\epsilon }^2\right)}^{}} \end{gathered}$$

(11)

where $\epsilon$ is a small value (on the order of 10⁻⁴ and less), and $g_a$, $g_b$ and $g_c$ are the functions of the gradients of variables in the cells surrounding cell i (i.e., cells a, b and c in Figure 1a) defined by:

$$g_a=\Vert \nabla W_a{\Vert }_2^2;g_b={‖\nabla W_b‖}_2^2;g_c={‖\nabla W_c‖}_2^2$$

(12)

In Equation (12), ${‖\nabla W‖}_2^2$ is the second norm of the unlimited gradient of a specified variable. By definition, the norm L₂ of a vector is the sum of the squares of its elements. Using the above equations, the limited gradient of a variable may be determined, and the data may be reconstructed at the cell boundaries.

The slope and friction source terms, G₁ and G₂, in Equation (3) may be defined as, [39]:

$$G_1=\left(\begin{array}{c}0\\ -gh{S}_{0x}\\ -gh{S}_{0y}\end{array}\right)\hspace{1em};G_2=\left(\begin{array}{c}0\\ -c_f{U}_x\left|U_x\right|\\ -c_f{U}_y\left|U_y\right|\end{array}\right)\hspace{1em};c_f^{}=\frac{g{n}^2}{h^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$$

(13)

where S_0x and S_0y are the bed slopes in the x- and y- directions, respectively, c_f is the bed coefficient of friction, and n is the Manning roughness. During the numerical computations, G₁ and G₂ are calculated based on the data obtained for the cell centers.

The Boussinesq assumption is the basis of all of the turbulence eddy viscosity models. It relates the Reynolds stresses to the mean velocity gradients via the eddy viscosity. Using this assumption in the averaged Reynolds stress models, the effects of the Reynolds stress in the shallow flows (G₃ in Equation (3)) may be written as:

$$G_3=\left(\begin{array}{c}0\\ \frac{\partial }{\partial x}\left({\upsilon }_{t}h\frac{\partial U_x}{\partial x}\right)+\frac{\partial }{\partial y}\left({\upsilon }_{t}h\frac{\partial U_x}{\partial y}\right)\\ \frac{\partial }{\partial x}\left({\upsilon }_{t}h\frac{\partial U_y}{\partial x}\right)+\frac{\partial }{\partial y}\left({\upsilon }_{t}h\frac{\partial U_y}{\partial y}\right)\end{array}\right)$$

(14)

where ${\upsilon }_t$ is the turbulence eddy viscosity. To discretize this term, a semi-implicit method may be applied. When the viscosity is large, it is required to discretize the diffusive term implicitly. This term can be divided in two parts:

$$G_3=G_{3,}{}_{\perp }+G_{3,}{}_{‖}$$

(15)

where $G_{3,}{}_{\perp }$ is the orthogonal viscosity, and $G_{3,}{}_{‖}$ is the non-orthogonal viscosity. For the momentum component in the x-direction, the two components of viscosity may be calculated from [39]:

$$G_{3,}{}_{\perp ,\hspace{0.17em}x}={\displaystyle \sum _{j\in K_i}{\mathsf{\Gamma }}_{D_{\perp }}}{U}_{x,j}-\frac{{\mathsf{\Gamma }}_{D_{\perp }}}{h_i}{q}_{x,i}$$

(16)

$$G_{3,}{}_{‖}={\displaystyle \sum _{j\in K_i}{\upsilon }_{t,}{}_{ij}}{h}_{ij}\frac{d_{ij}}{d_{\perp ,ij}}\left(U_{x,B}-U_{x,V}\right)\left({\tilde{\alpha }}_{x,ij}{\tilde{n}}_{x,ij}+{\tilde{\alpha }}_{y,ij}{\tilde{n}}_{y,ij}\right)$$

(17)

in which ${\mathsf{\Gamma }}_{D_{\perp }}={\upsilon }_{t,ij}{h}_{ij}\frac{\left|n_{ij}\right|}{d_{\perp ,ij}}$ is the orthogonal diffusion, $q_{x,i}$ is the unit discharge at a cell i in the x-direction, $h_{ij}$ and ${\upsilon }_{t,}{}_{ij}$ are the averages of the depth and turbulent eddy viscosity in cells i and j, $U_{x,B}$ and $U_{x,V}$ are velocities in the x-direction at points B and V in Figure 4, ${\tilde{\alpha }}_{ij}=({\tilde{\alpha }}_{x,ij},{\tilde{\alpha }}_{y,ij})$ is a unit vector perpendicular to the line that connects the centers of the cells i and j, and $d_{\perp ,ij}$ is the projection of the distance between the two cell centers i and j over a line perpendicular to the common face of the two cells. All of the variables in Equation (16) are evaluated at time t_n except the unit discharge $q_{x,i}$, which is calculated at t_n+1. Therefore, no additional system of equations must be solved to increase the computational cost. The turbulence terms in the y-direction may be calculated similarly.

In the above equations, the eddy viscosity term may be computed using any of the turbulence modeling theories. The details of methods are presented by Alamatian and Jafarzadeh [39].

4. Modeling of Dry Bed

Dam break simulation is possible in both dry and wet bed conditions. In 2012, Alamatian and Jafarzadeh [39] considered the initial computational domain as wet. In the dry bed simulation domain, the riverbed of the dam is initially considered as dry. In the wet bed modeling, by approaching the water depth as zero in every single cell ($h\approx 0$), the flux increases sharply and creates instability in numerical modeling. It is possible to consider a minimum initial water depth in computational domain cells to avoid a dry bed condition; however, the accuracy of this method for dam-break simulation is limited. Various methods for modeling dry and wet bed conditions have been introduced in different studies [22]. In some methods, the computational domain cells can consider both the dry and wet domain of the solution [6].

In this paper, the active method is used for wet-dry behavior [22]. This method of work is based on activating and deactivating each computational cell for wet and dry bed conditions, respectively. These zones may be added or removed from the computational domain, and computational cells would be directly considered as active (wet) and inactive (dry) in the solution matrix. In this paper, the cell is assumed dry when its water depth is less than 1.0 mm. When the cell is inactive (dry bed), the computational flux in the boundary of that cell is ignored. In this method, the bed elevation is updating for obtaining an exact balance at the wet-dry front between the bed slope and the hydrostatic pressure term for hydrostatic conditions. If the wet-dry front occurs between the cells Ci and Cj, the modified bed elevation at the front is defined as [22]:

$${\Delta }_{zb,ij}=\left\{{}_{zb,j-zb,i\hspace{0.17em}\hspace{1em}otherwise}^{h_i-h_j\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}if}\right.$$

(18)

This method is time-efficient since the computational domain decreases by drying the bed. When the computational domain contains shocks, this method is more accurate [15].

4.1. Numerical Model Verification

The numerical model results are compared with experimental data to verify the performance of the code. The experiments were carried out by Soares-Frazao in 2007 [40]. The length and width of the considered flume are 5.6 and 0.5 m, respectively. A sluice gate is placed at x = 2.39 with a reservoir depth of 0.111 m upstream. A symmetrical triangular barrier with 0.065 m height and ±0.14 slope on each side, at x = 4.45, is considered. The water depth of 0.02 m at the downstream side of the barrier is considered, and the other parts are assumed as dry (Figure 5).

By immediately opening the sluice gate, which simulates the breaking of the dam, three probes start measuring the water depth in different locations. In the numerical simulation, the no-slip condition is applied for all wall boundaries [22]. A total of 14,618 unstructured triangular cells are considered with a Manning roughness coefficient and Courant–Friedrichs–Lewy (CFL) number of 0.011 and 0.9, respectively. The water level fluctuation of G₁ and G₂ during the 10 and 30 s after the dam break is investigated to verify the numerical simulation. The results of the numerical modeling and the experimental results are shown in Figure 6 and Figure 7.

Conservation of Mass Verification

One of the crucial issues in every CFD numerical code is to check the continuity or conservation of mass. To verify the ability of the code to conserve the mass in a dry canal, a numerical wave tank (NWT) with a barrier is considered. The NWT is defined with a length and width of 4.5 and 0.5 m, respectively. A triangular barrier is modeled in the center of the canal with a base of 1.5 m, a height of 0.5 m and 0.5 m in width (Figure 8).

The computational domain is considered as 12120 unstructured triangular cells. The no-slip condition is applied for all wall boundaries [22]. The initial conditions are composed of a water column with a height of 2 m and a length of 1 m in which its width is equal to the canal’s width. The acceptable error of 10⁻⁸ is considered for the numerical code. Figure 9 shows the initial conditions of the test.

Figure 10 shows the flow passing through the barrier and the stability of water above the barrier at different time steps. The numerical test result is compared with the analytical method. In this regard, the volume of water and barrier volume are added together, and by dividing this volume with the solution domain, the water level after stability can be obtained.

Table 1. Analytical solution and water level after stability.

Properties	Water Volume (m3)	Barrier Volume (m3)	Total Volume of Barrier and Water (m3)	Domain Area (m2)	Water Level after Stability (m)
Values	2 × 1 × 0.5 = 1	[(0.5 × 1.5)/2] × 0.5 = 0.188	1 + 0.188 = 1.188	4.5 × 0.5	1.188/(4.5 × 0.5) = 0.53

presents the analytical solution and water level after stability.

Accordingly, the water level will be equal to 0.53 m along the canal after reaching the stable stage. Figure 11 shows numerical and analytical solution results in comparison. It is observed that numerical and analytical solution results are in good agreement, and the conservation of mass is preserved favorably.

5. Modeling of the Dam Break

The Golestan dam has been constructed on a river in the north of Iran. It has a length of 30.6 m at the bottom with a crest length of 130 m, it is 20 m in height, and its reservoir volume is 2.7 million cubic meters. The area of the watershed dam is 1,010,000 hectares. The Golestan dam is a historic barrier dating back 150 years that was built using traditional materials. This dam has been exposed to various earthquakes. The changing of the land use upstream of the dam has caused increasing water to enter it and increased the likelihood of failure. In this paper, the length of 5400 m of the river downstream of the dam has been modeled. The topography of the area, dam reservoir and pathway plan have been obtained. Figure 12 shows the river plan and land risks.

Two bridges, B1 and B2, are located at distances of 3900 m and 5320 m downstream of the dam, respectively. These two bridges are critical for the area, and the dam break could cause a serious problem for these bridges. Their properties are presented in

Table 2. Bridges Properties.

Bridge Details	B1 Bridge	B2 Bridge
Pier radius	0.6 m	0.6 m
Number of piers	48	16
Bridge height	6 m	7 m

.

Only bridge piers have been applied in the modeling, and due to the size of the computational domain, a slip condition has been used at the wall boundaries of the bridge piers [22]. At the downstream end, the dependent variables (h, qx and qy) were interpolated from the solution domain, which is reasonable due to the supercritical nature of the flow [39]. The Courant–Friedrichs–Lewy (CFL) number is considered as 0.8. Figure 13 shows a change of maximum flow depth with different numbers of computational cells at B1 and B2 bridges. According to this figure, 35,628 unstructured triangular cells were chosen.

Due to the sand-bed river without vegetation downstream of the dam, the value of the Manning roughness coefficient varies between 0.025 and 0.033 [41]. The maximum depths of flood flow versus Manning roughness are shown in Figure 14 for the points attached to bridge B1 and B2. It is observed that with increasing Manning roughness from 0.025 to 0.033, the maximum depth of the flood wave increases by about 11%. For reliability, Manning roughness is considered as 0.033.

6. Results and Discussion

The dam break is simulated using the above information. The analysis continued up to the point where water level approaches the elevation of the bridge. The travel time duration of the first wave that hits the bridges is calculated. Figure 15 shows flood zoning in different periods. In this figure, the river border and computational domain are distinguished by red lines. Flood zoning indicates that water depth decreases gradually in the reservoir after the dam breaks, and it takes 22 min and 30 s.

Figure 14 and Figure 15 show a closer look at the flood wave when it arrives at the B1 bridge at a time of 620 and 660 s, respectively. The flow depth distribution across the river and the bridge’s piers are seen in Figure 16 and Figure 17.

The velocity vectors and the contours of streamlines at the B1 bridge at time 1100 s are shown in Figure 18. According to this figure, the proper performance of the wall boundary condition is confirmed.

Flood flow hydrographs at bridges B1 and B2 are presented in Figure 19. It can be seen that the maximum flood discharge at the locations of bridges B1 and B2 is 120 (at time 16.5 min) and 110.2 (at time 28.5 min), respectively.

Figure 20 and Figure 21 present the water level at B1, B2 pier locations. The B1 deck bridge would suffer from flooding at 11 to 40 min after the break and encounter dangerous waves with greater heights than the deck bridge elevation. It is observed that the height of the waves reaches up to 4 m above deck bridge level, which could cause a disaster. According to Figure 20, it is observed that 11 min after the dam break, the height of the flood wave reaches the level of the bridge deck. Consequently, for this bridge, there are only 11 min after the breaking of the dam to apply necessary proceedings and evacuate the adjacent roads. Precautionary measures should be taken during the safe periods.

The results based on Figure 21 showed that the B2 bridge would be flooded from 21 to 50 min. The maximum wave height of 9.54 m will occur about 27 min after the breaking of the dam, which is 2.5 m greater than the bridge elevation. Therefore, for bridge B2, there is only 21 min available before the bridge is completely flooded, and road evacuation should be considered to prevent future risks.

7. Conclusions

Real-scale dam-break modeling using 3D numerical models has a high computational cost. It is very useful to use 2D models that can reduce the cost of calculations while being accurate. In this paper, a real-scale dam-break wave was simulated using the 2D finite volume Roe-TVD method. For this purpose, a numerical code was developed to solve the 2D depth average, shallow water equations on unstructured triangular cells considering turbulence terms and a dry bed front. In order to validate the code, at first, available experimental data were considered, the water level in the flume has been predicted and compared with the experimental data. Furthermore, a numerical tank was simulated to study the capability of the model in passing the flow over a barrier in the dry canal, conserving mass and to reach a steady flow case above the barrier level. Comparing the analytical and numerical solutions indicated that the conservation of mass is satisfied. After verifying the model, the real-scale dam break has been simulated, and the flow behavior from encountering the two bridges was analyzed along the pathway. The flood wave arrival time to the bridges, flooded area and the duration of flooding of the bridges were studied. The results of the dam break simulation showed that B1 and B2 bridges are at risk of flooding in the case of a dam break. Long waves affect the structures and can cause severe damage. Vehicle traffic should be banned at the moment of dam break. Furthermore, it is necessary to evacuate any vulnerable objects from these areas. The maximum required time for the evacuation of the two bridges is 12 and 21 min, respectively.