7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
CFD simulation of shear induced particle migration
Yezaz Ahmad, Madeti Girish Kumar and Anugrah Singh*
Department of Chemical Engineering
Indian Institute of Technology Guwahati
Guwahati 781039 INDIA
Email: anugrah@iitg.ernet.in
Keywords: Shear induced migration, Computational fluid dynamics, Suspension flow
Abstract
Flow of concentrated suspensions is encountered in various process industries. It is now well known that shear induced
migration is observed in solidliquid suspensions in which initially wellmixed particles in suspensions when subjected to
inhomogeneous shear migrate from region of high shear rate to the region of low shear rate giving a nonuniform concentration
distribution. Shear induced migration phenomena has important implications in manufacturing processes where the properties
of final products are greatly influenced by the degree of particle distribution. The models developed to explain this phenomenon
basically fall into two categories, namely, suspension balance model and diffusive flux model. These models approximate the
suspension rheology as a generalized Newtonian fluid with concentration dependent rheology. The models have been
successfully applied to explain migration behavior in simple geometries such as straight channels, axisymmetric pipe flow etc.
However, the industrial processes often involve flow in complex 3D geometries. In this work we provide a framework to use
these models in conjugation with CFD simulations so that they can be applied to any general flow situations. The model
equations are recast into a form which can fit to the standard convectiondiffusion equation which can be solved by any
commercial CFD solver. The simulations were validated with available experimental and theoretical results for 2D channel
flow and 3D flow through an asymmetric Tjunction bifurcation channel. In another application the volume of fluid (VOF)
model for multiphase flow is used in conjugation with diffusive flux model of shear induced migration to study free surface
flow of suspension. The results from simulations are in good agreement with the experimental data.
Introduction
Suspensions comprising of solid particles suspended in a
fluid medium are commonly used in various chemical and
biological processes. An understanding of their flow
behavior is the major requirement for manufacturing,
processing and transport of these materials. Knowledge of
the dynamics of such multiphase systems is thus of
paramount importance to understand the various interesting
phenomena exhibited by them. Shear induced migration of
particles is one such phenomenon which has received wide
attention. Recent studies in the past few years and
advancement in experimental and analytical techniques
have increased the understanding of the physics of
suspensions. However, to understand the behavior of
suspensions in real systems often have cases which do not
have analytical solutions. The experimental evidence
strongly supports the view that suspensions behave
macroscopically as nonNewtonian fluids whose
theological properties are influenced by a large number of
factors (Zarraga et al. 2000; Singh & Nott 2003).
GadalaMaria and Acrivos (1980) provided the first
evidence of shear induced migration where they observed
the steady decrease of the viscosity of a concentrated
suspension in a Couette viscometer. Later, Leighton and
Acrivos (1987) showed by further experiments that this was
due to migration of particles from high shear rate regions to
lower shear rate regions. They went on to propose a
theoretical model for the process which explained the
observation of decrease in viscosity. They proposed a
diffusion equation for particles which was solved in
conjunction with continuity and momentum equations for
the entire suspension. They also proposed that irreversible
inter particle interaction generated due to surface roughness
of particles was responsible for migration. Dynamic
simulations of the pressuredriven flow in a channel of a
nonBrownian suspension were conducted by Nott and
Brady (1994). They performed Stokesian Dynamic
Simulation (Brady and Bossis, 1988) for nonBrownian
suspension at zero Reynolds number and demonstrated that
particles did migrate from higher shear rate to low shear rate.
By demonstrating that the irreversible migration took place
even for perfect hard spheres, they discarded the notion that
surface roughness was not the sole criteria for shear induced
migration. They also put forward the 'Suspension Balance
Model' based on mass, momentum and energy balances for
the suspension by giving a suspension temperature to
account for the non local suspension behavior. The
argument of surfaceroughnessindependent particle
migration was further supported by the work of Brady and
Morris (1997). They concluded that the asymmetry in
particle phase distribution leads to the observed
nonNewtonian behavior in suspensions. They proposed
that the Brownian motion and shear forces create
asymmetry in particle phase distribution and leads to
anisotropic normal stresses. Further, Morris and Boulay
(1999) proposed that the migration of particles could be
explained by considering the effect of
normalstressdifferences created due to asymmetry in
particle phase distribution. They showed that by
considering the effect of normal stress differences in the
Suspension Balance Model of Nott and Brady (1994), the
results obtained agree well with the experimental results.
Attempts were also made to experimentally verify the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
normal stress differences in the viscometric flow of
suspensions. A systematic study on normal stresses in
Stokesian suspensions was done by Zarraga et al. (2000).
Singh and Nott (2003) presented experimental
measurements of normal stresses in sheared Stokesian
suspensions through more accurate experiments. They
observed that these stresses could come only if there is
anisotropy in the particle phase distribution or the shear
field thereby attributing the nonNewtonian behavior of
suspension to the anisotropy in particle phase distribution in
suspensions.
Both diffusive flux model and suspension balance model
has been applied to explain migration in Couette and
pressure driven flow in channels and tubes. A number of
studies have worked to extend the original diffusive flux
model of Phillips et al. (1992) to more generalized
applications (Zhang and Acrivos, 1994; Subia et al., 1998;
Rao et al., 2002; Fang et al., 2002) to general conditions.
The diffusive flux model has been extensively tested in
Couette flow (Altobelli et al., 1991; Tetlow et al., 1998),
and has been widely used in numerical simulations of
suspensions flow in various geometries (Zhang and Acrivos,
1994; Fang and PhanThein, 1995; Subia et al., 1998). The
diffusive flux model is simpler than the suspension balance
model, and can be applied to the numerical simulation of
general flows with ease.
In this work we have studied particle migration in
suspension flow in rectangular bifurcation geometry which
is motivated by the recent experimental work of Xi and
Shapley i,21',, Such flows in bifurcation channels have
relevant applications in biological and technological
processes, such as laminar flow of blood in human blood
vessels. Another application is 'Plasma skimming', a
separation technique in which the blood flows through a
network of precisely controlled microfluidic bifurcations
resulting into the separation of red blood cells and white
blood cells. It also has application in drug administration
at vessel junctions which could have significant
physiological consequences (Dedrick, 1988). Thus, in the
systems involving bifurcations it is important to find out
bulk suspension partitioning and the particle partitioning in
the daughter branches downstream of the channel inlet.
For example, in one study Krogh (1921) observed that
blood cells are not distributed in the same proportion as
volumetric blood flow and attributed this discrepancy to the
nonuniform upstream distribution of blood cells. Recently,
Roberts and Olbricht (2003) found that in a rectangular
channel bifurcation, the partitioning of particles is different
from the partitioning of fluid. In case of Yjunction
bifurcations, particles preferentially enter the downstream
branch and receive the greater volumetric flow rate. In an
attempt to understand this partitioning mathematical models
were established to compute the cell distribution at
bifurcations (SchmidSchonbein et al., 1980; Perkkio and
Keskinen 1983; Fenton et al., 1985). These models all
assume that red blood cells follow streamlines at the
junction and require knowledge of the velocity profile, the
cell distribution profile and the dividing streamline surface.
Rong and Carr (1990) for a bifurcation composed of circular
tube branches found that the dividing surface shape is a
function of the ratio of branch diameters and the fractional
flow splits at the junction. At low flow Reynolds number
and when all branches have equal diameters, the separating
surface is virtually flat. Most of the studies mentioned
above involve particles which are comparable in size to the
channel width or tube diameter (B/a < 4) and present in a
dilute concentration. Here, B represents the channel width
or tube radius and 'a' the particle radius. In contrast,
bifurcation flows of small particles (B/a > 20) are less
studied through experiments and numerical simulations as
noted by Xi and Shapley (2008). They have performed
experiments in rectangular bifurcation geometry for B/a =
35 and said that continuum models with effects of
shearinduced particle migration could predict the
partitioning profiles. In this work we provide numerical
validation of the experimental observations of Xi and
Shapley (2008) using the diffusive flux model.
The phenomenon of particle migration in
suspensions is not only confined to flow in confined
geometries. Timberlake and Morris (2005) had performed
experiments on gravity driven free surface flow of neutrally
buoyant suspension down an inclined channel of constant
width. Singh et al. (2006) has studied free surface velocity
of neutrally buoyant suspension of uniform spheres in an
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
inclined channel flow. However, all these experiment had
reported only velocity profile and no experiments yet has
estimated the particle concentration near the free surface. In
another application of CFD simulation of shear induced
migration we have studied the free surface flow of
suspension in a slot coating application. To examine the
particle migration in free surface flow of suspension, we
had used volume of fluid (VOF) model.
Nomenclature
g gravitational constant (ms 2)
P pressure (Nm 2)
V Velocity (ms1)
N Flux due to collision frequency
N Flux due to spatial variation in viscosity
K K Phillips model constant
a Particle radius (m)
V Velocity of fluid in the cell (ms1)
Vw Plate speed (ms1)
Greek letters
T Shear stress (Nm 2)
Local shear rate (s 1)
q7 Suspension viscosity (Pas)
S Particle fraction
77 Fluid viscosity (Pas)
ap Fluid's volume fraction in the cell
Diffusive Flux Model
The suspensions are considered as continuum fluid whose
viscosity depends on the solid fraction.
The steady state continuity and momentum conservation
equations for the suspension flow are given as,
V (pu) = 0, (1)
V.(puu)=Vp + V. (2)
In diffusive flux model, an additional diffusion equation is
coupled with this to predict the particle migration.
The steady species mass conservation equation is,
V (puo) =V Nt (3)
Where U is the bulk suspension velocity and 5 is the
particle volume fraction, p is the density of suspension,
p is the suspension pressure, Nt is total diffusive flux and
Y denotes the stress tensor which is modeled as
generalized Newtonian fluid which has the following
constitutive form,
Y = qr(o) E, (4)
where E is the rate of strain tensor. In the eq. (3), left hand
side is convective particle flux and right hand side as
diffusive flux. The diffusive flux encompasses all the
diffusion mechanisms, such as Brownian motion,
sedimentation, gradients in suspension viscosity and
curvature induced migration. Here we consider the original
model of Philips et al., (1992) where Nt constitutes only of
two mechanisms. The particle flux due to spatially varying
interacting frequency is given as,
N = pK a2(2 Vj + j VO). (5)
The viscosity gradientinduced migration flux is
N pK a2 02 I (V ) (6)
a o
Nt = N + N,.
Diffusion coefficients Kc and K are empirical constants
determined experimentally. As discussed by Tetlow et al.
(1998), better predictions can be made with the diffusive
flux models if the diffusivity coefficients, K, and K,
are chosen as functions of the local volume concentrations.
In a recent study, Ingber et al.(2008) incorporated two
modifications in the diffusive flux model. First, the
diffusion coefficients are modelled as linear functions of the
a Vyl
nonlinearity parameter e .
S+ YNL
SHowever, this
modification in itself leads to essentially cubic scaling of
migration on the particle radius which is somewhat larger
than one found in Couette experiment. To reduce the scaling
to match experimental results in the widegap Couette,
velocity slip boundary conditions are imposed at the wall.
This was not done in our numerical simulations rather the
constant values chosen by Philips et al., (1992) were used,
namely, K = 0.41 and K, = 0.62.
Numerical Scheme
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
To solve the above set of equations we have used a
commercial software package Fluent, which is based on the
FVM (finite volume method). The FVM method has been
successfully used by Fang and PhanThien (1995) for
diffusive flux model. Miller and Morris (2006) have used
FVM method for suspension balance model. The general
purpose FVM solvers such as Fluent, besides solving
continuity and momentum equation can also solve
additional scalar transport equation. The steady state
general scalar transport equation is given as,
v. (puo) = v.(ro) + S,. (8)
Where, F is the diffusion coefficient and S is the
source term.
F= pK a2 2 d (9)
i do
S,= V.(Nc). (10)
Note that the diffusion coefficient has units of kg and
m.s
source has units of kg to maintain dimensional
min .s
consistency. To the particle conservation equation which is
not in the standard form for Fluent, we set the
righthandside term of eq. (3) as sum of the diffusive flux
term and source term. These terms which are given below
are then modeled through user defined subroutines.
1 dq (11)
F= pK a2 2 12
1 do
S, V (Nc). (12)
In the pressure driven flows in channel or pipe the local
shear rate is zero at the centre, this makes the diffusivity in
our formulation to be zero there. This leads to numerical
instability and convergence problems as 0 > 0 To avoid
this we have used the nonlocal stress concept of Miller and
Morris (2006) where a small but constant nonlocal
contribution depending on particle size is added to the local
shear rate. This reflects the fact that the RMS of y is larger
than the mean shear rate, owing to fluctuations resulting
from finite size of the particles of radius a. Thus the
nonlocal contribution depending on the mean shear rate
and particle size is,
YNL = as (E) s 
U
In the above equation y = m and a (E)= 2 where
B
umax is the centerline velocity and e=a/B for a channel of half
width B. The finite YNL values satisfy YNL << (x) except
where y  0 giving the model the desired effect of
influencing results only near the centerline.
The governing equations (mass, momentum, and additional
scalar transport equation) are solved with segregated,
iterative multigrid solver. The SIMPLE scheme was
employed for the pressurevelocity coupling. The second
order accurate discretization schemes were used for all the
equations.
For free surface flow, we had used volume of fluid (VOF)
model. In this model, we solve the continuity and
momentum equation in each phase (i.e air and suspension)
separately. For suspension phase we solve equation (8)
along with continuity and momentum equations. Tracking
of the interface between the two phases is accomplished by
the solution of continuity equation for the volume fraction
of any of the phases which can be written as,
(ap) + V (avp) = 0 (14)
at
Where, ap is fluid's volume fraction in the cell, vp is the
velocity of fluid in the cell.
Validation of the numerical scheme
Before simulating the bifurcation problem, the numerical
method was validated by performing twodimensional
steady state simulation in a channel and the results were
compared with semianalytical solution given by Phillips et
al. (1992) and experimental data of Lyon and Leal (1998).
The length required to reach a steady configuration of a
suspension flow was calculated by formula given by Nott
and Brady (1994),
L 1 B L
B ,, 12 d(>) a
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The function d(,) accounts for the bulk particle
concentration dependence of the shearinduced coefficient
of diffusion. Leighton and Acrivos (1987) and later
Chapman (1990) have found that,
1 1
d(Q) =q (1 + e ) (16)
3 2
The boundary conditions for 2D channel are,
1. At inlet, average velocity and particle volume fraction
,,,, is specified.
2. At walls, no slip boundary condition is applied and
particle flux is taken to be zero.
3. At outlet, a fully developed outflow condition is
applied.
For the 2D simulations the relative particlesize (B/a) was
taken to be 18 which is same as in the experiments of Lyon
and Leal (1998) for channel flow. The computational
domain was discretized into 70,000 finite volume cells. We
have carried out simulations for the same flow rate as taken
in the experiments of Lyon and Leal (1998). Comparisons
of the CFD model predictions with experiments and
semianalytic solution for bulk concentrations 6,, = 0.3
and 06,, = 0.4 in planar Poiseuille flow is shown in Figs. 1
and 2 respectively. The velocity is normalized with
maximum velocity Vm, of a parabolic profile for
suspending Newtonian fluid with the same volumetric flow
rate. The velocity profile is in close agreement with the
semianalytical solution. The suspension velocity profiles
are blunted in comparison with Newtonian velocity profile.
The concentration profile also shows good agreement with
the experiments suggesting that the nonlocal formulation is
able to remove the sharp cusp observed for semianalytical
solution. The decrease in centerline value is due to addition
of nonlocal shear rate term. For the same size of particles,
as the bulk volume concentration increases more particles
accumulate near the midplane resulting in a larger local
concentration and suspension viscosity and therefore a
lower shear rate or a more blunted velocity profile. To check
the grid independence of solution we have carried out
simulations for two different grid sizes. Grid A had 1000
cells in flow direction and 70 cells in the gradient direction.
Grid B had 1300 cells in the flow direction and 80 cells in
the gradient directions. Simulations for both the grids are
identical which confirms grid independent solution. We
observe that near the channel wall our simulation predicts
little higher concentration. Very small grid size is expected
to remove this discrepancy. However, choosing very fine
grid size near the walls resulted into numerical divergence.
It is to be noted that near the wall the experimental data
show much lower concentration and this is not captured
even with the semianalytical solution. Proper choice of
KC and K7 in diffusion equation and providing
appropriate wall boundary condition for particle
concentration and velocity would reduce this difference
(Miller and Morris, 2006).
S0 '. 0o Experlmentaldata (LyonandLeal, 1998)
0 Semi analytcal siolution
0a 0 8 2D nurencal smulalon (gnid A)
o. 2D numeriencal simulation (gnd B)
06 0 06
o05 o 0'. ." .. .un>
Puc pfi(N d) ( Kr 0 0
o 0 05 1 05 0 05 1
y/B y/B
(a) (b)
Figure 1: Comparison of predicted velocity profile for
,0 = 0.3 and B / a = 18 with experimental data, with
semianalytical solution and with parabolic profile
0 rExpemea da (Ly an Lel, 1998)
 s~t~~ I '9]0
Figure 2: C comparison of predicted velocity profile for
= 04 and B / a= 18 with e ntal data, 9 th
semi nanaly soltical solution and with parabolic profile
Sspeion fa lolw in afld symmetric T u io
y/B 0 y/B
(a) (b)
Figure 2: Comparison of predicted velocity profile for
6,1 = 0.4 and B / a = 18 with experimental data, with
semianalytical solution and with parabolic profile
Suspension flow in asymmetric Tjunction
bifurcation channel
Asymmetric Tjunction bifurcation channel
represents a branching system in which flow divides into
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
two. In such geometries one would like to know the bulk
suspension partitioning and the particle partitioning in the
daughter branches. Figure 3 describes the computational
domain and finite volume mesh for the branched geometry
or asymmetric Tjunction bifurcation channel. The
dimensions of this geometry were similar to the one taken in
the experiments of Xi and Shapley i2,.'" i In the
experiments the length of the inlet section of the channel
was 3 cm which was connected to a rigid circular tube of
radius 0.95 cm and length 49.4 cm. They have reported that
the flow in the tube lead to particle migration to the tube
centre forming an inhomogeneous particle distribution at
the junction of tube and inlet section of the bifurcation
channel. In the simulations we have directly imposed the
velocity and concentration boundary condition which was
obtained from analytical solution of fully developed
suspension flow in the pipe [Phillips et al. (1992)]. The inlet
branch length and downstream branches lengths are
extended by the required steady state length estimated by eq.
(15). The inlet section had a length of 20 cm and the lengths
of downstream branches were 10 cm respectively. All other
dimensions and flow parameters taken were same as in the
experiments.
x
Figure 3: Computational geometry and finite volume grid
for the bifurcation channel. The inlet branch divides into
main branch and side branch at the bifurcation junction.
Velocity field
The contour planes of velocity magnitude are shown in
Figure 4. For comparison we have presented the contour
planes for Newtonian fluid flow and the suspension flow.
The figures 4a and 4d shows the front view in the xz plane
2 q
sl
at the centerline (y=A/2). Figures 4b and 4e are for the side
view (yz) plane in centerline of the main branch at x=1.5
mm. The corresponding side view for the side branch is
given in the figures 4c and 4f for the Newtonian fluid and
suspension flow respectively. It can be observed that the
Newtonian fluid shows nearly parabolic profile across the
narrow dimension (x) and the mostly flat profile along the
wider dimension (y). However for suspension flow we
observe blunt velocity profile in the inlet and main branch
for the front view. The side view for main and side branch
shows a different pattern marked by two peaks. The two
peaks are located near the walls and a valley in observed in
the centerline along the z axis. The maximum velocity in the
two downstream branches is peaked towards the outer walls.
Another interesting observation is that the Newtonian fluid
velocity profile in the inlet branch persists up to the center
of bifurcation where it divides into main and side branch
flow. This is indicated more clearly in figure 4b. The flow of
suspension presents a different behavior and we observe
that the deacceleration of the flow starts much before the
bifurcation (Figure 4d).
I I
( () (
Figure 4: Velocity contour planes for Newtonian fluid
(a,b,c) and suspensions (d,e,f) of average particle
concentration 0.4. Figures (a) and (d) shows the front view
in xz plane at the centerline (y=A/2). Figures (b) and (e) are
for the side view (yz) plane at x=1.5 mm for main branch
and figures (c) and (f) shows the side view (yz) plane for
the side branch at x = 8.5 mm. The flow rate for the
simulation with Newtonian fluid was 20.18 ml/min and for
the suspension flow it was 102.1 ml/min. The velocity unit
in the contour plot is in cm/s.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
To view the velocity profile across the channel we have
shown axial velocity map (xy plane) in the inlet, main and
side branch for suspension of 40% particle concentration.
The flow rate in the inlet branch was 102.1 ml/minute. The
slices shown in the figure are taken at the locations, z=7
mm for inlet branch (figure 5a) and z=7mm for the main and
side branches (Figures 5b and 5c). The axial velocity maps
show two symmetrical velocity peaks in the inlet branch.
The main and side branches also show these peaks which
are symmetric about the axis BB' but skewed about the axis
AA. The peaks in the velocity are closer to the outer wall
for both main and side branches. These observations are
similar to the experiments of Xi and Shapley 111 I',i Figure
6 shows comparison of the experimental and simulation
velocity profile in the spanwise directions along the
centerline of the channel. The velocity profiles clearly
exhibit the peakvalleypeak pattern as shown in contour
planes. In velocity magnitude, the inlet branch has the
greatest velocity while the main branch has slightly larger
velocity than the side branch. Near the channel walls, the
velocities are close to zero, due to noslip boundary. The
corresponding plot for the lateral direction is shown in
figure 7. Once again, the inlet branch exhibits greater
velocity magnitude than the two daughter branches.
Another interesting observation is that all the three branches
show different symmetric properties. The inlet branch has
velocity symmetry while the two branches slope in opposite
directions and thus the crossing of the two curves occurs.
Common to the two branches is that they both show
maximum velocity near the outer walls. We observe good
agreement with the experiments except for the inlet branch.
The reason for this difference is the inlet velocity and
particle concentration boundary conditions. The exact
velocity and concentration profile at the junction of circular
tube and rectangular channel inlet was not available to us.
Therefore we had taken the profiles for fully developed
flow in a pipe. It is to be noted that the profiles in the
experiments of Xi and Shapley was not fully developed due
to the short length of the tube. They have also reported that
the velocity measurements near the channel wall could be
erroneous due to air bubbles trapped. Overall there is good
qualitative agreement of the simulation and experimental
results.
'A
I
x A'
(b)
,A
,1 4
B'
S 2.
A 0
(c)
Figure 5: Axial velocity map (xy plane) in the inlet branch
for z=7 mm (a), main branch for z = 7 mm (b), and side
branch for z = 7mm (c).The average particle concentration
was 0.4 and the flow rate was 102.1 ml/min.
12
SInlet branch (CFD simulation)
t lain branch (CFD simulation)
10 A Side branch (CFD simulation)
0 Inlet branch (Expt.data); Xi and Shapley (2008)
o Main branch (Expt.data), Xi and Shapley (2008)
8 X Side branch (Expt.data); Xi and Shapley (2008)
o  
o Oo [o n \
> 4 ._El E
> 4
2
4
0
0 0.2 0.4 0.6 0.8
1 y/A
Figure 6: Comparison of the experimental (Xi and Shapley,
2008) and simulation velocity profile in the spanwise
directions along the centerline.
10
 lct blanch (CID smmidatlon)
Main branch (CFD simulation)
Side branch (CFD simulation)
O Mlet branch (Expl data), Xt and Shaplw (2008)
o Main branch (Expt.data), Xi and Shapley (2008)
x Side blanch (Expt data), Xi and Shapley (2008)
6
t o
SP . \
O 02 04 06 08
1 x/2B
Figure 7: Comparison of the experimental (Xi and Shapley,
2008) and simulation velocity profile in the lateral
directions along the centerline.
Particle concentration field
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Axial concentration slices were taken at different heights in
the inlet branch (slices at z = 10.5, 7.0, 3.5mm) and they
show no apparent variance in particle concentration. The
same is true of slices in the daughter branches (slices at z =
10.5, 7.0, 3.5mm). This confirms that fully developed
conditions were reached in the simulations. Slice at z = 3.5
mm in the inlet branch and slice at z = 7.0 mm in the main
and side branch are chosen to represent the concentration
slices. Xi and Shapley (2008) have also presented these
profiles at the same location. Figure 8 shows particle
concentration distribution maps obtained from numerical
simulations at these slices. In all three branches, the middle
region along the long dimension (AA) shows a higher
concentration than the sides. We also observe that in the
daughter branches the particle concentration is higher near
the inner walls where the velocity was lower compared to
the outer walls. The development from the initially uniform
particle distribution to this distribution pattern indicates
particle migration from the edges to the center occurring
along the flow loop in experiments and also in simulation.
The absolute values of particle concentration mainly fall in
the region of 0.350.50 for experiments and 0.350.56 for
simulation. The simulation results has higher concentration
region due to extension of the inlet branch which provides
more length for the migration of particles.
IA IA IA 05
r r 
B
L... 'A'
(a)
B' B B' )IB'
L
'A'
(b)
A'
(c)
Figure 8: Axial particle concentration map (xy plane) in the
inlet branch for z=7 mm (a), main branch for z = 7 mm (b),
and side branch for z = 7mm (c).The average particle
concentration was 0.4 and the flow rate was 16.1 ml/min.
Figure 9 presents the comparison of concentration
profiles along spanwise (y direction) at the above
mentioned slices. The simulation and experimental profiles
of the three branches show the same symmetric pattern. The
exceptionally low values on one side of the main branch are
due to image artifacts caused by air bubbles trapped in the
lowerleft channel corner in the experiments of Xi and
Shapley i,2'"II Figure 10 presents the comparison of
concentration profiles along the lateral (x direction) at the
above mentioned slices. The inlet branch profile is
essentially symmetric, while the two daughter branch
profiles are asymmetric in the x direction. In both daughter
branches, high particle concentration occurs near the inner
wall and the regions near the outer wall have low particle
concentration. In simulation side branch profile is virtually
flat as compared to experimental profile. But at walls it has
high concentration due to spread of high concentration
region that has been directed along the wall from inlet
branch. An increase in inlet branch length would reduce the
high concentration as some particles migrate towards other
side. Gradients in particle concentration are also responsible
for the noticeably asymmetric lateral concentration and
velocity profiles acquired from the two downstream
branches, even though the inlet channel has symmetric
lateral concentration and velocity profiles.
O6
05
4 ,.
cl I'7
0 Inlet branch (CFD simulation)
Main branch (CI) simulation)
Side blanch (CFD simulation)
0 1 0 Inlel branch (Expl dlata), Xi and Shapley (2008)
n Main branch (.xpt data). Xi and Shaploy (2008)
SSide branch (Expt data), Xi and Shapley (2008)
0 01 02 03 04 05 06 07 08 09 1
1 /A
Figure 9: Spanwise particle distribution profile along the
central lines (Y direction)
05 .
 
 
h let branch (CFD smnlation)
0 2  Main branch (CHI) simulation)
A Side branch (CFD shnulation)
0 Inlet branch (Fxpt data); Xi and Shaplev (2008)
0 1 Main branch (Expt data), X1 and Shapley (2008)
x Side branch (,xpt data); Xia nd Shaplev (2008)
0 02 0'4 0'6 08 1
1 x/2B
Figure 10: Spanwise particle distribution profile along the
central lines (X direction)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
As shown in Fig. 9, the inlet branch exhibits spatial
inhomogeneity in the particle concentration. At the center
maximum particle volume fractions for simulation is seen
and which has a value of =0.56 while the wall regions
have a minimum value of 0 =0.35. In the experiments by
Xi and Shapley (2008) the maximum concentration value
was =0.53 at the centre and minimum value =0.35
was observed near the wall. When the suspension leaves the
inlet branch, the suspension diverges into two parts along
the dividing streamline, which crosses the concentrated
center. The asymmetry is formed because the concentrated
inlet channel center flows into the inner wall region of the
daughter branches and the less concentrated region flows
into the outer wall regions of the daughter branches. In the
velocity profiles sharper gradients appear in the regions of
low concentration near the outer walls of the downstream
branches. It is obvious that the difference in the Newtonian
and suspension velocity profiles arises from the tendency of
the particles to redistribute into a nonuniform
configuration during flow. We observe in contour planes
that regions of high particle volume fraction and low shear
rate magnitude coincide. For example, near the short axis
centre (BB'), the concentration is high, while the shear rate
is low, while near the edges in the y direction, the
concentration is low and the shear rate is high. This
observation is consistent with the concept in numerical
model and does predict the same. The downstream
asymmetry of the lateral concentration profiles was also
detected and explained by Roberts and Olbricht (2006) in
their study of a Yjunction microfluidic bifurcation
suspension flow. They have explained that the suspension
near the inner walls of the downstream branches originates
from the particlerich central region of the upstream branch.
The suspension near the outside walls of the downstream
branches originates from the particlefree region near the
inlet branch walls. This happens when the particle size is
comparable to the channel width. The downstream inner
and outer regions inherit different parts of the upstream
branch, thus producing the particle centre profile skews in
the daughter channels. Our numerical study and
experiments of Xi and Shapley (2008) shows similar results
1
even for small particles (compared to the channel width) for
concentrated suspensions.
Particle migration in free surface flow
To study the particle migration in free surface flow
we have considered a simple slot coating problem as shown
in figure 11 where the suspension from the die falls to a
bottom plate moving with a speed of V, The distance
between the moving plate and the die exit is h. After
traveling some distance between the moving plate and die
the suspension encounters free surface flow as shown in
figure 12. The numerical parameters for this problem are
given in table 1.
For this problem, inlet velocity and particle
concentration is given as fully developed profile which was
obtained from semi analytical solution of Phillips model.
inlet
Slot exit
h ; . o. an .
Figure 11: Schematic of slot coating system
Table 1: Numerical parameters
Gap size of the slot (d) 0.3 mm
Gap size of the coating bed (h) 0.19 mm
Average inlet velocity 0.2 m/s
Contact angle 20
Surface tension 0.05 N/m
Liquid viscosity 0.05 Pas
Particle size (a) 8im
Particle loading 0.4
Plate speed (V,) 0.36 m/s
Velocity profile
First of all, we discuss about the velocity profiles as we
know that gradient in shear rate causes migration of
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
particles. Figure 13 shows the velocity profile at various
locations. We have analyzed velocity profiles at 5 different
locations, one at distance 2h inside exit of the slot, second at
the exit of slot, third at a distance 2h from the exit of the slot,
fourth at a distance 4h from the exit of the slot and the last
one is at a distance 6h from the exit of the slot.
Figure 12: Contours for phases. The red region is
suspension phase and yellow is air phase.
It is observed that the velocity profile at the free surface is
flatter while the velocity inside and at the exit of the slot has
its maximum at moving plate and minimum at the other side.
Thus we see that in the free surface region the driving force
of the migration i.e. the shear rate gradient becomes very
small.
Figure 14 shows plot of particle phase fraction
across the film thickness in the free surface region.
Interestingly, we had also observed a minute amount of air
that was trapped near the moving plate. The first two points
of phase fraction plot at 4h from the exit of slot shows that
the volume fraction is less than one which means that there
is the presence of air traps on the plate. This can also happen
in the real situations. It was observed that the air pocket was
localized to very small region and this didn't affect the
velocity profile significantly as seen in figure 13.
Figure 13: Velocity profile at five different locations
y/h
Figure 14: Plot of phase fraction at 2h and 4h from the exit
of the slot
Particle distribution
To analyze particle distribution, we consider the
contour map shown in figure 15. We have found that
particle distribution profile inside the slot is parabolic in
shape, but as soon as the suspension comes into contact with
the moving plate the peak in concentration has shifted
towards the moving plate side. After this the suspension
flows between a moving plate and a stationary wall for
some distance and we encounter a region of constant shear
rate throughout the gap. Due to this reason there is no
further change in particle distribution till free surface is
encountered. In the free surface region the shear rate is
maximum at the moving plate and zero at the free surface
since the free surface will have zero shear stress. This is also
clear from the velocity profiles shown in figure 13.
However as soon as the suspension comes to the free
surface region the film thickness decreases and the particle
concentration is more homogeneously distributed. From
application point of view it will be of interest to know the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
parameters which affect the particle distribution near the
moving plate. A good coating will require the particles to be
homogeneously distributed near the moving plate. The
speed of the plate is one such parameter which can have
influence on the film thickness as well as particle
distribution. The other parameters could be particle radius,
slot width. The studies of these parameters on particle
migration in slot coating is subject of our future studies.
5 .2801
5. 65cI
6.3 01
4.39e01
4.7OBll
a.TUe
3.14e 1
2.82e01
220eU
g.42612
6.28eO2
3.14 802
0.0l e'0
Figure 15: Contour map of particle concentration
Conclusions
We have used the diffusive flux model to carry out CFD
simulation of suspension flow in an asymmetric Tjunction
bifurcation channel. Through the conservative finite volume
method we are able to demonstrate the close agreement
between the simulation and experimental results of Xi and
Shapley 121" 11 I Qualitatively (perhaps quantitatively) the
effect of nonuniform concentration distribution including
the detailed concentration and velocity profiles in each
branch, the position of the dividing stream surface, and the
partitioning of the particles between the downstream
branches has been captured well by numerical simulation.
Similar to the experimental observations we also observed
in our simulations that the concentration distributions are
almost equally partitioned between the two downstream
branches. This indicates motion of particles across the
dividing stream surface occurs at the bifurcation section.
Because the particles are small compared to the channel
width, the time scale of direct particle interactions is too
long to produce a shift of the observed magnitude, instead
of that, enhanced spreading of high concentration (and
therefore high viscosity) regions of suspension toward the
side branch. Rearrangement of particles changes the
concentration and velocity profiles which have been well
captured in our numerical simulations. In another
application of CFD simulation using diffusive flux model
we have studied the particle migration in free surface flow.
The volume of fluid (VOF) model of multiphase flow was
coupled with the flow simulations to predict film thickness
as well as velocity and particle concentration profile in the
suspension film.
Acknowledgements
The authors acknowledge financial support from DST
through grant no. SR/CE/38/2006.
References
A. Krogh, Studies on the physiology of capillaries: II. The
reactions to local stimuli of the bloodvessels in the skin and
web of the frog, J. Physiol. (London) 55 (1921) 412422.
A. Singh, P. R. Nott, Experimental measurements of the
normal stresses in sheared Stokesian suspension, Journal of
Fluid Mechanics 490 (2003) 293320.
A. Singh, A. Nir, R. Semiat. Free surface flow of
concentrated suspension. Int. J. Multiphase Flow, 32 (2006)
775790.
B. W. Roberts and W. L. Olbricht, Flowinduced
particulate separations, AIChE J. 49 (2003) 28422849.
GadalaMaria and A. Acrivos, Shearinduced structure in a
concentrated suspension of solid spheres, Journal of
Rheology 24 (1980) 799.
B. D. Timberlake, J. F. Morris, Particle migration and
freesurface topography in inclined plane flow of a
suspension, J. Fluid Mech. 538 (2005) 309341.
B. M. Fenton, R. R. Carr, G. R. Cokelet, Nonuniform red
cell distribution in 20 to 100 micrometers bifurcations,
Microvasc. Res. 29 (1985) 103126.
B. K. Chapman, Shear induced migration phenomena in
concentrated suspensions," Ph.D. Thesis, University of
Notre Dame (1990).
B. W. Roberts and W. L. Olbricht, The distribution of freely
suspended particles at microfluidic bifurcations, AIChE J.
52 (2006) 199206.
C. Xi and N. C. Shapley, Flows of concentrated
suspensions through an asymmetric bifurcation, J. Rheol.
52 (2008) 625647.
D. T. Leighton Jr. and A. Acrivos, The shearinduced
migration of particles in concentrated suspension, J. Fluid
Mech. 181 (1987) 415439.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
F. W. Rong and R. R. Carr, Dye studies on flow through
branching tubes, Microvasc. Res. 39 (1990) 186202.
G. W. SchmidSchonbein, R. Skalak, S. Usami, S. Chien,
Cell distribution in capillary networks, Microvasc.Res. 19
(1980) 1844.
I. E. Zarraga, Hill, D. A., Leighton, D. T., The
characterization of the total stress of concentrated
suspensions of noncolloidal spheres in Newtonian fluids. J.
Rheol, 44 (2000) 185220.
J. F. Brady, G. Bossis, Stokesian dynamics, Annual Review
of Fluid Mechanics 20 (1988) 111157.
J. F. Brady, J. F. Morris, Microstructure of strongly sheared
suspensions and its impact on rheology and diffusion, J.
Fluid Mech., 348 (1997) 103 139.
J. F. Morris and F. Boulay, Curvilinear flows of noncolloidal
suspensions: the role of normal stresses, J. Rheol. 43 (1999)
12131237.
J. Perkkio and R. Keskinen, Hematocrit reduction in
bifurcations due to plasma skimming, Bull. Math. Biol. 45
(1983) 4150.
K. Zhang and A. Acrivos, Viscous resuspension in fully
developed laminar pipe flows, Int. J. Multiphase flow 20
(1994) 579591.
M. S. Ingber, A. L. Graham, L. A. Monday, Z. Fang, An
improved constitutive model for concentrated suspensions
accounting shearinduced particle migration rate
dependence on particle radius, Int. J. Multiphase Flow 35
(2009) 270276.
M. K. Lyon and L. G. Leal, An experimental study of the
motion of concentrated suspensions in twodimensional
channel flow. Part 1. Monodisperse system, J. Fluid Mech.
363 (1998) 2556.
N. Tetlow, A. L. Graham, M. S. Ingber, S. R. Subia, L. A.
Monday, Particle migration in a Couette apparatus:
experiment and modeling, J. Rheol. 42 (1998) 307327.
P R. Nott and J. F. Brady, Pressuredriven flow of
suspensions: simulation and theory, J. Fluid Mech. 275
(1994) 157199.
R. M. Miller and J. F. Morris, Normal stressdriven
migration and axial development in pressuredriven flow of
concentrated suspensions, J. NonNewt. Fluid Mech.135
(2006) 149165.
R. L. Dedrick, Arterial drug infusion: Pharmacokinetic
problems and pitfalls, J. Natl. Cancer Inst. 80 (1988) 8489.
R. J. Phillips, R. C. Armstrong, R. A. Brown, A. L. Graham,
J. R. Abbott, A constitutive equation for concentrated
suspension that accounts for shearinduced particle
migration, Phys. Fluids A 4 (1992) 3040.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
R. R. Rao, L. A. Monday, T. A. Baer, S. A. Altobelli, T. S.
Stephens, NMR measurements and simulations of particle
migration in nonNewtonian fluids, Chem. Eng. Comm. 89
(2002) 122.
S. R. Subia, M. S. Ingber, L. A. Monday, S. A. Altobelli, A.
L.Graham, Modeling of concentrated suspensions using a
continuum constitutive equation, J. Fluid Mech. 373 (1998)
193219.
S. A. Altobelli, R. C. Givler, E. Fukushima, Velocity and
concentration measurements of suspensions by nuclear
magnetic resonance imaging, J. Rheol. 35 (1991) 721734.
Z. Fang and N. PhanThien, Numerical simulation of
particle migration in concentrated suspensions by a finite
volume method, J. NonNewtonian Fluid Mech. 58 (1995)
6781.
Z. Fang, A. A. Mammoli, J. F. Brady, M. S. Ingber, L. A.
Monday, A. L. Graham, Flowaligned tensor models for
suspension flows, Int. J. Multiphase Flow, 28 (2002) 137 
166.
