Chapter 15

Free Surface Flow Simulation Using VOF Method

Mohammad Javad Ketabdari

Additional information is available at the end of the chapter

[Link]

Abstract
Viscous flow with moving free surface is an important phenomenon in nature which
has broad applications in engineering. For these flows, temporal and spatial position of
this moving free surface in unsteady or non‐uniform conditions is very complicated. In
this chapter, free surface simulation methods based on computational grid are
presented. Volume of fluid (VOF) is a powerful and the most prevailing method for
modeling two immiscible incompressible fluid‐fluid interfaces. Herein, the governing
equations of fluid flow including Navier‐Stokes coupled with VOF equation are
discussed and the most prominent VOF schemes hierarchically presented to the readers.
Meanwhile, Compressive Interface Capturing Scheme for Arbitrary Meshes (CIC‐
SAM), Higher Resolution Artificial Compressive (HiRAC), High Resolution Interface
Capturing (HRIC), Switching Technique for Advection and Capturing of Surfaces
(STACS), and some other newly proposed methods are introduced, and the accuracy
and time calculation of each method are evaluated. Moreover, surface tension modeling
and its discretization as one of the most demanding phenomena in the nature are
brought to the readers. Finally, two schemes of parametric study of interfaces are
discussed.

Keywords: viscous flow, free surface modeling, Eulerian approach, volume of fluid,
interface simulation

1. Introduction

