`ISRN Applied MathematicsVolume 2012 (2012), Article ID 521012, 16 pageshttp://dx.doi.org/10.5402/2012/521012`
Research Article

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

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

Received 4 December 2011; Accepted 19 January 2012

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 - [1720] 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.

Table 1: Coefficients for standard - turbulence model.

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.

Figure 1: Interface orientation in a free surface cell.

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.

Figure 2: Definition of velocity arrangement in new advection algorithm in MYV method.

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

Table 2: Velocity calculation for new arrangement in MYV method.

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.

Figure 3: Lid-driven cavity test: (a) boundary conditions; (b) induced eddies.

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.

Figure 4: Lid-driven cavity test; mesh independency for (a) Re = 100; (b) Re = 3200.

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.

Figure 5: Lid-driven cavity test, comparison between results of horizontal velocity of present model and those of Ghia et al. [27] for (a) Re = 1000; (b) Re = 10000.
Figure 6: Lid-driven cavity test, comparison between results of vertical velocity of present model and those of Ghia et al. [27] for (a) Re = 1000; (b) Re = 10000.
##### 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.

Figure 7: Comparison between new models result for sloshing problem and that of Nakayama and Washizu [28].

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.

Figure 8: Numerical results of advection test for validation of VOF model: (a) exact solution; (b) Hirt- Nichols; (c) YV; (d) MYV.
Figure 9: Final shape of circle after 2000 steps forward and then 2000 steps backward for different mesh sizes using YV (left) and MYV (right) methods: (a) and (b) ; (c) and (d) ; (e) and (f) .

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.

Figure 10: Results for shearing field using YV (left) and MYV(right) methods and : (a) and (b) after 1000 steps forward; (c) and (d) after 2000 steps forward; (e) and (f) after 2000 steps forward and then 2000 steps backward.
##### 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.

Figure 11: Dam breaking test, comparison of models’ results with that of Martin and Moyce.

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.

Table 3: Comparison of SSE and SAE errors of different VOF algorithms in Dam break test.

#### 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

1. B. Burnner and G. Tryggavson, “Direct Numerical Simulation of three-dimensional bubbly flow,” Physics of Fluids, vol. 111, pp. 1967–1969, 1999.
2. 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.
3. 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.
4. 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.
5. 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.
6. 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.
7. 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.
8. M. Rudman, “Volume-tracking methods for interfacial flow calculations,” International Journal for Numerical Methods in Fluids, vol. 24, no. 7, pp. 671–691, 1997.
9. E. S. Duff, Fluid flows aspects of solidification modeling: simulation of low pressure Die casting, Ph.D. thesis, University of Queensland, Brisbane, Australia, 1999.
10. W. J. Rider and D. B. Kothe, “Reconstructing volume tracking,” Journal of Computational Physics, vol. 141, no. 2, pp. 112–152, 1998.
11. 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.
12. 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.
13. 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.
14. 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.
15. 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.
16. R. Rafei, Numerical solution of incompressible 3D turbulent flow in a spiral channel, M.S. thesis, Amirkabir University of Technology, Tehran, Iran, 2004.
17. 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.
18. 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.
19. 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.
20. 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.
21. 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.
22. 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.
23. 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.
24. 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.
25. 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.
26. 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.
27. 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.
28. 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.
29. 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.