Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 9.4.1 - Interface advection with the TAU code and the CICSAM method on arbitrary grids
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00226
 Material Information
Title: 9.4.1 - Interface advection with the TAU code and the CICSAM method on arbitrary grids Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Gauer, M.
Hannemann, V.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: VOF
CICSAM
surface tension
interfacial modeling
multi-phase flow
 Notes
Abstract: The DLR computational fluid dynamics code TAU is enhanced with the volume-of-fluid method (VOF) for a modelling of fluid interfaces. Within Tau the compressible Reynolds-Averaged-Navier-Stokes equations are discretized by a finite volume technique. TAU works with both structured as well as unstructured grids. In the pre-processing step the primary grid is split into control volumes, which are again subdivided into tetrahedra. These tetrahedra are cut with a predefined interface to initialize the volume-of-fluid distribution. Thus, the algorithms must be implemented to work efficient with both types of discretizations. For the advection of the interface the CICSAM method (Compressive Interface Capturing Scheme for Arbitrary Meshes) is employed which is far easier to implement when the code is desired to work on unstructured as well as structured grids than geometrical-based methods like PLIC. A further advantage of CICSAM when applied in the TAU code is the fact that the cell information of the primary grid is not explicitly needed. The method is implemented to work with the given dual grid control volume information. The performance of the interface advection is evaluated with the transport of a spherical distribution along the edges of a square grid and a rectangular one transported through an oblique, homogenous velocity field. The test cases are conducted on both structured as well as unstructured grids.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00226
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 941-Gauer-ICMF2010.pdf

Full Text

Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010


Interface advection with the TAU code and
the CICSAM method on arbitrary grids


Markus Gauer and Klaus Hannemann

German Aerospace Center, DLR; Institute of Aerodynamics and
Flow Technology, Spacecraft Section
Bunsenstr. 10, 37073 Gdttingen; Germany
markus.gauer@dlr.de


Keywords: Volume-Of-Fluid, CICSAM, surface tension, interfacial modeling, multi-phase flow




Abstract

The DLR computational fluid dynamics code TAU is enhanced with the volume-of-fluid method (VOF) for a
modelling of fluid interfaces. Within Tau the compressible Reynolds-Averaged-Navier-Stokes equations are discretized by a
finite volume technique. TAU works with both structured as well as unstructured grids. In the pre-processing step the primary
grid is split into control volumes, which are again subdivided into tetrahedra. These tetrahedra are cut with a predefined
interface to initialize the volume-of-fluid distribution. Thus, the algorithms must be implemented to work efficient with both
types of discretizations. For the advection of the interface the CICSAM method (Compressive Interface Capturing Scheme
for Arbitrary Meshes) is employed which is far easier to implement when the code is desired to work on unstructured as well
as structured grids than geometrical-based methods like PLIC. A further advantage of CICSAM when applied in the TAU
code is the fact that the cell information of the primary grid is not explicitly needed. The method is implemented to work with
the given dual grid control volume information. The performance of the interface advection is evaluated with the transport of
a spherical distribution along the edges of a square grid and a rectangular one transported through an oblique, homogenous
velocity field. The test cases are conducted on both structured as well as unstructured grids.


1. Introduction

In the development process of cryogenic upper stages
numerical simulations of the flow behaviour play an
essential role. Experimental test runs are disproportionally
more expensive and thus cost-efficient and reliable
simulation tools are required. To obtain representative and
accurate results, it must to be ensured that the numerical
method is reliable and efficient in the computation of the
important flow regions. In particular, the flow phenomena
which might arise in propellant tanks and transfer lines
during different phases of a mission need to be solved
precisely.
The present work is embedded in the German research
cooperation upper stage which was initiated to coordinate
and perform research on advanced cryogenic upper stage
technologies, see Gerstmann et al. (2009). The consortium
is composed by leading industrial partners, Astrium Space
Transportation and MT-Aerospace, the University of
Bremen with the Center of Applied Space Technology and
Microgravity (ZARM) and Germany's national research
center for aeronautics and space, the DLR. The research
focuses on different selected core technologies and
amongst them is the enhancement of the DLR CFD code
TAU enabling numerical simulations of the flow in
cryogenic upper stages in the long term. The DLR TAU
code provides a numerical platform with a wide range of


