Large-Scale Modeling of Dam Break Induced Flows

Abstract

Nowadays, mathematical modeling of the dam break flooding flows has become very interesting and challenging problem for scientists and engineers. In order to analyze possible risks and safety, mapping of flooded areas corresponding to different scenarios and conditions for the destruction of dams and reservoirs is significant. The results of mathematical modeling of large-scale flows in areas with a complex topographic relief were presented in this chapter. The flow was numerically simulated by the volume of fluid (VOF) method-based solver interFoam of OpenFOAM package, which solving the Reynolds-averaged Navier-Stokes equations (RANS) with the k-ε turbulence model. The mathematical model adequacy is checked by comparing with experimental data. The efficiency of the applied technology is illustrated by the example of modeling the breaking of the dams of the Andijan (Uzbekistan) and Papan (near the Osh town, Kyrgyzstan) reservoirs.

1. Introduction

The mathematical modeling using the advanced applied packages of computational fluid dynamics is an efficient tool for predicting various man-made and natural phenomena. The water flood of the regions and settlements takes the first place among natural spontaneous cataclysms in terms of the repetitions, extension area and annual material damage.

In this chapter, the problem of the prediction of consequences of a large-scale man-made disaster caused by the dam break is posed by the example of the Papan (near the Osh town, Kyrgyzstan) and Andijan (Uzbekistan) reservoirs.

The adopted mathematical model is presented in Section 2. The general technology of the numerical solution of the adopted mathematical model, the initial and boundary conditions, the model discretization methods are presented in Section 3. Section 4 contains the results of verifying the package OpenFOAM by various dam break test cases. Finally, the examples of using the OpenFOAM package for modeling the breaks of the dams of the Papan and Andijan reservoirs are considered in Section 5.

2. The mathematical model of the dam break flows

The mass and momentum conservation laws for viscous incompressible liquid in the absence of mass forces lead to the following unsteady Navier-Stokes equations [1]:

$$\partial \left(\rho \overline{u_i}\right)/\partial x_i=0,$$E1

$$\frac{\partial }{\partial t}\left(\rho \overline{u_i}\right)+\frac{\partial }{\partial x_j}\left(\rho \overline{u_i}\overline{u_j}+\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial \overline{{\tau }_{\mathit{ij}}}}{\partial x_j},$$E2

where $\overline{u_i}$ are the mean velocity components, ρ is the density, $\overline{p}$ is the mean pressure, $\overline{{\tau }_{\mathit{ij}}}=\mu \left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)$ is the mean tensor of viscous stresses, μ is the dynamic viscosity. The averaging is done in time, and the prime denotes the fluctuation part of the velocity. In the presence of external forces, it is necessary to augment these equations by the corresponding terms.

To close the systems of Eqs. (1) and (2), it is necessary to introduce some turbulence model. Most turbulence models employed in practice are based on the concepts of turbulent viscosity and turbulent diffusion. For the flows of general form, the turbulent viscosity introduced by Boussinesq, which couples the Reynolds stresses with the mean flow gradients, may be written in the following form [1]:

$$-\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\mu }_t\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)-\frac{2}{3}{\mathit{\rho \delta }}_{\mathit{ij}}k.$$E3

The turbulence kinetic energy k and its dissipation rate ε are determined from the following transport equations [1]:

$$\frac{\partial \left(\mathit{\rho k}\right)}{\partial t}+\frac{\partial \left(\rho \overline{u_j}k\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}+P_k-\mathit{\rho \epsilon },$$E4

$$\frac{\partial \left(\mathit{\rho \epsilon }\right)}{\partial t}+\frac{\partial \left(\rho \overline{u_j}\epsilon \right)}{\partial x_j}=C_{\epsilon 1}{P}_k\frac{\epsilon }{k}-\rho C_{\epsilon 2}\frac{{\epsilon }^2}{k}+\frac{\partial }{\partial x_j}\left(\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j},$$E5

where $P_k={\mu }_t\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)\frac{\partial \overline{u_i}}{\partial x_j}$ is the rate of the generation of the turbulence kinetic energy by mean flow, and ${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon }$ is the turbulent viscosity.

The coefficients of the model have the following standard values [1]: $C_{\mu }=0.09,$$C_{\epsilon 1}=1.44,$$C_{\epsilon 2}=1.92,{\sigma }_k=1.0,{\sigma }_k=1.3$.

