ISRN Applied Mathematics

Volume 2012 (2012), Article ID 521012, 16 pages

http://dx.doi.org/10.5402/2012/521012

## A Novel Algorithm of Advection Procedure in Volume of Fluid Method to Model Free Surface Flows

^{1}Faculty of Marine Technology, Amirkabir University of Technology, 424 Hafez Avenue, P.O. Box 15875-4413, Tehran, Iran^{2}Department of Civil Engineering, Ferdowsi University, Mashhad, Iran

Received 4 December 2011; Accepted 19 January 2012

Academic Editor: H. Du

Copyright © 2012 M. J. Ketabdari and H. Saghi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

In this study, the developed procedure of advection in volume of fluid (VOF) method is presented for free surface modeling. The fluid is assumed to be incompressible and viscous and therefore, Navier-Stokes and continuity are considered as governing equations. Applying Youngs’ algorithm in staggered grids, it is assumed that fluid particles in the cell have the same velocity of the cell faces. Therefore, fluxes to neighboring cells are estimated based on cell face velocities. However, these particles can show different velocities between two adjacent cell faces. In developed model, the velocity in mass center of fluid cell is evaluated to calculate fluxes from cell faces. The performance of the model is evaluated using some alternative schemes such as translation, rotation, shear test, and dam break test. These tests showed that the developed procedure improves the results when using coarse grids. Therefore, the Modified Youngs-VOF (MYV) method is suggested as a new VOF algorithm which models the free surface problems more accurately.

#### 1. Introduction

In the numerical computations of free surface flows such as water waves and splashing droplets, accurate representation of the interface is very important. The volume of fluid (VOF) method is a convenient and powerful tool for modeling the free surface flows, where the fluid location is determined using related function. In this method, the VOF function is averaged over each computational cell. It is set as unity, zero, and a value between unity and zero for full fluid, empty, and free surface cells respectively. Using this function, the VOF method is capable of modeling flows with complex free surface geometries such as rising bubbles [1], the merging and fragmentation of the drop [2]. In addition, in comparison with other methods, yet it is remarkably economical in computational point of view. It is due to requiring only one array for storing the VOF function and a simple algorithm to advect the function during each computational time step. Several volume advection techniques have been developed with the aim of maintaining sharp interface. The more famous ones are the simplified line interface calculation (SLIC) method of Noh and Woodward [3] such as SOLA-VOF [4] or its corrected form as in NASA-VOF2D [5], the VOF method of Hirt and Nichols [6], and the method of Youngs [7, 8]. VOF advection algorithm can be classified based on free surfaces reconstructing technique in each cell and the method of computing boundary flux integration. Some VOF methods represent free surface interfaces as a line parallel to one of the grid coordinates which are referred as piecewise constant scheme. Some of them are the methods used by Nichols et al. [4], Hirt and Nichols [6], Torrey et al. [5], and Duff [9] where free surface interfaces are constructed in a stair-shaped profile. The alternative methods are known as piecewise linear schemes. They are developed by Rider and Kothe [10], Harvie and Fletcher [11, 12], Geuyffier et al. [13], and Scardovelli and Zaleski [14, 15]. In these methods, oriented free surface interface is in a direction perpendicular to the locally evaluated VOF gradient. These schemes are complex but more accurate than their piecewise constant counterparts associated with more computational costs. In this paper, a new advection model as modified Youngs’ VOF (MYV) method is presented.

#### 2. Governing Equations

The fluid is considered to be Newtonian and incompressible. Therefore, 2D continuity and Navier- Stokes equations (NSEs) are used as where is time, velocity vector, hydrodynamic pressure, kinematic fluid viscosity, fluid density, gradient operator, and is gravity based body force. In turbulent flow, the effect of turbulence can be considered using eddy viscosity models [16]. Researchers have used different models such as - and - [17–20] to model turbulent flow. In this paper, the standard two equation - model is used, where the first equation involves turbulence kinetic energy and represents the velocity scale. The second one takes into account turbulent dissipation rate and represents the length scale. The two-equation - turbulence model accounts for the effect of turbulence as follows: where Values of , , , and are selected as in Table 1.

Finally, the Reynolds-averaged Navier-Stokes equations (RANSEs) are used to model the turbulent flows as To model the free surface by VOF method, a step function of is used. This function is expressed as [3] Based on this equation, varies versus both time and space. Derivations of (and its continuous spectrum ) are used to determine the direction and curvature of the free surface. As direction, position and curvature of the free surface is determined, the surface tension can be computed. The continuous spectrum of is able to represented as Combining (2.5) and (2.6), the colour function is obtained as where, has a specific value in each scalar cell as Specific technique must be used to discrete (2.7).

#### 3. Solution Algorithm and Stability Criteria

In this study, two-step projection method is used to solve the time-dependent NSE. The governing equations are discretized on a Cartesian-staggered grid system [21]. In this method, the two-step advancement algorithms are used as where is the velocity field in old time level, intermediate velocity field, the new velocity field, convection term, diffusion term, and the body force. The present method is explicit and first order in time, whereas the order of the space discretisation depends on the scheme used for the convective terms. The velocity field must satisfy the continuity equation. Therefore taking the divergence of (3.2), the Poisson’s equation for the pressure is obtained: So, at first, the intermediate velocity field is obtained by (3.1). Then, the pressure distribution in the new time level is obtained using (3.3). In this step, the new velocity field is calculated using (3.2). Finite difference and three-diagonal matrix algorithm (TDMA) methods are used to discretize and solve the Poisson’s equation.

The position of the free surface is then updated using VOF method expressed by (2.7). In this procedure, proper selection of the time step () affects the accuracy of final results. Therefore, time step was selected based on two stability criteria [22, 23] as Courant and diffusion conditions, respectively as follows: Using the above-mentioned relationships, the time step must be calculated to consider the smaller one in the numerical simulation.

#### 4. Definition of the New Advection Method

In solution procedure of governing equations, it needs to decide where to store the velocity components. staggered and collocated grids can be used to evaluate this problem. On staggered grids, the velocity components are stored at the cell faces, and the scalar variables such as pressure are stored at the central nodes. However, on collocated grids, all parameters are defined at the same location at the central nodes. The staggered grids method gives more accurate pressure gradient estimation. However, collocated grids method is more simpler for solving the equations [24].

In Youngs’ VOF (Y-VOF) which is a volume of fluid method, free surface is tracked based on fluxes between two neighboring cells [25, 26]. In the staggered grids, the velocity components are defined in the cell faces, and fluxes are estimated using those velocity components. Flux translation is estimated as volume of fluid passed from the faces with constant velocities of faces. While the particles of fluid between two adjacent faces have different velocities. In this paper, a new advection method based on velocity in the mass gravity of fluid in a cell is introduced as MYV method. At first, estimation is made for the interface orientation . The interface within a cell is then approximated by a straight line segment with orientation as shown in Figure 1.

Free surface cuts the cell in such a way that the fractional fluid volume is given by . The geometry of the fluid resulting from this reconstruction is then used to determine the fluxes through any side on which the velocity is directed out of the cell. For example, flux from right cell face () for the cell shown in Figure 2(a) can be estimated as Similarly, it is assumed that fluid is passed from cell face with constant velocity equal to cell face velocity. So, in the developed model, fluxes are calculated based on horizontal and vertical velocity of mass center of fluid cells. Figure 2 shows new arrangement of velocity in MYV method.

The new equations of horizontal and vertical velocity ( and ) in the mass center are estimated using equations summarized in Table 2.

For example, fluxes from cell faces in (4.1) are estimated using these new velocities as

#### 5. Model Validation

To validate the modified model, a series of standard tests such as lid-driven cavity, sloshing problem, constant unidirectional velocity field, shear test, and dam break over a dry bed was carried out.

##### 5.1. Lid-Driven Cavity

Lid-driven cavity is the fluid flow in a rectangular container which moves tangentially to itself and parallel to one of the sidewalls. Due to the simplicity of the cavity geometry, applying a numerical method on this flow problem in terms of coding is quite easy and straightforward. Despite its simple geometry, the driven cavity flow retains a rich fluid flow physics manifested by multiple counter rotating recirculation regions on the corners of the cavity depending on the Reynolds number. The boundary conditions and induced eddies are shown in Figure 3.