applications.
The dynamics of liquid propellant is undesired but
unavoidable in propellant tanks of space vehicles, and
places high demands on both tank and feed system design.
It may lead to a lowering of the feed system performance,
and in a worst case scenario may even lead to a total
failure of the system chain. Violent sloshing in tanks can
cause a sudden decrease in tank pressure, which may lead
to mechanical instabilities. In a microgravity environment
in particular, sloshing may be initiated by very small
perturbations. Therefore, the accurate simulation of
interfacial flow phenomena is an important aspect. The
interface is characterized by a thin layer across which the
fluid properties are abruptly changing. This discontinuous
change yields a challenge in the accurate modelling of the
localized interface physics.
To the present day several numerical methods have been
developed to model these effects. Basically two different
approaches can be distinguished amongst these methods:
moving grid and fixed grid algorithms. With respect to the
flow phenomena which TAU will be required to model, e.g.
condensation, vaporisation, mass and heat transfer, the
fixed grid method delivers a suitable approach due to its
expandability to additional interface physics. Further,
regarding the simulation of violent sloshing and thus
collapsing interfaces, the applicability of moving grid
methods is restricted.






Paper No


One of the most popular and frequently used techniques
to capture interfaces is the Volume-Of-Fluid (VOF) method
developed by Hirt and Nichols (1981). It has been extended
and enhanced by several researchers, e.g. Youngs (1982),
Lafaurie et al. (1994), Rider and Kothe (1995). In the VOF
method, each computational cell includes the information of
the fluid volume f in this cell, initiated from a predefined
interface geometry. Only the interfacial cells have a value
0>f>1 representing the fractional volume of the fluid in this
cell, with the remaining cells having values f= 0 orf= 1,
depending on whether the cell is located outside or inside
the considered fluid.
An essential part of the VOF method is the computation
of the numerical fluxes. This can basically be done in two
ways, either geometrically or algebraically. The original
VOF method by Hirt is based on an algebraic reconstruction
and uses a donor-acceptor formulation with flux limiting
manipulations. The drawback of this method is a strong,
undesired numerical diffusion of the interface. Thus, the
overall performance and accuracy of VOF methods is very
dependent on the used interface reconstruction and
advection scheme.
More accurate results can be obtained by
geometric-based methods where the interface is explicitly
reconstructed after every step, e.g. Youngs (1982) and Rider
et al. (1998). Where earlier methods used simple line
interface approximations (SLIC), nowadays, more
sophisticated methods with a piecewise linear interface
calculation (PLIC) are widespread. However, the drawback
of these methods is a significant increase in algorithmic
complexity and computational time, in particular on
unstructured grids.
The VOF distribution is a step-like function, transported
with local flow velocities, therefore several researchers used
high-resolution schemes for the advection of the VOF
values, e.g. Rudman (1997), Xiao et al. (2003 and 2005).
Very promising results were presented with the THINC
scheme (tangent of hyperbola for interface capturing) by
Xiao et al.. Unfortunately this method is only applied and
designed for structured grids; no advantages to PLIC are
expected when this method might be extended for the
application on unstructured grids. Usually, methods for
structured grids have in common that the side lengths of the
cells, which are equal in all dimensions, are directly
incorporated in the interface reconstruction algorithm.
Thus, the extension to three dimensional-arbitrary meshes is
often not feasible.
Another alternative is the use of high resolution
differencing schemes like TVD methods, FCT schemes and
techniques where the Normalised Variable Diagram (NVD)
by Leonard (1991) is used. With this knowledge Ubbink
(1997) released a new promising scheme which he called
CICSAM (Compressive Interface Capturing Scheme for
Arbitrary Meshes). It applies the NVD concept and switches
between different differencing schemes yielding a bounded
scalar field. At the same time it preserves both the
smoothness and the sharp character of interfaces and is
therefore appropriate for sharp fluid interfaces. An
important feature of CISAM is the applicability on
three-dimensional arbitrary grids without the need of
geometrical reconstruction operations.
This paper outlines the first results of the
implementation of CICSAM in the TAU code. Test cases


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

have been computed where an initial rectangle or spherical
VOF distribution is transported through both the structured
and unstructured computational domain.