The method of determining the interface between two phases—water and air occupies special place at the modeling of the flow class under consideration. According to the main idea of the volume of fluid (VOF) method [1], one determines for each computational cell a scalar quantity, which represents the degree of filling the cell with one phase, for example, water. If this quantity is equal to 0, it is empty; if it is equal to 1, then it is filled completely. If its value lies between 0 and 1, then one can say, respectively, that this cell contains the free (interphase) boundary. In other words, the volume fraction of water α is determined as the ratio of the water volume in the cell to the total volume of the given cell. The quantity 1 − α represents, respectively, the volume fraction of the second phase—air in the given cell. At the initial moment of time, one specifies the distribution of the field of this quantity, and its further temporal and spatial evolutions are computed from the following transport equation [1]:

$$\frac{\partial \alpha }{\partial t}+\frac{\partial \left(\alpha \overline{u_i}\right)}{\partial x_i}=0.$$E6

The free boundary location is determined by the equation $\alpha \left(x,y,z,t\right)=0$. Therefore, the physical properties of the gas-liquid mixture are determined by averaging with the corresponding weight coefficient:

$$\rho ={\mathit{\alpha \rho }}_1+\left(1-\alpha \right){\rho }_2,\mu ={\mathit{\alpha \mu }}_1+\left(1-\alpha \right){\mu }_2.$$E7

Here, the subscripts 1 and 2 refer to the liquid and gaseous phases.

The essence of the VOF method implemented in the solver interFoam of the OpenFOAM package [2] lies in the fact that the interface between two phases is not computed explicitly, but is determined, to some extent, as a property of the field of the water volume fraction. Since the volume fraction values are between 0 and 1, the interphase boundary is not determined accurately; however, it occupies some region, where a sharp interphase boundary must exist in the proximity.

3. Modeling technology

3.1. Initial conditions

In the case of unsteady problem, it is necessary to specify for the initial values for all dependent variables. The values of all velocity components are equal to zero because according to the condition of the problem under study, there is no motion until the moment of time t = 0. The hydrodynamic pressure is also equal to zero since the used solver—interFoam calculates hydrodynamic pressure [2]. The turbulence kinetic energy and its dissipation rate have some small value, which ensures a good convergence of the numerical solution at the first integration steps. The initial distribution of the volume fraction α is nonuniform because not all the computational cells are filled with water.

3.2. Boundary conditions

The no-slip condition is specified at solid walls of the computational region, which gives the zero components of the velocity vector. The Neumann conditions are specified for the water volume fraction: and $\raisebox{1ex}{$\partial \alpha $}\!\left/ \!\raisebox{-1ex}{$\partial n$}\right.=0$; at all wall boundaries, the fixedFluxPressure boundary condition is applied to the pressure (hydrodynamic pressure) field, which adjusts the pressure gradient so that the boundary flux matches the velocity boundary condition for solvers that include body forces such as gravity and surface tension [2].

The boundary conditions for the turbulence kinetic energy k and its dissipation rate ε were specified with the aid of the technique of wall functions [1]. Systematic calculations performed in this work show that the minimum value of dimensionless distance y⁺ for all solid wall greater than 25, so we can use wall functions technique.

The influence of surface tension forces between the solid wall and the gas-liquid mixture were not taken into account.

The top boundary is free to the atmosphere so needs to permit both outflow and inflow according to the internal flow. That is why it is necessary to use a combination of boundary conditions for pressure and velocity that does this while maintaining stability.

3.3. Mesh generation and discretization of governing equations

BlockMesh utility was used for generation of mesh. The discretization of the computational domain is obtained by the control volume method [3]. The use of an upwind difference scheme for the convective and Gauss linear scheme for diffusion terms yields an acceptable accuracy of numerical computations.

The explicit Euler first-order method was used for the discretization of the unsteady term. Numerical solution of the unsteady equations coupled with the pressure was based on the PIMPLE method [2] with the number of correctors equal to 3.

The iterative solvers PCG and PBiCG—the methods of conjugate gradients and biconjugate gradients with preconditioning were used for solving the obtained system of linear algebraic equations. The procedures based on a simplified Cholesky’s incomplete factorization scheme DIC and on the simplified incomplete LU factorization DILU were used as preconditioners.

