Numerical Simulation of Turbulent Dam-Break Flow Using the Level Set Method
Ashraf Balabel*
Mechanical Engineering Dept., Taif University, Taif, Hawiyya, Kingdom of Saudi Arabia
Abstract
The present paper introduces numerical simulations for one of the most challenging problems of two-phase flows; namely, the unsteady movements of dam break flow considering the turbulence effects. A novel numerical method is developed and validated for solving such complicated problem by solving the Reynolds-Averaged Navier Stokes equations over a regular and structured two-dimensional computational grid using the control volume approach. The standard k-e turbulence model is applied for predicting the turbulence characteristics of the dam break flow. The transient evolution of the dam free surface is predicted by the level set method. The effects of the geometrical parameters of the initial dam shape and the density ration of the two phases on the dam front movement and dam topological changes are investigated. The obtained results showed a faster movement of the dam front in the downstream direction by increasing the dam height and the density ratio. Moreover, the topological changes of the dam free surface in later evolutions are also predicted.
Keywords
Dam Break Problem, Level Set Method, Numerical Simulation, Turbulence Modeling, Two-Phase Flow
Received: May 26, 2015
Accepted: June 13, 2015
Published online: July 15, 2015
@ 2015 The Authors. Published by American Institute of Science. This Open Access article is under the CC BY-NC license. http://creativecommons.org/licenses/by-nc/4.0/
1. Introduction
Many of industrial and engineering applications encounter a large number of two phase flows; e.g. atomization of liquid jets [Lefebvre (1989)], droplet dynamics [Fuster, Agbaglah, Josserand, Popinet and Zaleski (2009)], bubble flow [Sato (1975)]and other multiphase flow systems [Kolev (2007)].
Dam break problem is considered to be one of the most important two-phase applications in engineering and industrial fields, e.g. hydropower stations, marine hydrodynamics and coastal engineering. Consequently, dam break flow has been the subject of extensive experimental, theoretical and numerical investigations for more than a hundred year. Detailed literature review for such problem can be found in [Park, JinKim, Van (2012)].
Experimental studies on dam break problem have observed the shape of the water propagation front and its traveling distance in the horizontal direction. In particular, the representative experimental investigation, which is considered as a benchmark test for dam break problem, is referred to [Martin and Moyce (1952)]. This experiment was later repeated using another experimental techniques by [Koshizuka and Oka (1984), Stansby, Chegini and Barnes (1998)].
According to the experimental complexity of such problem, exact measurements of the interface shape are not available; however, some measuring data such as the reduction of the water column height and the horizontal traveling distance of the water front can be applied for the validation of the numerical methods through the quantitative comparison of the obtained results.
Moreover, the theoretical as well as the numerical treatments of the dam break problem have encountered some constraints related to the nature of the problem as it includes a transient, non-uniform free surface flow with large spatial and temporal gradients. Moreover, turbulence characteristics should be considered with the moving wave front which is driven by gravity developed turbulence. Therefore, most of the previous numerical treatments of such problem are usually simulated by neglecting turbulent stresses or using various empirical assumptions about the parameters of turbulence. Recently, there is a trend to replace the previous attempts with direct numerical simulations, large eddy simulations and advance turbulence modeling, see for more details [Oezgoekmen, Iliescu, Fischer, Srinivasan, and Duan (2007), Park, JinKim, Van (2012)].
According to the huge development in the numerical simulation of turbulent two-phase flows with wide range of length scales, carefully executed simulations in such context can virtually replace experiments [Eggers (1997)]. In general, the numerical predictions of turbulent dam break dynamics have been limited in accuracy partly by the performance of three key elements, viz.: development of the computational algorithm, interface tracking methods, and turbulence prediction models [Balabel (2013)].
A variety of numerical methods have been recently developed and validated to two-phase turbulent flow. However, an efficient and complete numerical method is not available. An extended review of numerical techniques applied for turbulent two-phase flow including adv/disadvantages can be found in [Shinjo and Umemura (2010)]. More recently, the author has developed a new numerical method, known as Interfacial Marker Level Set Method (IMLS), for predicting the complete dynamics of turbulent two-phase flows. This numerical method is validated and its accuracy is estimated through the performing of a wide range of industrial and engineering applications [Balabel (2011-2012-2013)]. This method is also further applied in the present paper and shortly explained in the following sections.
2. Computational Method
The physical configuration of the dam break problem is illustrated in Fig. 1. The width and the length of the computational domain are H and L, respectively. At the initial state, the water depth and width are ho and xo, respectively. The water column has the properties; density of 1000 kg/m3 and viscosity of 1.0´10-3 Ns/m2, while the surrounding air phase has the density of 1.0 kg/m3 and viscosity of 1.5´10-5 Ns/m2. The computational domain is extended in x- and y- directions and the computational structured grids are considered to be of uniform size in both directions.
The governing equations for the transient two-dimensional, isothermal and incompressible turbulent two-phase flow are first described in the next section. An explanation of the level set formulation is given as well. Following, the numerical method adopted for solving the appropriate governing equations with the associated boundary conditions are also discussed.
2.1. Reynolds-Averaged Navier-Stokes (RANS) Equations
The Reynolds form of the continuity and momentum equations for turbulent two-phase flow, called here RANS equations, at each point of the flow field including the gravitational field can be represented by the following equations:
(1)
(2)
Where the subscript a takes the values 1 and 2 and refers to each phase properties, i.e., the liquid and gas (l, g) phases, respectively. In the above system of equations, is the velocity vector, r is the density, p is the pressure, m is the molecular viscosity, is the strain rate tensor and is the turbulent stress tensor which are given as:
(3)
(4)
where is the Kronecker delta and are the average of the velocity fluctuations. The turbulent viscosity is defined as:
(5)
The turbulent kinetic energy k and its dissipation rate e can be estimated by solving the following equations:
(6)
(7)
The coefficients for the so-called STD k-e turbulence model are given as follows [Launder and Spalding (1974)]:
2.2. Boundary Conditions
Considering the computational domain, shown in figure 1, two kinds of boundary conditions are applied, the external domain and the interfacial boundary conditions. The numerical boundary conditions applied for external domain are based on the slip boundary conditions (ui. ni=0) for the bottom and sides of the computational domain, while zero stress boundary conditions are prescribed on the open upper boundary as explained by [Kaceniauskas (2008)]. The standard wall function for the turbulent boundary layer is employed so that both k and ε can be expressed as a function of the distance from the solid boundary and the shear velocity [Launder and Spalding (1974)].
In contrast to the previous two-phase boundary conditions [Brackbill, Lothe and Zemach (1992)], in which the surface tension effects are included in the momentum equations as an external body force, the present model sets the surface tension pressure explicitly in the interfacial boundary conditions formulation in the normal direction as:
(8)
In the present modelling, all the interfacial forces arising from the interfacial jump conditions are including in the pressure term in such a way:
(9)
The treatment of such terms in the adopted numerical algorithm can be seen in [Balabel (2012)]. The above system of equations and boundary conditions, Eqs.(1-9) are solved simultaneously to determine the topological changes of the free surface of the dam during its evolution. The level set method is applied for advection of the deformable interface as described in the following section.
2.3. Level Set Function
The level set method is a type of capturing methods where a defined phase function f, is smoothed over the entire computational domain. The level set function at any given point is taken as the signed normal distance from the interface with positive on the liquid phase (i.e. f>0), and negative on the gas phase (i.e. f<0). Consequently, the interface is implicitly defined as the zero level set of the level set function.
The transport equation of the level set function can be described by the following equation:
(10)
where is the velocity vector. The geometrical properties of the interface, (normal vectors and curvature), can be defined as:
(11)
The original work of level set method in two-phase numerical simulation is referred to [Sussman, Smereka and Osher (1994)]. A review of such works can be found in the cited review [Sethian and Smereka (2003)]. However, the application of the level set method in capturing the moving interfaces in turbulent flows is indeed very scarce.
2.4. IMLS Numerical Scheme
In the present paper, the so called Interfacial Marker Level Set method (IMSL) is applied for predicting the dam break problem in turbulent flow. This numerical method is developed and validated by the present author and applied for a wide range of industrial and engineering applications. More details about IMSL method can be found in [Balabel (2013)].
3. Validation
The physical phenomena, such as breaking wave, shockwaves, subcritical and supercritical flows, occurring in some cases of dam break problem includes it into the complex applications of free surface flows. The dam break problem introduces a great challenge of two-phase flow as it is a transient, non-uniform free surface flow with large spatial and temporal gradients [Yang, Lin, Jiang and Liu (2010)]. Moreover, including the turbulence in such case makes it much complex regarding its modeling from the numerical point of view. Recently, there are few attempts to simulate the dam break problem with LES as they allow a much deeper insight into the associated dynamics [Oezgoekmen, Iliescu, Fischer, Srinivasan and Duan (2007)].
In the present section, the present developed numerical method is validated by performing the dam break problem in turbulent regime and comparing the obtained numerical results with several sets of previous experimental measurements of [Stanby, Chegini, and Barnes (1998)] and the numerical results of [Shigematsu, Liu and Oda (2004)]. The dimensions of the considered problem are (L=2m, H=0.25m, ho=0.1m and Xo=1m).
Figure 2 shows the comparison of the obtained numerical results for the dimensionless displacement of the dam break in turbulent flow The dimensionless time and displacement are given by and x*=x/ho. First, a grid independent study is carried out using different grid resolutions. It can be shown that, starting from 200x50 grid points in x- and y-directions, the present numerical results are in better agreement with the experimental measurements than the previous numerical results of [Shigematsu, Liu and Oda (2004)].
4. Results and Discussion
Figure 3 illustrates the transient free surface evolution of the breaking dam along with the predicted turbulent kinetic energy distribution at different time instants. It can be noticed that, during the wave propagation of the dam break, turbulence is generated and it concentrates near the moving wave front. In general, it can be concluded that, the present numerical method can successfully predict the complex topological changes of the dam break problem in turbulent flow.
4.1. The Effect of Dam Height and Density Ratio on Dam Break
Figure 4 (a, b) shows the effects of dam height h*=ho/Xo and the density ratio r*=rwater/rair on the downstream movement of the dam front. As shown in the figure, the increasing of the dam height leads to a fast movement of the dam front that can reach to the other side of the computational domain in short time. The same effect can be also observed by increasing the density ratio of the two phases.
4.2. Prediction of Dam Free Surface at Later Evolutions
Figure 5 shows the topological changes of the dam free surface at later time of breaking for h*=0.2 and density ratio of 1000. It is clear that the present numerical method can successfully capture the complex topological changes of the turbulent free surface during the dam breaking. The formation of an entrapped air cavity is clearly visible. The treatment of such cavity during the numerical simulation is carried out in normal way without any numerical ad-hoc assumptions.
5. Conclusion
In the present paper, the dam break problem in turbulent field is numerically simulated. The numerical method is based on solving the RANS equations with a developed numerical method known as IMLS method. The tracking and capturing of the dam free surface is carried out by the level set method. The comparison of the obtained numerical results showed a good agreement with the previous experimental measurements better than previous numerical results from the previous publications. Generally, it can be concluded that the implementation of the developed numerical method in three dimensional two-phase turbulent applications can be straightforward and the extension of the numerical model to include a wide range of industrial and engineering applications could be easily task.
References