The paper is structured as follows. Section 2 gives a
brief overview of the DLR TAU code and an explanation of
its principal functionalities. The algorithm for the
initialisation of the VOF values on structured as well as on
unstructured grids is explained in Section 3. The governing
equations are listed in Section 4. The core equations of the
implemented CICSAM method are listed in Section 5
followed by the test results in Section 6. Finally, the
conclusions are given in Section 7.

Nomenclature


Bo
Ax
C
D
e
F
f
g
k
L
fi
n
q
R
S
u, v
V


Bond number
grid spacing in x-direction [m]
CFL number
diameter [m]
internal energy
normal vector of a cell face
Volume of fluid
gravitational acceleration, [m/s2]
thermal conductivity of the fluid, [W/m-K]
length, [m]
unit normal vector
heat flux, [W/m2]
radius [m]
rate of strain tensor
flow velocity [m/s]
volume [m3]


Greek letters
8 tensorial Kronecker delta
K curvature, [1/m]
p dynamic viscosity, [kg/m-s]
Q smoothing coefficient
p density, [kg/m3]
o surface tension, [N/m]
r viscous stress tensor


Subsripts
A Acceptor
CSF Continuum Surface Force
D Donor
f face
U Upwind


2. The DLR TAU Code

The compressible Reynolds-Averaged-Navier-Stokes
equations are discretized by a finite volume technique using
tetrahedra, pyramids, prisms and hexahedra. The TAU code
can be used for flows ranging from low subsonic to the
hypersonic flow regime. Applications range from the
computation of the outer aerodynamics of air and space
vehicles to the computation of nozzles or the simulation of
shock tube and wind tunnel flows and combustion
chambers. For time accurate computations, global






Paper No


time-stepping and the implicit dual time stepping method
are employed. With a preconditioning of the
time-dependant terms in the governing equations, an
acceleration of convergence is achieved when low speed
flows are simulated (Radespiel et al. (1995)). Further,
various turbulence models are available, including one- and
several two-equation eddy viscosity models.
An important feature of TAU with regard to the
initialisation of the VOF values is the fact that the initial
primary grid is broken down into a secondary or dual grid as
illustrated in Figure 1. The primary grid is composed of
polyedrical elements with triangular and quadrilateral
surfaces. The primary and the secondary grid share the same
points in physical space. However, the secondary grid
consists of control volumes surrounding each primary grid
point (in the illustrated example Pi-P4). For the illustrated
tetrahedron {P1, P2, P3, P41 the point Stet inside the
tetrahedron is used to subdivide the control volumes. In this
example, the control volumes around PI and P2 have the
triangles {Stet, Sr, Sel and {Stea,, Sei, S,;} in common. On
the boundary of the computational domain, the control
volumes of the secondary grid are closed with respect to the
boundary surfaces of the primary grid.


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

calculates the cut points {Cl, C2, C3} of the predefined
interface with the edges, shown in Figure 2c). These are
then needed for the computation of the cut volume {Pl, S2,
S5, C1, C2, C3}, illustrated as the shaded region in Figure
2d). The cut volumes Vcutetr of all tetrahedra are summed
and the VOF value for point P1 is f = (Vc,,et rr)/ Vdua cel
(please note that the aforementioned procedure is the same
for unstructured grids, even though Figure 2 illustrates
structured primary grid cells).
With the aforementioned algorithm and the benefits
inherited from the dual grid cell technique, the calculation
of the initial volume is very accurate also on coarser
meshes. The discrepancies to the analytical solution for the
reconstruction of a sphere are plotted in Table 1. For these
examples, the mesh size is fixed with 303 cells for all test
cases. The grid is structured and equidistant. To test the
accuracy of the initialisation algorithm, the radius of the
sphere differs amongst the test cases to obtain a coarser
discretisation. The sphere is centred in the cubical test
domain and the ratio D/Ax describes the ratio of the sphere
diameter to the grid spacing.


N'xrp


d) 4


and secondary grid points of a


3. Initialisation of the VOF values