Mesh sensitivity was analyzed for four mesh of 60 × 25 × 25, 90 × 40 × 40, 135 × 60 × 60 and 150 × 80 × 80, indicating that the mesh size was important only in the vicinity of the leading front location. Since of unimportant differences between the wave fronts using mesh 135 × 60 × 60 and mesh 150 × 80 × 80, a first mesh was adopted to reduce computational effort. In this case, the minimum value of dimensionless distance y⁺ was more than 15 for all coordinate axes.

More detailed information about the boundary and initial conditions, discretization techniques, and the solution of systems of algebraic equations may be found in the work [2].

4. The OpenFOAM package verification

The problem on the liquid column collapse in a horizontal duct of rectangular cross section [4] is the first test problem. At the initial moment of time, the rectangular column of a viscous incompressible fluid is at rest. The column starts collapsing under the gravity force.

The numerical results (solid line) are compared with experimental data (markers) of the work [4] in Figure 1, where a = 0.05715 m is the water column width in the x-axis direction, g = 9.81 m/s² is the gravitational acceleration.

Figure 1.

Comparison of numerical and experimental data. a—water column height change, b—the flow leading front location.

A good agreement between experiment data and results of numerical computation shows the reliability and accuracy of the present numerical modeling.

Figure 2 shows the configuration of the next model problem. The model represents a reservoir of hexahedral shape whose length is 3.22 m, height 2 m and width 1 m [5].

Figure 2.

Test problem configuration.

The water column of height H = 0.6 m, length 1.2 m and width 1 m lies in the left lower corner of the reservoir. To measure the pressure force of fluid flow on the right reservoir wall, the corresponding pressure probe was placed at point Р with coordinates x = 3.22 m and y = 0.16 m. Besides, the levels h₁ and h₂ of the water free surface were measured in two sections at x₁ = 2.725 m and x₂ = 2.228 m. The density of water was equal to 998.2 kg/m³, and that of air was equal to 1.225 kg/m³.

As initial configuration of the modeling with OpenFOAM, the water column in the left bottom part of the domain is at rest. When diaphragm is removed suddenly simulation is started, due to gravity water column drops and starts to move into the empty part of the domain. The total duration of modeling in time amounted to 2.5 s. In Figure 3, we have several snapshots of some stages of the simulation.

Figure 3.

The motion of water column at different time step.

The diaphragm is suddenly removed at the moment of time t = 0, and the water column runs under the gravity force into the right empty part of the reservoir. At the moment of time t ≈ 0.65 s, the water reaches the right wall and impinging on it under the inertia force, moves upward. The flow is thinned as it moves upward along the right wall and at the moment of time t = 1.3 s, when the gravity force exceeds the inertia force; there occurs a reverse flow of water with the formation of a typical bend of its surface. The reverse wave formed in such a way reaches the main flow, impinges onto it and forms the secondary wave, and so on (at t = 1.5 s and t = 1.8 s). After the moment of time t = 2.2 s, the water inertia drops significantly, and the further consideration of the motion is of no practical interest.

Figure 4 shows the numerical (solid line) and corresponding experimental (markers) data [5] on the water column height in sections with coordinates х₁ = 2.725 m and х₂ = 2.228 m. The coincidence between these data is satisfactory up to the moment of time t ≈ 1.5 s for section х₁ = 2.725 m (Figure 4a).

Figure 4.

Water column height in different sections. х₁ = 2.725 m (a) and x₂ = 2.228 m (b).

After this moment of time, the reverse wave moving oppositely to the main stream impinges onto the free surface, and this gives rise to certain inaccuracies both in numerical and experimental data. The inaccuracies of such a kind were also observed in the work [6], where the well-known commercial package FLUENT was used for numerical modeling. The results of this work are presented in Figure 5 in a form similar to Figure 4. The water column height (h₁, h₂) has been nondimensionalized here by quantity Н, and the time is presented in the dimensionless form τ = t(g/H)^0.5, where t is the physical time, g = 9.81 m/s² is the free-fall acceleration, Н = 0.6 m is the initial water column height. One can conclude from a comparison of Figures 4 and 5 that the free surface levels in two different sections have been predicted more accurately in the present work.

Figure 5.

Water column height from the work [6] at х = 2.725 m (a) and 2.228 m (b). 1—experiment [5], 2—computation of the work [6].

