Computers & Fluids 39 (2010) 709–721

Contents lists available at ScienceDirect

Computers & Fluids

journal homepage: www.elsevier .com/ locate /compfluid

Numerical simulations of particle migration in a viscoelastic fluidsubjected to shear flow

G. D’Avino a,*, T. Tuccillo a, P.L. Maffettone a, F. Greco b, M.A. Hulsen c

a Dipartimento di Ingegneria Chimica, Università di Napoli Federico II, Napoli I-80125, Italyb Istituto di Ricerche sulla Combustione, IRC-CNR, Naples, Italyc Department of Mechanical Engineering, Eindhoven University of Technology, Eindhoven 5600 MB, The Netherlands

a r t i c l e i n f o a b s t r a c t

Article history:Received 23 July 2009Received in revised form 1 October 2009Accepted 17 November 2009Available online 23 November 2009

Keywords:Particle migrationMicrofluidicViscoelasticityNumerical simulationsArbitrarian Lagrange Eulerian

0045-7930/$ - see front matter � 2009 Elsevier Ltd. Adoi:10.1016/j.compfluid.2009.11.005

* Corresponding author.E-mail addresses: gadavino@unina.it (G. D’Avi

(T. Tuccillo), pierluca.maffettone@unina.it (P.L. Ma(F. Greco), m.a.hulsen@tue.nl (M.A. Hulsen).

Particle migration is a relevant transport mechanism whenever suspensions flow in channels with gapsize comparable to particle dimensions (e.g. microfluidic devices). Several theoretical as well as experi-mental studies have been performed on this topic, showing that the occurring of this phenomenonand the migration direction are related to particle size, flow rate, and the nature of the suspending liquid.

In this work we perform a systematic analysis on the migration of a single particle in a sheared visco-elastic fluid through 2D finite element simulations in a Couette planar geometry. To focus on the effects ofviscoelasticity alone, inertia is neglected. The suspending medium is modeled as a Giesekus fluid.

An ALE particle mover is combined with a DEVSS/SUPG formulation with a log-representation of theconformation tensor giving stable and convergent results up to high flow rates. To optimize the compu-tational effort and reduce the remeshing and projection steps, needed as soon as the mesh becomes toodistorted, a ‘backprojection’ of the flow fields is performed, through which the particle in fact movesalong the cross-streamline direction only, and the mesh distortion is hence drastically reduced.

Our results, in agreement with recent experimental data, show a migration towards the closest walls,regardless of the fluid and geometrical parameters. The phenomenon is enhanced by the fluid elasticity,the shear thinning and strong confinements. The migration velocity trends show the existence of a mastercurve governing the particle dynamics in the whole channel. Three different regimes experienced by theparticle are recognized, related to the particle-wall distance. The existence of a unique migration behav-ior and its qualitative aspects do not change by varying the fluid parameters or the particle size.

� 2009 Elsevier Ltd. All rights reserved.

1. Introduction

Particle cross-streamline migration occurs in many practicalutilizations, i.e. fluid–solid separation where different sized parti-cles need to be separated, suspensions in microfluidic deviceswhere migration can lead to non-homogenous particle distribu-tions, cells in blood provoking clots obstructing the flow, removalof pollutants from a gaseous flow, etc. Recently, new applicationsexploiting the lateral particle motion are arising. Examples are iso-lation of human leukocytes [1] or focusing of red blood cells inmicrochannel flows for bio-sensing applications [2]. Due to itsimportance in applications dealing with solid-particle transportmechanisms, it has been extensively studied over many decades.

From an experimental point of view, the first work focused onparticle migration was performed by Segre and Silberberg [3,4].The authors studied the lift experienced by spheres in a dilute

ll rights reserved.

no), teresa.tuccillo@unina.itffettone), fgreco@irc.na.cnr.it

suspension in a Poiseuille flow at low Reynolds number. Theyfound the existence of an equilibrium height in the channel wherethe particles tend to migrate.

Many analytical theories were also derived in the past, givingexpressions for the lift force experienced by a particle in low Rey-nolds number Newtonian flows and simple geometries. An enlight-ening paper by Saffman [5] gave an expression for the potential liftforce on a sphere in an unbounded shear flow. Such a theory waslater generalized by Asmolov [6] and McLaughlin [7] removingsome constraints on the flow parameters.

The motion of a sphere in the presence of a flat wall in shearflow was studied by using a perturbative approach by Leightonand Acrivos [8], Cherukat and McLaughlin [9] and Krishnan andLeighton [10]. They found the critical conditions for Reynolds num-ber such that the sphere separates from the wall.

All the references cited above refer to a Newtonian suspendingfluid. The non-Newtonian (i.e. viscoelastic) nature of the suspend-ing medium is, however, relevant in many processes (polymermelts with fillers, rubbers, cosmetics, foods, etc.). In this case,few experimental data sets on migration are available.

710 G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721

Old data by Mason and co-workers [11–13] reported on migra-tion of spheres in viscoelastic fluids in Poiseuille and Couette flows.Only recently, a careful and complete analysis by Lormand andPhillips [14] reported on the migration of a sphere in a viscoelasticfluid in a Couette apparatus with small curvature. Their dataclearly show the tendency of the particle to migrate towards thewalls of the flow cell. The authors also investigated the influenceof particle dimension and external flow parameters on themigration.

On the theoretical side, and limiting to simple shear only, whichis the case of interest in our paper, we are aware of only two worksexplicitly addressing the viscoelasticity-induced cross-flow migra-tion, but at vanishing flow rates [15] or at modest flow rates in a 2Dapproximation [16]. At Re ¼ 0, Ho and Leal [15] predict the exis-tence of migration induced by normal stresses whenever there isa lateral variation of the shear rate in the undisturbed flow. Thefirst work showing migration in an inertialess viscoelastic fluid isthe pioneer paper of Huang et al. [16] where Direct Numerical Sim-ulations are performed. They reported a particle lateral motionwith direction depending on the presence/absence of a rate-depen-dent viscosity of the suspending liquid. In their numerical simula-tions, the influence on particle lift of fluid parameters, particle size/channel height ratio were also studied. Huang et al. [16] accountedfor shear thinning with an ad-hoc modification of the Oldroyd-Bmodel. They used an elastic–viscous-split-stress (EVSS) formula-tion to discretize the momentum balance and the solid–fluid cou-pling is treated by an Arbitrarian Lagrangian Eulerian (ALE) methodwith mesh velocities. As in Ho and Leal [15], the main conclusion ofHuang et al. [16] is that migration has to be ascribed to the pres-ence of normal stresses.

The purpose of the present work is to perform a systematicnumerical study on the migration of a particle suspended in a real-istic viscoelastic fluid under simple shear flow. The suspendingmedium is modeled as a Giesekus fluid [17], which is often capableof accurately describing experimental viscoelastic data. The studyis carried out by neglecting fluid and particle inertia. It is antici-pated here that our simulation results qualitative agree with theexperiments by Lormand and Phillips [14], i.e., migration is pre-dicted to occur towards the walls of the shear cell, although our re-sults are limited to 2D. The analysis is performed through 2D directnumerical simulations to take into account the non-linear natureof the problem.

