Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 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: VolumeOfFluid, CICSAM, surface tension, interfacial modeling, multiphase flow
Abstract
The DLR computational fluid dynamics code TAU is enhanced with the volumeoffluid method (VOF) for a
modelling of fluid interfaces. Within Tau the compressible ReynoldsAveragedNavierStokes equations are discretized by a
finite volume technique. TAU works with both structured as well as unstructured grids. In the preprocessing 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 volumeoffluid 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 geometricalbased 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 costefficient 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 MTAerospace, 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 VolumeOfFluid (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 donoracceptor 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
geometricbased 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 steplike function, transported
with local flow velocities, therefore several researchers used
highresolution 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 dimensionalarbitrary 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
threedimensional 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 30June 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 xdirection [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/mK]
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/ms]
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 ReynoldsAveragedNavierStokes
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
timestepping and the implicit dual time stepping method
are employed. With a preconditioning of the
timedependant 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 twoequation 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 PiP4). 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 30June 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 pointdata 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
ReynoldsAveragedNavierStokes 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 NavierStokes 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 2S6Vv (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 donoracceptor
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 30June 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 onedimensional 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
JAfU
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 ULTIMATEQUICKEST (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 30June 4, 2010
=1.5 and k = 2. These results are explained in Section
6.
The influence of the kvalue 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 UQscheme 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 timeaccurate with global timestepping.
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 30June 4, 2010
Results for the structured grid
Figure 7 illustrates the results for the structured grid
after one complete cycle with varying kvalues. 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 kvalue is increased in 0.5steps,
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 xand ydirection 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 threedimensional. 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 30June 4, 2010
vice versa. The shape after the transition is illustrated in
Figure 10 and is obviously solved without additional
smearing or parasitic breakups 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 ydirection.
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 taillike 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 30June 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
kvalue and the results for k=0.5, 1.0 and 1.5 are presented.
The kvalue 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 ULTIMEQUICKEST 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 sinelike curve. On both
grids, the lower left comer showed a slight taillike
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 201225 (1981)
Youngs, D.L., TimeDependant MultiMaterial Flow With
Large Fluid Distortion, Numerical Methods for Fluid
Dynamics, edited by Morton, K. W, Baines, M. J., pp.
273285 (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. 134147 (1994)
Rider, W.J., Kothe, D.B., Stretching and Tearing Interface
Tracking Methods, AIAA 951717 (1995)
Rider, W.J., Kothe, D.B., Reconstructing Volume Tracking,
Journal of Computational Physics, Vol. 141, pp. 112152
(1998)
Rudman, M., Volumetracking methods for interfacial flow
calculations, International Journal for Numerical Methods
in Fluids, Vol. 24, pp. 671691 (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. 187210
(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. 10231040 (2005)
Leonard, B.P, The ULTIMATE conservative difference
scheme applied to unsteady onedimensional advection,
Computer Methods in Applied Mechanics and Engineering,
Vol. 88, pp. 1774 (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 30June 4, 2010
Preconditioning Methods, DLR Forschungsbericht 9529,
DLR, 1995.
Brackbill, J.U., Kothe, D.B., Zemach, C., A Continuum
Method for Modeling Surface Tension, Journal of
Computational Physics, Vol. 100, pp. 335354 (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.61741 (1988)