The problem of determining the pressure P on the obstacle is very important at the solution of unsteady problems with free surface, in particular, at the interaction of forming waves with various obstacles. Figure 6 shows the numerical results (solid line) for the pressure of the fluid on the right wall at point P with coordinates (x = 3.2 m, y = 0.16 m) and the corresponding experimental data (markers). The exact pressure value at point Р cannot be measured because the pressure probes have a finite size—a circle about 90 mm in diameter.

Figure 6.

Pressure variation at point Р. Pressure at the center (a) and at the lower point (b) of the probe.

The numerical pressure at the pressure probe center (see Figure 6a) increases slowly with time, after the moment of time t = 1.5 s or after the second maximum the coincidence of experimental data with numerical results improves. The character of the numerical pressure at the lower point of the probe (see Figure 6b) agrees fairly well with the character of the variation of corresponding experimental data; however, the maximum values are slightly underestimated.

It is assumed at the numerical modeling that the diaphragm is removed suddenly that is its vertical velocity is infinite. On the other hand, there can be also different physical conditions, which are hard to take into account at the numerical modeling. A detailed analysis of the conditions for conducting the experiment shows that this velocity has a finite value. The verification experimental data under the same conditions give different results, which do not coincide with one another [5]. In addition, the above discrepancies between the computation and experiment after the moment of time t = 1.5 s can probably be explained by the two-dimensionality of the employed model. It is possible that the flow acquires a three-dimensional character at some points of the computational region.

Comparing the data of the present work (Figure 4) and the work [6] (Figure 5), one can assert that the numerical results of modeling the task under consideration, which were obtained with the aid of the open package OpenFOAM, are closer to the experimental data than the results obtained with the aid of the commercial package FLUENT.

The problem of the fluid column breakdown in a reservoir of rectangular shape with an obstacle [7] is the next, more complex test problem. The chosen coordinate system and the problem diagram without the geometric proportion preservation are shown in Figure 7.

Figure 7.

Configuration (on the left) and location of measurement pressure probes (on the right).

An open reservoir 3.22 m in length with cross section of 1 × 1 m² has been used in experiment. The reservoir was initially partitioned into two unequal parts by a vertical wall located in section х = 2 m. The water 0.55 m in height is located behind this wall, another reservoir part is empty. A container 40 cm in length with cross section of 16 × 16 cm² is located in this empty part of reservoir. The container left-face coordinate equals х = 0.67 m.

At the execution of experiment, the water column height and the fluid pressure on the container surface were measured. The location of measuring probes is shown in Figure 8. Four probes were used to measure the water column height: the one on the part filled with water—Н4, and the remaining ones Н1, Н2, and Н3 were located in the reservoir empty part. The coordinates of these probes are x = 0.5, 1.0, 1.5 and 2.66 m, respectively.

Figure 8.

Comparison of numerical (on the left) and experimental data [7] (on the right) at the moment of time t = 0.4 s.

The container was supplied with eight pressure probes: four probes on the exterior surface at points with coordinates z = 0.025, 0.063, 0.099 and 0.136 m, and the remaining four probes were located on the container upper side with coordinates x = 0.806, 0.769, 0.733 and 0.696 m. The probes on the exterior surface were located at the distance of 0.026 m to the left of the central line (y = 0), and the probes on the upper surface were located at the distance of 0.026 m of this line. The force exerted on the container on the water stream side was also measured in experiment.

The water was at rest up to the moment of time t = 0. At the moment of time t = 0, the separating wall was removed suddenly, and the water column ran into the reservoir empty part under the gravity force. The 180 × 60 × 80 computational grid was used, and the CPU time amounted to 6 s. The initial time step was 0.001 s, and it was varied further depending on the Courant number, which was equal to 0.85.

Figures 8 and 9 show the comparison of numerical and experimental data for the moments of time t = 0.4 and 0.6 s, respectively. The pictures of shooting done during the experiment are shown on the right. One can notice a fairly good visual coincidence of the numerical results with experimental data.

Figure 9.

Comparison of numerical (left) and experimental [7] (right) data at the moment of time t = 0.6 s.

The time of reaching the container by water flow both in experiment and at the numerical simulation is the same. Besides, the free surface shapes forming after the flow impact onto the container also coincide. One can note, however, that there are at the numerical modeling some imperfections of the free surface between water and the ambient medium—air.

Figure 10 shows the water flow height at two different points: in the reservoir and in the immediate proximity of the container. There is a fairly good agreement between them until the water returns from the back wall after the moment of time t ≈ 1.8 s. After that, the numerical data (solid lines) prove to be somewhat higher than the experimental ones (markers). At the moment of time of t ≈ 5 s, the secondary wave reaches the neighborhood of the probe Н2. This time is, however, equal to about 5.3 s at the numerical modeling. The general character of the variations of the numerical and experimental data nevertheless coincides.

Figure 10.

Water flow heights at points Н2 (left) and Н4 (right).

Figure 11 shows the moment of time t = 0.5 s when the wave reaches the container has been predicted with a good accuracy; however, the computed pressure value (solid lines) is slightly overestimated as compared to the experimental value (markers) (the left figure).

Figure 11.

Pressure at points Р2 (left) and Р7 (right).

The numerical values of the second pressure maximum at point P2 are, however, shifted in comparison with experiment to the right by 0.6 s, and at point P7, they are shifted by 0.5 s. As the experiment has shown, the moment of time when the flow reaches repeatedly the container (≈4.7 s) is seen in these figures fairly well. Besides, when comparing the numerical and experimental pressure values at points Р7 (the right figure) one can notice some differences. After t = 1.3 s, a small oscillation lasting for about 0.2 s takes place in numerical computations, which is not observed in experiment.

5. Dam break flooding flow modeling in real region

To illustrate the techniques of the application of numerical modeling of large-scale hydrodynamic computations, we consider the problem of computing the flood process in the areas near the dams of the Andijan and the Papan reservoirs (see Figure 12).

Figure 12.

Maps of the Andijan (left) and Papan (right) reservoirs.

It is to be emphasized here that the situation of a real breakthrough of the dam and the flood of the areas at the lower level is not modeled here but the fundamental possibility of using the above technology under the availability of necessary topography data is demonstrated. The topography data of Digital Terrain Elevation Data [8] were used in computations, which were converted subsequently into the STL format. The hexahedral background grid generated with the aid of the utilities blockMesh and snappyHexMesh of the OpenFOAM package was transformed into a three-dimensional surface, which is employed for modeling the flood process (Figure 13).

Figure 13.

Three-dimensional surface of the Andijan (left) and Papan (right) reservoirs areas.

For the Andijan reservoir, the computational field had the sizes 6000 × 4000 × 1500 m, the physical modeling time amounted to about 9 h for the 120 × 120 × 80 grid. Figures 14 and 15 show different stages of the flood in areas with real topology. The red corresponds to a pure water flow, and the blue corresponds to air flow (there is no water flow in blue regions). It is seen in Figure 14 that the leading front of water flow reaches during 240 s the lower boundary of the computational region passing a distance of about 6000 m, covers the most part of the area located downstream.

Figure 14.

Flow pattern for the Andijan reservoir.

Figure 15.

Flow pattern for the Papan reservoir.

The computational field for the Papan reservoir has the sizes 5000 × 5000 × 1300 m (see Figure 15).

The total computing time in the case of the 50 × 60 × 30 grid amounts to about 5 h. As is seen in Figure 15, after the moment of time t ≈ 200 s there forms a reverse flow (Figure 15dand c) and after the moment of time t = 260 s; it separates into two parts—one part is in the zone of reverse flows and the other continues its flow in the lower part of the river bed. The computations show that about 60% of the entire initial water volume remains in the zone of reverse flows.

6. Conclusions

The results of the mathematical modeling of complex hydrodynamic phenomena on the basis of unsteady three-dimensional Navier-Stokes equations describing the dynamics of a gas-liquid mixture with free boundary have been presented. The adequacy of the employed model has been verified by the example of the classical problems of computational fluid dynamics. Special attention has been paid to the accuracy of the computation of the water flow level and the gas-liquid flow pressure on the reservoir walls. The efficiency of the employed technology has been illustrated by the example of modeling the breaks of the dams of the Andijan (Uzbekistan) and Papan (near the Osh town, Kyrgyzstan) reservoirs. The developed technology is universal and can be used for the flood modeling for a real relief. It is shown that the relief features are a substantial factor.

Acknowledgments

The Dam break flooding flow modeling in real region part of this chapter was done within the framework of the Fulbright exchange program of the US State Department. Authors gratefully thanks to of Prof. Greg Olyphant from the Indiana University, Bloomington, USA for his supervision and to Sally L. Letsinger, PhD, research Hydrogeologist of the Center for Geospatial Data Analysis of Indian Geological and Water Survey of Indiana University, Bloomington, USA for helping to prepare all necessary topographic data.