A Novel algorithm of Advection procedure in Volume of Fluid Method to model free surface flows

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.


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 and is set as one and zero in full fluid and empty cells respectively.While between these values it presents the free surface.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 Nooh 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 according to free surfaces reconstructing technique in each cell and the method of computing boundary flux integration.Some VOF methods represents free surface interfaces as a line parallel to one of the grid co-ordinates which are referred as piecewise constant scheme.Some of them are the methods used by Nichols et al. [4], Hirt et al. [6], Torry 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], Harvies and Fletcher [11,12], Geuyffier et al. [13] and Scardovelli et al. [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 research, a new advection method in FCT and YV methods as Modified Flux corrected Transport (MFCT) and Modified Youngs' VOF (MYV) methods respectively are presented.

Governing equations
The fluid is considered to be Newtonian and incompressible.Therefore, 2D continuity and Navier-Stokes Equations (NSE) are used as: where, t is time, V velocity vector, p hydrodynamic pressure,  kinematic fluid viscosity and B is body force.In the turbulent flow, the effect of turbulence can be considered using eddy viscosity models [16].Researchers have used different models such as   k and w k  [17,18,19,20] to model turbulent flow.In this paper, the standard two equation where: C ,   and k  are selected as in Table 1.
Table 1: Coefficients for standard Finally, the RANSE is used to model 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, f varies versus both time and space.Derivations of f (and its continuous spectrum F) 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 f can be represented as: Combining Eqs 11 and 12, the colour function is obtained.
where, F has a specific value in each scalar cell as: Specific technique must be used to discrete Eq. 13.

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, n U is the velocity field in old time level,U ˆ intermediate velocity field, B 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 Eq. 16, the Poisson's equation for the pressure is obtained: So, at first, the intermediate velocity field is obtained by Eq. 15.Then, the pressure distribution in the new time level is obtained using Eq. 17.In this step, the new velocity field is calculated using Eq.16.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 Eq. 13.In this procedure, proper selection of t  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 in the numerical simulation.

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 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 Fig. 1.
Fig. 1: Interface orientation in a free surface (i,j ) cell Free surface cuts the cell in such a way that the fractional fluid volume is given by F(i,j).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 ( r F ) for cell shown in Fig. 2-a can be estimated as: Similarly, it is assumed that fluid is passed from cell face with constant velocity equal to cell  in the mass center are estimated using equations summarized in Table 2.
For example, fluxes from cell faces in Eq. 20 are estimated using these new velocities as:

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.

Lid-driven cavity
Lid-driven cavity is the fluid flow in a rectangular container which moves tangentially to    A sensitivity analysis is performed on mesh size.For example, horizontal velocities for Re=100 and 3200 and for different mesh sizes are shown in Fig. 4. As 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 Figs 5 and 6.
As can be seen in these figures, there is a good agreement between the results.

Sloshing problem
Sloshing of a liquid low amplitude 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   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 Fig. 7.The existence of a good agreement between results can be seen in this figure.

Simple Advection Test
Before attempting to couple the advection F 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.

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 shown in Fig. 9   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.

Shear test, t he Rudman vortex
The final advection test examined in this study employs a non-uniform vorticity field which stretch and shear free surface interfaces as fluid is translated through the computational domain.A 2D computational domain with dimensions of uniformly sized square cells is used in the shear test.The velocity field is specified as: Which A equals 1 for first N computational time step, and -1 for the second N time steps.The initial fluid geometry is a circle of radius  The results show that the present methods improve the results comparing with original ones for coarse grids.But they have the same accurate for fine grids.The results in three steps for NX=NY=100 and for YV and MYV methods are presented in Fig. 10 as a sample.

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 L=5.7 cm and height of 2L is assumed at the left 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 Fig. 11.In this figure X is the location of wave front in time t.It is seen that the developed models passes this test successfully as the calculated result well agree with experimental data.

Discussion and conclusion
In this research, modified Volume of Fluid method based on Youngs' VOF algorithm denoted by Modified Youngs' VOF (MYV) method is presented.In this method, 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 Liddriven 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 this fact that the velocity field is constant in this physical problem.Therefore velocities of mass center and cell faces are the same.In some other cases such as Shear test, modified model improves the results.However, they have the same accuracy in fine grids.In Dam break test, present method improves the results.Therefore, the MYV method is suggested as a new VOF algorithm which models the free surface problems more accurately.
where the first equation involves turbulence kinetic energy (k) and represents the velocity scale.The second one takes into account turbulent dissipation rate    and represents the length scale.The two-equation   k turbulence model accounts for the effect of turbulence as follows: face velocity.So, in the developed model, fluxes are calculated based on horizontal and vertical velocity of mass center of fluid cells.Fig. 2 shows new arrangement of velocity in MYV method.

Fig. 2 :
Fig. 2: Definition of velocity arrangement in new advection algor ithm in MYV method 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 straight forward.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 condition and induced eddies are shown in Fig.

Fig. 7 :
Fig.7: Comparison between new models result and that of Nakayama and Washizu[28] right-hand corner of the computational domain.Fig.8shows the fluid position computed using the Hirt and Nichols, YV and MYV methods, every 0.1sec until 0.7sec.The computational time step in these calculations number of 0.1.

Fig. 8 :
Fig. 8: Numerical results of advection test for validation of VOF model: a) Exact solution;

.
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 presented in Fig 9.

Fig. 9 :X
Fig. 9: Final shape of circle after 2000 steps forward and then 2000 steps backward for

Fig. 10 :
Fig. 10: Results for shearing field using YV and MYV methods respectively and

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

Table 2 :
Velocity calculation for new arrangement in M YV method Par.

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