In order to reconstruct the interface, the VOF
distribution needs to be initialised. The calculation of the
VOF values is included in the TAU preprocessing. The
VOF values are written in a point-data file which can be
evaluated at every time step.
To increase the accuracy of the VOF initialisation the
algorithm is run on the secondary grid rather than on the
coarser primary grid. The basic principle of the algorithm
can be exemplified with a spherical test case.
Around each point of the primary grid, a secondary grid
cell is constructed which is itself composed of hexahedra (a
quarter of the dual grid cell around P1 is shown in red in
Figure 2a). Each hexahedron is again subdivided into six
tetrahedra, one of which {Pl, S2, S4, S5} is depicted in
Figure 2b). For each tetrahedron, the initialisation algorithm


Figure 2: a) Primary and dual grid cells around point
P1, b) tetrahedron inside one hexahedron of the dual
grid cell, c) cut of the interface plane with the
tetrahedron, d) cut volume in the tetrahedron.


D/Ax Numerical deviation
of the volume [%]
6 1.07
12 0.173
18 0.0107
24 0.0017

Table 1: Numerical deviation of the sphere volume in
respect to the discretisation of the sphere.


4. Governing Equations

The DLR TAU code solves the compressible
Reynolds-Averaged-Navier-Stokes equations which are
discretized by a finite volume technique.


Figure 1: Primary
tetrahedron.





Paper No


Continuity Equation
For Newtonian fluids the conservation of mass is written
as

-+(Vpv) = 0, (1)
9t
with the density p, the time t and v as the local velocity.

Momentum Equation
The equation for conservation of momentum for a
compressible fluid can be written as
p+ V (vv) = -Vp + V r +pg +F (2)
8t
To account for forces due to surface tension, the source
term FCSF is employed. This term converts the interfacial
surface tension force into a respective volume force which
is then evaluated in the Navier-Stokes equation. The source
term comes into play when the curvature is nonzero, which
is the case for interfacial cells and their neighbours. This
reflects the approach of the Continuum Surface Force
method first introduced by Brackbill et al (1992).
The second term on the right hand side of Eq. (2),
V. r, represents the contribution of the viscous forces,
with the viscous stress tensor
2
r=/u 2S--6Vv (3)
3
The strain tensor S which describes the rate of
deformation of the fluid, is
1 T
S (Vv+ (Vv), (4)
2
where p is the dynamic viscosity, the superscript T refers
to the transpose and 6 is the tensorial Kronecker delta
68 (with = 1 fori = j and 6 = 0fori j ).

Energy Equation
The energy equation can be written in vector notation in
terms of the internal energy e in the following form
e-+V (vpe) =pV v+2 S2 (Vv) V q
t L 3
(5)
with the heat flux q = -kVT as defined by Fourier's law
for heat conduction, where k is the thermal conductivity of
the fluid.


5. Interface Advection

The quality of the solution stronlgy depends on the
method how the VOF values are evaluated on the cell faces
for the flux calculations, see Figure 3. The main focus is on
preserving the resolution of the interface, i.e. to avoid an
unphysical smearing of the interface from step to step over
a growing numbers of grid cells or even a deviation of the
interface shape due to numerical reasons.
The original VOF scheme used the donor-acceptor
philosophy to maintain the resolution of the interface. To
avoid a wrinkling of the interface, which would result
when downwind differencing is solely used (usage of the


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

acceptor cell value), the method switches to upwind
differencing (usage of the donor cell value) when the angle
between the interface and the direction of motion drops
below 45. However, this sudden switch in differencing
and the violation of local boundedness criteria results in an
unphysical deformation of the interface (Lafaurie et al.
(1994), Ubbink (1997)). This is where the improvements
of the CICSAM method may be found, namely in a
smoother switching between the varying differencing
methods under consideration of preserving the
boundedness of the interface.


pFIooVdMAO"


C

Sf~onar
amq


f c~P


Figure 3: Upwind, donor, acceptor and face
values for a one-dimensional control volume


Therefore, for its simplicity and for the ability to be
compliant with unstructured grids, the CICSAM approach
has been implemented in the DLR TAU code to calculate
the face volume fractions for the computation of the
numerical flux over the control volume faces. The working
order of the algorithm and the operations performed by the
solver are explained in the following.
The first step is the determination of the donor and
acceptor control volume according to the flow
direction which is obtained through the volumetric
flux. By doing so the VOF values of the donor cellfD ,
and the acceptor cell fA can be determined, see Figure
3. The VOF value will always be fluxed from the
donor to the acceptor cell via the control volume face,
which is in the following defined with the subscriptf
The upwind value is also required for the subsequent
operations with the normalized variables, which form
the basis for high resolution schemes. However, on
unstructured grids this value might not be available
and therefore the following approximation is used


f = min max([f, 2Vf, d],0),1}


