7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Characterization of droplet dynamics in a bifurcation channel
A. Carlson*, M. DoQuang* and G. Amberg*
Department of Mechanicics, Linne Flow center, The Royal Institute of Technology, Stockholm, Sweden
andreaca@mech.kth.se, minh@mech.kth.se and gustava@mech.kth.se.
Keywords: Droplets, twophase flow, unstable flow, Phase Field, three dimensional simulations
March 24, 2010
Abstract
We present here a phenomenological description of droplet dynamics in bifurcating channels that is based on
threedimensional numerical experiments. Droplet dynamics is investigated in a bifurcating channel, which has
symmetric outflow conditions in its daughter branches. For the parameter spaces explored here, we find two distinct
flow characteristics as the droplet interacts with the junction; splitting or nonsplitting. A flow map based on the
initial droplet size and the Capillary number characterizes the two flow regimes of splitting and nonsplitting droplets.
Close to the threshold between these regimes, the RayleighPlateau instability is identified as a driving mechanism for
the droplet breakup.
embolotherphy treatment.
Flow in bifurcating channels takes place in many
phenomena in Nature, as the transport of nutrition in
plant leaves (Katifori et al. (2010)) and blood flow
in veins and capillaries (Bull (2005)). Bifurcating
channels are also an unavoidable part in applications
in the medical technology, microsystem technology as
well as for transport of largescale quantities as oil and
gas from refineries and human waste in sanitary sewer
systems. A bifurcating channel consists of a primary
or parent channel that separates into two or more
secondary channels, also named daughter channels.
The flow as it meets the point of bifurcation, needs to
divide its mass flow into the daughter branches. The
present paper is devoted to the study of an unstable
flow phenomena at small scales as a droplet meets
the junction in a bifurcating channel, where special
emphasis is placed on understanding what governs the
splitting or nonsplitting behavior of the droplet in the
junction.
Droplets are an important part in microsystem tech
nology and medical technology. They have been used
in a wide range of application as both a passive and
active parameter. How these droplets behave in such
complex networks of channels can be crucial in order
to promote desired processes. Two important examples
of such processes is the mixing at microscale and
One avenue of droplet microfluidics is the use of
droplets as compartments for a desired physical phe
nomenon such as mixing, as shown by Garstecki,
Fischbach, and Whitesides (2005). Mixing in micro
scale is a wellknown obstacle, since it is merely
driven by molecular diffusion. One direction to follow
in order to foster mixing is emulsion, which can be
realized by design of dropletdroplet interaction. The
potential use of an emulsion technique is often limited
by the possibility to precisely control the droplets size
distribution.
Embolotheraphy is a medical treatment strategy that
exploits bubbles or droplets and their free surfaces (Bull
(2005)), where a bubble is introduced in order to occlude
the blood flow to certain parts of the tissue. This method
is primarily applied as a last treatment resort for cancer,
when all other conventional clinical methods have failed
to eradicate the tumor. It relies on a precise control
of the twophase flow in the bifurcating capillaries, as
one seek to occlude the flow by lodging the bubble in
one of the junctions daughter branches. This leads to
a starvation of the tumor by excluding its oxygen supply.
In the present paper we will present three dimen
sional simulation results of droplet dynamics in a
bifurcation channel, based on the Phase Field theory.
We show that two different flow regimes can be ob
Introduction
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
served as the droplet interacts with the junction; splitting
or nonsplitting. A distinct criterion is found to define
the threshold between the two possible outcomes. Close
to the threshold that determines splitting or nonsplitting
droplets, the RayleighPlateau instability is found to be
a driving parameter for droplet splitting.
1 Mathematical formulation
1.1 Phase Field theory
The phase field theory is based on the thermodynamical
consideration of the free energy of a binary system. The
two components are here assumed to be immiscible
and separated by a narrow diffuse interfacial region.
By considering that two immiscible components ac
tually mix over an interfacial region van der Waals
(1893) proposed the idea of a diffuse interface. The
composition profile of the interface can be seen as the
competition between the random molecular motion and
the molecular attraction.
Cahn and Hilliard (1958) derived the free energy
by making a multivariable Taylor expansion about the
free energy per molecule,
f = 3W'(C) + I VC12 (1)
2
following here the phase field derivation and notation
by Jacqmin (1999) where C is an order parameter. The
creation of an interface is established by the competi
tion between the bulk free energy /3((C), and the in
terfacial energy QIVC]2 due to composition variations.
a and 3 are constants that comes out directly from the
Taylor expansion, see Cahn and Hilliard (1958), and
are proportional to the surface tension coefficient a and
the interface width c, 3 ~ and a ~ ac. The free
energy functional and the phase field parameters a,
,3 and the mobility M control the interfacial dynam
ics and width. is taken as a doublewell function,
p(C) I((C + 1)2(C 1)2, which will give the two
equilibrium states at C=+l. Integration over the vol
ume of the systems gives the total free energy defined
by F J= fdV. The functional derivative of F with
respect to the order parameter C gives rise to the chemi
cal potential,
6F
C 
6C
Bt'(C) aV2C.
By minimizing the chemical potential with respect to C
we obtain the equilibrium profile to the interface, here
given in one dimension as Co(x) = tanh( ). c =
/is the mean field thickness and the surface tension
is defined by the integral,
f dj (d)2 d
Joo \ dx )
2 32
3
By taking into account the effects of the fluids motion a
convectiondiffusion equation is obtained, also referred
to as the CahnHilliard equation
+ u vc = v (MVy),
where the mobility M is considered as a constant. We
require that there is no flux of the chemical potential
through the boundaries of the domain, which is fulfilled
by the Neumann boundary condition
i9d
O 0, (5)
On
where n is the normal vector to the boundary. The sur
face free energy contribution is postulated as (Carlson,
DoQuang, and Amberg (2009))
Fwaii = [IsL + (asv TSL)g(C)] dA, (6)
where T( ) is the surface energy between the three dif
ferent phases; liquid (L), gas (V) and solid (S). g(C)
0.5+0.75C0.25C3 is a smooth function between zero
and one, and its derivative g'(C) will be nonzero only
at the diffuse contact line. In eq.(6) it is assumed that the
interface is at or close to equilibrium as it wets the solid
surface. 00 is the equilibrium contact angle, formed be
tween the tangent of the interface and the solid surface,
given by Young's Law; cos(0o) sv L. Through
the variational derivative of the total free energy of the
system, with respect to C, and integration by parts, the
natural boundary condition for the concentration at the
wall is obtained,
On
a^ n + a cos(Oo)g'(C) =0 (7)
governing the diffusively controlled wetting at local
equilibrium (Carlson et al. (2009)).
1.2 Governing equations
In addition to the formulation of the free energy, the
(2) continuity and Navier Stokes equations are defined, here
given in nondimensional form,
Vu 0,
S9
S+u VC
at
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
p + (u V) u VP
+ V (/t(Vu+(Vu)T))
Re (
CVo (10)
Cn Ca) (
These equations have been made dimensionless with,
L ,
x =Lx*, t = t*, P pU2P*,
U
S= U*, 40= 3 *T
a Ua 72=p
where Lc is a characteristic length scale and U is the
reference velocity. This results in four nondimensional
numbers appearing in eq.(9, 10),
2v2UeL 6
Pe=3 Cn=
3Mcr Lc
pUL 2JpU
Re pUL Ca = 2U
P 3ua
The Peclet (Pe) number expresses the ratio between
advection and diffusion. The Cahn (Cn) number ex
presses the ratio between the interface width and the
characteristic length scale. The Reynolds (Re) number
expresses the ratio between the inertia and the viscous
force. The Capillary (Ca) number expresses the ratio
between the viscous and the surface tension force.
1.3 Computational methodology
The governing equations have been solved numerically
with the finite element toolbox FemLego (Amberg
et al. (1999)). FemLego is an open source symbolic
tool for solving partial differential equations. Within
one single Maple worksheet the partial differential
equations, boundary conditions and numerical solvers
are all defined. It inherits parallel computational and
adaptive mesh refinement capabilities see DoQuang
et al. (2007) for details.
The CahnHilliard equation has been solved with a
type of precondition Conjugate Gradient (CG) solver,
described by Villanueva and Amberg (2006). The
multifluid NavierStokes equations have been solved
with a projection scheme, that was proposed by Guer
mond and Quartapelle (2000). The linear system for
the velocities and pressure have been, respectively,
solved with the General Minimum Residual method
and a CG solver. All variables are approximated with
linear basefunctions and the mesh consist of triangular
elements.
Figure 1: Illustration of the Couette flow and a descrip
tion of the minor (B) and major (L) axis of a
droplet at steadystate.
1.3.1 Numerical validation
The dynamic behavior of the surface tension force has
been verified by validating the numerical solution of the
deformation of a threedimensional droplet in a Couette
flow. Taylor (1934) derived analytically the deformation
(D) of a droplet in a shear flow in the limit of low Re
numbers (Re < 1), which he also verified in experi
ments,
19A + 16
D 6 6a
16A + 16
LB
L+B'
Here the viscosity ratio of the two liquids is A, L and
B are the major and minor axes, respectively. The
Ca number Ca = is defined by the shear rate
/, t the viscosity of the continuous liquid and d the
droplet diameter. The deformation parameter D can be
extracted from experiments and numerical simulations
by measuring the major and minor axis when the droplet
has reached steadystate. A twodimensional sketch of
the domain and a schematic description of the droplets
minor (B) and major (L) axis is given in fig.(1). Initially
the droplet is spherical in the simulations and it deforms
into an elliptic shape.
Table(l) gives the relative difference between the
analytical and numerical results for a droplet in a
Couette flow for various Ca numbers. All results are
found to be in good agreement with the analytical
prediction from Taylor given in eq.(13).
1.3.2 Geometrical description of the domain
Droplet dynamics is investigated in a threedimensional
Yjunction described by the twodimensional sketch
given in fig.(2), where the zdirection goes into the
plane of the figure. A velocity profile with a shape of
a paraboloid, with a nondimensional mean velocity
u = 1, has been prescribed at the inlet of the parent
channel, and a Neumann boundary condition is defined
for the pressure. A symmetric outlet condition is
given for the pressure (P = 0) at the upper and lower
daughter branch. The parent channel has a quadratic
cross section L2, L being the width of the channel,
(a)
Figure 2: Geometrical description of the numerical do
main.
and the daughter branches have a rectangular cross
section L LB where L is in the zdirection. The
daughter to parent channel area is LL = 0.75 and
0 is the bifurcation angle. The droplet has initially a
volume Vi and the nondimensional volume is defined
as V = The walls are hydrophobic, having an
affinity of the continuous phase, with an equilibrium
contact angle of 0= 180 degrees. This prevents wetting
of the dispersed phase on the channel surfaces. The
nondimensional numbers given in eq.(12) are defined
with the characteristic length scale chosen as the width
of the parent channel Lc = L and the reference velocity
has been chosen as the mean velocity at the inlet.
The channel has an extension in the xdirection of
about 10L. A nearly equidistant mesh has been em
ployed between the inlet and the straight section (in
xdirection) of the branches, with few skewed cells near
the point of the bifurcation. 50 node points discretize L.
In the straight section of the branches where they are
aligned with the xdirection, has a coarse mesh with
L~ 10 node points in order to reduce the computational
time. The extension of the branches ensures that the
droplet behavior is not influenced by the outlet bound
ary condition in the simulations. The mesh consists
of nearly 7.5 million tetrahedron elements and the
computational time on 1024 processors for a single case
was about 24 hours.
2 Results
Fig.(3, 4) displays the threedimensional droplet
isosurface for C=0 and the velocity field in the cross
section, taken at the centre line in the channel, with a
normal [0 ex, 0 ec, 1 ec] at three snapshots in time. In
these three simulations the droplets have the same initial
shape, as a liquid slug, with the aspect ratio L 2, a
volume V 0.67 and Cn = 0.06.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Ca 0.15 0.1 0.08 0.06
Deo,,,, 2.8% 1.1% 0.8% 1.0%
Table 1: The relative deviation between the analytical
and computed deformation parameter is given
for different Ca numbers. All other parameters
have been kept constant Re=0.01, Pe 103,
Ca=0.06 and A = 1 and the mesh is equidistant
with a spacing Ax = 0.02.
Fig.(3a) shows the droplet as it has deformed at
the tip of the junction, forming two symmetric liquid
fingers in the upper and lower branch of the channel.
The flow from the parent channel continues to deform
the droplet, and a liquid thread is formed at the bi
furcation. This thread is thinning in time, as the two
liquid fingers propagate into the daughter branches.
Finally it becomes so small that surface tension starts to
contract the thread radially, an instability resembling the
RayleighPlateau instability. After breakup two equally
sized droplets are formed in each branch, that rapidly
obtains their energy minimizing shape.
By just changing the magnitude of the surface ten
sion coefficient a very different flow behavior than
observed in fig.(3a) is obtained, see fig.(4). In fig.(4a)
the droplet has just deformed at the junction. Work
from the continuous phase is transferred into surface
energy through droplet deformation, making the droplet
decelerate at the bifurcation. Finally it obtains a quasi
steadystate condition as it sticks at the bifurcation. The
incipience of a numerical disturbance initiates droplet
migration into the lower branch, see fig.(4c). Finally
the droplet solely occupies the lower branch, leaving an
asymmetric distribution of both the continuous as well
as the disperse phase in the junction branches.
The droplet splitting or nonsplitting phenomenon
depends mainly on the relative dominance of the surface
tension force, represented in the Ca number. The
resulting droplet dynamics seem to also highly depend
on its initial configuration, identifying the droplet size
and the Ca number as two of the key parameters for
the definition of the multiphase flow characteristics
in the junction. Similar experimental observations
have been made in both T (Link et al. (2004)) and
Ajunctions (MenetrierDeremble and Tabeling (2006);
Eshpuniyani et al. (2005)). These parameters form a
nondimensional space, which is explored in numerical
experiments, describing the relationship between the
splitting and nonsplitting flow regimes as shown in
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
~~

 
(c)
Figure 3: Droplet shape with V=0.67 and Ca=3.0 10 3
at three snapshots in time; (a)=2.5, (b)=2.7,
(c)=3.0. Figure 4: Droplet shape with V=0.67 and Ca=1.8 103
at three snapshots in time; (a)=2.0, (b)=3.5,
5 (c)=4.0.
(c)
(b)


111*""
~
2~
x~o*~""" I:~rc=

17
      
10
log(Ca)
Figure 5: Flow map of splitting (filled markers) and
nonsplitting (hollow markers), where the full
line represents the threshold between split
ting and nonsplitting found as V 0.8 
0.5log(Ca).
fig.(5).
A distinct and welldefined criterion defines
here the splitting or nonsplitting characteristics
V 0.8 0.5log(Ca), indicated by the full line in
fig.(5). These findings identify that large droplets favor
splitting and that such a process is hard to obtain with
small droplets.
3 Conclusion
We have reported here threedimensional numerical
experiments based on phase field theory of droplet
dynamics in a bifurcating channel with symmetric
outflow conditions. Two distinct flow regimes are
identified as the droplets interact with the junction,
splitting and nonsplitting. In particular we show the
effects of the initial droplet size and Ca number on the
resulting twophase flow characteristics.
Droplets that split equally, produce a symmetric
distribution of both phases in the channel daughter
branches. Near the threshold between the two regimes,
we observe that the RayleighPlateau instability can be
a driving parameter for droplet division.
In the nonsplitting regime the droplet migrates
into one of the channel branches, leading to a strong
temporal asymmetric flow in the junction. These results
are illustrating that a symmetric placement of the droplet
in the parent channel is highly unstable. A more detailed
description of droplet dynamics in a bifurcating channel
K* t
Splitting
*
N^4
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
by the respective authors can be found in Carlson et al.
(2010).
References
G. Amberg, R. Tonhardt, and C. Winkler. Finite ele
ment simulations using symbolic computing. MATHE
MATICS AND COMPUTERS IN SIMULATION, 49(4
5):257274, Sep 1999. ISSN 03784754.
Joseph L Bull. Cardiovascular bubble dynamics. Crit
Rev Biomed Eng, 33(4):299346, 2005. ISSN 0278
940X (Print).
John W. Cahn and John E. Hilliard. Free energy of
a nonuniform system. i. interfacial free energy. THE
JOURNAL OF CHEMICAL PHYSICS, 28(2):258267,
1958. doi: 10.1063/1.1744102.
A. Carlson, M. DoQuang, and G. Amberg. Droplet
dynamics in a bifurcating channel. INTERNATIONAL
JOURNAL OF MULTIPHASE FLOW, 2010.
Andreas Carlson, Minh DoQuang, and Gustav Am
berg. Modeling of dynamic wetting far from equilib
rium. Physics of Fluids, 21(12):1217014, 12 2009.
M. DoQuang, W. Villanueva, I. SingerLoginova, and
G. Amberg. Parallel adaptive computation of some
timedependent materialsrelated microstructural prob
lems. BULLETIN OF THE POLISH ACADEMY OF
SCIENCESTECHNICAL SCIENCES, 55(2):229237,
Jun 2007. ISSN 02397528.
B. Eshpuniyani, J. B. Fowlkes, and J. L. Bull. A bench
top experimental model of bubble transport in multiple
arteriole bifurcations. INTERNATIONAL JOURNAL OF
HEAT AND FLUID FLOW, 26(6):865872, 2005.
P. Garstecki, M. A. Fischbach, and G. M. Whitesides.
Design for mixing using bubbles in branched microflu
idic channels. APPLIED PHYSICS LETTERS, 86(24):
8688,2005.
J. L. Guermond and L. Quartapelle. A projection fem for
variable density incompressible flows. JOURNAL OF
COMPUTATIONAL PHYSICS, 165(1):167188, Nov
2000. ISSN 00219991.
D. Jacqmin. Calculation of twophase navierstokes
flows using phasefield modeling. JOURNAL OF COM
PUTATIONAL PHYSICS, 155(1):96127, Oct 1999.
ISSN 00219991.
E. Katifori, G. J. Szollosi, and M. O. Magnasco. Damage
and fluctuations induce loops in optimal transport net
works. PHYSICAL REVIEW LETTERS, 104(4):048704,
NonSplitting
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Jan 2010. ISSN 00319007. doi: 10.1103/Phys
RevLett. 104.048704.
D. R. Link, S. L. Anna, D. A. Weitz, and H. A. Stone.
Geometrically mediated breakup of drops in microflu
idic devices. PHYSICAL REVIEW LETTERS, 92(5):
05450310545034, 2004.
L. MenetrierDeremble and P. Tabeling. Droplet breakup
in microfluidic junctions of arbitrary angles. PHYSICAL
REVIEW E, 74(3), 2006.
G. I. Taylor. The Formation of Emulsions in Definable
Fields of Flow. Royal Society of London Proceedings
Series A, 146:501523, October 1934.
J.D. van der Waals. The thermodynamic theory of cap
illarity under the i\ p'1,ihc ,i, of a continuous variation of
density. J.D. Verhandel. Konink. Akad. Weten. 1 (J. Stat.
Phys., 20, 1979 (in English)), 1893.
W. Villanueva and G. Amberg. Some generic capillary
driven flows. INTERNATIONAL JOURNAL OF MUL
TIPHASE FLOW, 32(9):10721086, Sep 2006. ISSN
03019322.
