ICMF 2010, Tampa, FL, May 30June 4, 2010
A simple and conservative algebraic VOF method for resolving interface in one
cell without geometrical calculations on structured and unstructured grids
Mingyu Sun
Center for Interdisciplinary Research, Tohoku University, Sendai 9808578, Japan
sNtIGcil Iohokli ,1c ip
Keywords: PLICVOF, PLIC/SN, algebraic VOF, conservation, unstructured, structured, adaptive grid
Abstract
A simple and exactly conservative method is proposed for capturing sharp interfaces on arbitrary grids. The classic
PLICVOF method resolves a sharp interface (e.g, Rider and Kothe 1998), but necessitates sophisticated geometrical
reconstructions that can be prohibitive for general grids. The idea of this work is to define the volume flux by
using linear algebraic formulas that approximate geometrical reconstructions of the PLICVOF. In order to control
the appearance of flotsams, the approximated volume flux of a small volume is maintained not less than that of
the geometric PLICVOF. It is not a method based on highorder finitedifference approximation to the advection
equation, which resolves a sharp but still smeared interface. The present method resolves the interface in one grid
cell as a typical PLICVOF without geometrical calculations. Since the method has approximated all geometrical
reconstruction and volume flux evaluation algebraically, its extension to arbitrary grids is straightforward by using the
finite volume method.
Introduction
The volumeoffluid (VOF) is one of the most widely
used methods for the numerical simulation of interfa
cial phenomena (Hirt & Nicholas,1981). In the VOF
method, the evolution of an interface is predicted gen
erally by solving the following advection equation
t + u V = 0, (1)
where u is the velocity. Function p is the volume frac
tion or color function at the discrete level. It is unity in
a cell filled with one phase, and is zero if the cell lies in
the other phase. There are two approaches to solve (1).
One approach is based on the interface geometry in
a cell, and then to update the volume fluxes in a La
grangian fashion. Excellent reviews on this approaches
have been given by Rider & Kothe (1998) and Scar
dovelli & Zaleski (1999). The majority of these geomet
ric VOF methods are based on two key procedures. One
is to reconstruct an interface in a cell, and the other to
advance the volume fraction in time and space. Both are
not trivial, because the volume fraction is discontinuous
across a sharply resolved interface. These procedures re
quires sophisticated geometrical calculations, and their
extension to general grid structures is even more diffi
cult. However, the accurate interface reconstruction with
a proper advection algorithm provides the best accuracy
in the VOF family.
An alternate approach is to discretize the advection
equation with a differencing scheme that guarantees
bounded volume fractions while preventing the smear
ing of interfaces over several mesh intervals. The well
known highresolution schemes such as TVD methods
are too diffusive for this problem. A rather successful
strategy is to introduce artificial compression by down
wind differencing (e.g., Ubbink and Issa 1999). This
approach relies much less on the cell geometry, and
can be fairly readily extended for general grid. How
ever, an interface is often spread in two and more cells.
The interface is practically represented by isocontours of
volume fractions, instead of reconstructed interfaces as
commonly adopted in the geometric VOF approach.
The present work, following the methodology of the
geometric VOF, tries to reconstruct and advect interfaces
using linear algebraic formulas that approximate geo
metrical reconstructions of the PLICVOF
Numerical methods
The advection equation (1) is solved by the finite vol
ume method. The volume flux is evaluated in the direc
tion normal to the grid line, and integrated in a unsplit
ICMF 2010, Tampa, FL, May 30June 4, 2010
manner. Given a flow field u, the advection equation is
rewritten as
Ot + V (up) = pV u,
4t + V (up) 0, (3)
for incompressible velocity field. Equation (3) is dis
cretized as,
2+' = k (un A kAt)j' (4)
where Pk is the volume fraction of phase k defined at the
grid interface, satisfying Ek =P 1. (2k is the volume
of phase k, 2k p = k2, satisfying
Q= 1: (5)
The normal velocity across grid interface u, equals to
Un = U S,
where s is the outward normal vector with a magnitude
of the length of grid line j. In short, (u, Pn At)j repre
sents the volume flux crossing face j in time At.
Without losing generality, the phase with the smaller
volume is denoted as slave volume, Q2, and the other
phase master volume, 2,. In the present method, for
maintaining the positivity or the boundedness of volume
fractions, the slave volume is integrated by (4),
2+ 1 = 2 5(u, ,pAt)j, (6)
where 4 is introduced to ensure the boundedness
2 1E [0, Q]. (7)
For the master volume,
2"+1 (2
"TO "T
S(uanpmAt)j,
where ,m = 1 (Ds. The limiter function 4) is defined
such that the outflow of the slave volume should be no
more than that inside, for cell i,
S for
1
fs > Q2
otherwise,
where fs is total outflow volume flux of the slave vol
ume,
fs= (unAt)j.
j,un>O
The limiter function at grid line is taken from the up
stream cell,
4i u, >0
Dj u, < 0'
The master volume is then updated using (5),
Q,91 Q91. (10)
With this saturated relation, it seems that the equation
(8) is not necessary for volume update. In fact, the small
and master phases are in general not the same fluid in
two neighboring cells, so both volume fluxes have to be
calculated in order to update the volume of the small
phase. Furthermore, the volume flux calculation in (8) is
required to define the numerical fluxes for other conser
vation laws.
It is obvious that the sum of (6) and (8) at the new
time step satisfies the saturated constraint (5), if
S(uaAt)j= 0, (11)
which is the finitevolume discretization of the
divergencefree velocity field. In other words,
divergencefree relation (11) is the sufficient condition
for maintaining the volume conservation of the present
method.
It is noted that the 4) limiter is different from the re
distribution algorithm adopted in some VOF algorithms.
The redistribution algorithm is often triggered when an
abnormal volume that is beyond [0, Q2] is found in the
solution, and then redistribute the volume to somehow
empirically chosen cells nearby. The D limiter is to ad
just the volume flux such that the abnormal volume will
not appear. This limiting method is simple and general
for any grid system, and it is also necessary for the fol
lowing approximate volume flux formulas to preserve
the boundedness of the phase volumes.
Before introducing the present algebraic VOF method
to evaluate the volume flux (unAt)j, we first re
view how a geometric VOF method works. The inter
face in a cell is not tracked explicitly, but reconstructed
and approximated by a simple geometry. The simple
line (piecewise constant) interface calculation (SLIC)
assumes the geometry is a line parallel to one of the
grid lines. Currently most widely used method is based
on the piecewise linear interface calculation (PLIC). A
historic review of the piecewise constant and linear re
constructions was given by Rider and Kothe (1998). In
PLIC, the surface normal vector, n, is required to con
struct the linear interface
n r dl = 0,
where d1 is the line constant. The normal vector is in
ferred from the spatial distribution of p. If p is a smooth
ICMF 2010, Tampa, FL, May 30June 4, 2010
function, the normal vector satisfies
n = V. (13)
In the present work, the normal vector is treated as an
independent variable that are solved separately follow
ing our previous work (Sun, 2010). Therefore in order
to uniquely reconstruct the line segment in the cell, one
needs to find the line constant di. Among many advec
tion scheme, a simple way is to define the volume or area
cut by the grid line j shifted upstream by u,At as the
volume flux (u,a At)j. These two procedures are geo
metrically complicated, especially in axisymmetric and
3D cases on unstructured grids. Therefore, in order to
simplify two procedures, instead of calculating the geo
metric relations precisely, this work tries to approximate
them using simple formulas. We need
1. a simple formula to approximate line constant di,
and
2. a simple formula to evaluate the area cut by the grid
line shifted upstream.
The distance function in the cell is bounded, d e
[dmin, dma,], where din and dax are two extreme
values in those defined at all vertexes of the cell,
dmin = min(dj), d.x. = max(dj).
It is clear that the volume fraction of the dark phase in
the figure satisfies
(d) dmax
1, d = din
The basic idea of this work is to make use of this prop
erty, and to construct a linear function to approximate
the line constant and the volume flux. The linear func
tion is unique,
di dmax s(dmax dmin),
for interface reconstruction. ~, is the volume fraction of
the slave phase. The tilde over p, indicates the value will
be modified to prevent the appearance of flotsams, to be
defined later. The linear function can be reformulated as
dmaxi d (15)
Ts = (15)
dmaxl dminl
for the calculation of volume flux, where dmint and
dma, are taken as the maximum and the minimum value
at two vertexes of the grid line ', as shown in Fig. lb.
Formulas (14) and (15) are exact for grid lines par
allel or perpendicular to the grid lines on a rectangular
grid. It is less accurate for other angles, and attains the
maximum error at 45. Denoting
sin20
1 + sin20'
where 0 is the angle between interface and grid line nor
mal vectors. We introduce a modification of 0s,
Ps, P r < s
Vo4's4, 0' > 4s
b. d_, rn
Figure 1: Line constant calculation (a) and volume flux
calculation (b) using algebraic approximation
Consider a general quadrilateral as shown in Fig. 1.
Given the surface normal n, we introduce a function, d,
representing the relative distance to the interface,
d = n r.
such that the volume flux of the slave phase is not
less than that of the exact geometric reconstruction for
square cells. The proof is neglected here.
Numerical results
The initial phase volumes intersected with interfaces are
exactly calculated. For example, the sum of all phase
volumes inside a circle differs from the exact area by the
SThe definition for the maximum and the minimum values in the eval
uation of volume fluxes is chosen, so that we may find a modifi
cation (16) to prevent this approximate VOF method from gen
erating flotsams. Other options may be possible, but still remain
unexploited, because of the difficulty in controlling flotsams.
ICMF 2010, Tampa, FL, May 30June 4, 2010
order of 0(1016). For all tests, the velocities at cell
faces are also exactly specified, to exclude any possible
error in the treatment of velocity. The CFL number is
taken as 1/8 if not specified. For drawing interfaces, lin
ear segments reconstructed in all interface cells are plot
ted. They are not contours of equal volume fractions.
Surface normal vectors in the interface cells and the cells
filled with the dark phase are also plotted in most fig
ures. The vectors in other cells are used in computation
as well, but not plotted for the sake of clarity. The ar
row indicates the direction of the normal vector, starting
from the cell center, and its length is proportional to the
magnitude of the vector.
The geometric error measure, L1 norm, is defined as
1 N
Cy= N = QL \0a
>2 ZJa
Ii (17)
where Ocal and exact are the calculated and exact vol
ume fractions in cell j with volume (2j. N is the total
number of cells used in the domain.
Translation of circular surfaces on square grid cell In
this test case, the computational domain is a square of
(0,1)2, with open or free boundary outside. Initial sur
face normal vectors for circles are specified the same as
those define in Sun (2010). The initial circle of 0.3 in
diameter centered at (0.25, 0.5) is translated in the con
stant velocity field of (u, v) = (1, 0) fort = 0.5. Fig.
2 shows the reconstructed interfaces at the last step, to
gether with associated surface normal vectors located at
the cell center. It is emphasized that, even for a low res
olution of d/Ax < 3, as shown in Figs. 2a, the particle
is reasonably advected and reconstructed.
Translation of hollow shapes The second test is the
advection of a hollow square and a hollow circle in an
oblique velocity field of (2, 1), which has been gener
ally tested previously (Rudmann 1997, Ubbink and Issa
1999). The side lengths of outer and inner interfaces of
the square is 0.8 and 0.4 respectively, and for the circle
the outer and inner radii are 0.4 and 0.2 respectively. The
shapes are initially centered at (1.2, 1.2) in a square
domain of (2, 2)2 with a 2002 mesh. The solutions
are recorded at t = 1.25. The final shapes are plot
ted in Fig. 3, and the corresponding error distributions
of volume fraction are shown in Fig. 4. The interfaces
are resolved sharply in one cell without smearing, which
has never been realized by the finite difference VOF ap
proach. The solution errors together with a few calcu
lated results in literature, are listed in Table 1. For the
consistency the literature, the error is slightly different
from (17),
Figure 2: A circular particle of d 0.3 after translating
from (0.25, 0.5) to (0.75, 0.5) in xdirection on different
meshes: a) d/Ax = 2.4, 8 x 8; b) d/Ax 4.8, 16 x 16.
c) d/Ax = 9.6, 32 x 32.
1 N
6 j=91 .exactJt j=l
i:I
ICMF 2010, Tampa, FL, May 30June 4, 2010
, '  i L I u ( . .1 u 1 n '
U u'I, ,' , i'
i~t
'4124
/'
"` FIE
Figure 3: Advection of a hollow square (a) and a hollow
circle (b) with an oblique uniform velocity (2,1) using a
2002 square grid.
The small time step is found to improve the accuracy of
the present algebraic method, and the accuracy is com
parable with other methods. A nottrivial improvement
is that all interfaces are resolved in one cell without pro
ducing flotsams.
Table 1: Geometric error for oblique translation. a er
rors obtained from Rudman (1997), and b obtained from
Ubbink and Issa (1999)
SLICa FCTVOFy CICSAMUb Present Present
(CFL=1/8) (CFL=1/16)
Square 1 32E1 1 63E8 3 97E2 6 83E2 6 28E2
Circle 9 18E2 3 99E2 2 84E2 3 83E2 2 14E2
The last example is the deformation of a circular sur
face of d 0.3 in a vortex velocity field. The circle is
initially centered at (0.5, 0.75). The timeresolved ve
locity field is specified as,
u = cos(7t/T)sin (7rz)sin(27y),
Figure 4: Contours for volume fraction error in advec
tion of a hollow square (a) and a hollow circle (b) with
an oblique uniform velocity (2,1) using a 200 2 square
grid.
v = cos(wt/T)sin (2(y)sin(2wx).
The test was taken from Rider and Kothe (1998). This
is a vortical velocity field, which will deform the cir
cle and promote topology changes. It is periodic with
a period of T. The circle will undergo the maximum
deformation at t T/2, and then return to its initial
state. After one period, the circle is supposed to re
cover its original shape. The sequential snapshots to
gether with the corresponding grids are plotted in Fig. 5.
A magnified view of the circle at the last step is shown
in Fig. 6. The solutionadaptive unstructured quadrilat
eral grid (Sun and Takayama 1999) is used. The back
ground coarse grid contains 620 cells, and the number
is increased to 1943 at the end by using a threelevel re
finement. The interface is shown by thick segments that
are exactly reconstructed from the volume fraction be
tween zero and one. Neither interface diffusion and nor
flotsams are seen. This can be more clearly seen from
I s1
La meI
ICMF 2010, Tampa, FL, May 30June 4, 2010
~u~40//i
3
fttra mlqt//tt
J
e. I I I I I I I I f. L++ + 
Figure 5: Results for the circular fluid body placed in
the timereversed singlevortex flow field on a solution
adaptive unstructured grid, T=2. (a) t=0.0; (b) t=2.0;
(c)t=0.4; (d)t=1.6; (e)t=0.8; (f)t=1.2. Figures in the right
column are supposed to be the same as those in the left.
the fact that fine grids are always clustered near the in
terface.
Concluding remarks
An algebraic VOF method has been proposed. The
method is based on a linear approximation to the PLIC
VOF method. Although its accuracy is not as good as ad
vanced geometric VOF methods (e.g. L6pez et al. 2005,
Pilliod & Puckett 2004), it is still comparable with those
finitedifference based methods. The real advantages of
this method are its simplicity and its easy implementa
tion on arbitrary grid structures. The method generates
no flotsams, without explicitly introducing debris sup
pressing tricks, and exactly conserves the phase volumes
with a divergencefree velocity field.
References
C.W. Hirt and B.D. Nichols, Volume of Fluid (VOF)
method for the dynamics of free boundaries, J. Comput.
Figure 6: The zoomed view of Fig. 5b. The interface is
shown by thick segments that are reconstructed exactly
for all cells with volume fraction between zero and one.
Neither interface diffusion and nor flotsams are seen.
Phys. 39 (1981) 201225.
J. L6pez, J. Hernmndez, P. G6mez and F. Faura, An im
proved PLICVOF method for tracking thin fluid struc
tures in incompressible twophase flows, J. Comput.
Phys. 208 (2005) 5174.
J.E. Pilliod Jr., E.G. Puckett, Secondorder accurate
volumeoffluid algorithms for tracking material inter
faces, J. Comput. Phys. 199 (211 14 465502.
W.J. Rider, D.B. Kothe, Reconstructing volume track
ing, J. Comput. Phys. 141 (1998) 112152.
M. Rudman, Volumetracking methods for interfacial
flow calculations, Intl. J. for Numerical Methods in Flu
ids, V24:671691, (1997).
R. Scardovelli, S. Zaleski, Direct numerical simulation
of freesurface and interfacial flow, Annu. Rev. Fluid
Mech. 31 (1999) 567603.
M. Sun, Volumetracking of subgrid particles, to appear
in Intl. J. for Numerical Methods in Fluids, (2010).
M. Sun, K. Takayama, Conservative smoothing on an
adaptive quadrilateral grid, J. Comput. Phys. 150 (1999)
143180.
O. Ubbink, R.I. Issa, A method for capturing sharp fluid
interfaces on arbitrary meshes, J. Comput. Phys. 153
(1999) 2650.
IEi
