Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
A VOF Method Based on the Unstructured Grid
Meng Huang, Lilong Wu, Bin Chen*
Xi'an Jiaotong University, Energy & Power Engineering, Thermal Power Engineering, State Key Laboratory of Multiphase
Flow in Power Engineering
No.28 Xianning West Road, Xi'an, 710049, China
chenbin@mail.xjtu.edu.cn
Keywords: VOF, SLIC, unstructured grid
Abstract
In this paper, an interface reconstruction technology similar with SLIC method on structured grid is developed for the
triangular unstructured grid. It is a totally Eulerian method. For a cell containing the interface, a straight line segment parallel
with one edge of the triangle cell is set to the interface according to the VOF function values in the three neighbor cells and
volume flux of the particular fluid is determined by the geometry relation between the interface line and volume flux on each
edge. Three classical cases are used to test this unstructured SLICVOF interface reconstruction method including advection
test of a hollowed square, Zalesak slotted disk rotation test as a test of scalar advection, and singlevortex shearing flow test,
which are effective to test interface capturing method. It was found that the interface transported back by the reversed velocity
field is comparative with SLICVOF on structure grid. However, further improvement is to be achieved.
Introduction
Two phase flow is a very common phenomenon in nature
and industry. When learning the flowing character of two
immiscible fluids, the interface locating becomes important.
Over the last several decades, mainly two groups of
approaches were proposed to compute the twophase flow
problems with moving interface: interface tracking (Front
tracking method) and interface capturing methods (Level set
and VOF methods). Each method has its own advantages as
well as drawbacks, and no method is applicable to the wide
range of possible twophase flow phenomena. In the front
tracking method (Unverdi & Tryggvason 1992, Tryggvason
et al. 2001), an Eulerian grid is used to solve the
NavierStokes equation for the two immiscible fluids and
the interface is explicitly represented by the Lagrangian
tracking and redistribution of a number of marker points.
This method can calculate curvature in high accuracy and
consequently preserve the sharpness of the interface.
However, it cannot easily handle significant topology
changes of the interface, for example breakup or
coalescence of bubbles or droplets.
Nevertheless, capturing method can deal with topological
change of the interface easily. Level set method (Osher &
Sethian 1988, Sethian & Smereka, 2003) uses distance
function to implicitly capture the interface efficiently. It has
advantages in computation of curvature but disadvantages in
mass conservation. As an alternative, VOF (VolumeOf
Fluid) method defines a volume fraction function F to track
the interface implicitly. Once obtain F in each cell by the
advection of flow field, the interface can be reconstructed
by different approximation schemes: SLICVOF (simple
line interface calculation) scheme (Noh & Woodward, 1976),
SOLAVOF (Hirt & Nichols 1981), PLIC (Piecewise Linear
Interface Calculation) scheme (Youngs, 1982), and other
higher order differencing schemes. Although the
threedimensional implementation is not so easy, VOF is
still a popular method to track interface because it can
preserve mass in a natural way.
The methods mentioned above originally resolve fluid
motion by a structured grid and the applications are
therefore restrained to simple domain. For geometrically
complex domain, bodyfitted grid is one choice but it is
usually limited to relatively simple geometries. Methods
based on unstructured grid with triangular (2D) and
tetrahedral (3D) cells are highly flexible to adapt the
domain boundaries and the computational expense could be
saved through the local refinement to a certain extent.
Recently, there is a growing interest to develop interface
tracking/capturing methods on unstructured grid. Gloth et al.
(2003) presented a front tracking method on unstructured
grid with a local level set formulation for the topological
treatment of fronts. Although the results of isolated
discontinuities have not been impaired, the concept of this
method is not being rigorously conservative. Unstructured
grid Level set methods proposed mostly based on finite
element methods, except Herrmann (2008) proposed a
balanced force refined level set grid method, which solves
the NavierStokes equations on an unstructured grid and
track the interface location by a refined level set method on
an auxiliary, highresolution, equidistant Cartesian grid.
As to the VOF method, Gao (1999) developed a finite
element VOF method for unstructured grid. Zhao et al.
(2002) discretized the VOF equation by a highorder
characteristicsbased finite volume scheme and a
matrixfree implicit dual timestepping scheme utilizing
Paper No
fully unstructured (triangular) meshes. Shahbazi et al.
(2003) developed a second order accurate a piecewise linear
remapping volume tracking algorithm on triangular meshes.
JI et. al (2005) validated the similar LagrangianEulerian
Remap (MLER) method by single vortex test and complex
deformation field test. This method is not easy to implement
because it contains a Lagrangian phase, a reconstruction
phase and a remapping phase. Yang & James (2006)
developed analytical relations for reconstructing piecewise
linear interfaces and volume fractions in triangular and
tetrahedral grids, i.e., PLIC method in triangular grids. Yang
et al. (2006) proposed an adaptive coupled level set and
volumeoffluid method for 2D problems on unstructured
triangular grids. Wang et al. (2009) also developed a
coupled level set and volumeoffluid method for sharp
interface simulation of plunging breaking waves, in which
the interface is reconstructed via a VOFPLIC scheme and
the Level set function is redistanced based on the
reconstructed interface.
To date, it will be of great theoretical importance and
practical value to develop efficient VOF method based on
the unstructured grid with high accuracy to capture the
interface in geometrically complex domain. In this paper,
we will develop an interface reconstruction technology on
unstructured grid similar with structured SLIC method. For
a cell containing the interface (also called a mixed cell), a
straight line segment parallel with one edge of the triangle
cell is set to the approximation interface according to the
volume fraction of the liquid in the three neighboring cells
and volume flux is determined by the geometry relation
between the interface line and the normal velocity on each
edge.
Nomenclature
f VOF function
F discretized volume fraction
AS, area of the mesh cell i (m2)
u velocity (ms1)
nk the normal unit vector on the kth edge
T period
Superscripts
e exact solution
n numerical simulation
Subscripts
r relative
g geometrical
m mass
Numerical Scheme
The VOF functionfis a continuous function defined in the
calculation domain and it denotes the exact distribution of
volume fraction for the two fluids. The volume fraction F,
an integrally averaged quantity, is defined as the portion of
the area (2D) or volume (3D) of the cell filled with the
specific phase fluid. Therefore, all cells can be classified
into three types: full with specific phase (F=l), empty with
specific phase (i.e., full with the other phase, F=0), and
interfacial cell (0
function f and the discretized volume fraction F are related
as follows
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
F, = f (x, y, t)ds (1)
AS,
where AS, is the area of the ith cell. For a given flow field u,
the evolution of VOF function f is governed by following
advection equation because this function must be conserved
as it is convected with the fluid
+ (uV) f = (2)
at
This equation can be recast in integration as follows:
f 3
+fds+ (ulkf)dl0, k =1,2,3 (3)
O1 AS k=1
where k denotes the each edge of the triangular cell i, I is the
edge length, and is the normal vector on edge k. Here
nonstaggered grid is used, i.e., all variables are stored in the
center of the triangular cell. Assuming, we can get the
following equation by the temporal discretization
3
ffS fd" fdS" + (f.ukdt)dl = 0 (4)
AS, AS, k=1
Substituting Equation (1) into Equation (4), we obtain
1 1 n 3 d
S F" ( f H dt) dl (5)
AS, k=1
As a result, the change of volume fraction F in current
triangular cell depends on the solution of the RHS terms,
which can be explained as the mass transportation rate
through each boundary. They can be solved if the mass
transportation between each cell are worked out.
After advection of volume fraction, interface reconstruction
techniques are essential to reconstruct an interface
approximately through the given volume fractions in every
grid cell inversely. Similar with SLICVOF, we use a
sequence of straight line segments aligned with the grid to
approximate the interface inside the triangular cell by means
of geometric method. Then the problem will be how to
locate the line segment in each cell.
Since the volume fraction of each grid is known at the
beginning of every time step, we could judge the position of
the line segment in a cell according to the volume fraction
of neighboring cells and that of itself. Consider an arbitrary
triangular grid cell of a given volume fraction F and assume
the volume fractions in its three neighboring cells satisfy
Fl>F2>F3, we can locate the interface and specific phase
location through the judgment of the following criteria:
K =(FF2)(F2 F3)= Fl2F2+F3 (6)
Typically, there will be two cases corresponding to K>0 or
K<0 shown in Figure 1, where the dash and solid lines
inside the central interface cell are the actual and
approximation interface, respectively. Shadow region is the
specific phase fluid, while the region filled with cross lines
is the other phase fluid. Figure l(a) shows the case K>0,
which means F1 is larger than F2 and F3, so the specific
phase inside the central cell should be near the F1 cell and
we can arrange the line segment (approximation interface)
parallel with the common boundary between the central and
F1 cell. Figure l(b) shows the other case K<0, which means
F3 is smaller than F1 and F2, this time the specific phase
inside the central cell should be away from the F3 cell and
we can arrange the line segment parallel with common
Paper No
boundary between the central and F3 cell. The exact
location of this line segment will be decided by the current
volume fraction F (area ratio) in the central cell.
4
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
'K
\\ <
{ \ \ / \ r '
I~~~
(a) (b)
Figure 1: Line segment approximation of the actual inter
face (a) K>0, (b) K<0
Once we know thf ^g gij of the line sgment for the
approximation of the actual interface, it is necessary to
compute volume fluxes by geometric calculation of the new
volume fraction F" 1 which can be obtained by solving
Equation (5) for each cell edge by edge in a geometric way.
As show in Figure 2, firstly we decompose the velocity
vector into three components perpendicular to each edge to
get Uk, and it is positive whenpoints to the outside of the
cell. We only calculate positivcvolume flux for current cef
so as to ensure that the volume flux of the specific phase
through every edge will be considered and none will be
reconsidered. For a cell fully filled with specific phase fluid
(F=1), the volume flux can be computed directly without the
consideration of the interface. However, for those cells
containing interfaces, it is necessary to take the interrelation
between volume flux and interface location into account.
Figure 3 illustrates four different situations by the ubiety of
the fluid distribution and the current edge 4. We intercept a
length of uk dt from edge 4 to get a net region (filled with
cross pattern) in each situation, then the area of the this
region will be the volume flux through edge 4 to outside. If
a cell is fully filled with the specific phase fluid, the volume
flux can be calculated in the way similar with Figure 2a. If
the cell is empty, nothing shall be done. Substituting the
calculated volume fluxes on each edge into Equation (5), the
new volume fraction F"+1 can be solved for the
reconstruction of interface in the next time step. Repeat this
process, the timedependent reconstruction of interface can
be accomplished for every time step.
Figure 2: Decomposition of velocity vector
(c) (d)
Figure 3: Different situation for the calculation of volume
flux in a triangular cell (a) case (b) case2 (c) case3 (d)
case4 F3 F2
Results and Discussion
As shown in Figure 4, three classical cases are tested to
confirm the validation of our method to capture the fluid
interface. Figure 4(a) is a hollowed square in advection flow.
Figure 4(b) is a Zalesak model in rotating flow and Figure
4(c) is a circle in shearing flow. Those test cases are usually
applied together to show a relatively complete estimate to
an interface capturing method. The first case is used to test
the capability of the method in advection flow; the second
one is to test that in rotating flow and the last one in
shearing flow.
The cases are described in the following three illustrations.
The calculation domain is a 1x1 square and the values off
denote the distribution of the fluids. We operate a reverse
process after a forward process with the same steps and
check the deformation by the comparison of the final and
the beginning fluid distribution, so that the effect of the
method can be estimated.
The velocity fields are described in Equations (79) for tests
a), b) and c) as follows
J =l
v = 1
S= 71(y 0.5)
Iv= 7(x0.5)
71cos 7(x0.5)sin7(y0.5)
7. sin7 (x 0.5)cos 7 (y 0.5)
Paper No
ftI
061
02 04 06 08
(a)
(a)
02 04 06 08
X
(b)
f=O
06
.0.5. 3)
02 f 1
02 04 06 08 1
x
(c)
Figure 4: three test cases. (a) advection test of a hollowed
square. (b) Zalesak slotted disk rotation test. (c) single
vortex shearing flow test
Two different interface tracking errors measurement criteria
Er and E, as well as the relative mass conservation error Em
are conducted to test the accuracy of our method (Zhang et
al. 2008). Er is the relative distortion defined as:
IZ \(T)f[,
Er= (10)
where fe (T) is the exact solution of volume fraction for
ith cell at the end of the test (t=T) and f" is the
computational result after n time steps, fo is the initial
solution. E, is the absolute difference or geometrical error
between the exact and the computed results:
E = (AS j (T) ") (11)
where AS, is the area of the it cell. Finally, the relative
mass conservation error is described as follows:
Zr(T) Zf
E = (12)
In order to compare our method with structured VOFSLIC
method and check the influence of grid quality on the result,
two different grids are used in this work: (1) unstructured
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
grid by Bubble Packing Method (BPM, K. Shimada 1995)
with 7600 cells; (2) 80x80 structured grid. The grids used in
this work are shown in Figure 5.
F" H;': ". . ..; '... ....:
,,;,1%, '"~I ;'%"
:". " .... ;. ;... i . :: ... '. I ...v .':
...... . .. ........
,. ... ..........
".*..  ..__ ..... ....
(a)
Figure 5: Grids (a) BPM (7600 cells)
(6400 cells)
(b)
(b) Structured grids
Figure 6 shows the advection test of a hollowed square with
velocity field (1, 1), respectively. In these figures, each
column shows initial state (0), state after half revolution
(T/2) and state after one revolution (T) for results on
different grid. It can be seen from Figure 6 that for
unstructured grid by BPM, the deformation of the edges is
as small as we expected. As to the result on the structured
grid, the four edges show very small deformation.
The interface tracking errors and the relative mass
conservation error are measured by Equations (1012) and
shown in Table 1. For the advection test of a hollowed
square, structure grid shows the smallest errors because the
interfaces are just parallel to the boundary of each cell. The
results of BPM are in the same order with structured grid. It
is due to the original phase distribution, grids and the
features of the approaches. Since the method based on
structured grid sets a segment parallel to the vertical or
horizontal direction and the signed phase is square shaped,
the deformation should be small. However, the segments set
and the grid edges are usually unparallel in the method
based on BPM grid.
C
a) BPM
b) Structured grid
Figure 6: Advection test of a hollowed square
Table 1: Error measurement of advection test
Grid mbe Er Eg Em
number
BPM 7600 7.84E02 6.37E03 5.17E03
Structured 80x80 2.63E02 1.97E03 1.73E03
Figure 8 shows the Zalesak slotted disk rotation test, which
rotates a slotted circle in a pure rotation velocity field and it
086
oel
Paper No
will return to its initial location after a full revolution of 27t
rotation. Each grid shows numerical diffusion and BPM
result looks best. The diffusion in outside boundary is more
serious for BPM grid, while the diffusion in inside boundary
is more serious for structured grid. As to the error
measurement results shown in Table 3, the interface tracking
errors by structured grid and BPM are in the same order,
while the relative mass conservation error by structured grid
is the smallest, which means the mass conversation ability
of structured grid is best.
b) Structured grid
Figure 8: Zalesak slotted disk rotation test
Table 2: Error measurement of rotation test
Cell
Grid number Er Eg Em
number
BPM 7600 5.26E02 2.34E02 8.61E03
Structured 80x80 5.04E02 2.26E02 2.65E03
Figure 9 shows the singlevortex shearing flow test, the
evolution of an initially circular fluid body placed in the
single vortex flow. The velocity field will stretch the fluid
body into a filament that spirals toward the vortex center
and then shrink the body back into the initial circle after a
full revolution of 2n evolution. This test shows numerical
diffusion by BPM is smallest because not only some
flotsams appear in the half and full revolution but also the
final shapes are apart from initial shape for structured gird.
It can also be seen from Table 4 that all the three errors by
BPM grid are the smallest since the cell shape matches the
rotating velocity field well, not only for the deformation but
also for the mass conservation ability.
a) BPM
e* 0
b) Structured grid
Figure 9: Singlevortex shearing flow test
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Table 3: Error rotation measurement of shearing flow test
Cell
Grid number Er Eg Em
number
BPM 7600 6.34E02 7.92E03 2.25E02
Structured 80x80 1.59E01 1.97E02 6.90E02
Conclusions
A new interface reconstructing method for VOF based on
triangular unstructured grids is presented in this paper. In
this method, line segment is set parallel with one edge of
certain triangular cell containing interface through the
judgment of a new criteria about the volume fractions in the
three neighboring cells of the current cell. Once the location
of the line segment approximation of the actual interface is
decided, the new volume fraction in next time step can be
computed through the geometrical calculation of volume
flux.
In this way, the timedependent reconstruction of interface
can be accomplished for every time step can the algorithm is
examined under three different evolutionary test cases: (a)
advection test of a hollowed square; b) Zalesak slotted disk
rotation test and c) singlevortex shearing flow test. In order
to compare our method with structured VOFSLIC method
and check the influence of grid quality on the result, 80 X 80
structured grid and unstructured grid by Bubble Packing
Method with 7600 cells are used in this work. The evolution
results and error measurements of these three cases on
different grids show that this method is generally capable to
deal with all kinds of simple flows with comparative
accuracy with similar method on structured grid. In shear
flow cases, the numerical diffusion and errors by BPM are
even smaller than that of structured grid best.
Acknowledgements
The authors would like to acknowledge the joint support of
Natural Science Foundation of China (No. 50976087),
Specialized Research Fund for the Doctoral Program of
Higher Education (No.20090201110001) and Program for
New Century Excellent Talents in University (NCET07
0661).
References
Blazek J.: Computational fluid dynamics: principles and
applications. Elsevier Science & Technology (2001)
Gao D. M.: A threedimensional hybrid finite
elementvolume tracking model for mould filling in casting
processes. Int. J. Numer. Meth. Fluids, 29: 877895 (1999)
Gloth O., Hanel D., Tran L.: Vilsmeier R., A front tracking
method on unstructured grids. Computers & Fluids, 32:
547570 (2003)
Herrmann M.: A balanced force refined level set grid
method for twophase flows on unstructured flow solver
grids. Journal of Computational Physics, 227: 26742706
(2008)
Hirt, C. W. and Nichols, B. D., Volume of fluid (VOF)
method for the dynamics of free boundaries, Journal of
a) BPM
Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Computational Physics, 39: 201225 (1981)
Noh W. F. and Woodward P.: SLIC (simple line interface
method), Lecture Notes in Physics, 24: 330 (1976)
Ji C., Wang Y, Wang J.: A Novel VOFType
VolumeTracking Method for FreeSurface Flows Based on
Unstructured Triangular Mesh. China Ocean Engineering,
19: 529538 (2005)
Osher S. and Sethian J. A.: Fronts Propagating with
CurvatureDependent Speed: Algorithms Based on
HamiltonJacobi Formulations. Journal of Computational
Physics, 79: 1249 (1988)
Sethian J. A.and Smereka P.: Level set methods for fluid
interfaces. Annual Review of Fluid Mechanics, 35: 341372
(2003)
Shahbazi K., Paraschivoiu M., Mostaghimi J.: Second order
accurate volume tracking based on remapping for triangular
meshes. Journal of Computational Physics, 188: 100122
(2003)
Shimada K., Gossard D. C.: Bubble Mesh: Automated
triangular meshing of nonmanifold geometry by sphere
packing, in: Proceedings of ACM Third Symposium on
Solid Modeling and Applications, (Salt Lake City, Utah,
USA) 409419 (1995)
Tryggvason G, Bunner B., Esmaeeli A., et al.: A
FrontTracking Method for the Computations of Multiphase
Flow. Journal of Computational Physics, 169: 708759
(2001)
Unverdi S. Q., Tryggvason G. A.: A fronttracking method
for viscous, incompressible, multifluid flows. Journal of
Computational Physics, 100: 2537 (1992)
Wang Z., Yang J., Koo B., Ster F.: A coupled level set and
volumeoffluid method for sharp interface simulation of
plunging breaking waves. International Journal of
Multiphase Flow, 35: 227246 (2009)
Yang X., James A. J., Lowengrub J., Zheng X., Cristini V:
An adaptive coupled levelset/volumeoffluid interface
capturing method for unstructured triangular grids. Journal
of Computational Physics, 217: 364394 2006)
Youngs D. L.: Timedependent multimaterial flow with
large fluid distortion. In: K.W. Morton and M.J. Baines,
Editors, Numerical Methods for Fluid Dynamics 24,
Academic Press, New York (1982)
Zhang Q., Liu P. L.F.: A new interface tracking method:
The polygonal area mapping method. Journal of
Computational Physics, 227:40634088 (2008)
Zhao Y, Tan H. H., and Zhang B.: A highresolution
characteristicsbased implicit dual timestepping VOF
method for free surface flow simulation on unstructured
grids. Journal of Computational Physics, 183: 233273
(2002)