The momentum balance is discretized through the Discrete-Elastic–Viscous-Split-Stress (DEVSS) method that is one of themost robust formulations currently available. The viscoelasticconstitutive equation is stabilized by implementing the Stream-line-Upwind-Petrov–Galerkin (SUPG) technique. Furthermore alog-conformation representation of the conformation tensor isused. Finally, an ALE particle mover [18] is adopted to handle the

Fig. 1. Schematic representation of the problem: a rectangular fluid domain ðXÞwith dimis located at the domain center.

particle motion. The numerical scheme implemented here leadsto stable and convergent results, and allows one to achieve sub-stantial flow rates (as compared with [16]). In order to efficientlymanage the particle motion, a trick is used whereby the particleonly moves along the y-direction (i.e., the migration direction).This is achieved through ‘backprojection’ of the computed fieldsalong the x-direction (the main flow direction) at any time step.In this way, remeshing due to ALE approach is only neededonce–twice per run, always preserving the accuracy of the solution.

The influence of flow rate as well as of particle dimension (ascompared to the gap size) on the migration velocity is also studied.

2. Governing equations

We consider a single, rigid, non-Brownian, inertialess, circularparticle (2D problem) moving in a channel filled by a viscoelasticfluid. The problem is schematized in Fig. 1: a particle with diame-ter Dp ¼ 2Rp, denoted by PðtÞ and boundary @PðtÞ, moves in a rect-angular domain, X, with dimensions L and H along x- and y-axisrespectively and external boundaries denoted by Ci ði ¼ 1; . . . ;4Þ.The Cartesian x and y coordinates are selected with the origin atthe center of the domain. On the upper and lower boundaries,equal and opposite velocities are imposed that, for an unfilled fluid,would generate the shear flow depicted on the right part of thesame figure.

The vector xp ¼ ðxp; ypÞ gives the position of the center of theparticle P. In order to evaluate particle rotation, an angular infor-mation, H ¼ Hk, is also associated with the particle, where k isthe unit vector in the direction normal to the x—y plane. The par-ticle moves according to the imposed flow and its rigid-body mo-tion is completely defined by the translational velocity, denotedby Up ¼ dxp=dt ¼ ðUp;VpÞ and angular velocity, x ¼ dH=dt ¼ xk.

The governing equations for the fluid domain, X� PðtÞ, neglect-ing inertia, can be stated as follows:

r � r ¼ 0 ð1Þr � u ¼ 0 ð2Þr ¼ �pI þ 2gsDþ s ð3Þ

Eqs. (1)–(3) are the equations for the momentum balance, the massbalance (continuity) and the expression for the total stress, respec-tively. In these equations u; r; p; I; D; gs, are the velocity vector,the stress tensor, the pressure, the 2 � 2 unity tensor, the rate-of-deformation tensor and the viscosity of a Newtonian ‘solvent’,respectively. The viscoelastic stress, s, is written as (for the consti-tutive model chosen, see below):

s ¼ gkðc � IÞ ð4Þ

ensions L� H with a single particle PðtÞ is considered. The origin of a Cartesian frame

G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721 711

where c is the ‘conformation tensor’, g is the polymer viscosity, andk is the relaxation time.

We will model the viscoelastic fluid with the Giesekus constitu-tive equation (for c):

k crþc � I þ aðc � IÞ2 ¼ 0 ð5Þ

where a is the so-called mobility parameter that modulates theshear thinning behavior. The symbol (r) denotes the upper-con-vected time derivative, defined as:

cr� @c@tþ u � rc � ðruÞT � c � c � ru ð6Þ

Notice that the zero shear viscosity for the total stress tensor r

is given by g0 ¼ gs þ g.The boundary and initial conditions are given by:

u ¼ Up þ x� ðx� xpÞ on @PðtÞ ð7Þu ¼ _cy; v ¼ 0 on C1 and C3 ð8Þuð�L=2; yÞ ¼ uðL=2; yÞ 8y 2 ½�H=2;H=2� ð9Þrð�L=2; yÞ ¼ rðL=2; yÞ 8y 2 ½�H=2;H=2� ð10Þcjt¼0 ¼ c0 ð11Þ

Eqs. (7) and (8) are the rigid-body motion on the particle boundaryand the shear flow conditions on the upper and lower fluid bound-aries, with _c the applied shear rate. Eqs. (9) and (10) state periodicboundary conditions on C2 and C4. Since inertia is neglected, no ini-tial condition of the velocity and pressure fields is required whereasthe initial conformation tensor condition is necessary (Eq. (11)). Inour simulations, we use a stress-free state, i.e. cjt¼0 ¼ I, as initialcondition over the whole fluid domain. Notice that, by consideringperiodicity along x-direction, no inflow section exists, therefore wedo not need to specify the conformation tensor on any boundary.

The equation governing the rigid-body motion (Eq. (7)) adds(for the 2D case) three additional unknowns, namely, the transla-tional and angular velocities of the particle. So, to obtain particlemotion, it is necessary to consider the balance equations for dragforce and torque acting on the particle boundary. Under theassumptions of absence of particle inertia, and of no ‘external’forces and torques (force- and torque-free particle), such balanceequations are given by:

F ¼Z@PðtÞ

r � nds ¼ 0 ð12Þ

T ¼Z@PðtÞðx� xpÞ � ðr � nÞds ¼ 0 ð13Þ

where F ¼ ðFx; FyÞ and T ¼ Tk are the total force and torque on theparticle boundary, respectively, and n is the outwardly directed unitnormal vector on @P.

The particle position and rotation are updated by integratingthe following kinematic equations:

dxp

dt¼ Up; xpjt¼0 ¼ xp;0 ð14Þ

dHdt¼ x; Hjt¼0 ¼ H0 ð15Þ

The analysis is carried out by making the above equationsdimensionless, using _c as the characteristic time scale and g _c asthe scale for the stress. Then, the Deborah number, De ¼ k _c, willappear in the equations. The Deborah number is simply a nondi-mensional measure of the intensity of the imposed flow rate.

3. Weak form and implementation

3.1. Weak form

The system of Eqs. (1)–(6) with initial and boundary conditions(7)–(11) and the balance Eqs. (12) and (13) form a well-posedproblem in the unknowns: p; u; c; Up; x. At any time step theproblem is solved and the flow fields, and rigid-body unknownsare evaluated. The kinematic Eqs. (14) and (15) are then integratedto update the particle position and rotation.

Such a system is solved by the finite element method. A log-conformation representation for the conformation tensor has beenused [19,20]. The original equation for the conformation tensor c,Eq. (5), is transformed to an equivalent equation for s ¼ logðcÞ:

_s ¼ @s@tþ u � rs ¼ gðruT ; sÞ ð16Þ

An expression for the function g for a Giesekus fluid can be found in[20]. Solving the equation for s instead of the equation for c leads toa substantial improvement of stability for high Deborah numbers.

The finite element method requires the weak formulation of theequations to be solved. Such a weak form for the X� PðtÞ domaincan be stated as follows: For t > 0, find u 2 U; p 2 P; s 2 S,G 2 G; Up 2 R2; x 2 R, k 2 L2ð@PðtÞÞ such that:Z