A sensitivity analysis is performed on the problem mesh size. For example, horizontal velocities for Re = 100 and 3200 and for different mesh sizes are shown in Figure 4. As can be seen in this figure, by increasing the mesh size, the results become independent to it.

The model is then performed for various Reynolds numbers at a range of 100 to 10000 and results were compared with those of Ghia et al. [27] which is regarded as the true solution of this problem. For example, the results for Re = 1000 and 10000, are presented in Figures 5 and 6. The existence of a good agreement between results can be seen in these figures.

##### 5.2. Sloshing Problem

Sloshing of a low-amplitude liquid wave under forced movement is another problem to test the interfacial flow solver. In this test, a rectangular tank with a width of 0.9 m and a water depth of 0.6 m was exposed to a horizontal periodic sway motion as . Therefore, was considered as exciting acceleration. The displacement of a node on the free surface in contact with right-hand sidewall was calculated by the model and compared with those of Nakayama and Washizu [28] in Figure 7. The existence of a good agreement between results can be seen in this figure.

##### 5.3. Simple Advection Test

Before attempting to couple the advection to solution of the momentum equations, a comparison of the volume-tracking schemes with analytical solution is performed to validate VOF solver. In this regard some standard tests are used.

###### 5.3.1. Constant, Unidirectional Velocity Field

The simplest test involves advection of a geometric shape in the computational domain. In this test, the geometric shape should remain intact, and total amount of the fluid within the region should be conserved. The test examined here is a hollow box being translated by a uniform constant velocity field. This is chosen to highlight the problem existing in some recent methods [13]. The 2D Cartesian region as shown in Figure 9 has dimension of and is composed equally sized cells. A square fluid blocks of dimensions moves with equal horizontal and vertical velocities of 1 m/s towards the top right-hand corner of the computational domain. Figure 8 shows the fluid position computed using the Hirt and Nichols, YV and MYV methods, every 0.1 sec until 0.7 sec. The computational time step in these calculations is leading to a courant number of 0.1.

It is evident in these figures that the geometric shape of the translated hollow square in modified model has not been improved comparing with original one. It is due to this fact that the velocity field is constant in this physical problem. Therefore, the velocity of mass center and cell faces are the same.

###### 5.3.2. Shear Test, the Rudman Vortex

The final advection test examined in this study employs a nonuniform vorticity field which stretch and shear free surface interfaces as fluid is translated through the computational domain. A 2D computational domain with dimensions of composed of uniformly sized square cells is used in the shear test. The velocity field is specified as: in which equals to 1 for first computational time step, and −1 for the second time steps. The initial fluid geometry is a circle of radius 0.2 m. Time step is selected by courant number of 0.25 based on the maximum velocity within the computational domain. A sensitivity analysis is performed on mesh sizes. The final shape of circle after 2000 steps forward and then 2000 steps backward for different mesh numbers are estimated using different VOF methods. The results are presented in Figure 9. The results show that the present methods improve the results comparing with original ones for coarse grids, but they have the same accuracy for fine grids. The results in three steps for and for YV and MYV methods are presented in Figure 10 as a sample.

##### 5.4. Dam Break over a Dry Bed

Another test problem used for free surface case is the collapse of water column over a dry bed. This problem was first studied and used as benchmark by the developers of SOLA-VOF (Nichols et al. [4]). It is a very useful benchmark providing extreme conditions to assess the numerical stability as well as the capability of the model to treat the free surface problem. In this test, a square computational domain with a length and height of 22.8 cm is set. A water column with the width of cm and height of is assumed at the left-hand side of the computational domain surrounded by walls with no-slip boundary condition. The spatial step sizes in the horizontal and vertical direction are . The comparison of the computed results using present methods with the experimental data given by Martin and Moyce [29] is shown in Figure 11. In this figure, is the location of wave front in time . It is seen that the developed models passes this test successfully as the calculated result are well agreed with experimental data.