Simulations of free surface flows have progressed rapidly over the last decade, and it is now
possible to simulate the motion of complicated waves and their interactions with structures
considering even deformable bubbles in turbulent flows. In the continuum mechanics, there are
two methods to express the motion in the environment. The first description is the Eulerian
approach. In this method, attention is paid to a special volume in the space. A mesh remains

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons
Attribution License ([Link] which permits unrestricted use, distribution,
and reproduction in any medium, provided the original work is properly cited.
366 Numerical Simulation - From Brain Imaging to Turbulent Flows

fixed in the Eulerian method and fluid regions change in shape and location on the mesh. It uses
a fixed grid system which is not transformed during the solution procedure. The fluid is studied
while passing this volume and continuously replaced in time. Therefore, this method is not
appropriate for formulation of basic equations of fluid movement. The Eulerian method has
some limitations. For example, when the portion of the perimeter to the area of a zone of fluid
is large, the error of this method is increased. In the Eulerian method, it is not possible to
decompose the equation on the boundaries with the same precision of inner region of fluid and
accordingly, the finer mesh should be used near the boundaries. Therefore, when the free surface
of a discontinuous region is modeled by this method, finer grid should be employed in order to
achieve more precise results, specifically if this surface has large deflections. This is crucial when
the portion of the area to the perimeter of a zone is low, for example on phase of a multiphase
fluid. In this case, using finer mesh could increase the portion of the number of the inner elements
to the boundary elements, which in turn, increases the precision of the numerical solution. The
main superiority of the Eulerian description is the possibility of modeling of complicated
surfaces. For example, the collapse of a column of a fluid could be modeled in the Eulerian grid
which is shown in Figure 1.

Figure 1. Fluid column in Eulerian grid: (a) before collapse and (b) collapsing flow.

Figure 2. A sample of Lagrangian grid in vertical direction.

 Free Surface Flow Simulation Using VOF Method 367
[Link]

In Lagrangian method, the flow field of the considered fluid is covered by a mesh moving with
the fluid. The fluid boundaries always coincide with the grid boundaries and the fluid inside
each cell of the grid always remains in that computational cell. Although this method is not
applicable to flows undergoing large distortions, where meshes can be twisted into unaccept‐
able shapes, but its advantage is the ease with which it handles free surfaces and interfaces,
which makes it applicable to a wide variety of problems. For example, the grid shown in
Figure 2 is Lagrangian in the vertical coordinate. For free surface problems, if the free surface
movement or the tangential acceleration gradient in the perpendicular direction to its surface
is not large, the Lagrangian method can be used to simulate free surfaces. The grid lines are
located on the free surface and move with it. Therefore, there is no need for any special
boundary condition in this location [1].

2. Governing equations

Governing equations for a compressible viscous fluid flow with no phase change are as follows:

¶r
+ Ñ × ( ru ) = 0 (1)
¶t

¶ ( ru )
+ Ñ × ( r uu ) + ÑP = r g + Fs + Ñ × éë m ( Ñu + Ñu T ) ùû (2)
¶t

In these equations, ρ , u, t , P , μ , g = (gx , gy ), and Fs are density, velocity vector, time, total
pressure, kinematic viscosity, gravity acceleration, and body forces, respectively. Body forces
include forces due to surface tension in the interface. Here, properties of a fluid such as density
and viscosity are included in the equations. However, it should be kept in mind that the
information changes from one fluid to another. Thus for mesh‐based numerical methods, new
properties based on fluid properties of both materials should be considered for Eq. (2) in the
cell containing the free surface, and the governing equations should be rewritten in the
following form:

¶
(ar1 ) + Ñ × (ar1u1 ) = 0 (3)
¶t

¶
¶t
( (1 - a ) r2 ) + Ñ × ( (1 - a ) r2u 2 ) = 0 (4)

¶ (ar1u1 )
+ Ñ × (ar1u1u1 ) + ÑP1 = ar1g + Fs + Ñ × éëam1 ( Ñu1 + Ñu1T ) ùû (5)
¶t
368 Numerical Simulation - From Brain Imaging to Turbulent Flows

¶ ( (1 - a ) r 2u 2 )
+ Ñ × ( (1 - a ) r 2u 2u 2 ) + ÑP2 = (1 - a ) r 2g + Fs +
¶t (6)
Ñ × éë(1 - a ) m 2 ( Ñu 2 + Ñu 2 T ) ùû

where indices 1 and 2 show first and second fluids properties, and α is a scalar phase indica‐
tor function which is defined as follows:

ì1 Control volume is filled only with phase 1

ï
a = í0 Control volume is filled only with phase 2 (7)
ï0 < a < 1 Interface present
î

This phase indicator function is the fluid property or volume fraction, which moves with it
and can be derived as follows:

¶a ¶U ia
+ =0 (8)
¶t ¶xi

This function can be used to calculate the fluid properties in each phase as a weight function.
In order to use a set of governing equations using the weight function, each fluid property
should be calculated based on the volume occupied by this fluid in the surface cell as ex‐
pressed in Eqs. (9) and (10) [2]:

r = a × r1 + (1 - a ) × r 2 (9)

m = a × m1 + (1 - a ) × m 2 (10)

Free surfaces considered here are those on which discontinuities exist in one or more varia‐
bles. This has been the challenge for researchers to omit or reduce this problem as much as
possible. The transient state as well as phenomena such as surface tension, changing of fluid
phase and Kelvin‐Helmholtz instability makes numerical simulation of such problems
cumbersome. It is expected that methods used to simulate interface of fluids have a number
of characteristics. These include mass conservation, simulating the interface as thin as possible,
 Free Surface Flow Simulation Using VOF Method 369
[Link]

being able to reproduce complicated topologies, generalization of expansion to 3D problems,

and being able to model surface phenomena and be computationally efficient.

3. Free surface modeling methods

There are different methods to simulate free surface flow, each of which has its own advan‐
tages and disadvantages:

3.1. Donor‐acceptor method

The main idea of donor‐acceptor approach is that the value of volume fraction in downwind
cell, the acceptor cell, is used for anticipation of transferring fluid in each time step. The
problem in this approach is that using downwind cell in calculations may lead to unreal
situations which are values out of zero and unity domain in surface cells. Figure 3a shows this
method with the first fluid with gray color and volume of fluid equals to unity. It could be
seen that using donor‐acceptor approach with downwind differencing scheme results in
values greater than unity in donor cell. It is because the second fluid in the acceptor cell is
greater than the value needed in the donor cell. Similarly in Figure 3b, using downwind
differencing scheme leads to negative values for volume of fluid, which is because the needed
fluid in acceptor cell is more than what is in the donor cell [3].

Figure 3. Schematic view of donor‐acceptor approach [4].

In order to be assured that volume of fluid is between zero and one, the amount of fluid or
volume of fluid in donor cell should be used to regulate the estimated fluid transferring
between two adjacent cells [5].

One drawback of donor‐acceptor method is that this method changes any finite gradient into
step, and consequently increases the slope of the surface model in the direction of flow. This
problem was alleviated by proposing a method to consider the slope of interface for flux
transferring in adjacent cells by Hirt and Nichols [6]. For this purpose, a donor‐acceptor
equation was proposed so that it could detect the direction of the flow in interface and then
define the upwind and downwind cells accordingly. Thereafter, this model was expanded
370 Numerical Simulation - From Brain Imaging to Turbulent Flows

for 3D domains by Torrey et al. [7]. The Surfer method is one version of volume of fluid which
deals with merging and fragmenting of interfaces in multiphase flows [8].

The volume of fluid method is one of the most popular methods for anticipation of interfa‐
ces, and many researches have been conducted based on this method including dam break,
Rayleigh‐Taylor instability, wave generation and bubble movement [6, 9–12]. This method was
modified in 2008 to get more accurate results by considering diagonal changes in fluxes of
adjacent cells for structured grid domains [13, 14].

3.2. The Hirt‐Nichols method

The volume of fluid (VOF) method was first proposed by Hirt and Nichols [6]. In this method,
similar to the SLIC method, free surfaces can be reconstructed based on parallel lines with
respect to one of the principal coordinates of the system. However, nine neighboring cells are
considered for flux changes and defining the normal vector in a desired cell. Then, free surface
is considered as either a horizontal or a vertical line in cell with respect to the relative normal
vector components. Figure 4 shows the actual free surface and what was simulated by Hirt‐
Nichols method.

Figure 4. Free surface (a) actual surface and (b) reconstructed surface based on Hirt‐Nichols method [6].

Figure 5. Hirt‐Nichols scheme (a) actual surface and (b) reconstructed surface.
 Free Surface Flow Simulation Using VOF Method 371
[Link]

Upwind fluxes are used for fluxes parallel to the reconstructed interface, while donor‐acceptor
fluxes are used for those fluxes normal to it. For instance according to Figure 5a, the inter‐
face in the cell (i, j ) is considered to have positive celerity direction with respect to x coordi‐
nate in the face i + 1 2 in donor‐acceptor method. Therefore, the reconstructed surface in the (i, j )
cell is vertical (Figure 5b), and this cell is considered a downwind cell for the cell (i + 1, j ),
(αi, j > αi+1, j ). According to donor‐acceptor method, transferred flux from the face (i + 1 2, j ) can
be calculated as follows:

{ ( )}
Fi + 1 2 , j = d y min éa i , jd x,U i + 1 2 , ja i +1, jd t + max 0,U i + 1 2 , j (1 - a i +1, j ) d t - (1 - a i , j ) d x ù
ë û (11)

where αi, j δx is the maximum fluid available for exiting the cell (i, j ); U α δt is the esti‐
i+ , j i+1, j
1

mation of downwind flux for the volume of the fluid, α ; U 1

i+ , j
(1 − αi+1, j )δt is the estimation of
2

downwind flux resulted from the convey of void portion of the cell (i, j ); and (1 − αi, j )δx is
the maximum void portion which can exit from the cell (i, j ).

The “min” operator has been designed to ensure the fluid leaving the cell (i, j ) is not more
than the calculated available fluid in it from the previous time step. As the fluid in a cell
transfers, so does the whole void space in the cell. Thus, the “max” operator has been designed
in order to assure that amount of void exits the cell is bounded by what was in it calculated
from the previous step. In this scheme, the combination of downwind and upwind fluxes has
been considered in such a way that not only the solution stability is guaranteed, but also avoids
the numerical diffusion.

3.3. Flux Corrected Transport (FCT) method

The FCT method is based on the idea to present a formulation which combines the upwind
and downwind fluxes. This formulation aimed to leave out upwind numerical diffusion and
instability of downwind scheme [15]. Idea of neighboring fluxes based on higher order
translate scheme was first proposed by Boris and Book [16] and then developed by Zalesak [17]
to multidimensional.

In this method calculations consist of some steps. First, an intermediate value of volume of
fluid, α *, must be defined based on a first‐order scheme. Figure 6 shows schematically the
solution for a 1D governing equation of fluid volume fraction for cell i as:

a i* = a in -
1
(
F L 1 - Fi -L 1 2
d x i+ 2
) (12)
372 Numerical Simulation - From Brain Imaging to Turbulent Flows

Figure 6. Three adjacent typical cells in FCT method.

where F L is the flux in the downwind cell, and F f is defined in the face f as:

Ff = U f d ta f (13)

Thereafter, an anti‐diffusive flux is needed to be defined (FL) in order to correct the diffusion
of the previous step. This is the difference between upwind and downwind fluxes as:

Fi +A 1 2 = Fi +H1 2 - Fi +L 1 2 (14)

To make this stable, a correction factor, q , is needed to modify the fluxes values. Finally, a value
for fluid fraction in next time step is defined as:

a in +1 = a i* -
(q i + 12 Fi +A 1 2 - qi - 1 2 Fi -A1 2 ) (15)
dx

3.4. Youngs’ method

This method was first proposed by Youngs in 1982 [18]. It was then developed by Rudman [19]
with more details. In this method, at first the slope of the interface position is estimated. Then,
the free surface is defined as a straight line with the slope of β in each cell of the numerical
domain. The position of this line segment in each cell is defined such that the area reconstruct‐
ed from the line and the perimeter of the cell is equal to the amount of volume of fluid, α . The
geometry of the polygon from this reconstruction is used to calculate the flux transferred from
the cell faces.
Assuming that αi, j is predefined in every cell, the first step is to calculate first‐order upwind
fluxes. Then, the Youngs exiting fluxes of every cell can be calculated by considering the values
in each cell. To do so, the angle β , between free surface and x‐coordinate, must be calculated.
Different methods can be used for calculation of β . One of them is first using the gradient
 Free Surface Flow Simulation Using VOF Method 373
[Link]

function of fluid volume fraction for defining unit normal vector of the free surface, and then
calculating β [20]. The method of defining the normal vector, however, can affect the accura‐
cy of the final results. This formulation for uniform grid is as follows:

Figure 7. Four possible positions for free surface in Youngs’ method [19].

a i +1, j +1 + 2a i +1, j + a i +1, j -1 - a i -1, j +1 - 2a i -1, j - a i -1, j -1

nix, j = (16)
dx

a i +1, j +1 + 2a i , j +1 + a i -1, j +1 - a i +1, j -1 - 2a i , j -1 - a i -1, j -1

niy, j = (17)
dy

Using components of the normal unit vector, the angle β can be calculated as:

æ -n ö
b = tan -1 ç x ÷
ç ny ÷ ( -p < b £ p ) (18)
è ø

The angle γ is also defined as:

ædx
è d y
ö
g = tan -1 ç tan b ÷
ø
(0 £ g £ p 2 ) (19)

It is possible to set 0 ≤ γ ≤ 90 by rotating the cell. Therefore, there are only four possible
positions for the free surface, which are depicted in Figure 7.

What is behind this conclusion is as follows:

374 Numerical Simulation - From Brain Imaging to Turbulent Flows

if g < p 4
if a £ 12 tan g
Case I
elseif a £ (1 - 12 tan g )
Case II
else
Case IV
else
if a £ 12 cot g
Case I
elseif a £ (1 - 12 cot g )
Case III
else
Case IV

Case I Case II

st 0 0

sr 2αtanγ
1
α + 2 tanγ

sb 2αcotγ 1

sl 0 1
α − 2 tanγ
if U t > 0 if U t δt ≤ (1 − sr )δy if U t δt ≤ (1 − sr )δy
Ft = 0 Ft = 0
else elseif U t δt ≤ (1 − sl )δy
U t δt − (1 − sr )δy
1 2cotγ
Ft = 2 Ft =
1
2 U t δt − (1 − sr )δy 2cotγ
else
F t = U t δtδx − (1 − α )δxδy

if U r > 0 if U r δt ≥ sbδx F r = U r δt (sr δy − 2 U r δttanγ )

1

F r = αδxδy
else
( )
F r = 2 U r δt 2 − U r δt s δx sr δy
1

b
 Free Surface Flow Simulation Using VOF Method 375
[Link]

Case I Case II

if U b > 0 if U bδt ≥ sr δy if U bδt ≤ sl δy

F b = αδxδy F b = U bδtδx
else elseif U bδt ≤ sr δy
(
F b = 2 U bδt 2 − U bδt s δy sbδx
1

r
) F b = U bδtδx −
1
2 (U bδt − sl δy )2cotγ
else
F b = αδxδy

if U l > 0 if U l δt ≤ (1 − sb)δx F l = U l δt (sl δy + 2 U l δttanγ )

1

Fl = 0
else
U l δt − (1 − sb)δx 2tanγ
1
Fl = 2

Case III Case IV

st 1
α − cotγ
2 1 − 2αcotγ
sr 1 1

sb 1
α + 2 cotγ 1

sl 0
1 − 2αtanγ
if U t > 0 F t = U t δt (st δx + 2 U t δtcotγ ) if U t δt ≥ (1 − sl )δy
1

F t = U t δtδx − (1 − α )δxδy
else
F t = U t δt (st δx + 2 U t δtcotγ )
1

if U r > 0 if U r δt ≤ st δx if U r δt ≤ st δx
F r = U r δtδy F r = U r δtδy
elseif U r δt ≤ sbδx else
F r = U r δtδy − 2 tanγ (U r δt − st δx )2
1
(U r δt − st δx )2tanγ
1
F r = U r δtδy − 2

else
F r = αδxδy

if U b > 0 F b = U bδt (sbδx − 2 U bδtcotγ ) if U bδt ≤ sl δy

1

F b = U bδtδx
else
(U bδt − sl δy )2cotγ
1
F b = U bδtδx − 2
376 Numerical Simulation - From Brain Imaging to Turbulent Flows

Case III Case IV

if U l > 0 if U l δt ≤ sbδx if U l δt ≥ (1 − st )δx

Fl = 0 F l = U l δtδy − (1 − α )δxδy
elseif U l δt ≤ st δx else
F l = U l δt (sl δy + 2 U l δttanγ )
1
U l δt − (1 − sb)δx 2tanγ
1
Fl = 2

else
F l = U l δtδy − (1 − α )δxδy

Table 1. Calculation of exiting flux in Youngs’ method.

Four side fractions (sl , sb, sr , st ) for up, right, down, and left faces can be calculated with the
selection of the free surface position in a cell. Thereafter, flow fluxes can be geometrically
computed for each face (F l , F b, F r , F t ) based on these side fractions. More details are present‐
ed in Table 1. In this table, positive value is set for velocities towards the outer edges of a cell,
and there is no flux calculation for negative velocities into the cell.

3.5. Piecewise Linear Interface Calculation (PLIC) method

To solve fluid volume transfer equation with FDM or FVM, diffusion error in interface
reconstruction occurs. This leads to poor modeling of free surfaces, specifically in the inter‐
face of two adjacent fluids with large density difference. PLIC is one of the methods to
reconstruct the interface between fluids with second‐order accuracy [20]. It can increase the
accuracy of transferred flux estimation and geometric fluid distribution in each cell. In this
method, unit normal vectors of the surface are calculated based on the volume fraction of fluid
using Youngs’ least square method as:

æ ö
r ç Ña i , j ÷
ni , j = - ç ÷ (20)
( Ña ) + ( Ña )
2 2
ç ÷
è x ,i , j y ,i , j ø

where (∇ α ) is defined as:

1é (21)
( Ña ) i , j = ( Ña )i - 12 , j - 12 + ( Ña )i - 12 , j + 12 + ( Ña )i + 12 , j + 12 + ( Ña )i + 12 , j - 12 ùû
4ë
 Free Surface Flow Simulation Using VOF Method 377
[Link]

1 æ a +1, - a i , j a i +1, j +1 - a i , j +1 ö r
( Ña ) i + = ç i j + ÷i +
2,
2 çè xi +1, j - xi , j
j + 12
xi +1, j +1 - xi , j +1 ÷ø
1

(22)
1 æ a i , j +1 - a i , j a i +1, j +1 - a i +1, j ö r
ç + ÷j
2 çè yi , j +1 - yi , j yi +1, j +1 - yi +1, j ÷ø

Figure 8. Different positioning of the interface for (0 ≤ θ ≤ π / 4).

where θ is the angle between the normal vector and the horizontal coordinate varies be‐
tween zero and 2π . For θ in the first one‐eighth space (0 ≤ θ ≤ π / 4), eight different conditions
are possible for the position of free surfaces as illustrated in Figure 8 [21]. All other situa‐
tions can be achieved with a mirror reflection of the first quarter with respect to the x and y
axes and bisectors between them. The exact position of the free surface is determined defin‐
ing surface unit normal vector using volume fraction of fluid in each cell. To do this, ex‐
treme values of αlim,1 and αlim,2 are determined as:

nmin
a lim,1 = , a lim,2 = 1 - a lim,1 (23)
2nmax

in which v (i + 1, j − 1 2) and nmax are calculated as:

378 Numerical Simulation - From Brain Imaging to Turbulent Flows

(
nmin = min nx , n y ) (
, nmax = max nx , n y ) (24)

αlim,1 and αlim,2 are shown in Figure 9.

When unit normal vector of a surface is defined, the true position of the interface can be easily
determined using volume of fluid in each cell.

Figure 9. Various positions of an interface in a cell.

3.6. Higher order differencing schemes

Another method to reconstruct the interface between two fluids is to discretize the convec‐
tion term using higher order differencing schemes or blended differencing scheme. The
accuracy of less/non‐diffusive schemes and compressive schemes was compared by Davis [22].
Less/non‐diffusive schemes prevent the interface profile from being diffused. Compressive
schemes not only prevent the interface from being diffused, but also omit any diffusion in the
neighboring of the interface. Thus, they are considered as powerful tools for thin interface
simulation.

Ghobadian [23] applied the higher order scheme proposed by Van Leer [24]. However, his
results showed that this scheme has poor ability in terms of removing diffusion. Therefore, he
proposed solutions for decreasing numerical diffusion. Other methods for omitting diffu‐
sion proposed by Pericleous and Chen [25] proved to be associated with interface diffusion.
Although first‐order upwind or downwind schemes lead to diffusion, higher order methods
result in numerical fluctuations in the interface. There are other methods for reducing the
interface as follows:

3.6.1. Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) scheme

The CICSAM scheme, presented by Ubbink, is a combined method to reduce the diffusion
problem in interface modeling. This method imposes some limitations on the fluid fraction
value. It is obvious that the value of a fluid in a cell should be constant in the absence of a
source. The CICSAM approach presents an equation for free surface volume fraction as:
 Free Surface Flow Simulation Using VOF Method 379
[Link]

a% f = g f a% f + (1 - g )a% f (25)
CBC
f UQ

where α̃ f CBC is an index for transport bound which defines a bound as follows:

ì ïì a% üï
ïmin í1, D ý when 0 £ a% D £ 1
a% f CBC
=í îï c f þï (26)
ï
îa% D when a% D £ 0, a% D >1

where α̃ D is defined based on the following equation:

a D - aU
a% D = (27)
a A - aU

ì ì 8c f a% D + (1 - c f )(6a% D + 3) ü
ïïmin í ,a% fCBC ý when 0£ a% D £1
a% f
UQ
=í î 8 þ (28)
ïîa% D when a% D < 0, a% D >1
ï

More details on determining α̃ f can be derived based on the combined schemes such as
Normalized Variable Diagram (NVD) and discussed by Ubbink and Issa [26]. Accordingly, a
local boundedness is defined for α̃ f UQ . In Eq. (25), γf is a weight factor which is defined as:

ì cos(2q f ) + 1 ü
g f = min íkg , 1ý (29)
î 2 þ

where kγ is a constant, usually set to unity, and θf is the angle between free surface normal
⇀
vector, (∇ α) D , and the vector d f such that it connects the center of the adjacent cells based on
Figure 10, and defined as:

r
( Ña ) D . d f
q f = arccos r (30)
( Ña ) D d f
380 Numerical Simulation - From Brain Imaging to Turbulent Flows

The CICSAM method well satisfies the bounds defined within it, and can be accurately
reconstruct the free surface. The basis of the method, however, is on the 1D equations and
linearization, which makes it less accurate for 3D modeling reconstructions.

Figure 10. Calculation of the upwind value for an arbitrary mesh.

3.6.2 THOR scheme

This scheme is based on the CICSAM and switches smoothly between the upper bound of the
universal limiter and ULTIMATE‐QUICK, a combination of the universal limiter and QUICK,
considering the angle between the interface and the direction of motion [27].

Analogous to CICSAM, this scheme is an algebraic advection scheme for the interface, which
is designed for the implicit time advancing algorithm. In this method αf is calculated using a
weighting factor βf as follows:

a f = (1 - b f )a D + b f a A (31)

where

a% f - a% D
bf = (32)
1 - a% D

α̃f and α̃D can be calculated by Eqs. (25) and (27), respectively.

3.6.3. Higher Resolution Artificial Compressive (HiRAC) scheme

HiRAC scheme is another modification of the CICSAM method [28]. This newly proposed
method tries to improve the computational efficiency and maintain the accuracy. In this
method, the weighing factor, γf , of Eq. (29) can be redefined as:
 Free Surface Flow Simulation Using VOF Method 381
[Link]

(
g f = min (q f ) ,1
m
) (33)

where θf was previously defined in Eq. (30). For m=2, the new formulation reduces to the
weighting function of Ubbink and Issa [26]. As m increases, the interpolation becomes more
biased towards the diffusive higher resolution scheme. It is shown that m=2 provides a good
balance between the compressive and diffusive higher resolution schemes.

3.6.4. High Resolution Interface Capturing (HRIC) scheme

This method is somehow similar to CICSAM, which benefits from a combined interpolation
scheme. In HRIC method, the difference between two upwind schemes is calculated based on
the normal vector angle of the free surface as [29]:

ìa% D a% D < 0
ï2a% 0 < a% D < 0.5
ï
a% **f = í D (34)
ï1 0.5 < a% D < 1
ïîa% D 1 < a% D

ìa% **f c f < 0.3

ï
ï 0.7 - c f
a% *f = ía% D + (a% **f - a% D ) 0.3 < c f < 0.7 (35)
ï 0.4
ïa% D 0.7 < c f
î

The portion of each of the two terms in the above equations can be defined as:

a% f = a% *f cosq + a% D (1- cosq ) (36)

In this way, αf is discretized based on the neighboring cells.

It should be noted that an improved scheme of HRIC, called Flux-Blending Interface-Capturing

Scheme (FBICS), has been recently proposed. In this method, analogous to CICSAM and HRIC,
382 Numerical Simulation - From Brain Imaging to Turbulent Flows

the difference between two upwind schemes is calculated based on the normal vector angle
of the free surface. Based on FBICS, Eq. (33) can be reformulated to obtain a more accurate
scheme as:

ì 3a% Da% D £ 1
ï 8
ï 11
ïa% D + < a% D < 3
a% *f , FBICS =í 4 8 4
(37)
ï 1 < a% D < 1
3
ï 4
ï a% Da% D £ 0 or a% D ³ 1
î

Some other modifications are also proposed by Tsui et al. [30].

3.6.5. Switching Technique for Advection and Capturing of Surfaces (STACS) method

One of the drawbacks of HRIC and CICSAM schemes is high Courant numbers. Both methods
lack a proper switching strategy to accurately model the interface when Courant number
increases. The Courant number, Cn, in HRIC method can be written as follows:

0.7 - Cn
a% f = a% *f + (a% *f - a% D ) (38)
0.7 - 0.3

which is suitable for Courant numbers between 0.3 and 0.7. For a Cn below 0.3 the scheme is
not modified, while for a Courant number above 0.7 the upwind scheme is used. This is true
for CICSAM when the Cn is equal to unity.

STACS method has been proposed to improve the accuracy and stability of the results
specifically in high Courant numbers by Darwish and Moukalled [31]. It uses an implicit
transient discretization, i.e. no transient bounding is applied, and in order to minimize the
stepping behavior of HRIC scheme, a modification is proposed. In this method, applying cos4θ
term is designed instead of cosθ in Eq. (36) as follows:

a% f = a% *f , sup cos 4q + a% *f , stoic (1 - cos 4q ) (39)

where α̃ * f , sup and α̃ * f , stoic are calculated as follows:

 Free Surface Flow Simulation Using VOF Method 383
[Link]

ì a% Da% D £ 0
ï1 1
ï + a% D 0 < a% D < 1
ì a% Da% D £ 0 ï2 2 2
ï ïï 3 3 1 5
a% *
= í10 < a% D < 1 a% *f , stoic = í + a% D < a% D <
(40)
f , sup
ï a% a% ³ 1 ï8 4 2 6
î D D ï 5
ï 1 < a% D < 1
ï 6
îï a% Da% D ³ 1

This enables a rapid but smooth switching strategy that works very well, especially where the
normal to the free‐surface face is not along the grid direction.

3.6.6. Inter‐gamma scheme

In this method, presented by Jasak and Weller [32], free surface compression is modeled using
additional compressive artificial terms.

¶g
+ Ñ × ( Ug ) + Ñ × ( U rg (1 - g ) ) = 0 (41)
¶t

where Ur is a velocity field for compressing the free surface. This artificial term is activated
only in the presence of free surface because of having the term γ (1 − γ ). The solution to this
equation is bounded from zero to unity with the inter‐gamma scheme. Eq. (40) can be rewritten
as:

¶g
¶t
( ) (
+ Ñ × f [g ] f (f , S ) + Ñ × frb [g ] f (f , S ) = 0
rb
) (42)

where ϕ = S ⋅ U f is the volume flux and ϕrb = (1 − γ ) f (−ϕ r ,S )ϕ r . Eq. (43) is used to calculate ∅r as:

ìï n*f üï
f r = K c n* max í 2 ý (43)
ïî S ïþ

