Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 16.4.3 - CFD simulation of shear induced particle migration
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00402
 Material Information
Title: 16.4.3 - CFD simulation of shear induced particle migration Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Ahmad, Y.
Kumar, M.G.
Singh, A.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: shear induced migration
CFD
suspension flow
 Notes
Abstract: Flow of concentrated suspensions is encountered in various process industries. It is now well known that shear induced migration is observed in solid-liquid suspensions in which initially well-mixed particles in suspensions when subjected to inhomogeneous shear migrate from region of high shear rate to the region of low shear rate giving a non-uniform 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 convection-diffusion 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 T-junction 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.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00402
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1643-Ahmad-ICMF2010.pdf

Full Text

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 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 solid-liquid suspensions in which initially well-mixed particles in suspensions when subjected to
inhomogeneous shear migrate from region of high shear rate to the region of low shear rate giving a non-uniform 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 convection-diffusion 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 T-junction 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 non-Newtonian fluids whose
theological properties are influenced by a large number of
factors (Zarraga et al. 2000; Singh & Nott 2003).
Gadala-Maria 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 pressure-driven flow in a channel of a
non-Brownian suspension were conducted by Nott and
Brady (1994). They performed Stokesian Dynamic
Simulation (Brady and Bossis, 1988) for non-Brownian
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 surface-roughness-independent 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
non-Newtonian 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
normal-stress-differences 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 30-June 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 non-Newtonian 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 Phan-Thein, 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
non-uniform 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 Y-junction
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 (Schmid-Schonbein 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
shear-induced 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 30-June 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 (ms-1)
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 (ms-1)
Vw Plate speed (ms-1)
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 gradient-induced 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 wide-gap 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 30-June 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 Phan-Thien (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
right-hand-side 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 non-local
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
non-local 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 pressure-velocity 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 two-dimensional
steady state simulation in a channel and the results were
compared with semi-analytical 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 30-June 4, 2010

The function d(,) accounts for the bulk particle
concentration dependence of the shear-induced 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 2-D 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 2-D simulations the relative particle-size (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

semi-analytic 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
semi-analytical 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 non-local formulation is
able to remove the sharp cusp observed for semi-analytical
solution. The decrease in centerline value is due to addition
of non-local 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 semi-analytical 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
semi-analytical 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
semi-analytical solution and with parabolic profile


Suspension flow in asymmetric T-junction

bifurcation channel

Asymmetric T-junction bifurcation channel
represents a branching system in which flow divides into


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 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 T-junction 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 x-z plane


2- q






sl








at the centerline (y=A/2). Figures 4b and 4e are for the side
view (y-z) 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 de-acceleration 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 x-z plane at the centerline (y=A/2). Figures (b) and (e) are
for the side view (y-z) plane at x=1.5 mm for main branch
and figures (c) and (f) shows the side view (y-z) 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 30-June 4, 2010



To view the velocity profile across the channel we have
shown axial velocity map (x-y 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 peak-valley-peak 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 no-slip 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 (x-y 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 30-June 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.35-0.50 for experiments and 0.35-0.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-' )--I-B'


L
'A'
(b)


A'

(c)


Figure 8: Axial particle concentration map (x-y 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

lower-left 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 30-June 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 non-uniform

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 (B-B'), 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 Y-junction microfluidic bifurcation

suspension flow. They have explained that the suspension

near the inner walls of the downstream branches originates

from the particle-rich central region of the upstream branch.

The suspension near the outside walls of the downstream

branches originates from the particle-free 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 30-June 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 30-June 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 .28-01
5. 65c-I
6.3 -01

4.39e-01
4.7OB-ll
a.TUe-
3.14e- 1

2.82e-01
220e-U


g.426-12
6.28e-O2
3.14 8-02
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 T-junction

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 non-uniform 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 blood-vessels in the skin and
web of the frog, J. Physiol. (London) 55 (1921) 412-422.

A. Singh, P. R. Nott, Experimental measurements of the
normal stresses in sheared Stokesian suspension, Journal of
Fluid Mechanics 490 (2003) 293-320.

A. Singh, A. Nir, R. Semiat. Free surface flow of
concentrated suspension. Int. J. Multiphase Flow, 32 (2006)
775-790.

B. W. Roberts and W. L. Olbricht, Flow-induced
particulate separations, AIChE J. 49 (2003) 2842-2849.

Gadala-Maria and A. Acrivos, Shear-induced structure in a
concentrated suspension of solid spheres, Journal of
Rheology 24 (1980) 799.

B. D. Timberlake, J. F. Morris, Particle migration and
free-surface topography in inclined plane flow of a
suspension, J. Fluid Mech. 538 (2005) 309-341.

B. M. Fenton, R. R. Carr, G. R. Cokelet, Nonuniform red
cell distribution in 20 to 100 micrometers bifurcations,
Microvasc. Res. 29 (1985) 103-126.

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) 199-206.

C. Xi and N. C. Shapley, Flows of concentrated
suspensions through an asymmetric bifurcation, J. Rheol.
52 (2008) 625-647.

D. T. Leighton Jr. and A. Acrivos, The shear-induced
migration of particles in concentrated suspension, J. Fluid
Mech. 181 (1987) 415-439.


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

F. W. Rong and R. R. Carr, Dye studies on flow through
branching tubes, Microvasc. Res. 39 (1990) 186-202.

G. W. Schmid-Schonbein, R. Skalak, S. Usami, S. Chien,
Cell distribution in capillary networks, Microvasc.Res. 19
(1980) 18-44.

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) 185-220.

J. F. Brady, G. Bossis, Stokesian dynamics, Annual Review
of Fluid Mechanics 20 (1988) 111-157.

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)
1213-1237.

J. Perkkio and R. Keskinen, Hematocrit reduction in
bifurcations due to plasma skimming, Bull. Math. Biol. 45
(1983) 41-50.

K. Zhang and A. Acrivos, Viscous resuspension in fully
developed laminar pipe flows, Int. J. Multiphase flow 20
(1994) 579-591.

M. S. Ingber, A. L. Graham, L. A. Monday, Z. Fang, An
improved constitutive model for concentrated suspensions
accounting shear-induced particle migration rate
dependence on particle radius, Int. J. Multiphase Flow 35
(2009) 270-276.

M. K. Lyon and L. G. Leal, An experimental study of the
motion of concentrated suspensions in two-dimensional
channel flow. Part 1. Monodisperse system, J. Fluid Mech.
363 (1998) 25-56.

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) 307-327.

P R. Nott and J. F. Brady, Pressure-driven flow of
suspensions: simulation and theory, J. Fluid Mech. 275
(1994) 157-199.

R. M. Miller and J. F. Morris, Normal stress-driven
migration and axial development in pressure-driven flow of
concentrated suspensions, J. Non-Newt. Fluid Mech.135
(2006) 149-165.

R. L. Dedrick, Arterial drug infusion: Pharmacokinetic
problems and pitfalls, J. Natl. Cancer Inst. 80 (1988) 84-89.

R. J. Phillips, R. C. Armstrong, R. A. Brown, A. L. Graham,
J. R. Abbott, A constitutive equation for concentrated
suspension that accounts for shear-induced particle
migration, Phys. Fluids A 4 (1992) 30-40.






7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 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 non-Newtonian fluids, Chem. Eng. Comm. 89
(2002) 1-22.

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)
193-219.

S. A. Altobelli, R. C. Givler, E. Fukushima, Velocity and
concentration measurements of suspensions by nuclear
magnetic resonance imaging, J. Rheol. 35 (1991) 721-734.


Z. Fang and N. Phan-Thien, Numerical simulation of
particle migration in concentrated suspensions by a finite
volume method, J. Non-Newtonian Fluid Mech. 58 (1995)
67-81.

Z. Fang, A. A. Mammoli, J. F. Brady, M. S. Ingber, L. A.
Monday, A. L. Graham, Flow-aligned tensor models for
suspension flows, Int. J. Multiphase Flow, 28 (2002) 137 -
166.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs