7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Parallel block adaptive mesh refinement for multiphase flows
D. Zuzio* and JL. Estivalezes*
Department of Models For Aerodynamics and Energy, ONERA Centre de Toulouse
2, av. Edouard Belin BP 4025 31055 Toulouse Cedex 4, France
davide.zuzio ~onera.fr and JeanLue.Estivalezes ~onera.fr
Keywords: DNS, Level Set, Ghost Fluid, AMR, Multigrid
Abstract
We present a two dimension parallelized adaptive algorithm for incompressible two phase flows. We use a coupled
LevelSet and GhostFluid method to track the interface between phases; the resolution of the Navier Stokes equations
is performed by a two steps projection method. We make use of the PARAMESH package to implement the adaptive
mesh and to get an efficient parallelization of our code. PARAMESH manages a blockbased type AMR: the
computational domain is recursively split into smaller subdomains where needed, until the desired resolution is
reached. Special attention is paid to the elliptic solver for the pressure equation: we have developed a BiCGStab
algorithm allowing the resolution of high density ratio flows over non uniform meshes; in order to speedup the
convergence we make use of a multigrid preconditioner based on a Fast Adaptive Composite algorithm. We refine
the mesh on the interface to the maximum allowed level, leaving the rest of the domain to be less resolved: given the
difficulty of solving the strong discontinuity lying on the interface, we want to show that a locally refined mesh can be
as effective as the corresponding uniform fine mesh. Some validation test cases are proposed, such as the dynamics of
static and rising bubble, and a RayleighTaylor type instability.
Introduction
In the vast topic of flows with interfaces, the study
of liquidgas interactions have a fundamental impor
tance, especially in combustion problems. The inter
facial instabilities, like the KelvinHelmholtz, produce
the droplets or sprays that afterwards participate to com
bustion. In particular where experiments are too diffi
cult or expensive to perform, the numerical simulation
becomes a powerful tool for predicting these physical
phenomena. The direct numerical simulation gives us a
"model free" approach to the Navier Stokes equations, at
the price, however, of a higher computational resources
demand. The strategy developed to improve the feasi
bility of extended domain DNS are the parallel computa
tion and the adaptative mesh. The PARAMESH package
(Peter MacNeice and Packer (2000)) is a set of libraries
developed to satisfy both approach: it is able to create an
adaptive mesh block type data structure, and to aggres
sively distribute the computation among the processor
in order to achieve load balance. In a two fluid config
uration the adaptive mesh is naturally advisable, as if
the precision of the interface discretization depends on
the resolution of the cells, it is often necessary to extend
the computational domain well beyond the interfacial re
gion of the flow, thus involving a lot of computational
effort spent on less interesting zones. Berger (1982)
and coworkers have developed the use of a hierarchy
of logically Cartesian grids and subgrids to cover the
computational domain. They allow for the subgrids to
overlap, have arbitrary shapes and be merged with other
subgrids at the same refinement level whenever appro
priate. This is a flexible and memoryefficient strat
egy, but the resulting code is complex and has proven
to be very difficult to parallelize. In a simplified variant,
the domain is divided in blocks which can be bisected
in each directions when needed (DeZeeuw and Powell
(1992)). This approach, followed by the PARAMESH
developers, intuitively produces meshes with lower re
finement efficiency, but the dedicated tree memory stor
age is lighter, and the cache managing as well as the
parallelization are easier. The resolution of an elliptic
equation for the pressure is an added difficulty for the
incompressible solvers, in particular when coupled with
the complex AMR mesh. The most attractive algorithms
are the Krilor subspace family methods, often coupled
with a preconditioner (Mark Sussman (1998), Teigland
and Eliassen (2001), Michael Oevermann (2000)). In
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
this work a preconditioned BiCGStab works in con
junction with an AMR adapted multigrid. The adap
tive mesh represents here an improvement over the al
gorithm of Couderc and Estivalezes (2003); it does not
involve radical changes of the original code, and retains
the Cartesian structured mesh discretization.
Governing equations
Two fluids are considered; in each of them the density is
given and constant, and the flow is given by the incom
pressible NavierStokes equations:
Numerical model
The variables are located into a staggered mesh, with
pressure and Level Set in the cell center and the veloc
ity components on the cell faces. The Level Set linear
equation (3) spatial derivatives are discretized by a fi
nite volume method and a fifth order Weno; the tempo
ral derivative is integrated by a third order Runge Kutta.
The momentum equation (2) is solved by a second or
der projection method: it works decoupling velocity and
pression by using first an explicit discretization for the
effects of convection in a predictor step, and then enforc
ing the compliance with a velocity divergence constraint
in a subsequent projection step:
*resolution of the momentum equation without pres
sure term
J 9 8x By
V u = 0
~u
+ u )u
 (V T +t f)