where n * is the free surface normal unit vector and K c is an adjustable coefficient, with an
appropriate value of 1.5, which defines the compression rate of free surface. The inter‐gamma
scheme is also similar to CICSAM and is based on a donor‐acceptor equation and normal‐
ized variable diagram, NVD. The formulation is as:
384 Numerical Simulation - From Brain Imaging to Turbulent Flows

ìa% D a% D < 0, a% D > 1

ï
a f = í-2a D + 3a D 0 < a% D < 1 2
% % 2
% (44)
ï1 1 <a %D < 1
î 2

Figure 11 shows the NVD for inter‐gamma scheme.

Figure 11. NVD for inter‐gamma scheme.

3.7. Integrated methods

As mentioned before, volume of fluid is among the most popular methods in free surface
modeling. Having in mind that this method is based on defining a discontinuous function, the
color function, there is not a unique form for free surface. Therefore, it is required to recon‐
struct the free surface using volume of fraction function. In one hand, VOF method satisfies
the conservation of mass while it is unable to calculate free surface parameters including
curvature radius and normal unit vector directly. On the other hand, in level set methods as
the distance function is smooth, the surface geometry can be easily calculated, while satisfy‐
ing the conservation of mass is very demanding. In order to resolve the problems of level set
methods, a number of different researches have been conducted. For example, higher order
schemes were proposed to improve the conservation of continuity equation by Peng et al. [33].
Adaptive mesh refinement techniques were also proposed to increase the accuracy of the local
mesh consistency. In 2009, an integrated method known as hybrid Particle Level Set (PLS) was
proposed to improve the accuracy of the results. However, the problem still remained in
relation with mass conservation.
In order to take the advantages of both methods and eliminate their disadvantages, integra‐
tion of volume of fluid and level set methods was proposed in a new scheme known as coupled
 Free Surface Flow Simulation Using VOF Method 385
[Link]

level set and volume-of-fluid (CLSVOF) method to model two‐phase incompressible flows by
Sussman and Pucket [34]. It should be noted that although accurate, this method cannot be
easily employed, because these two methods, VOF and level set, should be individually solved
and their effects need to be coupled based on the reconstructed interface.

4. Calculating surface tension

Defining the pressure difference inserted on the surface of two fluids with different densities
and tension stresses is one of the most demanding problems in fluid mechanics. One meth‐
od to do this is the Pressure Calculation based on the Interface Location (PCIL) method which
is presented here. Surface tension, that changes the value of variables in momentum equa‐
tions, imposes a discontinuity at the position of the interface between two fluids [21].

Stress from surface tension inserts a force upon the interface. The resultant force is perpen‐
dicular to the surface and its curvature is dependent on the geometry of the surface. Surface
tension can be considered in two ways. In the first approach, it is considered as a boundary
condition in the equations for the surface. This needs using an iterative method for true
approximation of pressure, which in result, increases the time and cost of calculation and
consequently makes it inefficient. In order to address this problem, some other methods have
been proposed in which the precise calculation of interface position is not necessary. In these
methods, the direct force of surface tension has been replaced with the body force in the
momentum equation. The Continuum Surface Force (CSF) method is a base method for
calculation of body forces of fluid surface tension [2]. The body forces can be considered to act
smoothly on a narrow strip of cells in interface zone. In this method the surface stresses are
replaced with the body forces which are calculated as:

Fs = s ò k nd ( x - x s )dS (45)
S

where σ , κ , n, δ (x − xs ) and Xs are surface tension force, curvature of the surface, normal unit
vector, Dirac delta function, and the position of a given base point on the interface S , respec‐
tively. This equation has been discretized for numerical methods of two‐phase fluids.

Another approach based on CSF method was proposed by Torrey et al. [7] called Contin‐
uum Surface Stress (CSS), in which body forces of CSF method were replaced by tension
tensors of surface tension based on the following equation:

Fs = sk nd s = -Ñ × T, T = -s ( I - n Ä n ) d s (46)

where I and T are unit tensor and tangential tension tensor of the interface respectively.

It should be mentioned that employing CSF or CSS methods has some drawbacks. For instance,
spurious velocities of the thinner fluid near the interface is one the reported problems.
386 Numerical Simulation - From Brain Imaging to Turbulent Flows

A number of researches have been performed in order to resolve the problem of spurious
velocities [35, 36]. In one approach, using virtual particles moving along with the surface could
improve the results [37]. One of the latest methods presented in this field is PCIL. This method
shows that having more precise border cells and calculating their associated pressure based
on the momentum equation can lead to significant reduction of bothersome flows near this
region. PCIL is a simple and efficient method of calculating free surfaces. The total pressure
on the left side of the cell can be calculated as (see Figure 12):

Figure 12. Cells on the interface in contact with left and right cells of the free surface [38].

lL Dy - lL
PL = P1L + P2 L = P1L H L + P2 L (1 - H L ) = P2 L + H L ( P1L - P2 L ) (47)
Dy Dy

where PL is the mean pressure on the left side of the given cell. In the same way, mean pressure
can be calculated for other faces of the cell as follows:

PR = P2 R + H R ( P1R - P2 R ) (48)

PT = P2T + H T ( P1T - P2T ) (49)

PB = P2 B + H B ( P1B - P2 B ) (50)

where H is a dimensionless number which shows the position of free surface in different
directions. For faces completely immersed in the main fluid, H is equal to 1, and for those
completely in the secondary fluid, H is zero. For other cases, H varies between zero and unity.
Note that for 3D models, it is just needed to replace edge faces with surfaces faces.
 Free Surface Flow Simulation Using VOF Method 387
[Link]

On the other hand, the change in pressure Ps in every point of the interface is calculated as:

Ps = P1 - P2 = sk (51)

Accordingly, the above equation can be reformulated for pressure in the Kth face of every
common cell as follows:

Pk = P2 k + H ksk (52)

where the second term introduces the normal force of the surface tension per unit area of the
interface. This can be presented in the vector form as:

Fs = H sk n (53)

where Fs is the surface tension force vector.

One of the most fundamental steps to perform surface tension calculations is defining the
curvature of the interface. Defining this curvature is not so demanding as long as the precise
position of the interface is known. However, using volumetric tracing methods and equiva‐
lent alternatives representing the interface position make the estimation of the curvature
cumbersome.

The method of volume of fluid presented by Hirt and Nichols [6] is one of the earliest methods
in this field. In this method, a curve y ( x ) is fitted to the nine neighboring cells of the inter‐
face. In this way, summation of the volume of fluid of three cells located in a column creates
a value for y . By fitting parabolic curve to the values of three neighboring columns and then
two times differentiating of these values, one can define the curvature of the surface. This
method, however, suffers from low accuracy and some limiting conditions. Another method
is presented by Chorin [39] in which a circle is defined based on trial and error in order to
satisfy the nine‐cell set in the best possible way. Thereafter, curvature of the surface can be
defined on the basis of the calculated circle. However the main problem of this method is its
dependency on the several times of trials and errors in order to calculate the best way of
estimating the circle, the thing that makes it practically inefficient. Ashgriz and Poo [40]
proposed another method in which a parabolic curve is fitted to the 9 or 25 neighboring cells.
This method is more accurate than the previous one; however, its applications are very limited
to specific cases.

The value of the surface curvature can be defined based on the distance function as:

Ñf
k = -Ñ × n = -Ñ × (54)
Ñf
388 Numerical Simulation - From Brain Imaging to Turbulent Flows

To discretize the above equation, it is required to first calculate the normal vector of the surface
based on Figure 13 and the following relations, and consequently estimate the curvature:

fi +1, j - fi , j fi +1, j - fi , j -1 + fi +1, j +1 - fi +1, j -1

divx ,i + 1 2 , j = 2 , div y ,i + 1
2, j = (55)
Dxi +1 + Dxi Dy j +1 + 2Dy j + Dy j -1

divx ,i + 1 2 , j
nx , i + 1 2 , j = , n i + 1 2 , j = divx2,i + 1 2 , j + div y2,i + 1 2 , j (56)
ni + 12 , j

nx , i + 1 2 , j - nx , i - 1 2 , j n y , i , j + 1 2 - nx , i , j - 1 2
ki, j = - - (57)
Dxi Dy j

Figure 13. Position of the normal vectors of the surface in the cell faces.

One of the advantages of this approach, the level‐set method, is using a distance function which
is smooth and uniform, so that it increases the simplicity of the calculation and accuracy of the
results.

The value of the body force from surface tension of the cell faces can be calculated in CSF
method as:

Ña
Fs = sk nd s = sk n (58)
[a ]
 Free Surface Flow Simulation Using VOF Method 389
[Link]

where δs is the delta Dirac function which is infinite on the interface and zero otherwise, and
its integration is equal to unity. The bar sign in Eq. (58) shows a smoothed (or filtered) value
of volume fraction. The bracket sign shows the difference between maximum and minimum
of volume of the fluid fraction. In CSF method, this function is estimated using | ∇ ᾱ | / α .
Some other references use the following equation to improve the results’ accuracy [41]:

Ña r
Fs = sk nd s = sk n (59)
[a ] [ r ]

in which the ratio of the densities is inserted in order to reduce spurious velocities of the thinner
fluid. The discretized version of Eq. (58) can be obtained as follows:

a i +1, j - a i , j ri + , j
1
Fx ,i + 1 2 , j = 2sk 2
(60)
Dxi +1 + Dxi r max - r min

Based on what was discussed for PCIL method, the following relation is adopted to the present
method:

a i +1, j - a i , j
Fx ,i + 1 2 , j = 2sk Hi+ 1
2, j (61)
Dxi +1 + Dxi

It can be seen that in this equation, the ratio of densities is replaced by the variable H.

5. Parametric method for calculation of curvature of free surfaces

This newly proposed method is based on two sub‐models, the Four-Point Method (FPM) and
the Three-Line Method (TLM). In the former sub‐model, a curve is fitted to the intersection of
the points of grid lines for central and two neighboring cells, while the latter fits a curve to the
free surface so that the distance between the curve and its linear interface approximation is
minimized [42].

5.1. The Four-Point Method (FPM)

In the four‐point method, free surface (as illustrated in Figure 14) is approximated using a
continuous function f ( ⋅ ) ∈ C 2( x1, x4) (a set of functions with continuous second derivation) so
that the distance between the function and the points is minimized according to Eq. (62), and
variations of curvature are bounded based on Eq. (63). In this case, the radius curvature can
be calculated as in Eq. (64).
390 Numerical Simulation - From Brain Imaging to Turbulent Flows

4
inf f (×)ÎC 2 ( x , x ) I ( f ( ×) ) = å f ( xi ) - yi (62)
1 4
i =1

ì k ( x ) - C1 £ e x Î ( x1 , x2 ]
ï
ïï k ( x ) - C2 £ e x Î ( x2 , x3 ]
S .t. í (63)
ï k ( x ) - C3 £ e x Î ( x3 , x4 ]
ï
ïîCi ³ 0 i = 1, 2,3

&&
f ( x)
k ( x) = (64)
1 + f& ( x)
2

Figure 14. Free surface modeling in FPM.

where ε is a small arbitrary given number.

In this method, the desired function is approximated using an n‐degree polynomial function
with unknown constant coefficients. Therefore, we have:

4
inf I ( Pn ( ×) ) = å Pn ( xi ) - yi (65)
i =1

&& ( x )
P
kn ( x) =
n
(66)
1 + P&n 2 ( x )

where κ is the curvature of the free surface.

 Free Surface Flow Simulation Using VOF Method 391
[Link]

It is supposed that Q is the set of f ( ⋅ ) ∈ C 2( x1, x4) such that Eq. (62) is feasible and Q(n) is the
set of Pn (.) such that Eq. (64) is feasible. Then, the following theorem proves that the se‐
quence of solutions for Eqs. (65) and (66) converges to the solution of Eqs. (62) and (64) as n
goes towards infinity.

Theorem 1: if η = in f Q I ( f (.)) and η (n ) = in f Qn I (Pn (.)) then η = limn→∞ η (n ).

Proof. It is obvious that Q (1) ⊂ Q(2) ⊂ ... ⊂ Q , then η (1) ≥ η (2) ≥ ... ≥ η . Therefore, {η (n )} is a non‐
∞

increasing and bounded sequence. Then it converges to a number called ξ . Set W = ∪ Q (n ) ;

n=1

therefore, infQn I (Pn ( ⋅ )) = ξ . Since W ⊂ Q, then ξ ≥ η . By the properties of infimum, for every ε > 0,
there exists f ( ⋅ ) ∈ Q such that:

h < I ( f ( ×) ) < h + e (67)

As we have f ( ⋅ ) ∈ C 2( x1, x4), there is a set of polynomials such that {Pn ( ⋅ )}, {Ṗ n ( ⋅ )}, {P̈ n ( ⋅ )} are
uniformly convergent. Therefore, there exists a natural number such that for every n ≥ N we
have:

Pn ( ×) - f n ( ×) ¥
<d (68)

P&n ( ×) - f&n ( ×) <d (69)

¥

&& ( ×) - &&
P f n ( ×) <d (70)
n ¥

Now, it is claimed that there is an N 1 ≥ N such that for i = 1, 2, 3, 4, the relation

