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.
- dam break flows
- topographic data
- k-ε turbulence model
- free surface
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 :
where are the mean velocity components,
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 :
The turbulence kinetic energy
where is the rate of the generation of the turbulence kinetic energy by mean flow, and is the turbulent viscosity.
The coefficients of the model have the following standard values : .
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 , 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
The free boundary location is determined by the equation . Therefore, the physical properties of the gas-liquid mixture are determined by averaging with the corresponding weight coefficient:
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  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
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 ; 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 .
The boundary conditions for the turbulence kinetic energy
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 . 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  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 .
4. The OpenFOAM package verification
The problem on the liquid column collapse in a horizontal duct of rectangular cross section  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  in Figure 1, where
A good agreement between experiment data and results of numerical computation shows the reliability and accuracy of the present numerical modeling.
The water column of height
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.
The diaphragm is suddenly removed at the moment of time
Figure 4 shows the numerical (solid line) and corresponding experimental (markers) data  on the water column height in sections with coordinates
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 , 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 (
The problem of determining the pressure
The numerical pressure at the pressure probe center (see Figure 6a) increases slowly with time, after the moment of time
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 . In addition, the above discrepancies between the computation and experiment after the moment of time
Comparing the data of the present work (Figure 4) and the work  (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  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.
An open reservoir 3.22 m in length with cross section of 1 × 1 m2 has been used in experiment. The reservoir was initially partitioned into two unequal parts by a vertical wall located in section
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
The container was supplied with eight pressure probes: four probes on the exterior surface at points with coordinates
The water was at rest up to the moment of time
Figures 8 and 9 show the comparison of numerical and experimental data for the moments of time
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
Figure 11 shows the moment of time
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
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).
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  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).
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.
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
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.
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.