X2gsDðvÞ : DðuÞdA�

ZXr � v pdAþ

ZX

aðrvÞT : rudA

�Z

XaðrvÞT : GT dA

þZ@PðtÞ½v � ðV þ v� ðx� xpÞÞ� � kds ¼ �

ZX

DðvÞ : sdA

ð17Þ

ZX

qr � udA ¼ 0 ð18Þ

ZX

H : G dA�Z

XH : ðruÞT dA ¼ 0 ð19Þ

ZXðS þ sðu� uÞ � rSÞ :

@s@tþ ðu� uÞ � rs� gðG; sÞ

� �dA ¼ 0 ð20Þ

Z@PðtÞ

l � ½u� ðUp þ x� ðx� xpÞÞ�ds ¼ 0 ð21Þ

s ¼ s0 at t ¼ 0 ð22Þ

for all v 2 U; q 2 P; S 2 S; H 2 G, V 2 R2; v 2 R and l 2 L2ð@PðtÞÞ,where U; P; S and G are suitable functional spaces. In Eq. (20) uis the velocity of the mesh nodes (see below for details).

To improve the numerical stability at high Deborah numbers, aDEVSS-G formulation is implemented [21,22] for the momentumbalance, Eq. (17), and the SUPG technique [23] is used for the con-stitutive relation, Eq. (20).

The s parameter in Eq. (20) is given by s ¼ bh=2Ue, where b is adimensionless constant, h is a typical size of the element and Ue is acharacteristic velocity for the element. In our simulations, we havechosen b ¼ 1 and for Ue we take the average of the magnitude ofthe velocities in all integration points. In addition, a in Eq. (17) ischosen as a ¼ g. We take the initial value s0 ¼ 0, correspondingto zero initial stress.

The rigid-body motion on the particle boundary is imposedthrough Lagrange multipliers k (see [24]). As a consequence, theparticle kinematic quantities are considered as additional un-knowns and are recovered by solving the full system of equations.

712 G. D’Avino et al. / Computers &

The particle motion is taken into account by using an ALE for-mulation [18], whereby at each time step the mesh nodes aremoved according to a mesh velocity u obtained by solving an extraequation. Notice that in general such node velocities are differentfrom the fluid velocities at the nodes. This can be done until themesh becomes too distorted. After that a new mesh is generatedand the solution on the old mesh is projected onto the new one.

The node mesh movement is taken into account in the discret-ized form of the governing equations by subtracting the meshvelocity to the fluid velocity in the convective terms. In this regard,note that, due to the fluid inertialess assumption, the convectiveterm only appears in the viscoelastic constitutive equation, so insuch a term (and in the test function) the velocity is actuallyu� u. For the latter velocity, we need an extra equation. FollowingHu et al. [18], in order to guarantee a smooth mesh motion, we usethe Laplace equation:

r � ð�ruÞ ¼ 0 ð23Þ

with boundary conditions:

u ¼ U þ x� ðx� xpÞ on @PðtÞ ð24Þu ¼ 0 on Ci; ði ¼ 1; . . . ;4Þ ð25Þ

In Eq. (23), the parameter � is taken equal to the inverse of thelocal element area in order to let the largest elements adsorb themost part of deformation. Eqs. (24) and (25) assure that the parti-cle boundary nodes move following the particle motion, whereasno movement is prescribed for domain boundary nodes. The weakform for Eq. (23) can be derived in a standard way.

The DEVSS-G/SUPG formulation with a log-representation ofthe conformation tensor and an ALE particle mover provides a veryefficient method in order to handle the system under investigation,giving satisfactory results for all the set of parameters consideredin this work.

3.2. Implementation

For the discretization of the weak form, we use triangular ele-ments with continuous quadratic interpolation ðP2Þ for the velocityu, linear continuous interpolation ðP1Þ for the pressure p, linearcontinuous interpolation ðP1Þ for velocity gradient G and linearcontinuous interpolation ðP1Þ for the log-conformation tensor s.

Regarding the time discretization, we decouple the momentumbalance and continuity equations from the constitutive one. Asshown in the weak form, the mesh velocity, u, appears in the con-stitutive equation only. This allows us to solve separately the La-place equation and the momentum balance.

Initially, the viscoelastic stress is set to zero in the whole do-main. Since we neglect inertia, the initial condition for the velocityis not necessary. So, we can solve Eqs. (17)–(19) and the constraintEq. (21) in order to get the distribution of the fluid velocity and therigid-body motion of the particles, at the initial time step. Then, atevery time step, we follow the following procedure:

Step 1: The particle position and rotation are updated. The newconfiguration is obtained by integrating the kinematicEq. (14). The explicit second-order Adams–Bashforthmethod is used:

xnþ1p ¼ xn

p þ Mt32

Unp �

12

Un�1p

� �ð26Þ

The first time step of the simulation is performed with an expli-cit Euler method:

xnþ1p ¼ xn

p þ Mt Unp ð27Þ

Step 2: The mesh nodes, xm, are updated according to:

Fluids 39 (2010) 709–721

xnþ1m ¼ xn

m þ Mt un ð28Þ

for the first time step and:

xnþ1m ¼ xn

m þ Mt32

un � 12

un�1� �

ð29Þ

when Eq. (26) is used to update particle position.

Step 3: The log-conformation tensor at the next time step, snþ1, isevaluated by integrating the constitutive Eq. (20). A sec-ond-order Crank–Nicolson/Adams–Bashforth scheme isused:

snþ1

Mtþ 1

2½ð2ðun � unÞ � ðun�1 � un�1Þ� � rsnþ1

¼ sn

Mt� 1

2ðun � unÞ � rsn þ 3

2gðGn; snÞ � 1

2gðGn�1; sn�1Þ ð30Þ

where g is the term which appears in the evolution equation of s(Eq. (16)). Notice that the velocity u� u in the convective terms isaccording to the ALE scheme. In the time integration scheme above,2un � un�1 and 2un � un�1 are used as velocity predictions for unþ1

and unþ1, respectively.

Step 4: The remaining unknowns ðu; p;G;Up;xÞnþ1 as well as theLagrange multipliers ðkÞ can be found by solving Eqs.(17)–(19) and (21) using the particle configuration andthe viscoelastic stress evaluated in the previous two steps:

Z

X2gsDðvÞ : Dðunþ1ÞdA�

ZXr � v pnþ1 dAþ

ZX

aðrvÞT : runþ1 dA

�Z

XaðrvÞT : ðGnþ1ÞT dA

þ v � V þ v� x� xnþ1p

� �� �; knþ1

D E@P

¼ �Z

XDðvÞ : snþ1 dA ð31Þ

ZX

qr � unþ1 dA ¼ 0 ð32Þ

ZX

H : Gnþ1 dA�Z

XH : ðrunþ1ÞT dA ¼ 0 ð33Þ

l;unþ1 � Unþ1p þ xnþ1 � x� xnþ1

p

� �� �D E@P¼ 0 ð34Þ

Step 5: Finally, the Laplace equation is solved:

r � ð�runþ1Þ ¼ 0 ð35Þ