| κn ( x ) − Ci | ≤ ε; ∀ x ∈ ( xi , xi+1 is true. Since otherwise for every n ≥ N , there is a
| κn ( x ) − Ci | ≤ ε; ∀ x ∈ ( xi , xi+1 . Therefore, limn→∞ | κn ( x ) − Ci | > ε which contradicts the assump‐
tion of f ( ⋅ ) ∈ Q . Thus, PN ( ⋅ ) ∈ Q ( N 1) ⊂ W ⊂ Q . Based on what was mentioned,
| I (P N ( ⋅ )) − I ( f ( ⋅ )) | < ε
1
or I (P N 1( ⋅ )) < I ( f ( ⋅ )) + ε < η + 2ε or η ≤ ξ + 2ε ; therefore, ξ = η , or
limn→∞ η (n ) = η .

Thus, our aim is to solve Eqs. (63) and (64), and one can write them in the following forms:

inf I ( Pn ( ×) )
392 Numerical Simulation - From Brain Imaging to Turbulent Flows

ì x2
ï ò k n ( x ) - C1 - e + k n ( x ) - C1 - e dx = 0
ï x1
ï x3
ï k x - C - e + k x - C - e dx = 0
n( ) n( )
S .t. í xò2
ï 2 2
(71)
ï x4
ï
ï ò k n ( x ) - C3 - e + k n ( x ) - C3 - e dx = 0
ï x3
ïîCi ³ 0 i = 1, 2,3, n = 1, 2,...

or

inf I ( Pn ( ×) )

ì h1
ï 2 éë E1n ( x1 ) + E1n ( x1 + h ) + ... + E1n ( x1 + m1h ) ùû = 0
ï
ï h2 é E x + E x + h + ... + E x + m h ù = 0
ï
S .t = í 2 ë 2 n 2
( ) 2n ( 2 ) 2n ( 2 2 )û
(72)
ï h3
ï éë E3n ( x3 ) + E3n ( x3 + h ) + ... + E3n ( x3 + m3h ) ùû = 0
ï2
ïîCi ³ 0 i = 1, 2,3. n = 1, 2,...

such that Ein ( x ) = ∥ κn ( x ) − Ci | − ε + ∥ κn ( x ) − Ci | − ε ∥ .

Now, the intervals x1, x2 , x2, x3 and x3, x4 are divided into three equal sections, m1, m2, m3,
x2 − x1 x3 − x2 x4 − x3
respectively. This means that h 1 = m1 , h2= m2 and h 3 = m3 . Thus, using a numerical
integration method such as trapezoidal rule, the problem above can be reformulated as:

4
inf å Pn ( xi ) - yi
i =1

ì h1
ï 2 éë E1n ( x1 ) + E1n ( x1 + h ) + ... + E1n ( x1 + m1h ) ûù = 0
ï
ï h2 é E x + E x + h + ... + E x + m h ù = 0
ï
S .t = í 2 ë 2 n 2
( ) 2n ( 2 ) 2n ( 2 2 )û
(73)
ï h3
ï éë E3n ( x3 ) + E3n ( x3 + h ) + ... + E3n ( x3 + m3h ) ùû = 0
ï2
ïîCi ³ 0 i = 1, 2,3. n = 1, 2,...
 Free Surface Flow Simulation Using VOF Method 393
[Link]

This is a nonlinear set of equations and can be easily solved using Matlab or Lingo software.

5.2. The Three-Line Method (TLM)

In this method as illustrated in Figure 15, the main goal is to find a function as f ( ⋅ ) ∈ C 2( x1, x4)
within x ∈ xj , x j+1 such that the distance between the function f ( x ) and the line L i ( x ) which
connects the given points ( xi , yi ) and ( xi+1, yi+1) is minimized for i = 1, 2, 3. Thus, variation of
curvature is bounded according to Eq. (74):

3
inf I ( f ( ×) ) = å ò
xi +1
f ( x ) - Li ( x ) dx
xi
i =1

ì k ( x ) - C1 £ e x Î ( x1 , x2 ]
ï
ïï k ( x ) - C2 £ e x Î ( x2 , x3 ]
S .t. í (74)
ï k ( x ) - C3 £ e x Î ( x3 , x4 ]
ï
ïîCi ³ 0 i = 1, 2,3

Figure 15. Free surface modeling by TLM.

Similar to what was discussed in the four‐point method, in this method the function f ( ⋅ ) can
be replaced with an n‐degree polynomial, Pn ( ⋅ ) as below:

3
inf I ( Pn ( ×) ) = å ò
xi +1
Pn ( x ) - Li ( x ) dx
xi
i =1

ì k n ( x ) - C1 £ e x Î [ x1 , x2 ]
ï
ïï k n ( x ) - C2 £ e x Î ( x2 , x3 ]
S .t. í (75)
ï k n ( x ) - C3 £ e x Î ( x3 , x4 ]
ï
ïîCi ³ 0 i = 1, 2,3, n = 1, 2,...
394 Numerical Simulation - From Brain Imaging to Turbulent Flows

Theorem 2: Sequence of the solution of Eq. (75) converges to the solution of Eq. (74).

Proof: The method of proof of this theorem is similar to the previous theorem. In the same
approach of FPM, the following problem is achieved:

æ Pn ( x1 ) - Li ( x1 ) + 2 Pn ( x1 + h ) - Li ( x1 + h ) + ... ö
3
hi ç ÷
inf I ( Pn ( ×) ) = å ç +2 Pn ( x1 + (mi - 1)h ) - Li ( x1 + (mi - 1)h ) + ÷
i =1 2 ç ÷
ç Pn ( x1 + mi h ) - Li ( x1 + mi h ) ÷
è ø

ì h1
ï 2 éë E1n ( x1 ) + E1n ( x1 + h ) + ... + E1n ( x1 + m1h ) ùû = 0
ï
ï h2 é E x + E x + h + ... + E x + m h ù = 0
ï
S .t = í 2 ë 2 n 2
( ) 2n ( 2 ) 2n ( 2 2 )û
(76)
ï h3
ï éë E3n ( x3 ) + E3n ( x3 + h ) + ... + E3n ( x3 + m3h ) ùû = 0
ï2
îïCi ³ 0 i = 1, 2,3. n = 1, 2,...

The steps of using the above equations are as follows:

Step 1: Read ε , and set n = 1.

Step 2: Solve Eq. (73) or (76) in the FPM or the TLM, respectively.

Step 3: If the previous step is infeasible, set n=n+1, and go to step 2, else set the value of target
function in I n .

Step 4: Set n=n+1, and solve Eq. (73) or (76) in the FPM or the TLM, respectively.

Then set the value of target function in I n .

Step 5: If | I n − I n−1 | > ε1 then go to step 4, else I n is the final answer.

6. Conclusions

In this chapter volume of fluid (VOF) scheme was introduced. This is one of the most effective
methods employed in the simulation of two fluid flows interfaces with dramatic changes in
density and viscosity. . These interfaces are represented implicitly by the values of a color
function which is the fluid volume fraction. The advantage of the method is its ability to deal
with arbitrarily shaped interfaces and to cope with large deformations, as well as interface
rupture and coalescence in a natural way. In VOF the mass is rigorously conserved, provid‐
ed the discretization is conservative. However, advecting the interface without diffusing,
dispersing, or wrinkling is a big issue. This can either be performed algebraically, in schemes
 Free Surface Flow Simulation Using VOF Method 395
[Link]

such as CICSAM or geometrically, in schemes such as PLIC. Herein, the viscous fluid
governing equations which are Navier‐Stokes coupled with VOF equation were presented.
Then the most popular VOF schemes such as donor‐acceptor, Hirt‐Nichols, FCT, Youngs, and
PLIC were explained. CICSAM, HiRAC, HRIC, STACS, and some other up‐to‐date pro‐
posed methods were introduced and the accuracy and time calculation of each method were
evaluated. Moreover, surface tension modeling and parametric study of interfaces were
discussed. The author hopes this brief presentation of the VOF method will be beneficial for
scientists and students in their further researches and will help them to massively and
continuously expand this very challenging field of fluid mechanics.

Author details

Mohammad Javad Ketabdari

Address all correspondence to: ketabdar@[Link]

Faculty of Marine Technology, Amirkabir University of Technology, Tehran, Iran

References

[1] Noh WF. CEL: A Time‐Dependent, Two Dimensional, Coupled Eulerian‐Lagrangian

Code. in Alder, B., Fernbach, S. and Rotenberg, M. (Eds.), Methods in Computational
Physics, Volume 3, Fundamental Methods in Hydrodynamics, Academic Press, New
York, NY 1964; 117–179.

[2] Brackbill JU, Kothe DB, Zemach C. A continuum method for modeling surface tension.
Journal of Computational Physics 1992; 100: 335–354.

[3] Ketabdari MJ, Saghi H. A novel algorithm of advection procedure in volume of fluid
method to model free surface flows. International Scholarly Research Network 2012;
10: 1–16.

[4] Ubbink O. Numerical prediction of two fluid systems with sharp interfaces. PhD Thesis,
University of London 1997.

[5] Ramshaw JD, Trapp JA. A numerical technique for low‐speed homogeneous two‐phase
flow with sharp interfaces. Journal of Computational Physics 1976; 21: 438–453.

[6] Hirt CW, Nichols, BD. Volume of Fluid (VOF) method for the dynamics of free
boundaries. Journal of Computational Physics 1981; 39: 201–225.

[7] Torrey MD, Mjolsness RC, Stein LR. NASA‐VOF3D: a three dimensional computer
program for incompressible flows with free surfaces. Technical Report LA‐11009‐MS,
Los Alamos National Laboratory 1987.
396 Numerical Simulation - From Brain Imaging to Turbulent Flows

[8] Lafaurie B, Nardone C, Scardovelli R, Zaleski S. Zanetti G. Modelling merging and

fragmentation in multiphase flows with SURFER. Journal of Computational Physics
1994; 113: 134–147.

[9] Heinrich P. Nonlinear numerical model of landslide‐generated water waves. Interna‐

tional Journal of Engineering Fluid Mechanics 1991; 4(4): 403–416.

[10] Tomiyama A, Zun I, Sou A, Sakaguchi T. Numerical analysis of bubble motion with
the VOF method. Nuclear Engineering and Design 1993; 141: 69–82.

[11] Tomiyama A, Sou A, Hisato M, Sakaguchi T. Numerical analysis of a single bubble by

VOF method. JSME International Journal Series B Fluids and Thermal Engineering
1993; 36(1): 51–56.

[12] Tomiyama A, Sou A, Sakaguchi T. Numerical simulation of a Taylor bubble in a

stagnant liquid inside a vertical pipe. JSME International Journal Series B Fluids and
Thermal Engineering 1996; 39(3): 517–524.

[13] Ketabdari MJ, Nobari MRH, Moradi Larmaei M. Simulation of waves group propaga‐
tion and breaking in coastal zone using a Navier‐Stokes solver with an improved VOF
free surface treatment. Applied Ocean Research 2008; 30: 130–143.

[14] Nobari MRH, Ketabdari MJ, Moradi Larmaei M. A modified volume of fluid advec‐
tion method for uniform Cartesian grids. Applied Mathematical Modelling 2009; 33(5):
2298–2310.

[15] Rezaei H, Ketabdari MJ. A modified flux corrected transport method for incompressi‐
ble free surface flow. Canadian Journal on Computing in Mathematics, Natural
Sciences, Engineering and Medicine 2012; 3(6): 213–219.

[16] Boris JP, Book DL. Flux‐corrected transport. I. SHASTA, a fluid transport algorithm
that works. Journal of Computational Physics 1973; 11: 38–69.

[17] Zalesak ST. Fully multidimensional flux‐corrected transport algorithms for fluids.
Journal of Computational Physics 1979; 31: 335–362.

[18] Youngs DL. Time‐Dependent Multi‐Material Flow with Large Fluid Distortion. in
Morton, K.W. and Baines, M.J. (Eds.), Numerical Methods in Fluid Dynamics 1982;
Academic Press, Reading, UK.

[19] Rudman M. Volume‐tracking methods for interfacial flow calculations. International

Journal for Numerical Methods in Fluids 1997; 24: 671–691.

[20] Rider WJ, Kothe DB. Reconstructing volume tracking. Journal of Computational
Physics 1998; 141: 112–152.

[21] Shirani E, Ashgriz N, Mostaghimi J. Interface pressure calculation based on conserva‐

tion of momentum for front capturing methods. Journal of Computational Physics 2005;
203: 154–175.
 Free Surface Flow Simulation Using VOF Method 397
[Link]

[22] Davis SF. Flux difference splittings and limiters for the resolution of contact disconti‐
nuities. Applied Mathematics and Computation 1994; 65: 3–18.

[23] Ghobadian A. Development of a method for numerical simulation of flows with

moving interfaces. National Power Report 1991; TEC/L/0077/M91.

[24] Van Leer B. Towards the ultimate conservative difference scheme. IV. A new ap‐
proach to numerical convection. Journal of Computational Physics 1977; 23: 276–299.

[25] Pericleous KA, Chan KS. The SEA Method for Free‐Surface Problems with Heat
Transfer and Change of Phase. in Crowe, C.T. (Ed.), Numerical Methods in Multi‐
phase Flows 1994, presented at the ASME fluids engineering division summer
meeting 1994; American Society of Mechanical Engineers, 227–236.

[26] Ubbink O, Issa RI. A method for capturing sharp fluid interfaces on arbitrary meshes.
Journal of Computational Physics 1999; 153: 26–50.

[27] Hogg PW, Gu XJ, Emerson DR. An implicit algorithm for capturing sharp fluid
interfaces in the volume of fluid advection method. European Conference on
Computational Fluid Dynamics 2006; ECCOMAS CFD, TU Delft, The Netherlands.

[28] Heyns JA, Malan AG, Harms TM, Oxtoby OF. Development of a compressive surface
capturing formulation for modelling free‐surface flow by using the volume‐of‐fluid
approach. International Journal for Numerical Methods in Fluids 2013; 71: 788–804.

[29] Muzaferija S, Perić M. Computation of Free Surface Flows Using Interface‐Tracking

and Interface‐Capturing Methods. in Mahrenholtz, O. and Markiewicz, M. (Eds.),
Nonlinear Water Wave Interaction, Computational Mechanics Publications, Sout‐
hampton 1998; p. Chap. 3.

[30] Tsui YY, Lin SW, Cheng TT, Wu TC. Flux‐blending schemes for interface capture in
two‐fluid flows. International Journal of Heat and Mass Transfer 2009; 52: 5547–5556.

[31] Darwish M., Moukalled F. Convective schemes for capturing interfaces of free‐surface
flows on unstructured grids. Numerical Heat Transfer, Part B: Fundamentals 2006; 49:
19–42

[32] Jasak H, Weller HG. Interface Tracking Capabilities of the Inter‐Gamma Differencing
Scheme. Imperial College, University of London, 1995; p. 9.

[33] Peng D, Merriman B, Osher S, Zhao H, Kang M. A PDE‐based fast local level set method.
Journal of Computational Physics 1999; 155: 410–438.

[34] Sussman M, Puckett EG. A coupled level set and volume‐of‐fluid method for comput‐
ing 3d and axisymmetric incompressible two‐phase flows. Journal of Computational
Physics 2000; 162: 301–337.

[35] Renardy Y, Renardy M. PROST: a parabolic reconstruction of surface tension for the
volume‐of‐fluid method. Journal of Computational Physics 2002; 183: 400–421.
398 Numerical Simulation - From Brain Imaging to Turbulent Flows

[36] Jamet D, Torres D, Brackbill JU. On the theory and computation of surface tension: the
elimination of parasitic currents through energy conservation in the second‐gradient
method. Journal of Computational Physics 2002; 182: 262–276.

[37] Popinet S, Zaleski S. A front‐tracking algorithm for accurate representation of surface

tension. International Journal for Numerical Methods in Fluids 1999; 30: 775–793.

[38] Bohluly A, Borghei SM, Saidi, MH. A new method in two phase flow modeling of a
non‐uniform grid. Scientia Iranica, Transaction B: Mechanical Engineering 2009; 16(5):
425–439.

[39] Chorin AJ. Curvature and solidification. Journal of Computational Physics 1985;
58:472–490.

[40] Ashgriz N, Poo JY. A computational method for determining curvatures. Journal of
Computational Physics 1989; 84: 483–491.

[41] Kothe DB, Rider WJ, Mosso SJ, Brock JS, Hochstein JI. Volume Tracking of Interfaces
Having Surface Tension in Two and Three Dimensions. AIAA 34th Aerospace Sciences
Meeting and Exhibit 1996; AIAA 96‐0859, Reno, NV, USA.

[42] Saghi H, Ketabdari MJ, Zamirian M. A novel algorithm based on parameterization

method for calculation of curvature of the free surface flows. Applied Mathematical
Modelling 2013; 37: 570–585.