To compare the results of different VOF algorithms, two error criteria as sum square error (SSE) and sum absolute error (SAE) are used. The results are summarized in Table 3.

#### 6. Discussion and Conclusion

In this paper, modified Youngs’ VOF algorithm denoted by MYV method is presented. In Youngs’ VOF algorithm, for staggered grids fluxes to neighboring cells are estimated based on cell face velocities. However in practice, these particles have a variable velocity between velocities of two adjacent cell faces. In the developed model, the velocity in mass center of fluid cell is estimated and used to calculate fluxes from cell faces. To validate the modified model, a serious of standard tests such as Lid-driven cavity, Sloshing problem, constant unidirectional velocity field, Shear test, and Dam break over a dry bed was carried out. The results showed that in some cases such as Constant, unidirectional velocity field, geometric shape of the translated hollow square in modified model has not been improved comparing with original one. It is due to the constant velocity field leading to equal velocities of mass center and cell faces in this physical problem. In some other cases such as Shear test, the modified model improves the results especially for coarse grids. In Dam break test, present method improves the results dramatically. Therefore, based on these results, the MYV method is suggested as a new VOF algorithm which models the free surface problems more accurately.

#### References

- B. Burnner and G. Tryggavson, “Direct Numerical Simulation of three-dimensional bubbly flow,”
*Physics of Fluids*, vol. 111, pp. 1967–1969, 1999. View at Google Scholar - B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti, “Modelling merging and fragmentation in multiphase flows with SURFER,”
*Journal of Computational Physics*, vol. 113, no. 1, pp. 134–147, 1994. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - W. F. Noh and P. Woodward, “SLIC (Simple line Interface Calculation),” in
*Lecture Notes in Physics*, A. I. van Dooren and M. J. Baines, Eds., vol. 59, pp. 273–285, Springer, New York, NY, USA, 1982. View at Google Scholar - B. D. Nichols, C. W. Hirt, and R. S. Hotchkiss, “SOLA-VOF: a solution algorithm for transient fluid flow with multiple free boundaries,” Tech. Rep. LA-8355, Los Alamos Scientific Laboratory, 1980. View at Google Scholar
- M. D. Torrey, L. D. Cloutman, and C. W. Hirt, “NASA-VOF2D: a computer program for incompressible flows with free surface,” Tech. Rep. LA-10462-MS, Los Alamos Scientific Laboratory, 1985. View at Google Scholar
- C. W. Hirt and B. D. Nichols, “Volume of fluid (VOF) method for the dynamics of free boundaries,”
*Journal of Computational Physics*, vol. 39, no. 1, pp. 201–225, 1981. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus - D. L. Youngs, “Time-depend multi-material flow with large fluid distortion,” in
*Numerical Methods for Fluid Dynamics*, K. W. Morton and M. J. Baines, Eds., pp. 273–285, Academic Press, New York, NY, USA, 1985. View at Google Scholar - M. Rudman, “Volume-tracking methods for interfacial flow calculations,”
*International Journal for Numerical Methods in Fluids*, vol. 24, no. 7, pp. 671–691, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - E. S. Duff,
*Fluid flows aspects of solidification modeling: simulation of low pressure Die casting*, Ph.D. thesis, University of Queensland, Brisbane, Australia, 1999. - W. J. Rider and D. B. Kothe, “Reconstructing volume tracking,”
*Journal of Computational Physics*, vol. 141, no. 2, pp. 112–152, 1998. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - D. J. E. Harvie and D. F. Fletcher, “A new volume of fluid advection algorithm: the defined donating region scheme,”
*International Journal for Numerical Methods in Fluids*, vol. 35, no. 2, pp. 151–172, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - D. J. E. Harvie and D. F. Fletcher, “A new volume of fluid advection algorithm: the stream scheme,”
*Journal of Computational Physics*, vol. 162, no. 1, pp. 1–32, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, and S. Zaleski, “Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows,”
*Journal of Computational Physics*, vol. 152, no. 2, pp. 423–456, 1999. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus - R. Scardovelli and S. Zaleski, “Analytical relations connecting linear interfaces and volume fractions in rectangular grids,”
*Journal of Computational Physics*, vol. 164, no. 1, pp. 228–237, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - R. Scardovelli and S. Zaleski, “Interface reconstruction with least-square fit and split Eulerian-Lagrangian advection,”
*International Journal for Numerical Methods in Fluids*, vol. 41, no. 3, pp. 251–274, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus - R. Rafei,
*Numerical solution of incompressible 3D turbulent flow in a spiral channel*, M.S. thesis, Amirkabir University of Technology, Tehran, Iran, 2004. - C. W. Li and Y. F. Zhang, “Simulation of free surface recirculating flows in semi-enclosed water bodies by a $k-w$ model,”
*Applied Mathematical Modelling*, vol. 22, no. 3, pp. 153–164, 1998. View at Publisher · View at Google Scholar · View at Scopus - H. Gao, H. Y. Gu, and L. J. Guo, “Numerical study of stratified oil-water two-phase turbulent flow in a horizontal tube,”
*International Journal of Heat and Mass Transfer*, vol. 46, no. 4, pp. 749–754, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus - B. Ren and Y. Wang, “Numerical simulation of random wave slamming on structures in the splash zone,”
*Ocean Engineering*, vol. 31, no. 5-6, pp. 547–560, 2004. View at Publisher · View at Google Scholar · View at Scopus - Y. M. Shen, C. O. Ng, and Y. H. Zheng, “Simulation of wave propagation over a submerged bar using the VOF method with a two-equation $k-\epsilon $ turbulence modeling,”
*Ocean Engineering*, vol. 31, no. 1, pp. 87–95, 2004. View at Publisher · View at Google Scholar · View at Scopus - K. Guizien and E. Barthélemy, “Short wave phase shifts by large free surface solitary waves: experiments and models,”
*Physics of Fluids*, vol. 13, no. 12, pp. 3624–3635, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus - D. P. Renouard, F. J. S. Santos, and A. M. Temperville, “Experimental study of the generation, damping, and reflexion of a solitary wave,”
*Dynamics of Atmospheres and Oceans*, vol. 9, no. 4, pp. 341–358, 1985. View at Publisher · View at Google Scholar · View at Scopus - M. J. Boussinesq, “Théorie de l’intumescence liquide, appelée onde solitaire ou de translation. se propageant dans un canal rectangulaire,”
*Comptes Rendus de l'Académie des Sciences*, pp. 755–759, 1871. View at Google Scholar - M. J. Ketabdari and H. Saghi, “Large eddy simulation of laminar and turbulent flow on collocated and staggered grids,”
*ISRN Mechanical Engineering*, vol. 2011, Article ID 809498, 13 pages, 2011. View at Publisher · View at Google Scholar - M. J. Ketabdari, M. R. H. Nobari, and M. Moradi Larmaei, “Simulation of waves group propagation and breaking in coastal zone using a Navier-Stokes solver with an improved VOF free surface treatment,”
*Applied Ocean Research*, vol. 30, no. 2, pp. 130–143, 2008. View at Publisher · View at Google Scholar · View at Scopus - M. R. H. Nobari, M. J. Ketabdari, and M. Moradi, “A modified volume of fluid advection method for uniform Cartesian grids,”
*Applied Mathematical Modelling*, vol. 33, no. 5, pp. 2298–2310, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - U. Ghia, K. N. Ghia, and C. T. Shin, “High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method,”
*Journal of Computational Physics*, vol. 48, no. 3, pp. 387–411, 1982. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus - T. Nakayama and K. Washizu, “Boundary element analysis of nonlinear sloshing problems,” in
*Developments in Boundary Element Method*, P. K. Bauerjee and S. Mukherjee, Eds., vol. 3, Elsevier Applied Science Publishers, New York, NY, USA, 1984. View at Google Scholar · View at Zentralblatt MATH - J. C. Martin and W. J. Moyce, “An experimental study of the collapse of liquid columns on a rigid horizontal plane,”
*Philosophical Transaction of the Royal Society of London*, vol. 244, pp. 312–324, 1952. View at Google Scholar