with boundary conditions:

u ¼ Unþ1 þ xnþ1 � x� xnþ1p

� �on @PðtÞ ð36Þ

u ¼ 0 on Ci; ði ¼ 1; . . . ;4Þ ð37Þ

and the mesh velocities are obtained.In Step 4, a sparse linear symmetric system needs to be solved

whereas Eq. (30) leads to an unsymmetric system. In both cases weuse the parallel direct solver PARDISO [25].

3.3. Node mesh updating

In principle the ALE method allows to use the starting mesh un-til it becomes too distorted. When this occurs, a new mesh cover-

G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721 713

ing the same domain of the old one is generated and the solution isprojected from the old on the new mesh. The projection step isgenerally done by interpolation: the element containing a newnode is searched through the old mesh and, by using the nodal ele-ment values of the last solution, the fields are interpolated to givenew nodal values. Such a procedure is time-consuming and leadsto interpolation errors. Therefore it should be reduced as muchas possible during the simulation.

Due to the particular geometrical system investigated in thiswork (a single circular particle in a channel), we may exploit thefact that the mesh velocities, u, can be chosen in an arbitraryway, and the remeshing step can be reduced to once–twice persimulation, or even be avoided.

Indeed, as will be shown in Section 5, the particle moves alongthe x-direction, while translating along the y-axis very slowly (par-ticle cross-streamline migration). This is equivalent to fix the par-ticle position along the x-axis and move the mesh nodes with an x-velocity given by the particle x-velocity, Up. However, due to theperiodicity along the x-axis, the new configuration (particle andmesh translated along the x-direction by a quantity related tothe particle x-velocity) is equivalent to the one of the previous timestep. It should be only recognized that the x-component of themesh velocity is ux ¼ Up. Of course the particle can still move alongthe y-axis.

Therefore, such a procedure is implemented by solving in theStep 5 only the y-component of Eq. (23) and setting u ¼ ðUp; uy)in Eq. (30). In the Step 2 the mesh nodes are moved along the y-direction only and the mesh distortion is limited along the y-axis.Our simulations show that at most once–twice remeshing and pro-jection steps are needed.

Finally, regarding the mesh quality, following Hu et al. [18], weuse two monitoring parameters:

f1 ¼ max16e6Nel

ðf e1 Þ and f 2 ¼ max

16e6Nel

f e2

� �ð38Þ

where Nel is the number of elements and:

f e1 ¼ log Ve=Ve

0

� ��� �� and f e2 ¼ log Se=Se

0

� ��� �� ð39Þ

with Ve and Ve0 the volume of the element e and its value in the

undeformed configuration, respectively; Se and Se0 are the aspect ra-

tios of the element e and its value in the undeformed mesh, respec-tively. The aspect ratio is defined as:

Se ¼ ðleÞ3

Ve ð40Þ

with le the maximum length of the sides of the element e.

Fig. 2. Typical mesh used in the simulations. The me

As soon as one of f1 or f2 is greater than 1.39 (i.e. element vol-ume or aspect ratio is larger than four times or smaller than 1/4of its original value) the mesh is considered too distorted andneeds to be regenerated.

4. Convergence and code validation

4.1. Convergence in space and time

In this Section, the results about the convergence in spaceand time of our method are presented. In what follows, we re-port a limited set of example cases in order to show the sensi-tivity of the method implemented to mesh resolution and timestep size. However it must be remarked that the convergenceis checked for any flow and geometrical parameter presentedin this work. In this Section as well as throughout the paper,non-dimensional quantities only are reported ( _c is used as char-acteristic time scale, and the gap size H as the characteristiclength).

In Fig. 2 a typical mesh used in the simulations is shown. Thechannel length (L) is chosen sufficiently large to avoid that the par-ticle can ‘feel’ itself through the periodic boundaries. Notice therefinement around the particle and between the particle and theclosest wall, where largest gradients are expected. The spatial con-vergence is checked by refining the mesh. In Table 1, some param-eters of the mesh used are reported. The mesh in Fig. 2 is the onelabeled as ‘M1’ (the coarsest one). The meshes in the table refer to aspecific particle size and starting position. However, the number ofelements is kept almost the same by changing the geometricalparameters.

In Fig. 3 the angular velocity, x, and the particle y-velocity, Vp,of a particle starting from yp;0 ¼ 0:3 are reported for the meshes la-beled as ‘M1’, ‘M2’, ‘M3’ and ‘M4’. The step size chosen is 0.005 andthe other parameters are: De ¼ 1:0; a ¼ 0:2, gs=g ¼ 0:1; Dp=

H ¼ 0:1.For any mesh considered here similar trends are found: the

migration velocity increases during the start-up and becomes al-most constant as soon as the stress develop. The convergence isachieved for the mesh ‘M2’, although the comparison with ‘M1’shows very small deviations (ffi1–2% for Vp and ffi0.5–1% for x).A similar situation is found for De ¼ 2:0 (not reported). However,we remark that the convergence severely depends on the fluidand geometrical parameters, and on the flow rate. In general, ourconvergence tests show that a finer mesh is needed as Deborah num-ber increases (for De ¼ 4:0 we need the mesh ‘M4’) and as a is re-

sh reported is the one labeled as ‘M1’ in Table 1.

Table 1Mesh parameters (for Dp=H ¼ 0:1 and yp;0 ¼ 0:3).

Mesh label #el. on the particle boundary #el. in the mesh

M1 40 4176M2 50 6096M3 60 8012M4 80 12,978

714 G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721

duced (with a! 0 the Giesekus model tends to the Maxwell model,which gives serious convergence problems, as it is well known). Onthe other hand, for different particle sizes, leading of course to a dif-ferent spatial distribution of triangles, about the same number ofelements are sufficient to achieve the convergence.

Finally, the influence of time step size, Dt, is investigated. Weconsider the mesh ‘M3’ and De ¼ 1:0. In Fig. 4, Vp and x by varyingDt are considered. Even for the largest Dt considered, the curvessuperimpose. However, for Dt ¼ 0:02 (dotted line), slight devia-tions are found in the start-up phase (see the inset on the rightof Fig. 4) that are more and more pronounced as the particle startsclose to the wall (such an effect has been reported in [26,27] too).On the other hand, by changing fluid and geometrical parameters,we do not find any different behavior with respect to the oneshown in Fig. 4. For these reasons, in all simulations in this work,Dt is chosen 0.015, giving a good compromise between simulationspeed and accuracy.