where f, is the predicted VOF upwind value, Vf, is
the volume fraction gradient in the donor cell and d is
the vector pointing from the donor to the acceptor cell
point. By doing so the dependency on acquiring an
upwind value has vanished. The value is bounded
with the physical bounds zero and unity.
With the predicted upwind value the normalized
variable f, for the donor cell is calculated, which
involves as well the donor and the acceptor cell value


f -fU
JA-fU





Paper No


The upper bound of the Convection Boundedness
Criterion (CBC) by Gaskell (1988) is used to calculate
the first normalised face value fj


{minG 1

SfD


when 0 < f < 1


when f, <0 or fD > 1


where c is the CFL number. The CBC is a very
compressive differencing scheme. Thus, this scheme
would be the most suitable when the jump is normal to
the flow direction. However, with decreasing angle
between the interface and the flow direction the CBC
may yield smooth gradients to be squeezed to a step
profile.
The second normalized face value is calculated with
the ULTIMATE-QUICKEST (UQ) scheme by
Leonard (1991). It is a more diffusive scheme
compared to CBC and is therefore more dominant the
closer the interface is aligned with the flow. The
normalized face value is evaluated as
I 1
I {min f 8c + (1- c)(6, + 3)
f = 8


when 0 jf < 1
(9)
when fI <0 or jD >(

Dependent on the angle between the interface and the
neighboring cell points the weighting factor 0 controls the contribution of the two differencing
schemes to the normalized face value calculated by


S = 4I, +(i- 7f) ,


The weighting factor y7 is calculated with the cosine
of the angle Of between the vector normal to the
interface Vf, and the connecting vector d of the donor
and acceptor cell with


O = arccos v J n


which is substituted in


S= min{cos(20)+12 1
[ 2 '


where k is a constant with k = 1 as the recommended
value (Ubbink (1997)). This value is introduced to
control the dominance of the different schemes in the
calculation of the normalized face value. However, in
some selected cases, dependent on grid topology and
interface shape, improved results were obtained with k


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

=1.5 and k = 2. These results are explained in Section
6.
The influence of the k-value choice shall be explained
by the graph plotted in Figure 4. The weighting factor
yf is plotted as a function of the angle Of with k=0.5,
1.0, 1.5, and 2.0. The angle Of is evaluated between the
vector normal to the interface Vf% and the vector
pointing from the donor to the acceptor cell point. For
small angles the interface is therefore oriented normal
to the flow direction, the sketch on the bottom left
applies. When the interface is rather oriented with the
flow, the angle is increasing and the sketch on the
bottom right applies. Since y) is bounded with a
maximum of 1, no contribution of the face value
calculated with the UQ-scheme results in the
normalized face value of Eq.(10). The increase of k
results therefore in a dominance of the face value
calculated with the CBC scheme. As illustrated in
Figure 4, fJ equals f), for the angles where y, (,) > 1,
which is the case where the cosine curves exceed the
dashed line at yf= 1.
Thus, the switching between the differencing schemes
for the curves with k=1.5 and 2.0 is not as smooth as it
is the case for k < 1.


FLOW
-7


FLOW


Figure 4: yf as a function of 9f for k=0.5, 1.0, 1.5 and 2.0.

The CICSAM weighting factor can now be computed
with the normalized face value f, as

f f f (13)


And finally the face value can be computed with


f =(l-,/) f( +,Off


for the calculation of the flux over the respective
control volume face.

6. Results and Discussion

The test cases which have been selected for the
validation of the interface advection algorithm are chosen
to test the performance of the algorithm in terms of
interface smearing and preservation. All test cases are





Paper No


performed on both structured as well as unstructured grids
with spherical and rectangular initial interface distributions.
The surface tension is set to zero to focus on the pure
advection of the interface. All computations are performed
inviscid and time-accurate with global time-stepping.
The performance of the CICSAM method was found to
be strongly dependent on the chosen CFL number. For
growing CFL numbers larger than 0.2 an increased
smearing of the interface was examined. For smaller
numbers no considerable improvement of the result can be
stated.

Test case 1: Circle transported through a rectangular
domain
In this test case the circle is initiated in the lower left
corer of the grid. The circle is transported along the edges
of the computational domain, which are normal to each
other, such that, after three direction changes, the circle is
transported back to its initial position, as illustrated in
Figure 5. The length of the grids' edges is ten times the
radius of the circle and the domain is discretized with
125x125 grid cells for the structured case.


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Results for the structured grid
Figure 7 illustrates the results for the structured grid
after one complete cycle with varying k-values. Figure 7a)
shows a zoomed view of the initial VOF distribution in the
lower left corer of the computational domain. From
Figure 7b) -d) the k-value is increased in 0.5-steps,
yielding the best result for k=1.5. The interface shape is
very well conserved and no smearing of the interface can
be examined.


0- 0


i 1

Go0


Figure 5: Circle transported along the
rectangular grid.


edges of the


For the unstructured mesh the grid cells are clustered in
the area where the circle is transported and coarsen to the
centre of the grid and in the vicinity of the edges. The
unstructured grid is illustrated in Figure 6. With the
coarsening of the grid cells the total number of cells can be
reduced to 4644 which speeds up the computations
compared to the structured case. The grid cells are coarser
between the edges of the circle's flow path to be able to
examine how the method is treating and preserving the
interface shape when transported from fine to coarse
arbitrary grid cells and vice versa.


Figure 6: Unstructured grid for test case 1.


A x
Figure 7: Results for test case 1 on the structured grid.
The three VOF isolines are plotted at VOF = 0.05, 0.5
and 0.95. a) Initialized VOF distribution, b) VOF
distribution with k=0.5, c) after one complete cycle
with k=l, d) with k=-.5.