With the following hypothesis:
T = pl+ D
D = p Vu + (Vu) )
where u [u, v]T represents the velocity vector field,
p the hydrodynamic pressure and p the density, p the
dynamic viscosity of the fluid, f represents the external
forces like gravity. The fluids are separated by a sharp
interface tracked via level set approach. It is implicitly
given by the zero level of a function (x, t), where 4
is defined as the signed distance to the interface, so that
any quantity a is
a (x, t) = Q(A) (x, t) if (x, t) > 0
al") (x, t) elsewhere
The evolution of the level set follows the linear advec
tion equation
+ u V = (3)
so that it is passively advected by the velocity field. The
jump conditions imposed on the interface are determined
by capillarity and viscosity:
p 8 u"
8 u"
v?~~j 1 =Ui~+ At" i
u"
+v"
1
2
P d2V Z
d **
By"
*resolution of a Poisson equation, given the diver
gence of the predicted velocities as right hand side
8J 1? 8p 8 18p
8x : p x y p y
} v
*correction of the intermediate velocities by the
pressure gradient
atn i
n+] *
U ~
i h~3 ,~3
n 1 atn
'U. 'U
"+1;3 i+ij+4

a/c (4)
0 (5)
[u]0
Here the notation [] (), ()A represents the jump
across the interface, n the normal vector pointing from
A to B, t the tangential one. a is the surface tension, k
the interface curvature.
The viscous terms are treated by a totally explicit sec
ond order centered scheme:
[p,] n [p VBu + (Vu!)t n
t [I Vru + (Vu)~ t
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
" 2 '3 % t13t~a~j
a22
d2Un d2Un)i ,3
dz2 dy2 1
d2VndZ2 d2Vndy2
12 r
U
i 2~3 1
,j + Uni t~jl
ay2
n 2vn vn ~ 1
ij+ 2 %~3 2
a22
vn 'U
a+ Lj+l 2Vn 2 n
2 il~j
ay2
Figure 2: Representation of a 4x4 block with 2 guard
cells in each direction.
than their parents. PARAMESH keeps also track of
physical boundaries on which particular boundary con
ditions are to be enforced, ensuring that child blocks in
herit this information when appropriate. The communi
cation between coarse and fine grids is assured by the
prolongation operation, that involves giving values at
the fine points overlapping coarser ones, and the recip
rocal one, the restriction. The interpolation algorithm
can be implemented by the user, or a default Lagrangian
polynomial interpolation is provided. The discretization
stencils are completed by a user defined set of ghost
cells around each block; Dirichlet boundary conditions
are exchanged among the blocks every synchronization
point: this structure allows the single grid discretization
to be applied to each block without modifications.
Conservation of mass is a key issue and to insure this
care must be taken at the interface between coarse and
fine patches. The fluxes along this interface are calcu
lated on both the coarse and fine patches; these fluxes
do not necessarily match because they are evaluated on
different grids, so that mass could be lost or gained. To
insure mass conservation the flux through a coarse cell
face is therefore corrected to equal the sum of the fluxes
through the corresponding fine cell faces, a procedure
called flux matching shown in figure 3.
Poisson SOlver
In order to solve the linear Poisson pressure equation
(7) within the AMR framework we have combined the
BiCGStab algorithm (VanDerVorst (1992)) with a Fast
Adaptive Composite (McCormick and Thomas (1986))
multigrid preconditioner. The BiCGStab is an evolu
tion of Krilor the subspace family methods CG, CGS
In this section the extension of the numerical method to
an adaptive mesh is explained. In the original Berger's
algorithm a certain number of cells are tagged for re
finement, then clustered to form new finer grids which
are superimposed to the coarser ones; they can over
lap and merge, but they have to be properly nested, so
that the refinement jumps are always one level. This
method is referred to as "patch based". In an alter
native formulation, known as "block based", the com
putational domain is split recursively into smaller ones
where needed, until the finest grid is generated (fig
ure 1). In the present work, the PARAMESH package
Figure 1: Representation of the two kind of AMR, left
blockbased, right patchbased.
(Peter MacNeice and Packer (2000),Olson and MacNe
ice (2005),Olson (2006)) has been chosen to be used
with the code. It is an opensource software package
of Fortran90 subroutines developed by Peter McNeice
and Kevin Olson. The package manages the creation
of a blocktype AMR: builds and maintains the tree
structure which tracks the spatial relationships between
blocks, distributes them among the available processors
and handles all interblock and interprocessor commu
nication. It distributes the blocks so that locality is max
imized and communications minimized. Each block is
an uniform cartesian mesh with fixed user defined di
mensions, shown in figure 2; typically 2D blocks are
4x4 or 8x8. The refinement ratio is fixed to r = 2,
that means refined cells are always two times smaller
Adaptive mesh refinement
I .
Figure 3: Representation of the 2D flux match
ing procedure. In black are represented the mesh
point; in grey the interpolated guardcells points.
Facarse =Fine + F/in /2
and BiCG. It is able to work with unsymmetrical linear
systems and does not require the explicit computation
of the system matrix transpose. When the Laplacian is
discretized on the adaptive mesh, the standard discretiza
tion stencil is applied to the interior points of each block,
given the updated interpolated values in the guardcells.
The global matrix is hence unsymmetrical due to the in
terpolation coefficients.
The multigrid preconditioning is widely em
ployed Coudere (2007),Sussman (2005),Vuik and
J.J.I.M vangan (2000) to speed up the convergence. A
natural way to implement a multigrid algorithm within
an AMR mesh is to utilize the coarse grid to perform
the basic two grid cycle. The AMR tree can be "inter
preted" in two different ways: a multigrid level may
coincide with a refinement level, so that the relaxation is
performed each level on a fraction of the computational
domain, or the level may be built with all the grid points
until the "local" level refinement is reached. The latter
has been chosen for some reason. The first is for parallel
efficiency purpose: in the PARAMESH framework the
repartition of the blocks among the processors is done
for the leaf blocks first, because the most of the time
is spent here. So, if a relaxation sweep is performed
on a fixed level of refinement, some processors only
will likely be at work. Instead, if it is done on the
entire domain at each level, it can benefit from a more
efficiently parallelized relaxation The second reason
is that PARAMESH multigrid routines are meant to
work in a global sense, so that it would not be possible
to easily perform "submultigrid" relaxation over a few
blocks only (like Sussman (2005)), and the FAC is
meant to work in a more global sense.
As for the advection equation, the coarsefine inter
face problem arises. The Dirichlet boundary condition
exchange between blocks makes coarse and fine solution
not sufficiently linked, since the pressure gradient are
computed on two different meshes. This causes a dis
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
continuity in the first derivative which is O(Axcourse),
that acts like a singular charge proportional to the first
derivative mismatch, and corrupts the solution prevent
mng the attainment of the fine grid precision or even
convergence. The solution proposed here is, as sug
gested by Teigland and Eliassen (2001), a Dirichlet and
Neumann boundary conditions coupling. Across the in
terface a preliminary Dirichlet conditions exchange be
tween blocks is allowed in the two ways; then the coarse
grid ghost cells values are corrected to make the first
derivative match the fine grid value, so that the coarse
grid actually receives an imposed Neumann boundary
values.
Numerical results
Static bubble
The simulation of a round static bubble is the first in
teresting test for the numerical method. Although very
simple in theory, it can give informations about the pre
cision of the interface treatment, in particular of the sur
face tension. A round bubble of a dense fluid is cen
tred in a square domain, surrounded by a lighter one;
there are no initial velocities and no gravity. In this case
the pressure jump is imposed by the interface bending
and the surface tension product, as given by the Laplace
equation
pint pext = a/R, if in 2D
Pint Pext = 2a/R, if in 3D
Any error resulting from the discrete computation of the
level set bending leads to an erroneous value of the pres
sure gradient between two adjacent cells: this non zero
pressure gradient then acts as a wrong correction for the
predicted velocity field. This the cause of the so called
spurious currents, non zero components of the velocity
that arise from the first iteration, and may in some cases
destroy the interface. A comparison between two dif
ferent test cases is proposed: in the first mesh conver
gence results from some uniform meshes are shown; in
the second one the same computation is performed with
an adaptive mesh localized on the interface. The param
eters of the simulation are taken from F. Coudere (2006):
L 4 ,cm, R 1 cm, pint 1000 kg.m3, pext=
1 kg.m3, Pint 103 Pa.s, Pext 10sPa.s,
a 0.1 N.m Only one iteration of level set reini
tialization is performed in order to reduce any artificial
displacement. In the first table the norm of he velocity is
calculated after one time step alone: here the predicted
velocities are equal to zero, so that the spurious currents
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Mesh
Uniform
Grid
16
32
64
128
256
16
32
64
128
lI .I
1.933e06
1.910e07
2.485e08
3.277e09
3.893e10
1.933e06
2.125e07
3.259e08
5.717e09
AMR
2562 9.451e10 2.60
Table 1: Spurious currents norm, first iteration. AMR
grid correspond to the finest level.
little more than in the uniform mesh; still it does not rise
dangerously, but it stabilizes as for the uniform mesh.
The convergence rates in table 2 are now near to 1.5.
Figure 4: Pressure field of a static round bubble after one
timestep captued with four grid levels
are determined by the corbure computation.
If the analytical bending is used as input, then the solu
tion is always precise to the roundoff error. As shown in
table 1, the convergence rate, given by
In (errg,)  In (err,)
In (2)
are, for the uniform mesh computation, close to 3:
this behaviour, justified by Oevermann and Berger
Michael Oevermann (2000), comes from the zero gradi
ent solution of the pressure on both the sides that nullify
the effects of the interface position relative to the grid.
In the adaptive mesh test a slight decrease of the code
performance, mainly for the finer resolutions, is visible.
It has to be kept in mind that there is no lower limit to
the refinement level, so that the lower resolution are kept
far from the interface as the resolution increases. The ef
fects of the coarser grids and the interpolation are visible
in the small spurious current increase. Still, the conver
gence rate never drops under 2.5.
The growth of the spurious currents are here checked
for a prolonged time of simulation: in this case the code
tries to achieve an equilibrium for those nonphysical ve
locities, with the generation of non zero vorticity near
the interface. An initial steep growth followed is fol
lowed by a stabilization in time. In the AMR compu
tation the initial error grows, as for the first iteration, a
Mesh
Uniform
AMR
Grid I
162 1.885E03
322 5.828E04 1.69
642 1.778E04 1.71
1282 3.879E05 2.20
2562 1.060E05 1.87
162 1.885E03
322 5.676E04 1.73
642 2.181E04 1.38
1282 5.545E05 1.98
2562 2.953E05 0.91
Table 2: Spurious currents norm, t = 1
correspond to the finest level.
s. AMR grid
Rising bubble
In this section the dynamics of a rising round bubble are
presented. The bubble is now immersed into a heavier
fluid and subjected to gravity. There are no analytical
solutions to compare; however, extensive research work
has been done in creating benchmarks for this problem,
mainly from Mark Sussman (2002) and Hysing et al.
(2008). In this sense, a convergence test for some com
puted quantities and a comparison with the Hysing's are
performed.
In the initial configuration there is a column of fluid 1
with no slip conditions imposed on the horizontal sur
faces and slip conditions on the vertical ones. In the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Time [s]
Figure 5: Evolution of the spurious currents for a long
time simulation, t = 1 .s
lower half of the [1x2] domain, centered at [0.5, 0.5] is
the circular bubble of radius [r = 0.25] composed of a
lighter fluid, all the velocities initialized to zero. As the
gravity applies its effect, the bubble starts rising, and at
the same time begins to change its shape. As in the Hys
ing's benchmarks, computations are performed for two
configurations.
Two dimensionless numbers can describe the regime
of the bubble rise. Given L 2rn and LT, ~
these are the Reynolds and the Eiitviis number:
Figure 6: Initial condition for the rising bubble test
This is an interesting quantity to compute, because it
does not only measure how the interface tracking algo
rithm behaves, but also it gives an idea of the quality of
the overall solution; this quantity checks also the reduc
tion of precision of the AMR grid far from the interface,
and its influence on the interface itself. The last quantity
is the circularity, an index of how much the shape of the
bubble differs from the circle:
pb
Here the numerator is the perimeter or circumference
of a circle with diameter cl,, which has an area equal to
that of a bubble with perimeter Pb: it always starts from
the maximum value of one, an then decreases as the bub
ble begin to be deformed. Some qualitative comparisons
between our results and Hysing's ones are presented in
figure 9; for the quantitative results, the results of the re
fined mesh are compared with a reference computation
made with an uniform fine mesh of [256x1024]. The dif
ference between the two values is calculated at the final
timestep: e rr = ahbs(qt=3 9,ef ,t=3).
The temporal evolutions of the bubble values follow
the benchmarks for both the configurations; however he
code could not track the circularity for the second test
case after t 1.5 .<, approximately when the bubble
forms the two sharp corners. In this location Hysing
could find, with a finite elements code, a small ligament
that subsequently gives birth to two small drops.
In the first simulation set the convergence towards the
reference solution is very good. The position of the cen
Eo =
Re =
In the first case the surface tension should be enough
to hold the bubble together, and the final shape should be
of an ellipsoidal regime, ; the density and viscosity ra
tio are equal to 10. The second test is more challenging,
as the density ratio reaches 1000 and the viscosity ratio
100; this bubble lies somewhere between the skirted and
dimpled ellipsoidalcap regimes indicating that break up
can possibly occur (Clift et al. (1978)). The evolution
of the two bubbles is tracked for three time units, during
which three quantities are measured. The first is the cen
troid of mass (.re, yc), of which, given the symmetry of
the problem, the y coordinate only is considered:
The second is the instantaneous rise velocity (u, v),
again y component only:
SI22 V cl#
foP cy
o "
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Test case
1
#1g #2 y Re
10 1 0.98 24.5 35
P1/# #/#
10 10
P1 #2
1000 100
2 1000 1
10 0.1 0.98 1.96 35 125 1000 100
Table 3: Physical parameters for the two test cases
Mesh value err
Ye 32 1.0894 4.03E3
64 1.0863 9.39E4
1.0856
1.0854
0.9313
0.9233
0.9206
0.9212
0.1986
0.1964
0.1955
0.5917
2.39E4
1.01E2
2.16E3
5.17E4
2.86E3
6.98E4
1.57E4
Vcise
Table 4: Computed values for the rising bubble, test 1.
Mesh value refers to the finest level.
(a) Solution test 1
(b) Solution test 2
Figure 7: Solution after t = 3 .s of the two configurations
of bubble, four levels of mesh up to 128x256
Mesh value
Ye 32 1.1619
64 1.1525
128 1.1477
ref 1.1434
32 0.8630
64 0.8532
128 0.8443
ref 0.8389
err
1.86E2
9.21E3
4.42E3
2.41E2
1.43E2
5.38E3
(b) Hysing, Re = 35, Eo
125
(a) Hysing, Re 35, Eo =
10
Figure 8: Zoom over the bubble shape, Hysing
ter of mass is almost the same for all the meshes, and
grows almost linearly in time; the convergence rate is
quadratic for this value. The circularity computation is
more difficult: the lowest mesh computed value seems
quite far from the good solution, mainly in the last half
of the simulation; still there are similar values of q. The
rise velocity seems a little bit underestimated, for the
lesser mesh, in the region around the maximum; the
terminal velocity is however well estimated for all the
grids, and q keeps quadratic.
The second test case is clearly more difficult. Again
the center of mass rises linearly in time; however, we
have linear convergence rate for its final position. The
0.2105 5.61E3
0.2125 3.59E3
Vcise
0.2144
0.2161
1.66e3
Table 5: Computed values for the rising bubble, test 2.
Mesh value refers to the finest level.
circularity tracking is lost for all the meshes after ap
proximately 1.75 .s. The instantaneous mean rise veloc
ity improves quite a lot with the finer meshes: with a
resolution of at least Ar: = 1/64 the second velocity
peak is found around t = 2 .s. For those values the con
vergence rate is linear or sublinear.
Time [s]
(b) Center of mass, test 2
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
r v 1/32
1/64
+ Hysing [1/320] 
1
~0.8 
0.6
r r1/32
v 1/64
+ Hysing [1/320]
Y,
.~.I
Time [s]
(a) Center of mass, test 1
0.8
AA 1/32
r v1/64
1/128
0.6 + Hysing [1/320]
4,
Z +i.~
. 1/128
+ Hysing [1/320]
1 2
Time [s]
(c)( l C I .. 1 I test 1
0.98 E
.
I,
. 0.96
Time [s]
(d)(l CI ...ni .. test 2
02 i h.*~*Ysys..r.~~;;;~~':;I
0.15 
0.05
0.0
0o
0.15
n .
 1/32
+ Hysing [1/320]
Time [s]
(e) Rise velocity test 1
Time [s]
(f) Rise velocity test 2
Figure 9: Comparison of the temporal evolution of three quantities
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Acknowledgements
This research project has been supported by a Marie
Curie Early Stage Research Training Fellowship of the
European Communitys Sixth Framework Programme
under contract number MESTCT2005020426.
The PARAMESH software used in this work was
developed at the NASA Goddard Space Flight Cen
ter and Drexel University under NASA's HPCC and
ESTO/CT projects and under grant NNGO4GP79G from
the NASA/AISR project.
References
M. J. Berger. Adaptive mesh refinement for hyperbolic
partial d!ifferentrial equations. PhD thesis, Stanford Uni
versity, 1982.
R. Clift, J. R. Grace, and M. E. Weber. Bubbles, drops
and particles. Academic Press, 1978.
F. Couderc and J. L. Estivalezes. Sharp interface captur
ing method based on levelset method for incompress
ible flow with surface tension. ASME Annual Summer
Meeting FEDSM200345151, 2003.
Frederic Couderc. Developpement d'un code de calcul
pour la simulation d'e'coulements de fluides non misci
bles Application ia la disintegration assiste'e d'un jet liq
uide par un courant gazeux. PhD thesis, ENSAE, 2007.
D. DeZeeuw and K. G. Powell. An adaptative cartesian
mesh for the euler equations. Journal of Computational
Physics, (245), 1992.
S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Bur
man, S. Ganesan, and L. Tobiska. Quantitative bench
mark computations of twodimensional bubble dynam
ics. MOXReport, (23), 2008.
John B. Bell Phillip Colella Louis H. Howell Michael
L. Welcome Mark Sussman, Ann S. Almgren. An adap
tive level set approach for incompressible twophase
flows. Journal of Computational Physics, 148, 1998.
Y. Yoshida Mark Sussman, M. Ohta. Numerical simula
tion of bubble motion by a coupled level set/vof method.
Proceedings of the 35th Annual Meeting of the Society
of Chemical Engineers, 2002.
S.F. McCormick and J. Thomas. The fast adaptive com
posite grid (fac) method for elliptic equations. Applied
Mathematics and Computation, (46), 1986.
Germany Marsha Berger Michael Oevermann. A pro
jection method for twophase incompressible flow with
surface tension and sharp interface resolution. Zuse In
stitute Berlin Report, (17), 2000.
RayleighTaylor instability
In this section the simulation of a RayleighTaylor insta
bility is presented. This phenomenon can be observed
by superimposing a heavy fluid over a lighter one; with
out any perturbation they keep in equilibrium. An initial
sinusoidal perturbation of a given wavelength is subject
to the destabilizing effect of gravity and the stabilizing
effect of surface tension: if the destabilizing effect pre
vails, then the perturbation grows, as the heavier fluid
moves down and deplaces the lighter one. The simu
lation is performed in a L, 1 m, L, 4 m do
main. The interface is placed horizontally in the mid
dle of the domain. The computational parameters are:
ph = .22o kg.m ,pi 0.1694 kg.m ,PA = t
3.13E 3 kg.m .sl ,o 0 N.m The initial per
turbation is y 0.05cos(kx) with k 2x7. Those pa
rameters are the same used by Popinet Popinet (2000):
into his work a comparison is made between an eule
rian VOF method and a lagrangian type interface track
ing. In this simulation the interface sustains strong non
linear deformations, the initial perturbation grows into
a "mushroom" shape until 0.8 seconds. At this time the
extremities begin to form two thin ligaments that quickly
shrink to the grid cell size. With insufficient resolution
the ligaments are broken and, subsequently, two drops
are detached, so that the computed velocity field deviates
quickly from the good solution. The precise lagrangian
solution shows us that until t = 0.9 s the ligaments are
kept. We have performed the first computation with the
same base resolution, 64 x 256 with four levels of AMR.
We can see that at the last timestep the ligaments are al
ready somewhat fragmented, even if the shape remains.
We have subsequently tried to improve the interface res
olution by adding an AMR level, reaching a fine reso
lution of 128 x 512: in this simulation the ligaments are
clearly maintained, showing us the improvement on of
the biphasic code by a local refinement over the inter
face.
Conclusions
A second order parallel adaptive projection method for
incompressible two phase flows has been developed.
The ability to perform high density ratio computations
has been demonstrated, in particular the ability of the
elliptic solver to converge for all the test cases. In this
paper it is proven that a local interface refinement AMR
simulations can nearly achieve the fine grid precision in
all the tests, as well as maintaining good convergence
rate. An extension to three dimension will soon be done,
as the benefit of the AMR would be even more evident.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 10: Results from the RayeighTaylor instability computations, five levels of mesh, finest grid equivalent to
64x2.56. t = 0.9 s
Figure 11: Results from the RayeighTaylor instability computations, four levels of mesh, finest grid equivalent to
128x512, t = 0.9 s
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Kevin Olson. Paramesh: A parallel adaptive grid tool.
in Parallel Computational Fluid Dynamics 2005: The
ory and Applications: Proceedings of the Parallel CFD
Conference, College Park, MD, U.S.A., eds. A. Deane,
A. Ecer, G. Brenner, D. Emerson, J. McDonough, J. Pe
riaux, N. Satofuka, and D. TromeurDervout (Elsevier),
2006.
Kevin Olson and Peter MacNeice. An over of the
paramesh amr software and some of its applications.
in Adaptive Mesh RefinementTheory and Applications,
Proceedings of the ( Ill. .G,. Workshop on Adaptive Mesh
Refinement Methods, Series: Lecture Notes in Computa
tional Science and Engineering, eds. 7 Plewa, 7 Linde,
and G. Weirs (Berlin: Springer), 41, 2005.
Clark Mobarry Rosalinda deFainchtein Peter MacNeice,
Kevin M. Olson and Charles Packer. Paramesh : A par
allel adaptive mesh refinement community toolkit. Com
puter Physics Communications, 126, 2000.
Sthephane Popinet. Stability' et formation de jets dans
les bulles cavitantes : De'veloppement dune mi'thode de
chaine de marqueurs adapti'e au traitement nume'rique
des equations de NavierStokes avec surfaces libres.
PhD thesis, Universit6 Pierre et Marie Curie Paris VI,
2000.
M. Sussman. A parallelized, adaptive algorithm for mul
tiphase flows in general geometries. Computers and
Structures, 2005.
R. Teigland and I. K. Eliassen. A multiblock/multilevel
mesh refinement procedure for efd computations. Inter
national journal for numerical methods in fluids, 2001.
H. A. VanDerVorst. Bicstab: a fast and \1n... .II1h con
verging variant of bicg for the solution of nonsymmetric
linear systems. SIAM Journal of Scientific Computing,
1992.
C. Vuik and P. Wesseling J.J.I.M vanKan. A black box
multigrid preconditioner for second order elliptic partial
differential equations. European Congress on Compu
tational Methods in Applied Science and Engineering,
2000.
