7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Simulating multiphase flows with multilength scales using moving mesh interface
tracking with adaptive meshes
Shaoping Quan
Large Scale Complex Systems, Institute of High Performance Computing, #1616 Connexis, Singapore 138632
Email: quansp@ihpc.astar.edu.sg
Keywords: Moving mesh interface tracking, mesh adaptation, mesh combination, mesh separation, droplets,
Arbitrary Lagrangian Eulerian method (ALE), boundary fitted mesh
Abstract
Dropletdroplet interaction and droplet creation are two of the most important phenomena in multiphase flows. The
length scales in the systems are usually much different, for example, when two drops are touching or a drop is about
to pinch off, the length scale of the thin film or the thin thread could be hundreds or thousands times smaller than
the drop radius. It is a challenge to numerically solve these problems with high fidelity, especially for the numerical
methods in which the fluids' properties are smoothed in a finite thickness interfacial region. In this paper, the moving
mesh interface tracking method (MMIT) with mesh adaptation by Quan and Schmidt (2007) is employed to address
this challenge. The interface is represented as surface triangle meshes and thus is zero thickness, therefore the jump
conditions across the interface are directly implemented without smoothing of fluids' properties. The interface moves
with fluid velocity in a Lagrangian fashion. To maintain good mesh quality and to resolve the different length scales as
well as to achieve computing efficiency, local and dynamics mesh adaptations, such as edge contraction, edge bisection
and edge swapping, are applied. To allow topological transitions, mesh combination and mesh separation (Quan et al.
2009a) are employed.
The interaction between a droplet and a substrate is simulated, and the thin film thickness is around 1.5% of the
diameter of the spherical droplet. To investigate the detailed physics of droplet necking and pinchoff, especially near
the moment just before the breakup, simulations on the relaxation and pinchoff of an elongated droplet are performed
with the smallest radius of the thin thread around 2% of the initial radius. The numerical scheme well predicts the
evolution pattern of the thin thread radius, and the results agree well with recent works by Burton et al. (2005) and
Thoroddsen et al. (2007).
Nomenclature
(You can add a nomenclature section with the tabbing
environment should you wish to do so.)
Roman symbols
CV control volume
CS control surfaces
p pressure
u fluid velocity
v moving mesh velocity
n unit normal vector of a control surface
Greek symbols
p density
p viscosity
a surface tension coefficient
K effective curvature
7 unit tangential vector of a surface
Subscripts
d droplet phase
s surrounding phase
Superscipts
T transpose
Others
I II difference across the interface
Introduction
Numerical simulation of multiphase flows is very chal
lenging not only because of the singularity of great fluid
property difference across a zerothickness interface but
also because of the different length scales involved in
the multiphase flow systems. The smallest length scales
could be thousands times smaller than the characteris
tic length scale for the overall systems, such as the ra
dius of a bubble or drop. These small length scales usu
ally appear in the droplet (bubble) breakup as the forma
tion of a neck region before the actual pinchoff and also
the coalescence of drops (bubbles) as a thin film formed
when the drops are so near to each other but before colli
sion. However, it is a challenge for numerical methods,
such as front tracking (Tryggvason et al. 2001), level set,
and volume of fluid (Scardovelli and Zaleski 1999), etc.
where the fluids' properties are usually smoothed in a in
terfacial region with a thickness of several cell sizes, to
compute these flows with high fidelity to investigate the
underlying physics in these small lengthscale regions.
However, the moving mesh interface tracking
(MMIT) method which is developed very recently by
Quan and Schmidt (2007) could be capable to deal with
this multilength scale problem in high fidelity within
the continuum assumption. In MMIT, the interface is
represented by cell edges (in 2D) or faces (in 3D), and
thus exact zerothickness. The interior meshes for the
two phases are connected by these interface elements
and formed to a single set of mesh to solve the governing
NavierStokes equations. As time progresses the inter
face elements move in a Lagrangian mode to naturally
conserve the volume (the mass) of each phase. This
kind of mesh configuration is also called boundary fit
ting mesh (Xu and Basaran 2007; Tsamopoulos et al.
2008), and the way of the motion of the mesh is also
named arbitrary LagrangianEulerian method (ALE) in
fluidsolid interactions (Feng et al. 1994). On the other
hand, mesh motion on the interface deteriorate the mesh
quality, and if no measures were taken, the simulation
will have great errors, or even worse, the simulation dis
continues. In order to maintain good mesh quality, a
number of mesh adaptation schemes, including 32 and
23 flips, edge bisection, and edge contraction, etc. are
implemented in MMIT and applied locally. These adap
tive mesh methods also help to achieve computing effi
ciency by distributing fine mesh in the region of inter
est and coarse mesh elsewhere. Actual breakup and the
coalescence of interface(s) are mimicked by the mesh
separation and mesh combination schemes (Quan et al.
2009a), respectively.
Governing equations, numerical methods
The fluids of multiphase flows are assumed incompress
ible and immiscible. The unsteady motion of the mul
tiphase flows is governed by the continuity and the
NavierStokes equations, here written in integral form
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and in a moving control volume scenario,
S I
cv
CV
cv
cv
C
cv
dv = jv nds,
cs
pdv + JJ p (u v) nds =
cs
cs
pudv + 1 pu (u v) nds
CS
pfdv pnds
Cs
+ ff(Vu+ Vu) nds.
cs
(1)
0, (2)
(3)
(4)
Equation (1) is a consistency requirement for the volume
of a moving and deforming control volume. In MMIT,
the interface moves with the fluid, i.e. u v, so the
jump in the fluid's densities across the interface presents
no numerical difficulties for computing the convection
terms in (2) and (4), as these terms identically equal zero
on the interfacial surfaces.
The jump condition for the normal stresses and the
continuity of shear stress across the interface read
[pH K + [2/ (Vu n) n.,
S((Vu + VU) n). = 0.
In the numerical method, the curvature of a surface is
calculated by a leastsquares parabola fitting in a lo
cal coordinate system; the jump in the normal viscous
stresses is implemented directly using a Taylor series ex
pansion for multivariable functions; the continuity con
dition of the stresses is satisfied by using a geometric
harmonic mean method, and this method is secondorder
accurate, because it calculates the exact shear stresses
for linear velocity distributions. The details and some
validations can be found in Dai and Schmidt (2005);
Quan and Schmidt (2007, 2006).
Mesh adaptation for small length scales
To obtain good mesh quality, a number of mesh adap
tation schemes are implemented and these shemes are
mainly based on the mesh quality metric (Knupp 2003),
local interface curvature, edge angles, dihedral angles,
averaged length of the initial configuration. However,
for cases where two flat interfaces, or an interface and
a solid wall are very near to each other forming a very
thin neck or a very thin film, the above mentioned mesh
adaptation criteria are not efficient for mesh refinement
in these smalllengthscale regions. Fig. 1 shows meshes
of a case where a droplet is approaching a substrate,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 1: Interior mesh of the droplet and the surround
ing phases for a droplet approaching a substrate. Blue:
droplet; black: the surrounding fluid. (a) overall view
and (b) enlarged view of the region between the droplet
and the substrate
where the mesh of the droplet phase is. For this simula
tion, all the mesh adaptation schemes based on the above
mesh criteria are employed. From Fig. 1 (a), it can be
seen that the mesh adaptation works well, as the mesh
near and inside the droplets is very fine and the mesh
far away from the droplet is coarse. However, from the
enlarged view of the region between the droplet and the
substrate, it is concluded that the mesh in this region is
too coarse to obtain reliable information for examine the
underlying detailed physics, as there are only two cells
in between the interface and the substrate for most of the
region. Apparently, to address this issue, extra length
criteria are needed.
From Fig. 1 (b), it is clearly that the curvature of the
interface cannot be used for the length scale to refine
the mesh, as the mesh on the interface is already fine
enough to capture the curvature changing. The averaged
length of the interior mesh in the suspending fluid can
neither be a good criteria, as the mesh in this region is
already much finer than other region. However, it is ob
served that unit normal of the interface and the substrate
Figure 2: Boundary condition for the Poisson equation
are changed dramatically over this small gap, if the nor
mal of the interface surfaces points inward to the droplet
while the normal of the substrate surfaces points away
from the simulation domain. So, we could solve a Pois
son equation with the surface normal of the interface and
outside boundaries as the boundary conditions (see Fig.
2, i.e.
V2yi 0 with i 1,2,3 (7)
4(interface) n(inter face)
boundaries) = n(boundaries),
where 0 is a vector and has three components in three
dimensions, and n is the unit normal vector. This equa
tion is solved by a conjugate gradient method, and the
solution is shown in Fig. 3. It is observed that the varia
tion of 4 away from the gap region is smooth, while the
change in 4 in the gap region is very sharp. There are
some edges on which the magnitude of the difference of
4 between two nodes is even larger than unity. There
fore, the magnitude of the difference of 4 between the
two nodes can be served as a criteria for edge bisection.
Fig. 4 (a) shows the mesh after the mesh bisections is
applied, and it can be seen that the mesh in the gap is
much finer with at least four cells in the interface nor
mal direction. The vectors of p are also displayed in
the figure, and the magnitudes of the difference of the
vectors between two nodes in the gap are also less than
unity. The distribution of 4 is more smooth. However,
it is noted that the quality of the meshes connected to
the interface is not good, as the edge on the interface is
at least two times longer than edges in the gap. To im
prove the mesh quality, interface edge bisections which
is based on the averaged length of interior edges in the
gap sharing the two nodes with the interfacial edge are
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0 0.2 0.4 0.6 0.8 1
Figure 3: (a) Contour of the solved normal vector mag
nitude; (b) vector distribution of the solution
applied, and the mesh is displayed in Fig. 4 (b). It is
clear that mesh uniformity is much better.
Results and discussion
To demonstrate the capability of the MMIT method in
capturing detailed physics at small length scales, two test
cases are simulated which include interaction between
droplet and the substrate and the relaxation and pinch
off of an elongated droplet.
Fig. 5 show the initial velocity distribution of a
droplet approaching a substrate, i.e. the boundary at the
bottom, in which the velocity of the droplet centroid (in
y direction) is perpendicular to the substrate. A vortex
ring is observed near the droplet which is denoted by
Figure 4: (a) mesh and vector of 4 after edge bisections
for the gap region; (b) mesh after edge bisections for the
interfacial edges and edges connected to the interface
the circle, and the center of the vortex on each side of
the droplet is outside of the droplet. The density ratio
(r) and viscosity ratio (A) between the droplet and the
surrounding fluids are 50, respectively. The Reynolds
number (Red) based on the droplet properties, droplet
initial diameter, and the droplet centroid velocity is 50,
and the Weber number (Wed) based on droplet is 25.
As the droplet approaches the substrate, the droplet will
deform, and the Fig. 6 shows the shape evolution for
the droplet at t =0, 0.25, 0.5, 7.5, 1.0, 1.25, 1.5, and
1.56, where t is nondimensionalized by 2ro/Uc with ro
being the radius of the initial spherical droplet and Uc
the initial centroid velocity of the drop. As time pro
gresses, the droplet first become a tear shape, and when
the droplet become very near to the substrate, the front
of the droplet become flattened. At t = 1.56, there is a
thin film formed between the droplet front and the sub
strate. To clearly view the thin film, an enlarged view
of the film with meshes for the droplet (blue) and the
surrounding flow (red) is plotted in Fig. 7. The small
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
est width of the thin film is around 1.5% of the diameter
of the spherical droplet, and there are at least three cells
in the thin film region in the width direction. This gives
confidence that reliable fluid dynamics in the film can be
obtained.
I. _ 
Figure 5: Initial velocity field for the headon droplet
substrate interaction. The circle inside the domain rep
resents the droplet and the black bold straight lines are
the walls with the one at the bottom mimicking the sub
strate. The time series from left to right and from the top
to bottom is 0, 0.25, 0.5, 7.5, 1.0, 1.25, 1.5, and 1.56.
0000
OOO
OQOO
Figure 6: Shape evolution for the headon droplet
substrate interaction
Fig. 8 shows the simulated shape evolution of an ini
tial elongated droplet in a suspending fluid with finite
viscosity. The density and viscosity ratio between the
droplet and the suspending fluids are 10 and 0.1, respec
tively. The Ohnesorge number (Ohd = d/Jpdgro)
based on the droplet properties is 0.037, and the length
ratio between the maximum length of the ligament to
the radius of the middle cylinder is 20. Initially, the
=
\\
Figure 7: Enlarged view of the thin film at t 1.56 with
meshes
COO
Oco
C^3
00
00
Figure 8: Shape evolution of the relaxation and pinch
off of an initially elongated droplet with length ratio
of 20, density ratio of 10, and viscosity ratio of 0.1.
The time sequence is 0, 435,3, 489.6, 509.8, 511.2, and
613.7.
fluids are at rest, however, due to the unbalanced sur
face tension forces, the ligament will contract as the sur
face tension forces serves as a restoring force. For this
case, there are two neck region formed and finally three
droplets are created with two primary droplets at the two
ends and one satellite drop in the middle. In the figure,
the time is nodimensionalized by t, rotd/Oa, where
ro is the radius of the initial cylinder at the middle sec
tion, Pd stands for the droplet viscosity, and a denotes
Figure 9: Enlarged view of mesh in the neck region of
the drop at t = 509.8
r*uI
U
tII I ..l l l .ll l l .ll II I
20 40 60
Figure 10: Neck radius evolution for the case in Fig.
8. Symobls: numerical results; line: curve fitting with a
slope of 0.89
the surface tension coefficient. However, the neck at
t 509.8 in Fig. 8 cannot be clearly observed as it
is so tiny. An enlarged view of the neck with meshes
is displayed in 9. The mesh in the neck region is fine
enough to capture reliable physics, and the mesh far
away from the neck is much coarser to achieve comput
ing efficiency. The smallest neck radius is around 2% of
the radius of the initial cylinder, however compared to
the radius of the primary droplet created after breakup;
it is less than 1%. Stone and Leal (1989) performed sim
ilar experimental work, and full threedimensional sim
ulations with parametric studied were report by Quan
et al. (2009b). The evolution of the smallest radius in
the neck region is displayed in Fig. 10. The time in the
figure is computed by t = tb t, where t, is the ac
tual simulation time and tb is the time for the breakup
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
instance, and the time here is nondimensionalized by
tc. Our simulation shows that the reduction of radius
follows a power law with respect to time, which agrees
with the previous published results (Burton et al. 2005;
Thoroddsen et al. 2007).
Conclusions
The capability of the moving mesh interface tracking
(MMIT) method in capturing detailed physics in tiny
length scales in the thin film region or the neck re
gion has been demonstrated. Furthermore, by using
the mesh adaptation schemes, especially the edge bi
section based on the distribution of the normal vectors
and the averaged length, MMIT is capable to distribute
mesh with enough resolutions in the regions with small
length scales, while much coarse mesh are in the region
away from the smalllengthscale regions. In such a way,
MMIT is robust in solving multiphase flows with differ
ent length scales (the ratio between the largest charac
teristic length and the smallest characteristic length can
be as large as 100) and at the same time the computing
efficiency is well achieved.
References
S. P. Quan and D. P. Schmidt. A moving mesh interface
tracking method for 3d incompressible twophase flows.
J. Comput. Phys., 221:761780, 2007.
S. P. Quan, J. Lou, and D. P. Schmidt. Modeling merg
ing and breakup in the moving mesh interface tracking
method for multiphase flow simulations. J. Comput.
Phys., 228:26602675, 2009a.
J. C. Burton, R. Waldrep, and P. Taborek. Scaling and
instabilities in bubble pinchoff. Phys. Rev. Lett., 94(18):
184502,2005.
S. T Thoroddsen, T G. Etoh, and K. Takehara. Exper
iments on bubble pinchoff. Phys. Fluids, 19:042101,
2007.
G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al
Rawahi, W. Tauber, J. Han, S. Nas, and Y. J. Jan. A
fronttracking method for the computations of multi
phase flow. J. Comput. Phys., 169:708759, 2001.
R. Scardovelli and S. Zaleski. Direct numerical simula
tion of freesurface and interfacial flow. Ann. Rev. Fluid
Mech., 31:567603, 1999.
Q. Xu and O. A Basaran. Computational analysis of
dropondemand drop formation. Phys. Fluids, 19:
102111,2007.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
J. Tsamopoulos, Y. Dimakopoulos, N. Chatzidai,
G. Karapetsas, and M. Pavlidis. Steady bubble rise and
deformation in Newtonian and viscoplastic fluids and
conditions for bubble entrapment. J Fluid Mech., 601:
123184,2008.
J. Feng, H. H. Hu, and D. D. Joseph. Direct simulation
of initial value problems for the motion of solid bodies in
a newtonian fluid. Part 1. Sedimentation. J. FluidMech.,
261:95134, 1994.
M. Z. Dai and D. P. Schmidt. Adaptive tetrahedral mesh
ing in freesurface flow. J. Comput. Phys., 208:228252,
2005.
S. P. Quan and D. P. Schmidt. Direct numerical study of
a liquid droplet impulsively accelerated by gaseous flow.
Phys. Fluids, 18:102103, 2006.
P. M. Knupp. Algebraic mesh quality metrics for un
structured initial meshes. Finite Elem. Anal. Des., 39:
217241, 2003.
H. A. Stone and L. G. Leal. Relaxation and breakup
of an initially extended drop in an otherwise quiescent
fluid. J. FluidMech., 198:399427, 1989.
S. P. Quan, D. P. Schmidt, J. S. Hua, and J. Lou. A nu
merical study of the relaxation and breakup of an elon
gated drop in a viscous liquid. J. FluidMech., 640:235
264, 2009b.