The test run with k=0.5 was stopped in the lower right
corer since the interface already diverged significantly
from the circular shape and showed the tendency to
converge to a rather cornered shape.
The result with k=- shows a reduced urge to converge
to a square shape. However, for Cartesian grids and the
propagation of spherical shapes in x-and y-direction the
shape tends to converge to an octahedron as illustrated in
Figure 7c). Therefore, test runs with higher values for k are
conducted to further increase the effect of the compressive
CBC scheme.
The result for k=1.5 is illustrated in Figure 7d). The
runs with this value delivered qualitatively the best results.
With k=-.5 the influence of the calculation of the
normalized face value conducted with the compressive
CBC scheme is increased yielding an improved
preservation of the circular shape. No further
improvements are seen for larger k values.
The discrepancies between the initial VOF distribution
and the final result for k=1.5 are illustrated in Figure 8.
The final VOF values (as shown in Figure 7d) are
subtracted from the initial values (Figure 7a). Figure 8
shows a 2D and 3D view of these discrepancies. The
rectangular field in the 2D view illustrates the area which
is plotted three-dimensional. Only half of the circle is
shown since the discrepancies are symmetrical. The solid
black lines illustrate the edges of the grid cells which


----------------LI






Paper No


emphasizes that discrepancies can only be noticed in a few
cells.


Figure 8: Visualization of VOF discrepancies
between the final and the initial VOF distribution
(k=1.5).

Results for the unstructured grid
Figure 9 shows the unstructured grid results in the same
order as listed before, the initialized distribution is shown
in Figure 9a) and b)-d) illustrate the results with k=0.5, 1.0
and 1.5. Again, the computation with k=0.5 delivered
unsatisfactory results regarding blur of the interface. The
interface is stretched over a larger number of grid cells.
On the unstructured grid the interface tends to converge
to a hexagon. Unfortunately, this effect can not be
attenuated when k is increased. Obviously, the solution is
highly grid dependent; the edges of the hexagon are clearly
aligned with the grid.


." .-.-.". *. -..


/ \< *'**.*;l/ ::*7; **:;


Figure 9: Results for test case 1 on the unstructured
grid. VOF isolines are plotted at VOF = 0.05, 0.5 and
0.95. a) Initialized VOF distribution, b) VOF
distribution after one complete turn with k=0.5, c)
with k=l, d) with k=1.5.

Further, the focus was also on the treatment of the
transition from the finer to the coarser grid spacing and


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