Fig. 3. Particle y-velocity ðVpÞ (left) and angular velocity ðxÞ (rig

Fig. 4. Particle y-velocity ðVpÞ (left) and angular velocity ðxÞ (r

4.2. Code validation

The method proposed (DEVSS/SUPG + ALE) is validated througha comparison with the fictitious domain method (FDM) [28]. In thelatter, a regular mesh is used covering the fluid as well as the soliddomain. The particle rigid-body motion is imposed inside the par-ticle (or on its boundary if inertia is neglected [29,24]) through La-grange multipliers. The details of the fictitious domain methodused in this work (weak form and implementation) are reportedin [24] with the only difference that here we used a weak imple-mentation of rigid-body constraints instead of collocation [30].We just give a brief description of the mesh used: since the particlemoves along the x-direction because of the external flow field andmigrates slowly (approaching the closest wall), the regular meshused is refined in a narrow channel where the particle is expectedto move. As soon as the particle, due to the migration along the y-direction, is going out the refined channel, a new mesh is generatedby translating the channel (that, of course, partially covers the re-fined region of the previous mesh). The solution is then projectedfrom the old to the new mesh. In such a way, we can concentrateelements where the particle is moving. Considering that the pro-jection is done at most 3–4 times in a simulation, the accuracy isnot compromised at all.

In Fig. 5 (left) the comparison between the two methods is doneby considering the particle y-position as function of time. A goodquantitative agreement is found for any particle starting position.

ht) as functions of time for the meshes reported in Table 1.

ight) as functions of time for different time step size ðDtÞ.

Fig. 5. y-Position of particle center ðypÞ (left) and particle y-velocity ðVpÞ (right) for ALE method (solid lines) and fictitious domain method (dashed line).

G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721 715

The FDM predictions are always slightly higher then the ALE ones.As indeed is evident in Fig. 5 where Vp versus time are reported(right), whereas a monotonic curve (initially decreasing and thenincreasing) is found for Vp in the ALE method, oscillations are evi-dent when the FDM is considered. Such oscillations are due to therelative motion of the particle boundary grid over the fixed grid. Ingeneral, the error is small but in the problem investigated Vp issmall as well and this results in strong oscillations around the ex-act value. Finally, considering that the oscillations are on averageabove the curve for ALE, the trajectories reported on the left ofFig. 5 can be justified. To decrease the velocity perturbations a fur-ther mesh refinement is needed and the dashed curves are ex-pected to tend to the solid one.

To close the comparison between ALE and FDM, it must be ta-ken into account the computational aspect as well. Our FDM sim-ulations requires a calculation time about 10 times higher thanthe ALE method. This is quite easy to understand if one considersthat the refinement around the particle is more effective when aboundary fitted mesh is implemented and an irregular, coarsemesh can be done far from the particle. On the other hand, theFDM mesh should satisfy some regularity far from the particle aswell, limiting its coarseness. In conclusion, the problem of a single

Fig. 6. y-Position of particle center ðypÞ as a function of time ðtÞ for different starting posshaded area is the channel region unaccessible to the particle due to its finite size.

particle in a rectangular box is more effectively tackled by using aboundary fitted method. Instead, for a many particle system, wherefrequent remeshing and projection steps are required, along withthe difficulties related to the mesh generation itself, it would bepreferable to use a fictitious domain method, maybe together witha proper algorithm deforming and refining the regular mesh alongthe particle boundaries.

5. Simulation results

In this Section our simulation results are presented and dis-cussed. As stated in Section 4, the time step size Dt ¼ 0:015 pro-vides time convergent solutions for any case investigated in thiswork. On the contrary, the spatial convergence is therefore checkedfor all the situations investigated.

In Fig. 6, the y-position of the particle center yp as function of timeis reported for different starting positions yp;0. The parameters cho-sen are: De = 1.0, a ¼ 0:2; gs=g ¼ 0:1 and Dp=H ¼ 0:1. Due to thesymmetry with respect to the centerline (at y ¼ 0), the trajectoriesare symmetric as well, and only one-half domain ðy P 0Þ is consid-ered. The shaded region in the upper part delimits the channel zone

itions ðyp;0Þ. The other parameters are: De = 1.0, a = 0.2, gs=g ¼ 0:1; Dp=H ¼ 0:1. The

Fig. 7. Particle migration velocities ðVpÞ as a function of time ðtÞ for different starting positions ðyp;0Þ. The other parameters are: De = 1.0, a = 0.2, gs=g ¼ 0:1; Dp=H ¼ 0:1.

716 G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721

where the particle cannot access because of its finite size: when theparticle center is at y ¼ 0:45, the particle touches the wall.

The trends in Fig. 6 show that, from any starting position, theparticle moves towards the closest wall (the upper one). Of course,due to the symmetry mentioned above, a particle initially locatedwith the center exactly at the channel centerline will stay thereand only rotate. Notice also that, close to the upper wallðyp > 0:3Þ, the trajectories are steeper than those in the channelcore, i.e., the particle migration is faster near to the walls. Whatis found here is completely different from the case of shearing inan inertialess Newtonian suspending fluid where no migration oc-curs, independently of the initial position. Finally, note that thecurves do not touch the shaded region because our numericalmethod breaks down as soon as the particle is about 1=5Dp fromthe wall. The analysis of the situation of a particle touching thewall is beyond the scope of this work and requires different tech-niques of analysis. Although a direct quantitative comparison with

Fig. 8. Particle angular velocities ðxÞ as a function of time (t) for different starting po

Huang et al. [16] results cannot be carried out due to the differentconstitutive model used, our simulations qualitative agree withthose of their work (of course the case of an inertialess and shearthinning fluid is chosen).

Particle y-velocities (migration velocities), Vp, are reported inFig. 7 and the angular velocities, x, are plotted in Fig. 8. In both fig-ures, after an initial transient due to the development of viscoelas-tic stresses, the kinematic quantities change slowly in time forparticle near to the channel center, and at a slightly faster pacefor particles nearer to the walls. In other words, in the system un-der investigation, no steady state is ever achieved and the particlemotion is always transient in nature. Exceptions are a particle lo-cated on the channel centerline and a particle that reaches thewall. The outwardly directed migration can then be looked at asan instability of the dynamics at the centerplane: when the particlecenter of mass is midway between the walls, no transverse motionoccurs, but any disturbance in the vertical position triggers the

sitions ðyp;0Þ. The other parameters are: De = 1.0, a = 0.2, gs=g ¼ 0:1; Dp=H ¼ 0:1.

Fig. 9. Particle migration velocities ðVpÞ as a function of particle y-positions ðypÞ for different starting positions ðyp;0Þ. The other parameters are: De = 1.0, a = 0.2,gs=g ¼ 0:1; Dp=H ¼ 0:1. The same symbology of Fig. 6 is used. The grey line is the function reported in Eq. (41). The shaded area is the channel region unaccessible to theparticle due to its finite size.

G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721 717

migration. A very similar instability was very recently experimen-tally observed by Sullivan et al. [31] in pressure-driven flow of abubble suspended in a viscoelastic liquid in microchannels.

The dynamics of the particle so far described is determined bythe three state variables position, migration velocity, and angularvelocity, reported as a function of time in Figs. 6–8 respectively.It is instructive to analyze such dynamics directly in the ‘phasespace’ of the particle, i.e., the space formed by the state variablesonly. This is done in Figs. 9 and 10 where Vp and x, respectively,are reported versus the vertical position. For the sake of clarity,the fast initial dynamics is not included in Figs. 9 and 10. The sym-bols in Figs. 9 and 10 refer to different starting positions, with thesame symbology of Fig. 6. Notice that symmetry with respect tomidplane dictates the parity (oddity) of the function xðyÞ ðVpðyÞÞ,and for this reason only the upper half channel is shown in Figs.9 and 10. Concerning the migration velocity shown in Fig. 9, wemay note that: (i) all the data collapse onto a single curve, regard-less of the starting position, (ii) in the channel core the migration

Fig. 10. Particle angular velocities ðxÞ as a function of particle y-positions ðypÞ atdifferent starting positions ðyp;0Þ. The other parameters are: De = 1.0, a = 0.2,gs=g ¼ 0:1; Dp=H ¼ 0:1. The same symbology of Fig. 6 is used. The grey line is thefunction reported in Eq. (42). The shaded area is the channel region unaccessible tothe particle due to its finite size.

velocity scales linearly with position (see later), (iii) by approach-ing the wall a faster migration is observed and (iv) the curve passesthrough a maximum to abruptly decrease very near to the wall. Inshort, three different dynamics can be distinguished, pertaining tothree different zones of the channel.

Concerning the angular velocity x in Fig. 10, again data collapseon a single curve (when disregarding the initial transients), and dif-ferent dynamics are found in different channel regions. In particu-lar, notice the steep decrease of the angular velocity whenapproaching the wall.

The data reported in Figs. 9 and 10 are well described by the fol-lowing functional forms:

Vp ¼ m1yp þm5y5p ð41Þ

x ¼ a0 þ a2y2p þ a6y6 ð42Þ

Some remarks: (i) the constant term in Eq. (41) is missing be-cause of the symmetry around the channel centerlineðyp ¼ 0) Vp ¼ 0Þ, (ii) the linear term in Eq. (41) governs the dy-namic in the core regime and (iii) the constant term in Eq. (42) isthe angular velocity of a particle at the center plane under confinedshear flow (compare with [27]).

It should be pointed out that the functions (41) and (42) do notcome from any theory or physical consideration but are just thesimplest polynomial fits of the simulation data. Indeed, the fittingpolynomials will be different with different choices of Dp=H; Deand a (see below), not only in terms of coefficients but even inpower laws. Independently from those geometrical and fluidparameters, however, the functional forms in the core regionðy=H� 1Þ are found to be universal, i.e., it is always m0 = 0,m1 – 0, a0 – 0; a2 – 0. Notice that for fitting the migration veloc-ity we excluded the last four points because Eq. (41) cannot beused to predict the decreasing trend very close to the wall.

The functional form, Eq. (41), can be integrated to recover theparticle trajectory over the channel:

dyp

dt¼ Vp ¼ m1 yp þm5 y5

p ) yp

¼yp;0 m1=4

1 expðm1tÞ½m1 þ y4

p;0 m5ð1� expð4m1tÞÞ�1=4 ð43Þ

Fig. 11. Master curve given by Eq. (43) (grey line) and the simulated trajectories reported in Fig. 6 (symbols). The parameters are: De = 1.0, a = 0.2, gs=g ¼ 0:1; Dp=H ¼ 0:1. Arepresentation of the particle is reported as well. The shaded area is the channel region unaccessible to the particle due to its finite size. In the figure is also reported aschematic representation of the particle, for Dp=H ¼ 0:1, to show its dimension compared to the channel size.

Fig. 12. Particle migration velocities ðVpÞ as a function of particle y-positions ðypÞfor different Dp=H values. The other parameters are: De = 1.0, a = 0.2, gs=g ¼ 0:1.The grey lines are the functions reported in Eq. (41) with the parameters calculatedby fitting the simulation data. The shaded area is the channel region unaccessible tothe particle due to its finite size (depending on the particle size).

718 G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721

with initial condition t ¼ 0) yp ¼ yp;0. The meaning of m1 coeffi-cient is apparent from Eq. (43): indeed, it represents the growth rateof the instability from the center line [31]. Eq. (43) is the mastercurve governing the particle trajectories in the whole channel (ex-cluded the ‘terminal’ part). In other words, the behavior in the coreportion of the channel as well as in the not-too-close-to-wallportion is completely defined by Eq. (43), regardless of the startingposition. Notice that the core regime is governed by an exponential-like trajectory, as follows by setting m5 ¼ 0 in Eq. (43).

In Fig. 11, Eq. (43) is reported with a grey line. The symbols referto the curves calculated through numerical simulations reported inFig. 6. The overlap for any starting position yp;0 confirms the valid-ity of Eq. (43) in the whole channel (except the terminal region): inconclusion, each individual trajectory, just translated in time, ispart of an unique behavior.

It seems to appropriate to emphasize once more that such arelationship, Eq. (43), is not valid in the initial transient behavior,due to the ‘rapid’ development of fluid stresses. However, as evi-dent in Fig. 7 or 8, 3–4 relaxation times are sufficient to extinguishthe initial start-up and entering the ‘slow’ dynamics.

The just discussed master curves depend on Dp=H ratio and fluidparameters (De and a). In the following, the effects of geometricaland fluid parameters are separately investigated.

5.1. Effect of particle size

The influence of the particle size is investigated by varying Dp=Hratio (dimensionless diameter) at De ¼ 1 and a ¼ 0:2. In Fig. 12 theparticle y-velocity Vp is reported as a function of particle y-positionyp for four different Dp=H values. The full circles refer toDp=H ¼ 0:10 and are the same data as in Fig. 9. The three differentregimes identified in the previous section, the core, the intermedi-ate and the close-wall regimes, exist even by varying the confine-ment. We notice that a stronger confinement (higher Dp=H) leadsto a faster migration in the whole channel, in particular, the insta-bility growth rate monotonically increases by increasing thedimensionless diameter, following a quadratic law. For whateverdimensionless diameter, however, simulation data are still verywell fitted by the functional form in Eq. (41). The fitting curves(with the two parameters m1 and m5 evaluated by a least square

method) are reported as grey lines in the figure. It is apparent thatthe curves well describe the particle dynamic except very near tothe wall.

The validity of Eq. (41) for any Dp=H ratio investigated impliesthat the particle trajectories are in the form of Eq. (43). Such curvesare reported in Fig. 13 (for yp;0 ¼ 0:05) where the faster migrationat stronger confinement is evident.

Coming back to the three dynamical regimes mentioned above,Fig. 14 summarizes their sizes, ðHÞ, as the dimensionless diameteris varied. The separatrix between the core and intermediate re-gimes is identified by those data in Fig. 12 showing a deviationof 5% from the linear trend. The separatrix between the intermedi-ate and the close-wall regimes is identified by 5% deviations of datafrom the grey curves in the same figure. Notice that the core re-gime covers the most part of the channel, especially as the confine-ment is weak, whereas the close-wall regime becomes larger andlarger as the confinement is stronger.

Fig. 13. y-Position of particle center ðypÞ as a function of time ðtÞ for different Dp=Hratio. The other parameters are: De = 1.0, a = 0.2, gs=g ¼ 0:1. The shaded area is thechannel region unaccessible to the particle due to its finite size (depending on theparticle size).

Fig. 14. Dimensionless channel size ðH=HÞ versus Dp=H. The curves delimit thecore, the intermediate and the close-wall regime. The shaded area is the channelregion unaccessible to the particle due to its finite size (depending on the particlesize). The other parameters are: De = 1.0, a = 0.2, gs=g ¼ 0:1.

Fig. 15. Particle migration velocities ðVpÞ as a function of starting positions ðyp;0Þ fordifferent Deborah numbers ðDeÞ. The other parameters are: a = 0.2, gs=g ¼ 0:1;Dp=H ¼ 0:1.

Fig. 16. y-Position of particle center ðypÞ as a function of time ðtÞ for differentDeborah numbers ðDeÞ. The other parameters are: a = 0.2, gs=g ¼ 0:1; Dp=H ¼ 0:1.

G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721 719

5.2. Effect of fluid rheology

The effect of Deborah number is now investigated, at a fixedDp=H ¼ 0:1. In Fig. 15 the migration velocities as a function ofthe particle y-position are reported for different De values. Theblack circles are the same data as in Fig. 9. It is immediately evi-dent that the overall migration behavior is still found, i.e. we dis-tinguish the linear core regime, a faster intermediate regime anda close-wall regime. Furthermore, the transitions among thosethree dynamics seems to be independent from Deborah number(for example look at the maxima of the curves).

An important effect of the De-value is the non-monotonicbehavior of the instability growth rate visible in Fig. 15. Indeed,the slope of the migration velocity curve in the core regime first in-creases with De. The curve for De ¼ 4:0 (black squares), however,lies in between those for De ¼ 1:0 (black circles) and De ¼ 2:0(white squares). As a consequence, beyond some critical De-num-ber (larger than unity), a particle moving in the central channel re-gion migrates increasingly slower as De increases.

In the intermediate regime, on the other hand, the migrationvelocity increases monotonically with De. This latter behavior im-plies that the functional form in Eq. (41) cannot be used to describethe data in Fig. 15. Indeed, although the linear term is always pres-ent, the exponents of the nonlinear terms vary with De.

The just described dependence on De of the migration velocityis perhaps more evident by looking directly at particle trajectories,as shown in Fig. 16.

The effect of the Deborah number on the particle angular veloc-ity is reported in Fig. 17 where x as a function of the particle posi-tion is reported for different De-values (the Newtonian case,De ¼ 0, is reported as well). A decreasing, monotonic trend is evi-dent as De increases. The decreasing of x with De is somehow sim-ilar to that observed by D’Avino et al. [27] for the case of a spherewith its center at the midplane of the shear cell.

Another interesting aspect is the effect of the shear thinning onthe particle motion. The Giesekus model, Eq. (5), predicts the shearthinning phenomenon in shearing flows modulated by the param-eter a: The thinning is more and more pronounced as a increasesand, for a ¼ 0, the upper convected Maxwell model is recovered(no shear thinning). The influence of thinning on particle migrationis then investigated by varying a. In Fig. 18 the particle trajectories

Fig. 17. Particle angular velocities ðxÞ as a function of particle y-position ðypÞ atdifferent Deborah numbers ðDeÞ. The other parameters are: a = 0.2,gs=g ¼ 0:1; Dp=H ¼ 0:1. The hexagons refer to the Newtonian case.

720 G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721

for three different values of a are reported at fixed De ¼ 1. Byincreasing a, the particle moves faster towards the wall, i.e. themigration velocity is higher. Thus, a particle suspended in a Max-well fluid is expected to experience the slowest migration rate.

Finally, we checked the validity of the behaviors just reportedfor a different Dp=H. By choosing Dp=H ¼ 0:2, we found the sameeffect of De and a on the particle migration.

5.3. Particle translational x-velocity

We finally consider the last state variable in a 2D one-particlesystem, namely, the particle x-velocity, Up. Indeed, due to theexternal imposed flow, the particle translates along the x-directionfor both Newtonian and viscoelastic suspending fluid, with a veloc-ity that is in general different from the one ðux;unfilled ¼ _cyÞ experi-enced by the unfilled fluid at the same vertical position as theparticle center.

In Fig. 19 the relative particle x-velocity, Up � ux;unfilled, usuallytermed the slip velocity, is reported as function of particle positionfor Dp=H ¼ 0:1 for both the Newtonian and a Giesekus fluid withDe ¼ 2:0; a ¼ 0:2. In the whole channel, and for both fluids, theparticle moves faster than the unfilled fluid, i.e., the particle ‘leads’

Fig. 18. y-Position of particle center ðypÞ as a function of time ðtÞ for differentGiesekus constitutive parameter ðaÞ. The other parameters are: De = 1.0,gs=g ¼ 0:1; Dp=H ¼ 0:1.

the fluid. For the Newtonian case, such a result was analyticallydemonstrated long ago [15].

The dashed curve in Fig. 19 is the Newtonian prediction. Due tothe absence of the migration phenomenon, a particle in a Newtonianfluid keeps its horizontal velocity during the motion. In other words,the dashed curve only gives the locus of the translational velocitiesfor different vertical positions. On the other hand, a particle in a vis-coelastic fluid migrates while horizontally translating, hence its hor-izontal velocity changes in time: the solid line in Fig. 19 is travelledin time. As a matter of fact, an initial transient behavior character-ized by a fast change in Up is found at start-up (the almost verticallines departing from the Newtonian curve). The start-up is clearlyvisible in the inset where Up versus time is reported for two differentstarting positions. The dash-dotted line is the Newtonian predictionthat is of course constant in time. In the non-Newtonian case, afterthe extinction of the initial fast transients, the horizontal velocityreaches a master curve from whatever initial starting position. Evenwith this state variable, therefore, a slow dynamics is attained, sim-ilarly to what observed for the migration velocity. The viscoelasticcurve in Fig. 19 is higher than the corresponding Newtonian one ina large central part of the channel, whereas the trend is inverted asthe particle is close to the wall.

6. Conclusions

The migration of a 2D single, circular particle in a sheared visco-elastic liquid is studied through direct numerical simulations. Inorder to pick the influence of the viscoelasticity of the suspendingfluid, particle and fluid inertia are neglected.

An ALE particle mover is combined with a DEVSS-G/SUPG for-mulation together with a log-representation of the conformationtensor, in order to easily manage the particle motion and guaranteenumerical convergence at finite Deborah numbers. Finally, the ri-gid-body motion on the particle boundary is imposed through La-grange multipliers.

Particle trajectories as well as translational and angular veloci-ties are reported for different fluid parameters and particle dimen-sions. Our simulations show that the particle migrates towards theclosest wall regardless of its starting position, fluid and geometricalparameters, i.e. no stable equilibrium position in the channel isfound. Indeed, the cross-flow migration scenario is that of an insta-bility of the centerplane dynamics. Migration towards the walls isindeed observed in recent, accurate experimental data by Lormandand Phillips [14].

The main finding of this work is that, for a given set of geomet-rical and fluid parameters, the migration velocity behavior in thewhole channel can be described by a single curve. As a conse-quence, regardless of the starting position (and after extinction ofthe initial transients due to the development of stresses), the par-ticle trajectories collapse on a master curve.

By analyzing the migration velocity curves, three regions in thechannel are readily identified, characterized by different particledynamics. When the particle is close to the channel centerline,the migration velocity is found to always be linear with the particlevertical position. The slope of this linear portion of the migrationvelocity curve is the growth rate of the instability from the center-plane. Approaching the wall, a faster migration sets in, whereas,quite close to the wall, the migration velocity abruptly decreases.

The existence of these three regimes is confirmed by varyingthe particle dimensions and the fluid rheology. As expected, astronger confinement leads to a reduction of the linear regionand to a faster migration towards the wall. The instability growthrate is found to depend non-monotonically on Deborah number.

The different dynamics depending on the particle positionsthrough the channel and on the particle size can be exploited in

Fig. 19. Particle slip velocity ðUp � ux;unfilledÞ as a function of position ðypÞ for different starting positions (solid lines). The other parameters are: De = 2.0, a = 0.2,gs=g ¼ 0:1; Dp=H ¼ 0:1. The dashed line refers to a Newtonian suspending fluid. In the insets Up � ux;unfilled is reported as function of time for yp;0 ¼ 0:2 and yp;0 ¼ 0:4.

G. D’Avino et al. / Computers & Fluids 39 (2010) 709–721 721

applications where the migration plays a role. For an example, influid–solid separation devices, the inlet section could be conve-niently chosen such that a small particle experiences the linearslow dynamics whereas a larger particle is in the fast regime, opti-mizing the separation process.

Regarding the particle angular velocity, a decreasing trend isfound as the particle approaches the wall and as the Deborah num-ber increases. Finally, the particle horizontal velocity results show,for both Newtonian and viscoelastic fluid, an higher velocity thanthat of the unfilled liquid at the same vertical position, i.e., the par-ticle always ‘leads’ the fluid.

References

[1] Skora J, Kordecka A, Wlosiak P, Madeja Z, Korohoda W. Separation methods forisolation of human polymorphonuclear leukocytes affect their motile activity.Eur J Cell Biol 2009;88:531.

[2] Kim YW, Yoo JY. Three-dimensional focusing of red blood cells in microchannelflows for bio-sensing applications. Biosens Bioelectron 2009;24:3677.

[3] Segre G, Silberberg A. Radial Poiseuille flow of suspensions. Nature1961;189:209.

[4] Segre G, Silberberg A. Behaviour of macroscopic rigid spheres in Poiseuilleflow. Part 2. Experimental results and interpretation. J Fluid Mech1962;14:136–57.

[5] Saffman PG. The lift on a small sphere in a slow shear flow. J Fluid Mech1966;22:385–400.

[6] Asmolov ES. Dynamics of a spherical particle in a laminar boundary layer. FluidDyn 1990;25:886–90.

[7] McLaughlin JB. Inertial migration of a small sphere in linear shear flows. J FluidMech 1991;224:261–74.

[8] Leighton DT, Acrivos A. The lift on a small sphere touching a plane in thepresence of a simple shear flow. Z Angew Math Phys 1985;36:174–8.

[9] Cherukat P, McLaughlin JB. The inertial lift on a rigid sphere in a linear shearflow field near a flat wall. J Fluid Mech 1994;263:1–18.

[10] Krishnan GP, Leighton DT. Inertial lift on a moving sphere in contact with aplane wall in a shear flow. Phys Fluid 1995;7:2538–45.

[11] Gauthier F, Goldsmith HL, Mason SG. Particle motions in non-Newtonianmedia. Part I: Couette flow. Rheol Acta 1971;10:344–64.

[12] Gauthier F, Goldsmith HL, Mason SG. Particle motions in non-Newtonianmedia. Part II: Poiseuille flow. J Rheol 1971;15:297–330.

[13] Batram E, Goldsmith HL, Mason SG. Particle motions in non-Newtonian media.Part III: further observations in elastoviscous fluids. Rheol Acta1975;14:776–82.

[14] Lormand BM, Phillips RJ. Sphere migration in oscillatory Couette flow of aviscoelastic fluid. J Rheol 2004;48:551–70.

[15] Ho BP, Leal LG. Migration of rigid spheres in a two-dimensional unidirectionalshear flow of a second-order fluid. J Fluid Mech 1976;76:783–99.

[16] Huang PY, Feng J, Hu HH, Joseph DD. Direct simulation of the motion of solidparticles in Couette and Poiseuille flows of viscoelastic fluids. J Fluid Mech1997;343:73–94.

[17] Larson RG. Constitutive equations for polymer melts and solutions.Butterworth-Heinemann 1988.

[18] Hu HH, Patankar NA, Zhu MY. Direct numerical simulations of fluid–solidsystems using the arbitrary Lagrangian–Eulerian technique. J Comp Phys2001;169:427–62.

[19] Fattal R, Kupferman R. Constitutive laws for the matrix-logarithm of theconformation tensor. J Non-Newton Fluid Mech 2004;123:281–5.

[20] Hulsen MA, Fattal R, Kupferman R. Flow of viscoelastic fluids past a cylinder athigh Weissenberg number: stabilized simulations using matrix logarithms. JNon-Newton Fluid Mech 2005;127:27–39.

[21] Guenette R, Fortin M. A new mixed finite element method for computingviscoelastic flows. J Non-Newton Fluid Mech 1995;60:27–52.

[22] Bogaerds ACB, Grillet AM, Peters GWM, Baaijens FPT. Stability analysis ofpolymer shear flows using the extended pom–pom constitutive equations. JNon-Newton Fluid Mech 2002;108:187–208.

[23] Brooks AN, Hughes TJR. Streamline upwind/Petrov–Galerkin formulations forconvection dominated flows with particular emphasis on the incompressibleNavier–Stokes equations. Comput Methods Appl Mech Eng 1982;32:199–259.

[24] D’Avino G, Maffettone PL, Hulsen MA, Peters GWM. Numerical simulation ofplanar elongational flow of concentrated rigid particle suspensions in aviscoelastic fluid. J Non-Newton Fluid Mech 2007;150:65–79.

[25] Schenk O, Gartner K. Solving unsymmetric sparse systems of linear equationswith PARDISO. Future Gener Comput Syst 2004;20:475–87.

[26] D’Avino G, Hulsen MA, Snijkers F, Vermant J, Greco F, Maffettone PL. Rotationof a sphere in a viscoelastic liquid subjected to shear flow. Part I: Simulationresults. J Rheol 2008;52:1331–46.

[27] D’Avino G, Cicale G, Hulsen MA, Greco F, Maffettone PL. Effects of confinementon the motion of a single sphere in a sheared viscoelastic liquid. J Non-NewtonFluid Mech 2009;157:101–7.

[28] Glowinski R, Pan T-W, Joseph DD, Hesla TI. A distributed Lagrangianmultipliers/fictitious domain method for particulate flows. Int J MultiphaseFlow 1999;25:755–94.

[29] Hwang WR, Hulsen MA, Meijer HEH. Direct simulation of particle suspensionsin sliding bi-periodic frames. J Comput Phys 2004;194:742–72.

[30] D’Avino G, Hulsen MA. A comparison between a collocation and weakimplementation of the rigid-body motion constraint on a particle surface. IntJ Numer Methods Fluids 2009. doi:10.1002/fld.2185.

[31] Sullivan MT, Moore K, Stone HA. Transverse instability of bubbles inviscoelastic channel flows. Phys Rev Lett 2008;101:244503.