Abstract
A numerical model based on Navier-Stokes equation is developed to simulate the interface evolution of flooding flows. The two-dimensional fluid domain is discretised by structured rectangular elements according to finite volume method (FVM). The interface between air and liquid is captured through compressive interface capturing scheme for arbitrary meshes (CICSAM) based on the idea of volume of fluid (VOF). semiimplicit method for pressure linked equations (SIMPLE) scheme is used for the pressure-velocity coupling. A second order upwind discretization scheme is applied for the momentum equations. Both laminar flow model and turbulent flow model have been studied and the results have been compared. Previous experiments and other numerical solutions are employed to verify the present results on a single flooding liquid body. Then the simulation is extended to two colliding flooding liquid bodies. The impacting force of the flooding flow on an obstacle has been also analyzed. The present results show a favourable agreement with those by previous simulations and experiments.
1. Introduction
The interaction of flooding flows with structures is a classic and important problem in many engineering applications, such as the green water loading [1, 2], dam break flowing over an obstacle [3, 4], and the wave breaking problem [5]. The unsteady loads caused by the flooding flow on fixed or floating structures have unnegligible influence on the design of marine structures [6] and marine vessels [7, 8]. In addition, some natural disasters featuring as flooding flow, such as tsunami, lava flow, and landslide [9], have attracted much of our attention because of their huge damage to the environment and economy. Thus, it is of great interest to study the behavior of flooding flow.
Among various problems related to flooding flow, the dam break problem with consequent wall impact is widely used to benchmark various numerical techniques that tend to simulate interfacial flows and impact problems [10]. Koshizuka [11] experimented on a liquid column collapsing and hitting an obstacle. Zhou et al. [12] performed a dam break flow experiment in a water tank to analyze the impact on a rigid wall. Besides the experimental studies, numerical simulation has become a more efficient and economical way in the investigation of flooding flow with the development of high performance computers. Among these numerical simulations, the potential theory [13–15] has been used for decades. A more popular way is to employ the incompressible N-S equations without neglecting the fluid viscosity. Generally, these viscous models are more suitable for handling strongly nonlinear problems. Scardovelli and Zaleski [16] provided a comprehensive review of the viscous analysis on the interface capture of multiphase flow problems. Abdolmaleki et al. [17] simulated the impact flow on a vertical wall resulting from a dam break problem.
The modelling of interface evolution is both critical and challenging in the numerical simulation of flooding flow because of the complex air-liquid interaction [18, 19]. There are mainly three numerical schemes to track the free surface of the flooding flow which are boundary element method, smoothed particle hydrodynamics (SPH) method, and volume of fluid (VOF) method [8]. Using the dynamic and kinematic conditions on free surface, the boundary element method usually employed a time stepping scheme to track the free surface based on nonviscosity assumption [20]. The SPH method offers a variety of advantages for fluid modelling, particularly for the splashing droplets problem [21, 22]. Colagrossi and Landrini [23] used SPH method to predict the dynamic behaviour of flooding flow. Proposed by Hirt and Nichols [24], the VOF method has become popular for calculating interfacial flows. Panahi et al. [25] analyzed the numerical simulation of floating or submerged body motions based on a VOF-fractional step coupling. Kleefsman et al. [26] numerically investigated a dam break problem and water entry problem. Hänsch et al. [27] simulated the two-phase flows with multiscale interfacial structures, which is a dam break model with an obstacle.
In the present study, a numerical tool based on the Navier-Stokes equations is developed to simulate the viscous flooding flow in a water tank. The paper is divided into five sections. In Section 1, the recent developments of the numerical methods for the dam breaking problem are introduced, which mainly include smoothed particle hydrodynamics method and volume of fluid method. In Section 2, the main methodology and the numerical scheme are presented in detail, especially the discretization of the fluid domain and CICSAM used to capture the air-liquid interface. In Sections 3, 4, and 5, we, respectively, studied the dam breaking problem of a single fluid body, the collision between two liquid bodies, and the impact of a flooding flow on an obstacle. A detailed analysis on the interface evolution is presented and the impact force of the flooding flow is investigated.
2. Methodology
2.1. Governing Equations
The Navier-Stokes equation and continuity equation on viscous flows can be described as follows [24]:where and are the velocity components in directions, respectively, is the mixture density of the fluids, is the pressure, and is the gravitational acceleration.
By applying the gauss theorems, the integration of (1) over each cell with respect time can be presented aswhere the subscript represents the centre of the control volume, is the centre of the cell boundaries, is the volume of the control volume, is the volumetric flux at the face of control volume, and is the area vector of the control volume face.
Three schemes, explicit scheme, the Crank-Nicolson scheme, and Euler implicit scheme can be used for the temporal discretisation of the N-S equations. In this paper, Euler implicit scheme is chosen and (2) of the velocity in direction is discretised asSimilarly, we can discretise (2) of the velocity in direction.
In order to keep the numerical stability, the normalized variable diagram (NVD) method is applied when handling the convection item. In this paper, a second-order upwind discretization scheme is used to calculate the values on the control volume faces. The algebraic equation obtained ultimately for each variable in each control volume is described as follows:where denotes the scalar quantity of general variable, denotes the coefficient of the linear equations, and the subscripts , , , and denote the upper, bottom, left, and right boundaries of the cell.
For the interfacial flow, the water and air are assumed to be one type of fluid with different densities. Thus, we can use a single set of (1) to describe the entire flow field and interface terms. The mixture density can be described asin which and are the density of water and air, respectively, is the fluid volume fraction which is set to 1 in the water region, and 0 is in the air region. In the vicinity of the interface, is between 0 and 1.
Similarly, the mixture viscosity of the fluids can be described asin which and are the viscosity of water and air, respectively.
The conservative form of the scalar convection equation for the volume fraction is as follows:
The volume fraction of each phase is solved in all computational cells. In each time step, should satisfy the governing equations (1). In each cell, only knowing the value of is not enough to determine the local interface position and direction. As a step function, the volume fraction also often causes numerical diffusion when using finite difference method to obtain variable derivatives. In addition, the amount of flow convected over a cell during a time step should be less than the amount available in a doner cell. The computational grid should be fine enough with respect to the maximum flow velocity, also for a distinct gas-liquid interface. Thus, the inclusion of water and air with different densities will greatly complicate the numerical method.
As the real liquid has viscosity, the RANS (Reynolds-averaged Navier-Stokes) equations are used to solve the mean flow velocity and turbulent model is used for the closure of RANS equations [28]. In standard model, is the turbulent kinetic energy and is turbulent dissipation rate. Then the turbulent kinetic viscosity can be presented aswhere is empirical constant. Taking and as basic unknown variables, the corresponding transport equations arewhere and are the turbulent kinetic energy caused by mean velocity gradient and buoyancy, respectively. , , and are the empirical constant, and are the Prandtl Numbers of and . where is the component of the flow velocity parallel to the gravitational vector and is the component of the flow velocity perpendicular to the gravitational vector. .
2.2. Free Surface Capture
For the cells in the vicinity of the air-liquid interface, the volume fraction undergoes a step change, which is a challenge to simulate the surface. In present paper, we adopt the CICSAM (compressive interface capturing scheme for arbitrary meshes) method developed by Ubbink and Issa [29], which is capable of predicting a well-defined interface. In this method, the convection of fluid is simulated by weighting the available fluid in the donor cell with a weighting factor based on the face Courant number and the cell Courant number . The face fraction values are predicted with linear interpolation:where subscripts and represent the value at the centre and border of cell, respectively. Subscripts and indicate the donor cell and the acceptor cell, respectively, is the amount of element grids, is the volume of the cell, is the flux out of the donor cell, is the CICSAM weighting factors, and upwind form and downwind form are both used and switched flexibly to get smooth or sharp interface [30].
For transient calculations, the initial velocity and density fields can be specified according to the specific test cases. Also the implicit body force formulation can be used in conjunction with the VOF method to improve the convergence of the solution by accounting for the partial equilibrium of the pressure gradient and body forces in the momentum equations. In solving the Navier-Stokes and continuity equations, one needs physical properties density and viscosity distribution in the computational domain. The N-S equations are solved in every cell with fluid containing.
To increase the stability and the accuracy of the present simulation, Euler implicit scheme is applied for the temporal discretisation, which can guarantee the stability of iteration process even with a relatively large time step. Meanwhile, the courant numbers and are adjusted to be small enough to ensure the accuracy of the simulation. The second order procedure for solving N-S equations can also increase the accuracy, as the central difference finite volume method is used for viscous terms. Also a second-order upwind discretization scheme is used to calculate the values on the control volume faces. In handling convection item, the normalized variable diagram (NVD) method also stabilizes the computation effectively.
3. The Interface Evolution of a Single Flooding Body in a Square Tank
In this section, a dam break experiment [12] involving complex free surface evolution is used to verify present simulations, which is performed in a tank measuring 3.22 × 1 × 1.8 m. Two probes A and B located at m and m on the tank bottom are used to measure the water heights. According to the experimental configuration, a 2D numerical model is established, as shown in Figure 1. A water column located on the left side of the tank is contained by the tank boundaries and a vertical board. Once the board is lifted, the water with an initial height and width flows freely.
3.1. Meshing Strategy and Boundary Conditions
After some meshing tests, 24804 quadrilateral structured grids which have been proven enough in the simulation [31] are finally used to describe the computational domain. The maximum grid size is m and the minimum grid size is m. Generally, denser grids should be distributed in the vicinity of the evolving interface between the air and water. While calculating the velocity and pressure fields, the fractional step method of Kim and Choi [32] is applied. It should be assured that, at each time step, the fluid particles are evolving within the local grids. For the maximum water velocity approximated as below, a steady time step is considered proper:
The nonslip wall condition, which requires the fluid to stick to the wall, is imposed on the boundaries except the top one on which an opening boundary is set. The relative pressure is set as 0 Pa. The air at standard atmospheric condition is used as the primary phase while fresh water at 20 degrees Celsius is used as the secondary phase.
For each time step, the residual of continuity equation is achieved at a high level of within to minimize the effect of numerical diffusion. A maximum iteration number of 100 is set to save the computer resources.
3.2. Verification of the Present Results
Figure 2 presents the snapshots of the flow at different time instants. The time is nondimensionalized as . The comparison in Figure 2 shows that present simulation of interfacial flow agrees quite well with those through SPH method by Gao et al. [31], especially for the results by the laminar model.
![(a) Reference [31]](https://static.hindawi.com/articles/sv/volume-2015/794069/figures/794069.fig.002a.jpg)
(a) Reference [31]
(b) Laminar model
(c) model
In Figure 2, the water column is allowed to flow at . As time progresses, the flow front impacts on the right vertical wall at about . Then an upward water jet is formed instantly and the jet goes up until the gravitational effect counteracts the upward momentum. Due to the combination of the backflow and oncoming momentum, the flow overturns away from the boundary. As the flow head retouches the free surface (around ), a secondary jet is created by the impact at about . Comparing the interface snapshots in Figures 2(a), 2(b), and 2(c), we can see that the laminar model gives more similar results with those by SPH than the turbulent model. At , the interface by laminar model is smoother than that by turbulent model. At , we can find that both the liquid fronts by SPH and laminar model almost reach the right boundary, while there is an obvious distance between the right boundary and the liquid front by the turbulent model. At , we can find a thinner liquid layer climbing up along the right boundary. However, at the same time instant, the liquid front by turbulent model is relatively rounder and shorter.
The computed water heights , at the monitoring points A and B during the whole time stepping process compared to the experimental results in Figures 3 and 4. The time history of the water front is given in Figure 5. The dimensionless displacement is given by . From Figures 3, 4, and 5, we can see the present simulations are in good agreement with those by experiments [3] especially for the laminar model.
4. Impact of Two Colliding Flows
To evaluate the splashing problem, a case of two colliding flows is studied in this section. The computational domain shown in Figure 6 is discretised by 55728 quadrilateral grids. The simulation tank size is 1.8 m × 6.44 m. The splashing jet caused by the collision between two collapsing liquid bodies with the dimension shown in Figure 6 is simulated and compared with that by two symmetric bodies with the dimension of 0.6 m in height and 1.2 m in length, as shown in Figure 7.
(a) Symmetric liquid bodies
(b) Asymmetric liquid bodies
From Figure 7(a), we can see that, for the collision of two symmetric liquid bodies, the interface evolves symmetrically with respect to the vertical centre line of the computational domain. At , the two liquid bodies have started collapsing. As time goes on, a splashing jet forms and then falls at both sides due to the gravitational effect. When the falling liquid bodies reconnect with the upper surface of the shallow water, two cavities are formed and two secondary jets appear at the laterals. The water splashing for two asymmetrical water bodies is shown in Figure 7(b), from which we can see when the liquid fronts collide, liquid jet splashes inclined to smaller liquid body. Then the jet slightly touches the wall and falls down on the horizontal interface. At , the left jet reconnects with the horizontal interface and a subsequent jet occurs.
5. Impact of Flooding Flow on an Obstacle
5.1. Simulation Setup
A more interesting case of dam breaking occurs when a small obstacle is located in the way of the water front [11]. The case has also been used by other authors for testing their interface tracking or interface capturing method [9, 27].
The dimensions of the computational domain, liquid body, and the flap are shown in Figure 8. The minimum grid size is chosen to be and the time increment at each step to be . An opening condition is imposed to the top boundary of computational domain with a relative pressure of 0 Pa. Nonslip wall condition is imposed on the other part of the boundary, which requires the fluid to stick to the wall. In the experiments some water spilled over the tank. To guarantee the mass conservation, the computation domain is set to be higher than the actual tank in the experiments.
5.2. The Dam-Breaking Flow without Discharging Board
In Figure 9, five snapshots of the simulation at different stages are shown together with the corresponding images of the video from the experiment [27]. In the experiments, the initial liquid body is actually blocked by a board. However, the period of the board lifting is neglected. Meanwhile, the three-dimensional water tank is simplified as a two-dimensional computational domain as in Figure 8. Thus, we can see an obvious disagreement between the experimental and the computational result at . The disagreement can also be found in Figure 9 that there is a thin liquid layer in the experimental results at . At , it is found that the liquid front gradually detaches from the splashing jet and several liquid droplets are formed. However, the strong splitting of the front of the splashing jet in the experimental results was not reproduced in the numerical simulation.
(a)
(b)
(c)
(d)
(e)
As time goes on to , the jet has impinged on the wall and entrapped the air beneath it. Meanwhile, it is found in the snapshot of both experiment and simulation that a small liquid tongue forms above the small flap. At this time instant, a thin uprushing liquid layer is formed on the wall. Due to the gravitational effect, the layer then begins to fall downwards. At , the liquid tongue continues to develop. The middle part of the liquid sheet becomes thinner and is lifted by the movement of the entrapped air. If the simulation has continued, the air cavity may burst from the location where the sheet is the thinnest. Combining the uprushing momentum and the gravitational effect, a liquid block was formed above the impinging point at . It is noted that some upper part of the simulation domain has been treated to be invisible for a better comparison.
The comparison of pressure at the monitor point between present simulation and that by Hänsch et al. [27] is presented in Figure 10. From the present simulation, we can see that when the liquid front begins to touch the flap, the pressure obtains a sudden increase and the first peak appears at about . Then the monitored pressure begins to fall back until the second peak is at about due to the collision between the liquid front and the right wall. The third small peak at about corresponds to the collision of the liquid tongue on the bottom wall (see Figure 9(e)). Generally, the present simulation agrees very well with the works done by Hänsch et al. [27], despite the time lag of the pressure peaks and the difference of the second peak amplitude.
5.3. The Dam-Breaking with a Discharging Board
In many engineering applications, the flooding flow needs to be controlled by a discharging board. The opening size of the flooding flow is controlled by the height , which is defined as the vertical distance between the lowest point of the board and the tank bottom. The interfaces are at different time instants while , 0.095 m, and 0.048 m, respectively, are presented in Figure 11. Form this figure, we can see that, under the three heights, all the liquid fronts have started to move at , among which the front at is expected to be the first one to touch the flap. Despite the velocity difference of the fronts, there is no much difference between the interface under the three heights at . At , the longest splashing jet was generated by the flooding flow as . At , a large droplet was separated from the fronts as and 0.098 m. The splashing jet as is the first one among the three to impinge on the wall while the splashing jet as has the most curved front. At time goes on to 0.5 s, both the splashing jets at and 0.098 have impinged on the wall while there is an obvious time lag at .
(a)
(b)
(c)
The torque about the upper end and the horizontal force of the board caused by the flooding flow are shown in Figures 12(a) and 12(b), respectively. For the initial state at , the force and torque on the board mainly resulted from the hydrostatic effect due to the gravitation. As the liquid starts to move, the hydrodynamic effect becomes stronger and leads to the first peak at about . Then the force and torque begin to decrease with time. From Figure 12, we can see that both the force and torque tend to zero at about and for and , respectively, which are the time instants as the liquid detaches from the board. For , the detaching time is longer than 0.7 s. Once the splashing jets hit on the right boundary, the suddenly disturbed air will cause a small fluctuation of the force of torque, which can be seen from the second small peaks of the three curves in Figure 12.
(a) Force
(b) Moment
The momentum about the lower end of the flap and the horizontal force on the flap caused by the flooding flow are shown in Figure 13. In this figure, there are generally three major peaks on the time history curves of the force and torque despite that the third peak at lags too much and is not shown in the figure. We also can find that as becomes smaller, the first peak appears earlier. This means that the increase of the board blockage is actually accelerating the front of the flooding flow before it touches the flap. However, the time interval between the first and second peaks is increasing as the height decreases, which means the front velocity of the flooding flow obtains larger horizontal velocity.
(a) Force
(b) Moment
6. Conclusion
In the present paper, the interface evolution of flooding flow and the subsequent liquid impact was simulated using a Navier-Stokes solver with a VOF-based interface capture scheme. The application of laminar and turbulence models were discussed. Although there is some difference between the results by laminar model and the turbulent model, both can give reasonable results of the dam breaking problem. In general, the comparison proves that the present simulation has captured the main characteristics of the interface evolution. The computational results coincide favourably with the experimental results by Koshizuka [11]. Satisfactory results were also obtained for collision of two liquid bodies and the impact of the flooding flow on an obstacle. The present numerical technique may be applied for prediction of the hydrodynamics of structures with liquid impact.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The author is grateful for the support of the National Natural Science Foundation of China (Grant nos. 11302057 and 11102048) and the Research Funds for State Key Laboratory of Ocean Engineering in Shanghai Jiao Tong University (Grant no. 1310).