vice versa. The shape after the transition is illustrated in
Figure 10 and is obviously solved without additional
smearing or parasitic break-ups of the interface. The
interface compactness is preserved; however, the grid
dependency is again very strong and may be improved
with dynamic grid adaptation where the VOF gradient
serves as the adaptation criterion.















Figure 10: VOF distribution after the transition from
the fine to the coarse discretization.


Test case 2: Rectangle in an oblique velocity field

This test case transports an initial distribution which
has a rectangular shape from the lower left in the upper
right comer. The transport is achieved through an oblique
velocity field with u=v, where u is the velocity in x and v
the velocity in y-direction.




2









Figure 11: Rectangle transported in an oblique velocity
field, with u=v.


Structured grid
Figure 12 shows the results for two different grids, one
with 1252 Figure 12a) and b), and the other with 2502 grid
cells, visualized in Figure 12c) and d). All runs are
conducted with k=-. The rectangle is resolved within 20 or
40 grid cells, respectively.
In both cases the interface is very well conserved. The
interface seems to be stretched for the coarser case due to
the larger grid spacing. No parasitic artefacts can be
examined, except in the lower left comer of the rectangle,
where, in a tiny area, VOF values are tail-like dragged.
This effect is highlighted in Figure 13. This plot also
emphasizes that three more peaks can be examined in the
remaining comers of the rectangle. These peaks can be
ascribed to the initialization of the VOF values on the finer
dual grid. As explained in Section 3, the dual grid control


,~~~g~-~:
.-:ir-.
i :
.~~


z
\






Paper No


volumes are again subdivided in tetrahedra, where the
interface is initially geometrically reconstructed. Since the
flux balance is computed via control volume faces this
finer resolution can, as a matter of fact, not be preserved.


x x
Figure 12: Results for test case 1 on the structured
grid. a) Initialized VOF distribution in the lower left
corner of the 1252 grid, b) final VOF values in the
upper right corner, c) Initialized VOF distribution of
the 2502 grid, d) final VOF distribution.


7,5 6.5
X


.. .. .. . . ... . . .
.. .

.H . . . . . c


.. ..... .. .. ..
iii ..........., ............. i::: C'~
I::1:1::11: ::11^^ ^ :::::1^ :::: C'O


...... ....n
.~ I -t4 I -~~ ..~~~~~ Cl


Figure 13: 2D visualization of the discrepancies
between the final numerically solved VOF distribution
and a physically correct distribution on the 2502 grid.



Unstructured grid
The unstructured grid with 6266 grid cells is illustrated
in Figure 14 as well as the initial and the final VOF
distribution. The interface shape is not as well preserved as
it was examined in the structured grid computations and is
again significantly grid dependent. The upper edge
develops a wavy shape aligned with the grid cell directions.
The upper and lower corer of the rectangle in respect to
the flow direction tends to round off.
However, the interface is not smeared and the thickness
is clearly preserved. The dimensions of the square are
maintained even though the adapted grid cell corridor is


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

becoming narrow to the end of the flow path.
Improvements are expected in combination with the
dynamic grid adaptation.


Figure 14: a) Unstructured grid with the initialization in
the lower left corner, b) initial VOF distribution on the
unstructured grid and b) final distribution.



7. Conclusions

The CICSAM method was implemented in the DLR
TAU code to improve the advection of the VOF values on
both structured and unstructured grids.
Two varying test scenarios are chosen to investigate the
performance of the method. In the first test case, a sphere
is transported along the edges of a rectangle. After one
complete cycle, the final VOF distribution would
analytically and visually be congruent with the initial
distribution, thus, easing the evaluation of results.
The computations were performed with a varying
k-value and the results for k=0.5, 1.0 and 1.5 are presented.
The k-value is introduced in the original CICSAM method
as a weighting factor to control the dominance of the two
different schemes for the calculation of the normalized
face value, the ULTIME-QUICKEST and CBC. The best
result on the unstructured grid was achieved with k=1.5.
The result for k=1.0 showed the trend to converge to an
octahedron. The interface and the spherical shape were
preserved and remained within two grid cells. On the
unstructured grid the sphere converged to a hexagon and
lost its spherical shape. It was not possible to improve the
result with varying k. The results showed a strong
dependence on grid alignment, although the interface
experienced no blurring when transported from fine to
coarse grid discretization and vice versa.
The second test case is the advection of an initial
square VOF distribution through an oblique velocity field,
i.e. 45deg flow. Very satisfying results with few
discrepancies are again achieved on the structured grid. A
similar behaviour as mentioned before is examined on the
unstructured grid. The upper edge of the square deviated


