7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Interfacial Magnetophoresis in Magnetic Fluids
P. Yecko, A.D. Trubatch and W.K. Leet
Department of Mathematical Sciences, Montclair State University, Montclair, NJ 07043, USA
t Advanced Photon Source, Argonne National Laboratory Argonne, IL 60439, USA
philip.yecko@montclair.edu and david.trubatch@montclair.edu
Keywords: Magnetic fluid, bubbles, magnetophoresis, interfacial flow, ferrofluid
Abstract
We report on two examples of magnetophoretically driven microfluidic flows. First, we examine the aggregation of
bubbles into chains within a magnetic fluid, or ferrofluid. Bubbles are examined both directly, using Xray phase
contrast imaging, and computationally, using a volume of fluid (VOF) simulation code. High resolution Xray images
in the bulk of waterbased ferrofluid (EMG607/707) reveal that gas bubbles with diameters of a few 100 pm readily
attract one another along the field direction, forming linear chains of two or more bubbles. For the simulations,
magnetic stresses are captured using an interfacial force that is derived from a multigrid solution of Maxwell's
equations. A novel multiple color function scheme for the two phase flow allows nearby bubbles to come into contact
without coalescing, as would normally occur in a VOF model. Measurement of trajectories in experiments shows that
attraction is driven by the magnetophoretic force resulting from the induced fields of the bubbles. Bubble trajectories
in both simulations and experiments also shows evidence of unsteady drag effects. To address magnetic interfacial
forces more generally, we examine a second example: instability driven by magnetic interfacial stress in layered
channel flow in a microchannel. Linear theory makes precise predictions of modal growth rates for uniform fields
applied parallel to and normal to the channel walls, indicating that modest fields may be used efficient to control
processes such as pumping and droplet generation in channels a few 100 pm wide.
Introduction
Fluids with suspended nanoscale magnetic particles, re
ferred to as magnetic fluids or ferrofluids, posses a range
of scientifically interesting and technologically desirable
properties, most notably that their properties and flows
can be modified by applied magnetic fields (Holm and
Weis 2005; Rosensweig 1985). Well known applica
tions of ferrofluids include liquid seals in devices such
as harddisk drives and crystal growing furnaces, and as
coolants in loudspeakers (Rinaldi et al. 2005). Recently,
ferrofluids have begun to be utilized in microscale and
nanoscale self assembly (Ganguly et al. 2005a; Yellen
et al. 2005), microfluidic pumping (Hartshorne et al.
2004), eye surgery (Afkhami et al. 2008), drug de
livery (Gangulya et al. 2005b; Alexiou et al. 2001),
chemotherapy (Alexiou et al. 2005), magnetic cell sepa
ration (Zborowski and Chalmers 2005) and as magnetic
MRI contrast agents (Pankhurst et al. 2003). Ferroflu
ids are thus increasingly utilized at small scales, and in
configurations involving an interface with a gas or an
other liquid; several other potential applications are re
viewed in (Rinaldi et al. 2005). Relative easy manipula
bility makes magnetic fluid a very promising alternative
in the advancement of microfluidic technologies, such
as precisely tuned pumping or droplet formation(Ozen
et al. 2006; Lin et al. 2004), where magnetic fluid can
be used to drive a juxtaposed nonmagnetic fluid. Mag
netic control is particularly attractive in light of the chal
lenges in controlling fluids using electric fields (which
may require high voltages), thermal forcing (which may
require high temperatures) or complex channel geome
tries (which may be technically difficult).
A robust and flexible numerical model capable of ac
curately simulating interracial ferrofluid flows would be
an efficient and valuable alternative to experimentation
in the discovery and design of new applications. But
so far, direct numerical simulation (DNS) of interfacial
ferrofluid flows has been seldom performed. In most
cases, only static fluid configurations have been sim
ulated, in which the dynamic Navier Stokes equations
are not solved and which employ numerical methods
that cannot be easily extended to dynamical problems
(Lavrova et al. 2006). The DNS code that we apply
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and validate here is flexible and robust; in the absence
of magnetic effects it has proved to accurately model a
broad range of interfacial flows, including rising bubbles
(Gueyffier et al. 1999) and primary atomization (Boeck
et al. 2007). Bubbles are an important instance of an in
terfacial flow which are ubiquitous in many industrial
and natural examples, often exerting a ,igililk.ii im
pact on the transport properties of the flow. Broadly
speaking, a "bubble" may also refer to a droplet of a
nonmagnetic fluid immersed in ferrofluid.
Direct observations within ferrofluids are severely
quashed by the opacity of the fluid to visible light. To
circumvent this difficulty, most of the previous experi
mental work using optical techniques has been limited
to thin films (tens of pm) or highly diluted samples.
Specifically, thin film studies have been performed to in
vestigate the formation of clusters (Hayes 1975; Horng
et al. 2001; Jones 1985; Taketomi et al. 1991) and the
(2dimensional) shape of bubbles (Bashtovoi et al. 2005;
Drenckhan et al. 2003). However, such thin sample en
vironments are subject to wall effects (Wang et al. 2003;
Wiedenmann and Heinemann 2005) and may not be rep
resentative of a 3dimensional system. In particular, the
bubbles in these studies had diameters larger than the
thin sample constraints and are thus not "free" but at
tached to the walls of the sample cells.
Xrays have been used previously to image bubbles in
bulk ferrofluids, but, in those studies (Jeyadevan et al.
1999; Nakatsuka et al. 1999), the resolution has been
insufficient to quantify interface dynamics and shapes.
Indeed, to the best of our knowledge, there has not been
any previously reported highresolution measurement of
bubble shapes or dynamics in a bulk ferrofluid sample.
More generally, there is a lack of experimental data re
garding multiphase ferrofluid systems. This lack of ex
perimental data is ,igmi.lk.ll because most of the new
applications being discussed for ferrofluids, such as bio
medicine, involve multiphase configurations.
Here, we describe observations obtained by syn
chrotron xray phasecontrast imaging. By this method,
we measured the shape and dynamics of small (< 1
mm diameter) bubbles in a tube (9.5 mm outer diam
eter and 6.5 mm inner diameter) filled with ferrofluid.
We used a waterbased suspension of magnetite parti
cles with mean diameter of 10nm (FerroTec EMG607
or EMG707) and bubbles of both air and water vapor.
The phaseeffects dramatically enhance the image con
trast, which enables one to track the bubbles and their
shapes. Bubbles with diameters as small as 10 pm can
be seen and easily quantified. In particular, the forma
tion of bubble chains aligned parallel to the applied uni
form field can be seen (Fig. 1). This is the first time
that the bubble dynamics within a bulk ferrofluid sample
have been tracked with pmlevel accuracy.
Figure 1: A four bubble chain (boxed) in EMG607,
bubble diameters ~ 0.1 mm.
Acceleration between two bubbles lined up along
the applied field will result from the attractive magne
tophoretic force arising from the the nonuniform exter
nal fields induced by each bubble. To compute this force,
a solution for the field (or potential) is needed.
Our observations have shown bubble Reynolds num
bers up to 0(1) are common during aggregation events.
In this regime, it is unreasonable to rely entirely on an
unsteady Stokes flow simulation, such as a boundary
integral method. Rather, we are required to solve the
Navier Stokes equations for interfacial flow of magnetic
fluid. The numerical modeling of incompressible flu
ids with interfaces is a broad problem with many ap
proaches (see Scardovelli and Zaleski 1999; Tryggva
son et al. to be published 2010, for a review). One of the
most successful and robust approaches is to treat a multi
component fluid as a single fluid with abrupt changes in
density and viscosity and use an auxiliary function to
capture the interface. The two most common interface
capturing methods are the Level Set method, in which
the interface is the zero contour of a smooth function
such as the signed distance from the interface, and the
Volume of Fluid (VOF) method, in which a color (or
phase) function gives the fraction of each computational
cell that is occupied by the reference phase. Levelset
methods naturally allow precise computation of the ge
ometric quantities, such as the normal and curvature of
the interface, yet tend to conserve volume (mass) poorly.
VOF methods, on the other hand, conserve volume more
precisely, but generally involve estimates of the interface
curvature. The inaccurate curvatures estimates of VOF
methods are one of the sources of the socalled spuri
ous currents (Scardovelli and Zaleski 1999) observed in
VOF simulations.
All singlefluid methods, including VOF, Level Set
but also others, such as the CIP and phase field meth
ods, are prone to inaccuracies, one source of which is
the need to solve the pressure Poisson equation with a
coefficient which is discontinuous at the interface, where
the density jumps in value (Tryggvason et al. to be pub
lished 2010). This problem worsens for large density
jumps, such as for air and water. While the problem
may be mitigated by using a more robust elliptic solver,
this approach is computationally very expensive. An al
ternative is to abandon the single fluid approach in favor
of socalled sharp interface methods, such as the Im
mersed Interface or Ghost Fluid methods, but at the cost
of a considerably more complex numerical scheme. For
problems that are also physically complex, such as mag
netic fluids, increased code complexity is highly unde
sirable. Instead, we adopt a common approach, and use
an artificially large gas density, preserving both speed
and accuracy with little sacrifice, as the density ratio
plays a minor role in many problems involving bubbles
(Tryggvason et al. to be published 2010).
In a previous work (Korlie et al. 2008), we described
a volumeoffluid (VOF) numerical method for the sim
ulation of dynamic twophase flows of magnetic flu
ids, including rising bubbles and falling droplets. Aside
from the jump in density and viscosity generally found
at fluid interfaces, in a magnetic fluid there will usually
be a jump in magnetic properties. The influence of the
magnetic field on the fluid dynamics originates from the
Maxwell stress and can be reduced to a normal stress
acting at the interface. This stress depends on the mag
netic field at the interface and on the relative magnitudes
of the magnetic permeabilities of the two fluids.
In the current work, we have expanded our previous
VOF method to include multiple color functions as in
dicators of separate bubbles. By computing the inter
facial stresses from distinct color functions, we avoid
the merger of bubbles that occurs automatically in VOF
methods (when interfaces of distinct bubbles approach
one another). We are therefore able to simulate aggrega
tion of bubbles. Although this method is not directly mo
tivated by molecular and physical properties of the fluid,
it nevertheless plays the role of an effective disjoining
pressure and has been used, previously, in a combined
VOF/LS framework to effectively study the impact of a
droplet onto a flat interface (Coyajee et al. 2006).
We simulate and measure the dynamics of aggrega
tion using two dimensional computational models. Sim
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
ulation results are compared to theoretical predictions
stemming from a simple balance of drag and dipole in
duced magnetophoretic forces. The theory does not ac
count for the nonlinearity of the magnetic material, non
zero Reynolds number effects, nor the deformability of
the bubbles, all of which are likely to be important in
both simulations and experiments. Nevertheless, the
ory does provide the common benchmark of a simple
scaling relation prediction for both the 2D and 3D cases
that should be valid over some part of the aggregation
path. In numerical simulations we have generally used
the same physical parameters as in experiment, but there
are differences due to the computational constraints, as
described above, and because some quantities (e.g. the
density inside a bubble) are uncertain and not easily
measured. Our goal is to explain the observed bubble
behavior in terms of the hydrodynamic and Maxwell
stresses and to demonstrate that VOFbased numerical
methods are able to efficiently and accurately model ag
gregation. Although the experiments are, naturally, three
dimensional, the value of twodimensional simulations
is that they are both efficient and they can provide valu
able insight into the dynamics where analytic solutions
are intractable and experimental measurements are diffi
cult.
In the second part of this paper, we study the stabil
ity of a sheared ferrofluid interface in order to provide a
more complete understanding of multiphase ferrohydro
dynamic flow. Low Reynolds number flow in a straight
channel is examined due to its relevance to microflu
idic applications and to our experiments. Our desire is
to know which fields are most effective at destabiliza
tion, as evidenced by the growth rate of linear instability,
and at what scales this instability will develop, as evi
denced by the wavenumber of fastest modal growth. We
also plan to use the stability predictions of this work to
validate our numerical simulation code, currently being
used to model bubbles in magnetic fluids (Korlie et al.
2008). The magnetic flow stability problem was stud
ied more completely in the recent paper (Yecko 2010),
henceforth Y10. In that work stable nonmagnetic base
flows are considered, the Reynolds and Weber numbers
are fixed (Re We 1 ), and only exactly parallel or
exactly normal field orientations are considered, while
a fully nonlinear magnetic material model based on the
Langevin description is employed.
1 Bubble Magnetophoresis
Mathematical Background
We consider nonmagnetic bubbles in a magnetic fluid
with permeability p /1 po and susceptibility X
p//po 1, where the material is linear (i.e.x indepen
dent of H) and the magnetization, M = MoxH is par
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
T~I
YL 2 (a)
y
n(:
m, n, r, Pr
(2)
x t2 P2 U (y))
x,t)
y=0
y=L1
f2)
SX(H)
(i)
Hp No
Figure 2: Sketch of the base steady flow and field con
figuration.
allel to the H field. Each fluid phase has constant den
sity, p, and dynamical viscosity, r. The gives a modified
NavierStokes equation
S( +u. Vu)
Vp+pg + lVu1H H Vt
(1)
where g is gravity and we assume incompressible media,
such that
V u 0. (2)
In (1), the pressure p represents only the hydrodynamic
pressure.
At the interface, some quantities are required to be
continuous, for example the velocity field, while others,
such as the pressure, have jumps as a result of discon
tinuities or singularities in the mathematical representa
tion of physical quantities, such as viscosity. Specific
jump conditions at the interface result from a force bal
ance of discontinuous and singular terms.
To model the system for simulation by a VOF method,
we treat the fluid and bubbles as a single fluid with spa
tially varying density, viscosity and magnetic perme
ability. The phases are distinguished by the value of
a characteristic function (x, t), which is equal to one
in the reference or primary phase, and zero in the sec
ondary phase. If there is no phase change, the charac
teristic function is passively advected by the flow and
its time evolution is described by a standard advection
equation, d/dt = 0. The spatial distribution of the
density p in the single fluid approach is then defined as
P = Pi + (1 7) p2. A similar expression holds for
the viscosity and the susceptibility.
Magnetic forces are captured in the rightmost term of
(1). While this term vanishes in the bulk of each fluid,
where p is constant, on the interface, where p jumps
in value, V/ = Sn, where 6 is a deltafunction (the
"derivative" of a step function) and n is the normal to
the interface (Melcher 1963; Rosensweig 1985).
The effects of the magnetic field reduce to a singular
normal force acting on the interfaces, just as for the force
of interfacial tension. We are then able to describe a two
fluid system with a dynamic interface using the single
equation
p( + u Vu) Vp+ pg +V uD (3)
+acOTs n + FM6s n,
where: Dij (duj/dxi + dui/dxj) is the rateof
strain tensor; 8s is a distribution identifying the inter
face; a is the interfacial tension coefficient; K is the in
terface curvature; and TM is the magnitude of the mag
netic interfacial force. In terms of the jump in magnetic
normal stress at the interface,
M = n n (/HHH H2)], (4)
where the square brackets denote the jump across the
interface and H = HI.
One can use the continuity of normal B and tangential
H to find a more convenient expression for the magnetic
interfacial force,
M = ',) [H +/HJ (5)
where the subscripts i and o refer to inner and outer,
while n and t denote normal and tangential components.
This expression clearly shows that the difference / 
po determines the sign of the magnetic interfacial force,
while the ratio pI/p~, determines the relative importance
of the normal field component.
In the absence of free currents, the magneto
quasistatic Maxwell's equations are
VxH 0 V B 0,
where B /o(H + M), are more conveniently posed
in terms of a magnetic potential, p:
V. (1 + x)V 0, (6)
where H Vp. On boundaries and interfaces it is re
quired that [n B] 0 and [n x H] 0 Equation (6)
is also solved in a single domain where X varies discon
tinuously from the value Xi to Xo at the interface. In this
way, the conditions that normal B and tangential H re
main continuous are naturally captured by the solution.
Magnetophoretic Attraction between Bubbles
Because there are no tangential interfacial stresses in
a magnetic fluid, three types of interactions are expected
to occur in a pair configuration: (i) magnetophoretic ef
fects, due primarily to the nonuniform field induced by
the companion bubble; (ii) drag, including that depend
ing on the flow induced by the motion; and (iii) hydro
dynamic effects due to film drainage at small separa
tions. Drainage of the fluid separating two approaching
interfaces is a lowRe phenomenon that controls the fi
nal contact. Drainage is not of primary concern here; its
effects will be noted but not examined in detail.
Computing the magnetophoretic force between two
objects is a complex problem which has, as yet, no
simple closed form result even in the simplest case 
when the objects are identical, spherical and aligned
either along or normal to the imposed uniform field.
In discussing this problem, we refer extensively to the
problem of the dielectrophoretic force between dielec
tric objects, which has received considerably more at
tention. For perfect dielectrics and linear magnetic ma
terials these two problems are exact analogs upon iden
tifying the fields E, D with H, B and the dielectric per
mitivity E with the magnetic permeability p. In our ex
periments, the largest bubble exposed to maximum field
exhibited a deformation D ~ 0.04. Bubbles undergo
ing aggregation in our experiments have D ~ 10 2 and
remain nearly spherical. We therefore choose to restrict
our attention to spherical objects from here on.
A spherical bubble immersed in magnetic fluid and
placed in a uniform field will modify the field, adding a
dipole contribution. This solution can be found in elec
tromagnetism texts (i.e. Stratton 2007). In the case of a
sphere of radius R centered at the origin with uniform
imposed field H, Hay, the potential 40 outside the
sphere is given by
co(r, 0) = H rcos 0 Ha cos ,
3 + 2X r2
from which we can compute the radial field component
H,r dO/Or, or:
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
order contributions. Thus, the expression (7) is not ex
act for the force between two bubbles in an external
uniform field, since each individual bubble imposes a
nonuniform field on its companion. When the field non
uniformity is not too large, however, (7) will give ap
proximately correct results. Computing the net force in
the case of the bubble pair requires a summation over
all the active multiple pairs for each bubble, a signif
icant computational challenge. In solving the Maxwell
equations and computing the normal stresses in our VOF
simulations, we are implicitly performing exactly this
calculation in the twodimensional case. Later, in Fig
ure 6, we will show the magnetic potential 0 as com
puted by the multigrid algorithm within our simulation
code. The field is found directly from this potential and
magnetophoretic force is realized as the normal stress of
(5).
To estimate the dynamics of the observed bubble ag
gregation, we assume that one bubble's nonuniform
field acts as a driving field Hd on the companion bubble.
We further assume that only viscous stresses oppose the
motion, to obtain the balance
67rlRU = FMAP,
where Stokes drag appears on the left hand side and U is
the velocity of the bubble. Hence, in threedimensions,
the bubbles should approach with velocity
R2
U3d OC 4
,4
while, in two dimensions,
H, = H cos 0 (1
For an ellipsoid, the X factors are modified but the 0 and
r dependence is the same. In two dimensions (for an
infinite cylinder) the equivalent expressions are:
o02(r, 0) = H r cos 0 Ha X cCs 0,
2+X r
and
H,2 Hcs 2+ r2 x )
A nonuniform field Hd drives a magnetophoretic
force on a spherical dipole of radius R given by
FMAP = 27po(1 + x)R3 VH2, (7)
(for a derivation, see e.g. Jones 1995). A uniform
field induces a dipole around a spherical bubble; a non
uniform field will also induce quadrupolar and all higher
Computational method
Equations (2)(4) are solved numerically by including
the magnetic interfacial force (5) and a multigrid relax
ation algorithm to solve (6) within an existing VOF code
(Korlie et al. 2008) derived from SURFER (Lafaurie et al.
1994). In VOF methods, the interface is captured by us
ing a color function C, with values in [0, 1], that is the
discrete version of the characteristic function, 7. In par
ticular, C'i 1 when the computational cell (i,j) is
fully occupied by the reference phase, C" 0 when
only the secondary phase is present and 0 < Cj < 1
when the interface cuts the cell where the value of the
color function represents the fraction of the cell occu
pied by the reference phase. The C function is used
to track the identity of the fluid, to reconstruct the in
terface (here with a piecewiselinear interface calcula
tion, VOFPLIC) and to compute its normal and curva
ture at any computational point. In terms of the color
2X R3
3 + 2 r3 )
R
U2d o r3
function C, the material properties in the single fluid ap
proach are defined as follows: p = p1C + (1 C)p2,
ri = lC + (1 C)12 and x = XC + (1 C)X2.
In our experimental observations, when smallscale
bubbles come into contact, they do not coalesce, most
likely a direct result of the large capillary forces found
at such small radii. However, in VOF models, coales
cence occurs automatically whenever two interfaces first
occupy the same computational grid box. We have cho
sen to address coalescence suppression as a numerical
issue, rather than as a physical one. The use of a single
color function does not allow more than one interface to
be represented inside a single computational cell. Be
cause of the local linear interface reconstruction, for a
good approximation its radius of curvature should be at
least four times the grid spacing. Furthermore, when
different interface lines approach each other they affect
the reconstruction when they are within the same sten
cil, which is used to compute the normal vector, at a
distance of two to three grid spacings. In a standard
twodimensional VOF simulation, the approaching in
terface lines will be locally highly distorted and auto
matically merge in a few time steps. We avoid this nu
merical merging between separate bubbles by introduc
ing multiple color functions, C, one for each distinct
bubble. The advection of the different color functions,
the reconstruction of the associated interfaces and the
computation of the surface tension forces are performed
separately for each bubble. For all other aspects of the
computation, including the multigrid solution of the in
compressibility condition, the separate color functions
are treated as a single function. Interface coalescence
is thus largely suppressed without introducing an artifi
cial suppression mechanism. In exceptional cases, merg
ing/overlapping still occurs, mainly due to the linear in
terface reconstruction, but this has happened rarely and
does not influence the results in an appreciable manner.
To our knowledge, multiple color functions have been
used previously only with a combined VOF/LS tech
nique (Coyajee et al. 2006).
The solution of (6) is made difficult by the fact that the
coefficient, 1 + X, experiences a jump at the interface.
We use a multigrid algorithm, as described in (Korlie
et al. 2008). The magnetoquasistatic Maxwell system is
solved using Neumann boundary conditions on the po
tential at the top and bottom of the domain correspond
ing to a uniform vertical field while, on the horizontal
walls, periodic boundary conditions are applied.
Experimental procedure
The experiments were conducted at the XOR 32ID
beamline at the Advanced Photon Source in Argonne
National Laboratory. The ferrofluid setup consisted of
a 9.6 mm outer diameter (6.5 mm inner diameter) plexi
glass tube held fixed between the poles of an electromag
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
net. The tube axis was vertical and the applied field was
horizontal. Probe measurements confirmed that the field
in the middle of the poles, where the measurements were
done, was fairly uniform to about 1 mT/mm. The field
can be varied by changing the applied current to the elec
tromagnet; the maximum field (at 6 amps) was measured
to be 0.175 T (H = 1.4 x 105 A/m). 25 keV xrays were
used, with a sampledetector distance of about 0.5 m for
good phasecontrast. A 100 pm thick yttrium aluminum
garnet scintillator was used to convert the transmitted x
rays into visible light that was then imaged onto a video
rate (30 frames per second) CCD camera (Cohu, USA)
using a 5X microscope objective with the accompany
ing tube lens. The CCD camera exposure rate was set to
1 ms so as to minimize bubble motion related blurring.
The field of view was about 1.3 mm horizontal x 0.9
mm vertical and the spatial resolution of the setup was
estimated to be 23 pm. The video was recorded onto
miniDV tapes and subsequently digitized using commer
cial software (Apple iMovie).
Two waterbased ferrofluids were used, EMG607 and
EMG707, having cationic and anionic surfactant coat
ings respectively; both were obtained from Ferrotec Inc.
The following physical characteristics were provided by
Ferrotec: density p 1.1 g/cc, viscosity v 5 m Pa s,
initial susceptibility Xo 1.7 (for EMG707, Xo 1.5),
saturation magnetization Mo 11 mT and average par
ticle diameter d 10 nm. These values were also ap
plied in the numerical models. We also examined a di
lute solution (1:4 dilution of EFHl:mineral oil), also
from Ferrotec.
Initially, we attempted to generate small bubbles with
a syringe, but in order to prevent the syringe needle from
introducing nonuniformities in the magnetic field, the
needle had to be placed far from the region of interest.
However, this configuration created the problem of di
recting bubbles created far away into the uniform field
region because the field gradients around the uniform
field region tend to expel the nonmagnetic bubble. The
result was that we were unsuccessful in introducing bub
bles created from a syringe needle outside the uniform
field region into the uniformfield region.
Instead, we relied on the high intensity of the xray
beam to create bubbles via ionization. In the presence of
an imposed field, the agglomerations of magnetic parti
cles serve as nucleation sites for the bubbles. The advan
tage of this process is that the bubbles are created within
the uniform field and exactly in the location where we
are imaging. The disadvantage is that there is little con
trol over the number of bubbles or their locations within
the xray beam. In Fig. 1, the magnetic particle ag
glomerations are clearly visible as dark bands aligned
parallel to the imposed magnetic field. Such macro
chains are well known and were seen in early experi
ments by Hayes (1975) and Kreuger (1980) and were
first explained theoretically even earlier by de Gennes
and Pincus (1970).
2 Sheared Microchannel Interfaces
To more directly examine the effect of magnetic forces
on the dynamics of an interface, we examine a very gen
eral flow stability problem. This flow consists of two
viscous layers in a channel, separated by an initially flat
horizontal interface and acted on by a horizontal pres
sure gradient and either a normal or a parallel magnetic
field, as depicted schematically in Fig. 2. The fluids have
constant densities, Pi and P2 and viscosities, p1 and p2
and occupy depths L1 and L2, and interfacial tension a
acts on the initially interface, whose location is given by
y =r(x, t).
The magnetization M of the magnetic fluid is given
by the Langevin constitutive relation:
M Ms (coth
1) H
E/ H,
where 3 = 4 MH and Ms Mdyv is the sat
uration magnetization; Md is the domain magnetization
of particles of radius a and Yv is their bulk volume frac
tion in the ferrofluid. For strong fields we find M 
MsH/H. For weak fields, the relation M = IH is a
good approximation, which we have used above in ex
amining the bubbles. The quantity X, is the initial sus
ceptibility and is given by
Ms, 47ra3 poM3n v
XI = (11)
3H 9 kT
XI is clearly constant for isothermal flows, as we assume
here. Without loss of generality, in place of (10) we may
write M = X(H)H, where
(H) = (cothi~ ) (12)
The magnetoquasistatic Maxwell equations now ap
pear in the form
V. ([1+ x()]V) = 0. (13)
Steady solutions
A uniform magnetic field H (H,, H,) is applied to
the flow; only purely parallel (P) fields Hp (Ha, 0),
and purely normal (N) fields HN = (0, Ha) are consid
ered, where Ha is a constant. The steady unperturbed
flow U ) (U(J)(y), 0) is quadratic in y and is iden
tical to that of layered channel flow of two ordinary vis
cous fluids; the functional form is readily found in the
literature (Hooper 1989; Yecko 2010).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Perturbation equations
Perturbing the steady field and flow solution results
in a set of equations for the perturbation velocity, pres
sure, and magnetic potential, which are neither de
rived not reproduced here, but can be found in Yecko
(2010). Each perturbation quantity, y, has the usual
form y(x,y,t) = (y' Wt), leading to an Orr
Sommerfeld type equation in each layer, augmented by a
magnetic potential perturbation, interface matching con
ditions on y = and boundary conditions on the walls.
The total stress tensor is the sum of fluid and Maxwell
contributions: ij = Eij + Mij, where Eij is the ordi
nary fluid stress (see Y10 for details). A key property
determining the stability is the jump in normal stress
component at the interface, which is balanced by surface
tension: [7 .. ](2) [7 ,... \ 1) a(2l/ zx2 +
dl/9z2). We adopt the same form of the magnetic stress
tensor as in Y10:
Mij = (1 + x)H;Hj 6ij (Po MdH + t H2
(14)
which leads to a net force only on an interface where
magnetic properties change (Rosensweig 1985; Oden
bach 2009; Yecko 2010).
Nondimensionalizing leads to the following control
parameters: density ratio r := p2/pi, the Weber num
ber We = PUL1 and a magnetic fluid parameter
Ma P ug2, and the (layerj) Reynolds number based
on the interface velocity: Re(j) = pj UoL
The linear stability problem is then solved for eigen
values, c, and eigenvectors, which are computed numer
ically using the code developed in (Yecko 2008, 2010);
instability occurs for cw > 0, where cw := Im (w).
3 Results: Bubble Aggregation
In a sequence of experiments, water vapor bubbles cre
ated by Xray induced bubble nucleation were observed
and their subsequent dynamics recorded. Most bubbles
that we tracked had diameters in the range of 50 to 200
pm and readily aggregated with neighboring bubbles.
In this work, we use the term "aggregation" to refer to
the motion of two initially separated bubbles coming to
gether to form a single bubble pair, as depicted in Fig.3
and Fig.4. We never observed a bubble pair merge to
form a single larger bubble. Aggregative motion oc
curred approximately along the field lines, as the magne
tophoretic force is expected to be attractive in this direc
tion, but deviations were present due to the various ini
tial locations of the two bubbles. Repulsion of two bub
bles whose line of centers is perpendicular to the field
was also observed, but these events are easily missed
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 3: Aggregated water vapor bubbles (boxed) in EMG607 ferrofluid, frames extracted from 30 fps digitized
video; each frame has actual height 533 pm and width 316 pm; frames are 0.2 sec apart. The large overlapping bubble
in the lower right comer of the frames is in another plane.
due to the lack of a clear, easily distinguished end state.
The presence of magnetic nanoparticle chains appears to
strongly suppress motion laterally to the chains, but the
aggregation proceeds lengthwise along the chains.
Figure 3 shows an example of aggregation between
two bubbles. The Reynolds number in the figure se
quence is Re ~ 0(10 2), although the value clearly
increases in the final moments before contact. As is typ
ical in most of our experiments, additional bubbles are
present in the frame. This is because we have no di
rect control over the number and location of the bubbles
generated by the xray beam. Many of the bubbles that
appear in the image are located at different depths and
even though they may appear, in projection, to intersect,
in reality they are not in contact. This is due in part to the
projection of bubbles at different depths, normal to the
image plane. In spite of the crowded environment, we
were able to identify several aggregation events clearly.
Once in contact, as in the last frame of Figure 3, the
pair persisted with its axis parallel to the field and of
ten accumulated other bubbles, forming a longer chain.
We observed chains of up to seven bubbles. Note that
in the experimental images displayed here, the field is
oriented vertically, with a slight skew due to equipment
misalignment. In the laboratory, the field was applied
horizontally, but we have rotated the display to match
numerical simulations. There is no risk of confusion
since buoyant motion is negligible, as we know from the
small (gravitational) Bond number.
Next we present direct numerical simulation results of
bubble aggregation and then we analyze the dynamics
as a consequence of magnetophoretic force. In particu
lar, we examine pairs of ferrofluid bubbles aligned such
that the separation of their centers is parallel to the ap
plied uniform field H,. Without such initial alignment,
torques will also act on the bubbles and the resulting dy
namics becomes much more complex. Simulations were
performed in 2D, but the results can be easily general
ized to 3D geometry, as we do below.
The following conditions are fixed: fluid density
pp 103kg m 3, bubble density p, 102kg m 3,
fluid viscosity Ip 10 3Pa s, bubble viscosity p/
10 4 Pa s, interfacial tension a 0.07 N m 1. The
uniform field IH, I 20 kA m 1 and the bubble sep
aration are parallel and along the yaxis. The bubble
radii are both set equal to rp 100/m and the initial
separation of their centers is Lo = 2rp. Note that with
the above value of the magnetic field the magnetic bond
number BoM z 1 when X 3. In the following results,
BOM will be varied by changing X.
Figure 4 shows a typical aggregation of a pair of bub
bles in the case X = 3. The dynamics of the physical
process can be divided into four consecutive steps: slight
prolate deformation, magnetophoretic acceleration, hy
drodynamically damped magnetophoretic motion, and
a final deceleration mediated by "drainage" of the fluid
film separating the adjacent bubbles. A zoom of the con
tact region is shown in Figure 5, where the effectiveness
of the multiple color function approach can be inspected.
The potential solution from which the magnetophoretic
force originates is presented in Fig.6; this p field was
found using multigrid relaxation on the magnetoquasi
static Maxwell equations and boundary conditions.
As can be estimated from the figure, Re ~ 0(0.1) for
the attractive motion of this case. At this Reynolds num
ber the deviations from Stokes flow may not be negligi
ble. Bubble motion under acceleration may also be mod
ified by Bassett (or history) forces arising from boundary
layer development on the interface and by added mass
effects due to the displacement of ambient fluid (Crowe
et al. 1997). Stresses may also arise due to a nonuniform
ambient velocity field (Faxen forces) and influence the
dynamics of bubble aggregation.
In a laboratory experiment where many bubbles are
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
t=8
t=12
t=16
t=18
Figure 4: Simulation of two bubble aggregation in ferrofluid: color function field time sequence; time units are
10 6 s and length units are 10 5 m. Initial separation is 1 bubble diameter.
t=0 2
t=18
Figure 5: Detail of the contact region between two
bubbles; (left): as it occurs at t 18 in the simulation of
figure 3; and (right): as it occurs in experiments; frame
width is 285 pm.
present, field nonuniformities due to the surrounding
bubbly medium are prevalent. Nevertheless, we have
seen in both simulations and experiments that aggrega
tion dynamics is in good agreement with the above scal
ings. In Figure 7(a) we compare a sequence of aggrega
tion simulations with the prediction (9) where X is varied
to produce three distinct cases. Except for the onset of
motion (rightmost points) and the final approach (left
most points) there is good agreement with the expected
r 3 dependence. As we have shown previously (Kor
lie et al. 2008), bubbles under acceleration show devia
tions in their trajectories that can be explained by added
mass effects. In comparing our numerical results to the
simple scaling relation (9), we stress that there are sys
tematic deviations from an exact r 3 dependence which
are expected in a full NavierStokes based simulation in
05 0 05
05 0 05
Figure 6: Magnetic potential function, 0, correspond
ing to droplet experiment of figure 3, as found by multi
grid relaxation; in the presence of an imposed uniform
field and a neighboring bubble, the field of each bubble
is not an exact dipole, but approximately so; the conse
quent gradients drive the magnetophoretic attraction.
t=0.2
t=4
2
1
0
1
2
1
1
0 1
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
which hydrodynamic effects, including added mass, are
naturally captured.
In Fig.7(b) we show the equivalent behavior as ob
served in experiments. The data in the figure include
the aggregation depicted in Fig. 3 and two similar events
at different field strength and for different size bubbles.
After onset and before final approach, a clear r 4 de
pendence is apparent, consistent with the prediction of
(8) and with the analogous twodimensional result of
the simulations. Although our imaging technique has
allowed us to visualize bubbles in ferrofluid in unprece
dented detail, it remains an experimental challenge to
further resolve the temporal dynamics for an isolated
bubble pair. Until these challenges are surpassed, DNS
remains a viable tool to examine the dynamics of mag
netophoresis in magnetic fluids.
4 Results: Interfacial Instability
Of greatest interest to the control of flows in small de
vices are the parameters required to ensure stability or
instability of the interface, including the applied mag
netic field (Ha), its orientation (P or N), the relative layer
thickness (n), and the relative magnetic permeabilities
(/r). Because ferrofluids exhibit a wide range of vis
cosities, which may also vary with the applied field, the
effect of viscosity ratio (m) is also examined. Several
other control parameters are held fixed in the stability
calculations, namely: Re 1 We 1 Fr 0
and r = 0.9. Control parameters are also given refer
ence values at which each is fixed when not being var
ied; these are: n 2, m 1.5, and p, 1/3. Note
that the permeability ratio pr is a control parameter and
refers only to the relative initial permeability values. For
a lower ferrofluid layer, the initial permeability ratio is
related to the initial susceptibility as p 1/1 + Xi,
while for an upper ferrofluid layer, p = 1 +X. The ref
erence value = 1/3 is therefore equivalent to XI = 2
in a bottom layer of ferrofluid.
The base nonmagnetic flow is that of two superposed
fluids in a channel, first studied by Yih (1967) who found
that a longwave mode, dependent on the viscosity strati
fication, can be unstable even as Re 0. This socalled
interfacial (or Yih) mode was was found in later stability
studies (Renardy 1987; Hooper 1989; Yiantsios and Hig
gins 1988) to exhibit a stable regime, sometimes called
the "thin layer effect" since in one manifestation long
wave stability is ensured by choosing m and n such that
the less viscous layer is thin enough. In the cases exam
ined here, the base nonmagnetic flow is always main
tained within this longwave stability regime, while sur
face tension stabilizes the shorter waves. Thus, any in
stability must be due to magnetic effects.
Figure 8a depicts the wavenumber of the fastest grow
Figure 7: Velocity of bubble pair aggregation, U (in
pm/sec) as a function of centroid separation, r (in pm),
as found in simulations (a) and experiments (b), includ
ing best fit lines. The fits were performed only using the
central data points (filled symbols): startup and con
tact transients are displayed using open symbols. In
(a) the simulations are 2D and an r 3 law is expected;
the plotted best fit lines have slopes (top to bottom) of
2.8, 2.9, 3.1 and 3.2 with corresponding R2 fit
values (again top to bottom) of 0.953, 0.956, 0.954 and
0.973. In (b) an r 4 dependence is expected; the fit
values are 4, 4.6 and 4.2 for cases A,B,C, respec
tively with corresponding R2 fit values of 0.915, 0.980
and 0.872.
X2
A X=3
+ X=4
X=5
035 04 045 05
log r
055 06
ing instability as a function of Ha at fixed layer thickness
ratio n = 2, comparing a purely parallel (P, top data)
field to a purely normally (N, bottom data) oriented one.
Note that for the P field, the finite wavelength mode is
stabilized while in the N case the growth rate remains
S0. A numerical calculation of the neutral stability
curves and growth rates when the magnetic field is ap
plied normal (N) to the interface is presented in Fig. 8b
for a range of applied field strength and layer thickness
ratios. The plotted results have also been maximized
over the wavenumber range, a = [0 5]. That the
N configuration is "more destabilizing" than the parallel
(P) field orientation, as can be seen by comparison with
the P neutral curve (dotted) and by the reduced growth
rates, approximately onehalf of the equivalent N values.
It is noteworthy that for n < 1, corresponding to a thick
bottom layer of ferrofluid, P fields do not destabilize the
interface even at the largest Ha presented. Moreover, for
both P and N orientations, instability properties become
insensitive to the relative layer thickness for n > 5. Be
haviors such as this are directly relevant to the design of
ferrofluid flows having a set of desired stability proper
ties.
In related work on leaky dielectric fluids (Uguz et al.
2008; Uguz & Aubry 2008) the stability of microfluidic
channel flows was examined as a function of the rela
tive permeability (and conductivity) of the fluids. Those
works presented regimes of "destabilization," defined
as regions in parameter space where the (electric) field
increases the growth rates of the interfacial mode even
if the growth rate remains negative. Destabilization in
Uguz et al. (2008); Uguz & Aubry (2008) was mea
sured by examining the Maxwell stress contributions to
the interfacial stress balance in the limit Re 0, in
cluding the contributions of the base field solution. Li
et al. (2007) computed a similar quantity, but by us
ing longwave asymptotics. An equivalent result in the
present work is not possible because the base solution
and the perturbations depend on the nonlinear formula
(10). Instead, the following calculation was performed:
the growth rate of the most unstable mode was com
puted for zero magnetic field, Ha 0, and for a small
parallel magnetic field, both calculations done over a
range of parameters. The difference of these two growth
rates can be interpreted as a measure of the field de
stabilization (at zero field) and is depicted in Fig. 9a in
n m parameter space. Equivalent calculations per
formed around different values of the imposed field (i.e.
around H = 1, 2, 3,...) each give a different result as
a result of (10). Nevertheless, for comparison and in an
attempt to quantify this property, we have also computed
the neutral curve in n m space (Fig. 9b) for Ha = 5,
chosen arbitrarily. While the two results in Fig. 11 are
qualitatively similar, it is also clear that they are not sim
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
14
1 2
08
oa
max
06
04
02
0 v 5 H10 15 20
H
a
1 (b) N
15 2 25 3 35 4 45 5 55 6 65
n
Figure 8: (a): Wavenumber of maximum instability
growthrate for Parallel field (upper blue data, triangles)
and Normal field (lower red data, circles) at n = 2,
m = 1.5 and = 1/3; for parallel field there is no
instability below Ha 5; (b): Growth rates and neu
tral stability curves (heavy lines) in Ha n space for
m = 1.5, = 1/3 and field applied Normal to the in
terface; the heavy dotted line in (b) corresponds to the P
neutral curve, for comparison.
ilar enough that "de,i.iIlii,.iiiun i predictions would be
useful for nonlinear magnetic material. To accurately
predict instability of magnetic fluid flows, it is necessary
to perform an exact eigenvalue calculation, accounting
for the nonlinearity of the magnetic material.
Conclusions
We have presented an experimental and numerical anal
ysis of micron scale bubbles in a ferrofluid, visualized in
high detail by using a new Xray phase contrast method
and simulated using a magnetic VOF method. Under
applied magnetic fields up to 0.2 T, such small bubbles
are only slightly deformed and remain nearly spherical.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
15 2 25 3 35
Figure 9: Stability measures in n m space: (a) shows
regions of stabilization (S) and destabilization (D) for
parallel field and = 1/3; (b) shows the growth rates
and neutral curve (heavy line) for parallel field of H,
5.
Nevertheless, our experimental technique was capable
of resolving enough detail to distinguish between the
slightly different susceptibilities of the two ferrofluids,
EMG607 and EMG707.
Our VOFtype, NavierStokesbased code is not ide
ally suited to compute static equilibrium solutions of
slightly deformed single bubbles. Moreover, the surface
tension is large at these scales and the ParkerYoungs
method of computing curvature produces errors that ac
cumulate over time to form spurious velocities which al
ter the equilibrium shape. In a separate work, a high
order height function based curvature method is being
implemented to eliminate these errors.
The dynamics of multiple bubbles is dominated
by magnetophoresis, originating in the approximately
dipole fields induced by nearby bubbles, and by hydro
dynamic forces due to their motion, mainly drag. Drag
limited magnetophoretic dipole attraction has been ver
ified by comparing numerical and experimental results,
and by comparing both to the predicted scaling of the
aggregative velocity. The deviations from the scaling
prediction seen in the simulations are due to hydrody
namic effects, such as added mass, that are captured by
the DNS model but obscured in the data. Although we
have restricted our attention to gas bubbles, the results
here should apply also to droplets of a second liquid or
to solid particles immersed in ferrofluid.
Subsequent experimental and numerical efforts will
also examine the field dependent anisotropic viscosity
induced by the presence of the particle chains that are
clearly visible in the experimental images.
The linear stability problem for interfacial flow in
a microchannel, where one fluid is a ferrofluid having
a nonlinear magnetic permeability, indicates that both
parallel and normal oriented magnetic fields can pro
duce instability. Parallel fields generally lead to smaller
growthrate instabilities at smaller wavelength. Notably,
stability results were largely insensitive to whether the
ferrofluid was the top or bottom layer, but sensitive to
the layer thickness and viscosity ratios. Stability re
sults clearly identify choices of ferrofluid (initial) sus
ceptibilities and layer thicknesses where stability is as
sured as well as other choices where strong instability
is guaranteed. We can estimate that choosing channel
sizes (1 mm) and flow rates (1 pl/sec) consistent with
those in the experiments in (Ozen et al. 2006; Song et al.
2007; Hatch et al. 2001), places us in the correct regime
Re 0 (1). Upon applying a B field magnitude of about
0.1 T (corresponding to fields H ~ 0(105) A/m and
typical of a rare earth permanent magnet) and assum
ing typical ferrofluid susceptibility Xi = 2 (correspond
ing to = 1/3), we expect growthrates ci ~ 0(1)
are possible, especially for larger layer thickness ratios,
such as n z 5. Ten efolding times of growth then oc
cur in ten (nondimensional) distance units downstream,
equivalent to two channel widths or about 2 mm. In
even smaller channels the same result applies as long as
flow rates are reduced so that Re 0 (1). Such rapid in
stability growth has great potential to generate droplets
from a fluid stream, although a more conclusive predic
tion of droplet generation demands a full nonlinear cal
culation.
Aubry and collaborators (Aubry & Singh 2008;
Aubry et al. 2008) have shown that a twofluid EHD
interface can be exploited as a means to direct self
assembly of micro and nanoparticles. This possibil
ity is of great value to micro and nano manufacturing,
most notably chip fabrication. A similar technique had
also been proposed in the context of ferrofluids (Yellen
et al. 2005) where the need for a complete understand
ing, including flow, has been emphasized in Song et al.
(2007).
Acknowledgements
Use of the Advanced Photon Source at Argonne Na
tional Laboratory was supported by the U. S. De
partment of Energy, Office of Science, Office of Ba
sic Energy Sciences, under Contract No. DEAC02
06CH11357
References
Afkhami S., Renardy Y., Renardy M., Riffle J.S., and
St. Pierre T., Fieldinduced motion of ferrofluid droplets
through immiscible viscous media, JFM 610, pp. 363
380, 2008
Alexiou C., Lubbe A.S. and Bergemann C., Clinical ap
plications of magnetic drug targeting, J. Surgical Res.,
95 (2), pp. 200206, 2001
Alexiou C., Jurgons R., Schmid R., Hilpert A., Berge
mann C., Parak F., and Iro H., In vitro and in vivo
investigations of targeted chemotherapy with magnetic
nanoparticles, J. Magnetism Magnetic Mat. 293, pp.
389393,2005
Aubry N., Singh P., Janjua M. and Nudurupati S., Micro
and nanoparticles selfassembly for virtually defectfree,
adjustable monolayers, PNAS 105, pp.3711, 2008
Aubry N. and Singh P, Physics underlying controlled
selfassembly of micro and nanoparticles at a twofluid
interface using an electric field, Phys. Rev. E 77, 056302,
2008
Bashtovoi V., Kovalev M. and Reks A., Instabilities of
bubbles and droplets flows in magnetic fluids, Journal of
Magnetism and Magnetic Materials, 289, pp. 350352,
2005
Boeck T., Li J., L6pezPag6s E., Yecko P and Zaleski S.,
Ligament formation in sheared liquidgas layers, Theo
retical and Computational Fluid Dynamics, 21, pp. 59
76, 2007
Coyajee E.R.A., Delfos R., Slot H.J. and Boersma B.J.,
Interface Using a LevelSet/VolumeOfFluid Method
with multiple Interface Marker Functions, chapter 1.
Number VI in Direct and LargeEddy Simulation.
Springer, 2006
Crowe C.T., Sommerfeld M. and Tsuji Y, Multiphase
Flows with Droplets and Particles. CRC Press, 1997
de Gennes PG. and Pincus PA., Pair correlations in a
ferromagnetic colloid, Zeitschrift Physik B Cond. Mat.
11 (3), pp. 189198, 1970
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Drenckhan W., Elias F., Hutzler S., Weaire D., Janiaud
E. and Bacri J.C., Bubble size control and measurement
in the generation of ferrofluid foams, J. Appl. Phys., 93,
pp. 1007810083,2003
Ganguly R., Zellmer B. and Puri I.K., Fieldinduced
selfassembled ferrofluid aggregation in pulsatile flow,
Phys. Fluids 17 (9), 097104, 2005
Gangulya R., Gaind A.P., Senb S. and Puri I.K., Ana
lyzing ferrofluid transport for magnetic drug targeting,
J. Magnetism Magnetic Mat. 289, pp. 331334, 2005
Gueyffier D., Nadim A., Li J., Scardovelli R. and Zaleski
S., Volume of fluid interface tracking with smoothed sur
face stress methods for threedimensional flows, J. Com
put. Phys. 152, pp. 423456, 1999
Hartshorne H., Backhouse C.J. and Lee W.E., Ferrofluid
based microchip pump and valve, Sensors and Actuators
B,99,pp.592600, 2004
Hatch A., Kamholtz A.E., Holman G., Yager P and
Bohringer K.F., A ferrofluidic magnetic micropump, J.
Microelectromech. Sys. 10, pp. 215221, 2001
C. F. Hayes, Observation of association in a ferromag
netic colloid, J. Colloid Interf. Sci., 52 (2): 239243,
1975
Holm C. and Weis J.J., The structure of ferrofluids: a
staus report, Curr. Op. Colloid Interf. Sci. 10 (34), pp.
133140, 2005
Hooper A.P, The stability of two superposed viscous
fluids in a channel, Phys. Fluids A 1 (7), pp. 11331142,
1989
Horng H.E., Hong, C.Y., Yang S.Y. and Yang H.C.,
Novel properties and applications in magnetic fluids, J.
Phys. Chem. Solids 62, pp. 17491764, 2001
Jeyadevan B., Torigoe T., Nakatsuka K., Nakatani I., Fu
jita T. and Oka H., Xray image analysis and ultramicro
scope study of magnetic fluids used as working liquid in
heat transfer experiments, J. Appl. Phys. 85, pp. 5968
5970, 1999
Jones G.A., Aggregation of waterbased magnetic liq
uids observed with the polarizing microscope, J. Phys.
DAppl. Phys. 18, pp. 12811290, 1985
Jones T.B., Electromechanics of Particles. Cambridge
University Press, 1995
Korlie M.S., Mukherjee A., Nita B., Stevens J.G.,
Trubatch A.D. and Yecko P., Modelling bubbles and
droplets in magnetic fluids, J. Physics: Cond. Mat. 20
(20), 204143, 2008
Kreuger D., Review of agglomeration in ferrofluids,
IEEE Trans. Magnetics 16 (2), pp. 251253, 1980
Lafaurie B., Nardone C., Scardovelli R., Zaleski S. and
Zanetti G., Modelling merging and fragmentation in
multiphase flows with Surfer, J. Comput. Phys. 113, pp.
134147, 1994
Lavrova O., Matthies G., Mitkova T., Polevikov V. and
Tobiska L., Numerical treatment of free surface prob
lems in ferrohydrodynamics, J. Physics: Cond. Mat. 18
(38), pp. S2657S2669, 2006
Li F., Ozen O., Aubry N., Papageorgiou D.T. and
Petropoulos P.G., Linear stability of a twofluid interface
for electrohydrodynamic mixing in a channel, J. Fluid
Mech. 583, pp. 347377, 2007
Lin H., Storey B.D., Oddy M.H., Chen C.H. and Santi
ago J.G., Instability of electrokinetic microchannel flows
with conductivity gradients, Phys. Fluids 16, pp. 1922
1935, 2004
Melcher J.R., Fieldcoupled surface waves,: A compar
ative study of surfacecoupled electrohydrodynamic and
magnetohydrodynamic systems. (M.I.T. Press, Cam
bridge, MA, 1963)
K. Nakatsuka, B. Jeyadevan, Y. Akagami, T. Torigoe,
and S. Asari, Visual observation of the effect of mag
netic field on moving air and vapor bubbles in a mag
netic fluid, J. Magnetism Magnetic Mat. 201, pp. 256
259, 1999
Odenbach S., editor, Colloidal Magnetic Fluids, Vol
ume 763 of Lecture Notes in Physics. (SpringerVerlag,
Berlin, 2009)
Ozena O., Aubry N., Papageorgiou D.T. and Petropolous
P.G., Monodisperse drop formation in square mi
crochannels, Phys. Rev. Lett. 96, pp. 144501, 2006
Pankhurst Q.A., Connolly J., Jones S.K. and Dobson J.,
Applications of magnetic nanoparticles in biomedicine,
J. Physics D 36, pp. R167, 2003
Renardy Y., The thinlayer effect and interfacial stability
in a twolayer couette flow with similar liquids, Phys.
Fluids 30 (6), pp. 16271637, 1987
Rinaldi C., Chaves A., Elborai S., He X.T. and Zahn M.,
Magnetic fluid rheology and flows, Curr. Op. Colloid In
terface Sci. 10, pp. 141157, 2005
Rosensweig R.E., Ferrohydrodynamics, Cambridge,
Cambridge, 1985
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Scardovelli R. and Zaleski S., Direct numerical simu
lation of freesurface interfacial flow, Ann. Rev. Fluid
Mech. 31, pp. 567603, 1999
Song W.B., Ding Z., Son C. and Ziaie B., Aqueous mi
crodrop manipulation and mixing using ferrofluid dy
namics, Appl. Phys. Lett. 90, pp. 092501092503, 2007
Stratton J.A., Electromagnetic Theory, Adams Press,
2007
Taketomi S., Takahashi H., Inaba N. and Miyajima H.,
Experimental and theoretical investigations on agglom
eration of magnetic colloidal particles in magnetic flu
ids, J. Phys. Soc. Japan 60, pp. 16891707, 1991
Tryggvason G., Scardovelli R. and Zaleski S., Di
rect Numerical Simulations of GasLiquid Multiphase
Flows. Cambridge, Cambridge, to be published 2010
Uguz A.K. and Aubry N., Quantifying the linear stabil
ity of a flowing electrified twofluid layer in a channel
for fast electric times for normal and parallel electric
fields, Phys. Fluids 20, 092103, 2008
Uguz A.K., Ozen 0., and Aubry N., Electric field effect
on a twofluid interface instability in channel flow for
fast electric times, Phys. Fluids 20, 031702, 2008
Wang Z., Holm C. and Muller H.W., Boundary condition
effects in the simulation study of equilibrium properties
of magnetic dipolar fluids, J. Chem. Phys. 119, pp. 379
387, 2003
Wiedenmann A. and Heinemann A., Fieldinduced or
dering phenomena in ferrofluids observed by small
angle neutron scattering, J. Magnetism Magnetic Mate
rials 289, pp. 5861, 2005
Yecko P., Disturbance growth in twofluid channel flow:
The role of capillarity, Intnl. J. Multiphase Flow 34 (3),
pp.272282,2008
Yecko P., Effect of normal and parallel magnetic fields
on the stability of interfacial flows of magnetic fluids in
channels, Phys. Fluids 22, 022103, 2010
Yellen B.B., Hovorka 0. and Friedman G., Arranging
matter by magnetic nanoparticle assemblers, PNAS 102
(25), pp. 88608864, 2005
Yiantsios S.G. and Higgins B.G., Linear stability of
plane poiseuille flow of two superposed fluids, Phys.
Fluids 31, pp. 32253238, 1988
Yih C.S., Instability due to viscosity sII.iiilik.ii.ni J.
Fluid Mech. 27, pp. 337352, 1967
Zborowski M. and Chalmers J.J., Magnetic cell sorting,
Meth. Molec. Bio. 295, pp. 291, 2005