S----- --- -- -- -



.. . . . .. . . .. . .. . .


- -- -- -- --
- - - -

- --
- --- -
- -- -
- ---- -
- -

r


';`~ ~ :; :;,::::::
::: ':';:
:~~ ::~~

~ :'i:
:::::. .~~i ;;:;
~: : ''''

i::::: ~;~~ ~_~:-:
:.~~
'';" :'."
':':~:~i~:i~~ ,iiji~~l,, ~~ i~:


::~1

iii: ::i: iiii:i::iIj

ii;I.i.:iii::. i::i'i1.::iiii; i:1 i
it i.:
ii::ii i : iii ii:i


9.5 10






Paper No


slightly from a straight line to a sine-like curve. On both
grids, the lower left comer showed a slight tail-like
blurring in flow direction.
The CICSAM method proved to be a simple, easy to
implement method with few algebraic operations. It
showed very good results on structured grids in the frame
of the chosen test cases. It is aimed to improve the results
on unstructured grids with a dynamic adaptation where the
VOF gradient is serving as the adaptation criterion.


References

Gerstmann, J., van Foreest, A., Dittrich, L., Gauer, M.,
Manfletti, C., Hiihne, C., da Costa, R., Staiger, A., German
Research Cooperation Upper Stage, CEAS 2009 European
Air and Space Conference, Manchester, UK (2009)

Hirt, C.W., Nichols, B.D., Volume of fluid (VOF) method
for the dynamics of free boundaries, Journal of
Computational Physics, Vol. 39, pp 201-225 (1981)

Youngs, D.L., Time-Dependant Multi-Material Flow With
Large Fluid Distortion, Numerical Methods for Fluid
Dynamics, edited by Morton, K. W, Baines, M. J., pp.
273-285 (1982)

Lafaurie, B., Nardone, C., Scardovelli, R., Zaleski, S.,
Zanetti, G, Modelling Merging and Fragmentation in
Multiphase Flows with SURFER, Journal of
Computational Physics, Vol. 113, pp. 134-147 (1994)

Rider, W.J., Kothe, D.B., Stretching and Tearing Interface
Tracking Methods, AIAA 95-1717 (1995)

Rider, W.J., Kothe, D.B., Reconstructing Volume Tracking,
Journal of Computational Physics, Vol. 141, pp. 112-152
(1998)

Rudman, M., Volume-tracking methods for interfacial flow
calculations, International Journal for Numerical Methods
in Fluids, Vol. 24, pp. 671-691 (1997)

Xiao, F. Ikebata, A., An efficient method for capturing free
boundaries in multi fluid simulations, International Journal
for Numerical Methods in Fluids, Vol. 42, pp. 187-210
(2003)

Xiao, F., Honma, Y, Kono, T., A simple algebraic interface
capturing scheme using hyperbolic tangent function,
International Journal for Numerical Methods in Fluids, Vol.
48, pp. 1023-1040 (2005)

Leonard, B.P, The ULTIMATE conservative difference
scheme applied to unsteady one-dimensional advection,
Computer Methods in Applied Mechanics and Engineering,
Vol. 88, pp. 17-74 (1991)

Ubbink, O., Numerical Prediction of Two Fluid Systems
with Sharp Interfaces, Ph.D. thesis, University of London
(1997)

Radespiel, R., Turkel, E., Kroll, N., Assessment of


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Preconditioning Methods, DLR Forschungsbericht 95-29,
DLR, 1995.

Brackbill, J.U., Kothe, D.B., Zemach, C., A Continuum
Method for Modeling Surface Tension, Journal of
Computational Physics, Vol. 100, pp. 335-354 (1992)

Gaskell, H., Lau, A.K.C., Curvature compensated
convective transport: SMART a new boundedness
preserving transport algorithm, International Journal of
Numerical Methods for Fluids, Vol. 8 pp.617-41 (1988)




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs