Inertial focusing of finite-size particles in microchannels
Evgeny S. Asmolov, Alexander L. Dubov, Tatiana V. Nizkaya, Jens Harting, Olga I. Vinogradova
aa r X i v : . [ phy s i c s . f l u - dyn ] O c t Under consideration for publication in J. Fluid Mech. Inertial focusing of finite-size particles inmicrochannels
By Evgeny S. Asmolov , † , Alexander L. Dubov , Tatiana V. Nizkaya , JensHarting , , , and Olga I. Vinogradova , , ‡ A.N. Frumkin Institute of Physical Chemistry and Electrochemistry,Russian Academy of Sciences, 31 Leninsky Prospect, 119071 Moscow, Russia Institute of Mechanics, M. V. Lomonosov Moscow State University, 119991 Moscow, Russia Helmholtz Institute Erlangen-N¨urnberg for Renewable Energy, Forschungszentrum J¨ulich,F¨urther Str. 248, 90429 N¨urnberg, Germany Department of Applied Physics, Eindhoven University of Technology, PO box 513, 5600MB Eindhoven,The Netherlands Faculty of Science and Technology, University of Twente, 7500 AE Enschede, The Netherlands Department of Physics, M. V. Lomonosov Moscow State University, 119991 Moscow, Russia DWI - Leibniz Institute for Interactive Materials, Forckenbeckstr. 50, 52056 Aachen, Germany(Received Received: date / Accepted: date)
At finite Reynolds numbers, Re, particles migrate across laminar flow streamlines to their equi-librium positions in microchannels. This migration is attributed to a lift force, and the balancebetween this lift and gravity determines the location of particles in channels. Here we demon-strate that velocity of finite-size particles located near a channel wall differs significantly fromthat of an undisturbed flow, and that their equilibrium position depends on this, referred to asslip velocity, difference. We then present theoretical arguments, which allow us to generalize ex-pressions for a lift force, originally suggested for some limiting cases and Re ≪
1, to finite-sizeparticles in a channel flow at Re ≤
20. Our theoretical model, validated by lattice Boltzmannsimulations, provides considerable insight into inertial migration of finite-size particles in mi-crochannel and suggests some novel microfluidic approaches to separate them by size or densityat a moderate Re.
1. Introduction
Microfluidic systems have been shown to be very useful for continuous manipulation and sep-aration of microparticles with increased control and sensitivity, which is important for a widerange of applications in chemistry, biology, and medicine. Traditional microfluidic techniques ofparticle manipulation rely on low Reynolds number laminar flow. Under these conditions, whenno external forces are applied, particles follow fluid streamlines. Contrary to this, particles mi-grate across streamlines to some stationary positions in microchannels when inertial aspects ofthe flow become significant. This migration is attributed to inertial lift forces, which are cur-rently successfully used in microfluidic systems to focus and separate particles of different sizescontinuously, at high flow rate, and without external forces (Di Carlo et al. et al. † Email address for correspondence: [email protected]‡ Email address for correspondence: [email protected]
E. S. Asmolov et al lift forces on particles in confined flows. We mention below what we believe are the most relevantcontributions.Inertial lift forces on neutrally buoyant particles have been originally reported for macroscopicchannels (Segr´e & Silberberg 1962). This pioneering work has concluded that particles focus toa narrow annulus at radial position 0 . et al. (2016) for recent reviews). In these microfluidic applications the inertial lift has beenbalanced by the Dean force due to a secondary rotational flow caused by inertia of the fluid itself,which can be generated in curved channels (Bhagat et al. et al. et al. et al. et al. et al. et al. et al. et al. (2004) have shown that the Segr´e-Silberbergannulus for neutrally-buoyant particles shifts toward the wall as Re increases and toward the pipecenter as particle size increases. At large Re ≥ et al. (2017) have found that the inner annulus is not a trueequilibrium position, but a transient zone, and that in a long enough pipe all particles accumulatewithin the Segr´e-Silberberg annulus. It has also been found that equilibrium positions of slightlynon neutrally-buoyant particles in a horizontal pipe are shifted toward a pipe bottom (Matas et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. (2009) where a pipe flow has been addressed. The approachcan be applied when the particle Reynolds number, Re p = a G / ν , where a is the particle radius, G is the characteristic shear rate and ν is the kinematic viscosity, is small. If so, to the leadingorder in Re p , the disturbance flow is governed by the Stokes equations, and a spherical particleexperiences a drag and a torque, but no lift. The Stokeslet disturbance originates from the parti-cle translational motion relative to the fluid and is proportional to the slip velocity V ′ s = V ′ − U ′ ,where V ′ and U ′ are forward velocities of the particle and of the undisturbed flow at the particlecenter. The stresslet is induced by free rotation of the sphere in the shear flow and is proportionalto G . The lift force has then been deduced from the solution of the next-order equations whichaccounts a non-linear coupling between the two disturbances (Vasseur & Cox 1976): F ′ l = ρ a (cid:0) c l a G + c l aGV ′ s + c l V ′ s (cid:1) , (1.1)where ρ is the fluid density. The coefficients c li ( i = , ,
2) generally depend on several di- nertial focusing of finite-size particles in microchannels z / a , H / a , V ′ s / U ′ m , and on the channel Reynolds number, Re = U ′ m H / ν , where z is the distance to the closest wall, H is the channel thickness, and U ′ m is themaximum velocity of the channel flow. Solutions for c l have been obtained in some limitingcases only, and no general analytical equations have still been proposed for finite-sized particlesin the channel. Thus, Vasseur & Cox (1976) have calculated the coefficients c VCl , c VCl , c VCl forpointlike particles at small channel Reynolds numbers, Re ≪
1, which depend on z / H only andare applicable when z ≫ a . Cherukat & Mclaughlin (1994) have later evaluated the coefficients c CMli ( z / a ) for finite-size particles near a single wall in a linear shear flow assuming that z ∼ a andproposed simple fits for them. However, it remains unclear if and how earlier theoretical resultsfor pointlike particles at Re ≪ z and a finite Re of a microfluidic channel.According to Equation (1.1) the contribution of the slip velocity to the lift forces dominateswhen V ′ s ≫ Ga . Since the slip velocity is induced by external forces, such as gravity, it is be-lieved that it impacts a hydrodynamic lift only in the case of non-neutrally buoyant particles. Forneutrally buoyant particles with equal to ρ density, the slip velocity is normally considered to benegligibly small (Ho & Leal 1974; Hood et al. z ≫ a , but hydrodynamic interac-tions at finite distances z ∼ a can induce a finite slip, V ′ s ∼ Ga , so that all terms in Equation (1.1)become comparable (Cherukat & Mclaughlin 1994). The variation of the slip velocity of neu-trally buoyant particles in a thin near-wall layer can impact the lift force, but we are unaware ofany previous work that has addressed this question.The purpose of this introduction has been to show that, in spite of its importance for inertialmicrofluidics, the lift forces of finite-size particles in a bounded geometry of a microchannelstill remain poorly understood. In particular, there is still a lack of simple analytical formulasquantifying the lift, as well as of general solutions valid in the large range of parameters typicalfor real microfluidic devices. Given the current upsurge of interest in the inertial hydrodynamicphenomena and their applications to separation of particles in microfluidic devices it would seemtimely to provide a more satisfactory theory of a hydrodynamic lift in a microchannel and alsoto bring some of modern simulation techniques to bear on this problem. In this paper we presentsome results of a study of a migration of finite-size particles at moderate channel Reynoldsnumbers, Re ∼
10, with the special focus on the role of the slip velocity in the hydrodynamic lift.Our paper is arranged as follows. In § G m a . Toaccess the validity of the proposed theory we use a simulation method described in §
3, and thenumerical results are presented in §
4. We conclude in §
2. Theory
In this section we propose an analytical expression for the lift force on neutrally buoyant andslightly non-neutrally buoyant particles of radius a , which translate parallel to a channel wall.Our expression is valid for a / H ≪ z from the channel wall.We consider a pressure-driven flow in a flat inclined microchannel of thickness H . An incli-nation angle α ≥ x is parallel to the E. S. Asmolov et al F IGURE
1. Sketch of a migration of a particle of radius a to an equilibrium position in a pressure-drivenflow. The locus of this position is determined by the balance between lift, F ′ l , and gravity, F ′ g , forces. channel wall, and the normal to the wall coordinate is denoted by z . The geometry is shown inFigure 1. The undisturbed velocity profile in such a channel is given by U ′ ( z ) = U ′ m z ( − z / H ) / H . (2.1)Let us now introduce a dimensionless slip velocity V s = V ′ s / ( aG m ) , where G m = U ′ m / H is themaximum shear rate at the channel wall. We can then rewrite Equation (1.1) as F ′ l = ρ a G m c l , (2.2)with the lift coefficient c l = c l + c l V s + c l V s , (2.3)which depends on the slip velocity, V s , which in turns can be determined from the Stokes equa-tions (a zero-order solution). Therefore, to construct a general expression for a lift force actingon finite size particles in the channel it is necessary to estimate V s as a function of z .We begin by studying the classical case of neutrally buoyant (i.e. force- and a torque-free)particles with a density ρ p equal to that of liquid, ρ . The expression for V s in a linear shear flownear a single wall has been derived before (Goldman et al. V s are given in Appendix A,Equations (A 2)-(A 4). We first note that depending on z / a one can distinguish between twodifferent regimes of behavior of V s . In the central part of the channel, i.e. when z / a ≫
1, the slipcontribution to the lift decays as ( a / z ) (Wakiya et al. z / a − ≪
1, the slip velocityvaries very rapidly with z / a (Goldman et al. V nbs = − + . . − .
200 log ( z / a − ) . (2.4)As a side note we should like to mention here that a logarithmic singularity in Equation (2.4)implies that in the near-wall region the lift coefficient, Equation (2.3), cannot be fitted by anypower law ( a / z ) n as it has been previously suggested (Di Carlo et al. et al. et al. z = a , theslip velocity is largest, V s = −
1. In this limiting case the lift coefficient also takes its maximumvalue, c KLl ≃ .
257 (Krishnan & Leighton Jr. 1995). Far from the wall, the slip velocity is much nertial focusing of finite-size particles in microchannels c l ≃ c l . Therefore, when a ≪ z ≪ H , thevalue of c l in Equation (2.3) is equal to c CVl | z / H → = π / ≃ . c l varies significantlyin the vicinity of the wall due to a finite slip.We now remark that the Stokeslet contribution (the second and the third terms in Equa-tion (2.3)) is finite for z ∼ a only and vanishes in the central part of the channel. Within theclose proximity to the wall we may neglect the corrections to the slip and the lift of order a / H due to parabolic flow (Pasol et al. c CMli . The stresslet contribution to the lift (first term in Equation (2.3)) is finitefor any z . Close to the wall, the effect of particle size for this term is negligible as the coefficient c CMl ( z / a ) is nearly constant (Cherukat & Mclaughlin 1994). So we may describe the stressletcontribution by the coefficient c VCl obtained by Vasseur & Cox (1976). This enables us to con-struct the following formula for the lift coefficient: c l = c VCl ( z / H ) + γ c CMl ( z / a ) V s + c CMl ( z / a ) V s , (2.5)where γ = G ( z ) / G m = − z / H ≤ c VCl , Equation (A 9) to calculate c CMl , and Equation (A 10) for c CMl . Note that in the second term of Equation (2.5) we have introduced a correction factor γ ,which takes into account the variation of G in the second term of Equation (1.1) and ensures thelift to remain zero at the channel centerline.We recall, that Equation (2.5) is asymptotically valid for any z when a / H ≪ ≪ ≤
20. By this reason constructed for Re ≪ F ′ g , which in dimensionless form can beexpressed as F g = F ′ g ρ a G m = π g aG m ∆ρ , (2.6)where ∆ρ = ( ρ p − ρ ) / ρ . The gravity influences both the particle migration and equilibriumposition when F g = O ( ) . It also induces an additional slip velocity which is of the order of theStokes settling velocity, V St = F ′ g πµ a G m = Re p F g π , (2.7)where µ is the dynamic viscosity. The effect of this velocity on the lift is comparable to F nbl when V St = O ( ) , i.e., at large gravity, F g ∼ π Re − p ≫
1, and is very important for vertical ornearly vertical channels. For horizontal channels, the slip velocity is equal to that of a neutrallybuoyant sphere since F x =
0. Equation (2.5) can also be applied in this case since the slip velocityremains small far from walls. Equilibrium positions of particles, z eq , can then be deduced fromthe balance between the lift and the gravity, c l ( z eq ) = F g . (2.8) E. S. Asmolov et al
Equation(2.8) may have two, one or no stable equilibrium points depending on F g , and the sensi-tivity of equilibrium positions to the value of a or ∆ρ is defined by the value ∂ c l / ∂ z . Thus, whenthe derivative is small, small variations in F g will lead to a significant shift in focusing positions.We finally note that the range of possible z eq can be tuned by the choice of U ′ m .
3. Simulation method
In this section, we present our simulation method and justify the choice of parameters.For our computer experiment, we chose a scheme based on the lattice Boltzmann method (Benzi et al. et al. et al. z = z = δ , so that in all simulations H = δ , andtwo periodic boundaries with N x = N y = δ , where δ is the lattice spacing. Spherical particlesof radii a = δ − δ are implemented as moving no-slip boundaries (Ladd & Verberg 2001;Janoschek et al. et al. − ∇ p . We use a 3D, 19 veloc-ity, single relaxation time implementation of the lattice Boltzmann method, where the relaxationtime τ is kept to 1 throughout this paper. Different flow rates are obtained by changing the fluidforcing. We use two channel Reynolds numbers, Re = . .
6. To simulate the migrationin an inclined channel we apply the gravity force directed at an angle α relative to the z -axis atthe center of the particle. In our simulations the values of dimensionless F g vary from 0 (neutrallybuoyant particle) to 13 . x − and z − componentsof the particle velocity to find the dimensionless slip, V s = ( V ′ x − U ′ ( z )) / ( aG m ) , and migrationvelocities, V m = V ′ z / ( aG m ) . To suppress the fluctuations arising from the discretization artifactswe average the velocities over approximately 4000 timesteps. The error does not exceed 3% forthe particles of a = a . The lift force can then be found from thesecalculations, by assuming that the particle motion is quasi-stationary. The lift is balanced by the z − component of the drag, F ′ l = − F ′ dz . Following Dubov et al. (2014) we use an expression F ′ dz ≈ − πµ aV ′ m f z ( z / H , a / H ) , (3.1) f z = + az − a + aH − a − z , (3.2)where the second and the third terms are corrections to the Stokes drag due to hydrodynamicinteractions with two channel walls. In what follows c l = π V m f z Re − p . (3.3)Second method to calculate the lift (and to check the validity of the first approach) uses thebalance of the lift and the gravity forces described by Equation (2.8). By varying the gravityforce F g one can, therefore, comprehend the whole range of equilibrium positions within thechannel to obtain c l ( z ) . The advantage of such an approach is that it does not require the particlemotion to be quasi-stationary. However, the disadvantage of this method is that the convergenceto equilibrium can be slow in the central zones of the channel, where the slope of c l ( z ) is small.Therefore, we use this computational strategy only in the near-wall region. nertial focusing of finite-size particles in microchannels z/H -4-20246 V nb m × -3 (a) V nb s (b) F IGURE
2. (a) Dimensionless migration velocity computed as a function of z / H for particles of a = δ (symbols). The location of the particle in a contact with the wall, z = a , is shown by a vertical dotted line.Dashed curve plots theoretical predictions for pointlike particles. Solid curve shows a polynomial fit ofsimulation data. (b) Dimensionless slip velocities computed for the same particles (symbols). Solid curveplots the slip velocity in a linear shear flow near a single wall. Dashed line plots the Faxen correction.Vertical dotted line indicates the location of z = a .
4. Results and discussion
In this section, we present the lattice Boltzmann simulation results and compare them withtheoretical predictions. 4.1.
Neutrally buoyant particles
We start with neutrally buoyant particles and first calculate their migration V nbm and the slip V nbs velocities as a function of z / H . Figure 2 plots simulation data obtained for particles of radius a = δ . Here we show only a half of the channel since the curves are antisymmetric with re-spect to the channel axis z = H /
2. These results demonstrate that migration velocity differssignificantly from the velocity c l Re p / ( π ) , where Re p = a G m / ν , predicted theoretically forpointlike particles (Vasseur & Cox 1976). We also see that the equilibrium position, V nbm =
0, offinite-size particles is shifted towards the channel axis compared to that of pointlike particles,which is obviously due to their interactions with the wall resulting in a finite slip velocity. In-deed, Figure 2(b) demonstrates that computed V nbs grows rapidly near the wall being close tothe theoretical predictions for a linear shear flow near a single wall (Goldman et al. et al. (1967), the computed slip velocity does not van-ish in the central part of the channel. Its value is roughly twice larger than the Faxen correction4 U ′ m a / ( H ) (Happel & Brenner 1965). Note that a similar difference has been obtained in sim-ulations of the migration of finite-size particles based on the Force Coupling Method (Loisel et al. c l for particles of a = δ and 8 δ . The lift coefficient has been obtained fromthe migration velocity and from the force balance as specified above, and simulations have beenmade for two moderate Reynolds numbers, Re = . .
6. As we discussed above, if Re ≤
20 a potential dependence of c l on Re could be ruled out a priori , and this is indeed confirmed byour simulations. Therefore, below we provide a detailed comparison of our simulation data withasymptotic solutions obtained for Re ≪
1, which should be applicable for finite moderate Re.Figure 3 also includes theoretical predictions by Vasseur & Cox (1976) and curves calculatedwith Equation (2.5). One can see that simulation results show strong discrepancy from point-particle approximation, especially in the near-wall region, where hydrodynamic interactions are
E. S. Asmolov et al z/H c l F IGURE
3. Lift coefficient, c l , for neutrally buoyant particles of a = δ (circles) and 8 δ (squares) obtainedfrom the migration velocity at Re = . . a = δ and 8 δ , dashed curve plots predictions for point–like particles. Vertical dotted lines show z = a . Black symbols show c l obtained for non-neutrally buoyantparticles of a = δ from the force balance at Re = .
6. The inset plots c l in the central part of the channel. a/H z eq / H F IGURE
4. Equilibrium positions for neutrally buoyant finite-sized (gray circles) and point-like (whitecircle) particles. Dashed curve shows predictions of Equation (2.5). significant. This discrepancy increases with the size of particles. We can, however, conclude thatpredictions of our Equation (2.5) are generally in good agreement with simulation results. Thus,for smaller particles, of a = δ , Equation (2.5) perfectly fits the simulation data in the near-wall region, where the theory for point-like particles fails. Simulation results slightly deviatefrom predictions of Equation (2.5) near the equilibrium positions and in the central part of thechannel. For bigger particles, of a = δ , these deviations are more pronounced. We emphasize,however, that they are still much smaller than from the point-particle theory by Vasseur & Cox(1976).To examine a significance of the particle size in more detail, we plot in Figure 4(a) computedequilibrium position, z eq / H , as a function of a / H . We recall that the lift c nbl ( z ) is antisymmetricwith respect to the midplane of the channel axis, so that neutrally buoyant particles have a secondequilibrium position at H − z eq . In a point-particle approximation z eq / H ≃ .
19 (Vasseur & Cox1976). We see that for finite-size particles z eq / H is always larger, and increases with the particlesize. Note that the increase in z eq / H is nearly linear when a / H ≤ .
1. Also included in Figure 4 nertial focusing of finite-size particles in microchannels z/H -10-8-6-4-202 V m × -3 F IGURE
5. Migration velocity of non-neutrally buoyant particles in a horizontal channel. Symbols showsimulation data. Solid curve is calculation with Equation (4.1) using data for neutrally buoyant particles. are predictions of Equation (2.5). One can conclude that the theory correctly predicts the trendobserved in simulations, but slightly deviates from the simulation data. A possible explanationfor this discrepancy could be effects of parabolic flow (which are of the order of O ( a / H ) ) onthe slip velocity and the stresslet (see Yahiaoui & Feuillebois 2010; Hood et al. Non-neutrally buoyant particles
We now turn to the particle migration under both inertial lift and gravity forces.4.2.1.
Horizontal channel
Let us start with the investigation of migration of particles in a most relevant experimentallycase of a horizontal channel ( α = ◦ ).We first fix a weak gravity force, F g = . a = δ in a horizontal channel. Simulation results are plotted in Figure 5. We see that V m ( z ) is no longer antisymmetric, as it has been in the case of neutrally buoyant particles. Themigration velocity can be calculated as V m = V nbm − V St / f z , (4.1)where we use a fit for V nbm computed for neutrally-buoyant particles (see Figure 2(a)). The agree-ment between simulation data and calculations using Equation (4.1) is excellent, which confirmsthat Equation (2.5) remains valid in the case of slightly non-neutrally buoyant particles. We re-mark that due to gravity V m is shifted downwards relative to V nbm ( z ) shown in Figure 2. As aresult, with the taken value of F g the second equilibrium position disappeared.We recall that this type of simulations allows one to find values of c l ( z ) in the vicinity of thewall by varying F g . We have included these force balance results in Figure 3 and can concludethat they agree very well with data obtained by using another computational method and forneutrally-buoyant particles. This suggests again that above results could be used at moderateReynolds numbers, Re ≤
20, since in this case the lift coefficient does not depend on Re.Figure 6(a) shows z eq / a computed at several F g . It can be seen that when the gravity forceis getting larger, the equilibrium positions decrease rapidly. This trend can be used to separateparticles even when ∆ρ is very small. To illustrate this we now fix Re = .
3, inject particles of a = δ close to the bottom of the channel, z = . a , and simulate their trajectories at different F g . In Figure 6(b) we plot trajectories of particles, z / a , as a function of xG m a ν . The data show0 E. S. Asmolov et al g z eq / a ⋅ aG m / ν z / a F IGURE
6. (a) Equilibrium positions of non-neutrally buoyant particles ( a = δ ) in a horizontal channel;(b) trajectories of the same particles released at z = . a computed at different F g = .
268 (squares),0 .
804 (diamonds), 1 .
340 (triangles), 2 .
681 (circles), 9 .
383 (diamonds). that if F g is large enough, particles sediment to the wall. However, when F g is relatively small,particles follow different and divergent trajectories, by approaching their equilibrium positions.We stress that at a given F g and a / H trajectories, shown Figure 6(b), remain the same for anyRe ≤
20 (see Appendix B). Therefore, even in the case of very small ∆ρ , one can always tunethe value of Re to induce the required for separation difference in F g . For example, we haveto separate particles of a = µ m and different ∆ρ in the channel of H = µ m. If we choseRe = .
3, the separation length L = xG m a ν of Figure 6(b) will be ca. 3 . ∆ρ with Equation (2.6), we can immediately conclude that trajectories plotted in Figure 6(b) fromtop to bottom correspond to ∆ρ = . . .
037 and 0 . Inclined channel
When F g is large enough, it can also influences the slip velocity, and therefore, change the liftitself. This effect is especially important for vertical channels. Note that due to the linearity ofthe Stokes equations, which govern a disturbance flow at small particle Reynolds numbers, wecan decouple the contributions of the particle-wall interaction and of the gravity force into theslip velocity: V s = V nbs + ∆ V s sin α , (4.2)where ∆ V s = V St / f x is the gravity-induced slip velocity for a vertical channel ( α = ◦ ) and f x ( z / H , a / H ) is the correction to the drag for a particle translating parallel to the channel walls.The slip and the migration velocities of particles of a = δ in a vertical channel computed byusing several values of F g are shown in Figures 7(a) and 8(a). Note that the slip velocity, V s , growswith F g since the Stokes velocity, V St , is linearly proportional to F g (see Eq.(2.7)). We now usesimulation data presented in Figures 2(a) and 7(a) to compute ∆ V s , and then ∆ V s / F g . The resultsfor ∆ V s / F g are shown in Figure 7(b), and we see that all data collapse into a single curve, whichconfirms the validity of Equation (4.2). Figure 7(b) also shows that ∆ V s / F g is nearly constant inthe central region of the channel, being smaller than V St , but the deviations from V St grow whenparticles approach the wall. These results again illustrate that hydrodynamic interactions with thewalls significantly affect motion of particles in the channel.We recall that the variation of the slip velocity caused by gravity is small for slightly non- nertial focusing of finite-size particles in microchannels z/H -0.2-0.15-0.1-0.0500.050.1 V s (a) z/H ∆ V s / F g (b) F IGURE
7. Slip velocities (a) and ∆ V s / F g (b) computed for non-neutrally buoyant particles of a = δ ina vertical channel. The data sets correspond to F g = .
475 (circles), 6 .
956 (triangles), 10 .
44 (squares) and13 .
91 (diamonds). Dashed line shows V St / F g , vertical dotted lines plot z = a .F IGURE
8. Migration velocities (a) and ∆ V m / F g (b) computed for non-neutrally buoyant particles of a = δ in a vertical channel. The data sets correspond to F g = .
475 (circles), 6 .
956 (triangles), 10 .
44 (squares)and 13 .
91 (diamonds). Vertical dotted lines plot z = a . Solid curve shows a polynomial fit of data. neutrally buoyant particles (see Figure 7), so that Eq.(2.5) can be linearized with respect to ∆ V s : c l ≃ c nbl + ∆ V s ∂ c l ( V nbs ) ∂ V s , (4.3)where c nbl = c l ( V nbs ) is the lift coefficient for neutrally buoyant particles. By using Eq.(3.3) wecan then calculate the migration velocity V m = V nbm + ∆ V m = V nbm + ∆ V s ∂ c l ( V nbs ) ∂ V s Re p π f z . (4.4)The computed migration velocity is shown in Figure 8(a). We see that it decreases with F g , andthe equilibrium position shifts towards the wall, which is since ∆ V s / F g is positive while ∂ c l / ∂ V s is negative.We can now evaluate ∆ V m / F g by using simulation data presented in Figures 2 and 8(a), andthese results are presented in Figure 8(b). As one can see, the data collapse into a single curve,thus confirming the validity of our linearization, Equation (4.4).2 E. S. Asmolov et al α , ° z eq / H F IGURE
9. Equilibrium positions z eq / H for a = δ and F g = . V m =
0, where V m is calculated with Eq. (4.6). Finally, we briefly discuss the case of an arbitrary inclination angle α , where the z -componentof the force can be written as F z = c l ( V s ) + F g cos α . (4.5)By using Eqs.(4.2), (4.4) and (4.5), we can express the migration velocity as V m = V nbm + ∆ V m sin α + F g cos α Re p π f z , (4.6)where ∆ V m is evaluated for a vertical channel (see Figure 8(b)). The equilibrium positions canbe found by using a condition V m =
0, where V m is calculated with Eq.(4.6). The results of thesecalculations made at a fixed F g = .
475 and different α are plotted in Figure 9 together with directsimulation data, and one can see that they practically coincide. Our results show that in a verticalchannel two stable equilibrium positions coexist. They are symmetric relative to the midplaneand are located close to walls. Another, third equilibrium position has a locus at the midplane,but is unstable. A similar result has been obtained earlier (Vasseur & Cox 1976; Asmolov 1999).If we slightly reduce α both stable equilibrium positions become shifted towards the lower walldue to gravity as well seen in Figure 9. These two positions coexist only for α ≥ . ◦ . Onreducing α further the upper equilibrium position disappears, and only one, a lower, equilibriumposition remains. This obviously indicates that the inertial lift cannot balance gravity anymore.We note that this remaining single equilibrium position becomes insensitive to the inclinationangle when α ≤ ◦ .
5. Concluding remarks
In this paper we have studied the inertial migration of finite-size particles in a plane channelflow at moderate Reynolds numbers, Re ≤
20. We have shown that the slip velocity, V s , which isfinite even for neutrally buoyant particles, contributes to the lift and determines the equilibriumpositions in the channel. We have proposed an expression for the lift which generalizes theo-ries, originally applied for some cases of limited guidance, to finite-size particles in a channelflow. When the size of particle turns to zero, our formula recovers known expression of a point-particle approximation (Vasseur & Cox 1976). For particles close to the walls we recover earlierpredictions for finite-size particles in a linear shear flow (Cherukat & Mclaughlin 1994). Ourtheoretical model, which is probably the simplest realistic model for a lift in the channel that onemight contemplate, provides considerable insight into inertial migration of finite-size particles in nertial focusing of finite-size particles in microchannels ≤
20, thatequilibrium positions of heavy particles in a horizontal channel can be accurately determined byusing data for the neutrally buoyant case, and more.Several of our theoretical predictions could be tested in experiment. In particular, we haveshown that particles of a very small density contrast should follow divergent trajectories, sothat channel flows with low Reynolds numbers Re ∼ ≥ et al. et al. et al. et al. et al. et al. Appendix A. Fits for the slip velocity and the lift coefficients
In this Appendix we summarize known results for the slip velocity and the lift coefficientsfor finite-size particles in a linear shear flow near a single wall and for point-like particles in aPoiseuille flow. The velocity of a freely translating and rotating particle in a linear shear flow isgiven by (Goldman et al. V nb ′ x = U ′ ( z ) h , (A 1)where h is the correction function which depends on z / a only. We use (A 1) to estimate the slipvelocity in channel flow, i.e., we neglect the effects due to parabolic flow, so that V nbs = z ( H − z )( h − ) aH . (A 2)4 E. S. Asmolov et al
The correction factor fitting the results by Goldman et al. (1967) in the near-wall region reads(Reschiglian et al. h = . b − ( . b + ) ζ − − . − . b + . b − ζ < , (A 3)where ζ = z / a and b = log ( ζ − ) . Note that here we have reformulated the original equa-tion (Reschiglian et al. et al. (1967), h = − ζ − + ζ − − ζ − − ζ − − ζ − + ζ − − ζ − − ζ − at ζ ≥ . (A 4)Vasseur & Cox (1976) have obtained the lift force on a particle in a channel flow at Re ≪ c VCl = c VCl + Ha c
VCl V s + c VCl V s , (A 5)where coefficients c VCl , c VCl , c VCl depend on z / H only. Later Feuillebois (2004) has proposed asimple fitting expression: c VCl = . ( z / H − . ) − . ( z / H − . ) . (A 6)The expression for a lift coefficient of a finite-size particle in a linear shear flow near a singlewall has been suggested by Cherukat & Mclaughlin (1994) c CMl = c CMl + c CMl V s + c CMl V s , (A 7)where the coefficients c CMl , c CMl , c CMl depend on ζ only: c CMl = . + . ζ − − . ζ − + . ζ − , (A 8) c CMl = − . ζ − . − . ζ − + . ζ − , (A 9) c CMl = . + . ζ − − . ζ − + . ζ − . (A 10) Appendix B. Governing equations for trajectories of particles
In this Appendix, we derive equations which govern particle trajectories. The components ofthe particle velocity can be written as dxdt = V ′ x = U ′ ( z ) h = G m z ( − z / H ) h , (B 1) dzdt = V ′ m = F ′ l − F ′ g πµ a f z = ( c l − F g ) a G m πν f z . (B 2)The last equality indicates that the migration time, i.e., the time required for a particle to migrateat distance of the order of its radius a , is equal to ν / ( G m a ) = ( G m Re ) − ( H / a ) . Since theright-hand sides of (B 1) and (B 2) do not explicitly include time, one can formulate an equationgoverning the particle trajectory as dzdx = a G m πν c l − F g f z z ( − z / H ) h . (B 3)Let us now turn to dimensionless coordinates ζ and ξ = xG m a / ν . Equation (B 3) can then be nertial focusing of finite-size particles in microchannels d ζ d ξ = π c l − F g f z h ζ ( − ζ a / H ) . (B 4)We stress that Equation (B 4) does not depend on Re. Indeed, f z and h are dimensionless func-tions of ζ and a / H only, and the lift coefficient, c l , is also not sensitive to the Reynolds numberwhen Re ≤
20. This implies that at given F g and a / H trajectories satisfying Equation (B 4) areuniversal, i.e. remain the same at any Re ≤ R E F E R E N C E SA
SMOLOV , E. S. 1999 The inertial lift on a spherical particle in a plane Poiseuille flow at large channelReynolds number.
J. Fluid Mech. , 63–87.A
SMOLOV , E. S., D
UBOV , A. L., N
IZKAYA , T. V., K
UEHNE , A. J. C. & V
INOGRADOVA , O. I. 2015Principles of transverse flow fractionation of microparticles in superhydrophobic channels.
Lab. Chip (13), 2835–2841.B ENZI , R., S
UCCI , S. & V
ERGASSOLA , M. 1992 The lattice Boltzmann equation: theory and applications.
Physics Reports , 145.B
HAGAT , A.A.S., K
UNTAEGOWDANAHALLI , S.S. & P
APAUTSKY , I. 2008 Continuous particle separationin spiral microchannels using dean flows and differential migration.
Lab. Chip (11), 1906–1914.C HERUKAT , P. & M
CLAUGHLIN , J. B. 1994 The inertial lift on a rigid sphere in a linear shear flow fieldnear a flat wall.
J. Fluid Mech. , 1–18.C
HOI , Y.-S., S EO , K.-W. & L EE , S.-J. 2011 Lateral and cross-lateral focusing of spherical particles in asquare microchannel. Lab. Chip (3), 460–465.C HUN , B. & L
ADD , A. J. C. 2006 Inertial migration of neutrally buoyant particles in a square duct: Aninvestigation of multiple equilibrium positions.
Phys. Fluids , 031704.C OX , R. G. & H SU , S. K. 1977 The lateral migration of solid particles in a laminar flow near a plane. Int.J. Multiph. Flow , 201–222.D AVIS , A. M. J., K
EZIRIAN , M. T. & B
RENNER , H. 1994 On the Stokes-Einstein model of surfacediffusion along solid surfaces: Slip boundary conditions.
J. Colloid Interface Sci. (1), 129–140.D I C ARLO , D., E DD , J. F., H UMPHRY , K. J., S
TONE , H. A. & T
ONER , M. 2009 Particle segregation anddynamics in confined flows.
Phys. Rev. Lett. (9), 094503.D I C ARLO , D., I
RIMIA , D., T
OMPKINS , R.G. & T
ONER , M. 2007 Continuous inertial focusing, ordering,and separation of particles in microchannels.
Proc. Natl. Acad. Sci. U.S.A. (48), 18892–18897.D
UBOV , A. L., S
CHMIESCHEK , S., A
SMOLOV , E. S., H
ARTING , J. & V
INOGRADOVA , O. I. 2014Lattice-Boltzmann simulations of the drag force on a sphere approaching a superhydrophobic stripedplane.
J. Chem. Phys. (3), 034707.D
UTZ , S., H
AYDEN , M.E. & H ¨
AFELI , U.O. 2017 Fractionation of magnetic microspheres in a microflu-idic spiral: Interplay between magnetic and hydrodynamic forces.
PLOS ONE (1), e0169919.F EUILLEBOIS , F. 2004 Perturbation problems at low Reynolds number.
Lecture Notes-AMAS .F EUILLEBOIS , F., B
AZANT , M. Z. & V
INOGRADOVA , O. I. 2010 Transverse flow in thin superhydropho-bic channels.
Phys. Rev. E , 055301(R).G OLDMAN , A. J., C OX , R. G. & B RENNER , H. 1967 Slow viscous motion of a sphere parallel to a planewall - II Couette flow.
Chem. Eng. Sci. , 653–660.H APPEL , J. & B
RENNER , H. 1965 Low Reynolds number hydrodynamics with special applications toparticulate media.
Prentice-Hall .H ARTING , J., F
RIJTERS , S., R
AMAIOLI , M., R
OBINSON , M., W
OLF , D. E. & L
UDING , S. 2014 Recentadvances in the simulation of particle-laden flows.
Eur. Phys. J. Spec. Topics , 2253–2267.H O , B.P. & L EAL , L.G. 1974 Inertial migration of rigid spheres in two-dimensional unidirectional flows.
J. Fluid Mech. , 365–400.H OOD , K., K
AHKESHANI , S., D I C ARLO , D. & R
OPER , M. 2016 Direct measurement of particle inertialmigration in rectangular microchannels.
Lab. Chip , 2840–2850.H OOD , K., L EE , S. & R OPER , M. 2015 Inertial migration of a rigid sphere in three-dimensional Poiseuilleflow.
J. Fluid Mech. , 452–479.J
ANOSCHEK , F. 2013 Mesoscopic simulation of blood and general suspensions in flow.
Eindhoven Univer-sity of Technology (Ph. D. thesis) . E. S. Asmolov et al J ANOSCHEK , F., T
OSCHI , F. & H
ARTING , J. 2010 Simplified particulate model for coarse-grained hemo-dynamics simulations.
Phys. Rev. E , 056710.K ILIMNIK , A., M AO , W. & A LEXEEV , A. 2011 Inertial migration of deformable capsules in channel flow.
Phys. of Fluids (12), 123302.K RISHNAN , G. P. & L
EIGHTON J R ., D. T. 1995 Inertial lift on a moving sphere in contact with a planewall in a shear flow. Phys. Fluids (11), 2538–2545.K UNERT , C., H
ARTING , J. & V
INOGRADOVA , O. I. 2010 Random-roughness hydrodynamic boundaryconditions.
Phys. Rev. Lett. (1), 016001.L
ADD , A. J. C. & V
ERBERG , R. 2001 Lattice-Boltzmann simulations of particle-fluid suspensions.
J. Stat.Phys. (5), 1191.L IU , C., H U , G., J IANG , X. & S UN , J. 2015 Inertial focusing of spherical particles in rectangular mi-crochannels over a wide range of Reynolds numbers. Lab. Chip (4), 1168–1177.L IU , C., X UE , C., S UN , J. & H U , G. 2016 A generalized formula for inertial lift on a sphere in microchan-nels. Lab. Chip (5), 884–892.L OISEL , V., A
BBAS , M., M
ASBERNAT , O. & C
LIMENT , E. 2015 Inertia-driven particle migration andmixing in a wall-bounded laminar suspension flow.
Phys. Fluids (12), 123304.M ARTEL , J. M. & T
ONER , M. 2014 Inertial focusing in microfluidics.
Annu. Rev. Biomed. Eng. , 371–396.M ATAS , J.-P., M
ORRIS , J. F. & G
UAZZELLI , E. 2004 Inertial migration of rigid spherical particles inPoiseuille flow.
J. Fluid Mech. , 171–195.M
ATAS , J.-P., M
ORRIS , J. F. & G
UAZZELLI , E. 2009 Lateral force on a rigid sphere in large-inertialaminar pipe flow.
J. Fluid Mech. , 59–67.M
IURA , K., I
TANO , T. & S
UGIHARA -S EKI , M. 2014 Inertial migration of neutrally buoyant spheres in apressure-driven flow through square channels.
J. Fluid Mech. , 320–330.M
ORITA , Y., I
TANO , T. & S
UGIHARA -S EKI , M. 2017 Equilibrium radial positions of neutrally buoyantspherical particles over the circular cross-section in Poiseuille flow.
J. Fluid Mech. , 750–767.N
ETO , C., E
VANS , D., B
ONACCURSO , E., B
UTT , H. J. & C
RAIG , V. J. 2005 Boundary slip in newtonianliquids: a review of experimental studies.
Rep. Prog. Phys. , 2859–2897.P ASOL , L., S
ELLIER , A. & F
EUILLEBOIS , F. 2006 A sphere in a second degree polynomial creeping flowparallel to a wall.
QJ Mech. Appl. Math. (4), 587–614.P IMPONI , D., C
HINAPPI , M., G
UALTIERI , P. & C
ASCIOLA , C. M. 2014 Mobility tensor of a spheremoving on a superhydrophobic wall: application to particle separation.
Microfluidics Nanofluidics ,571–585.R ESCHIGLIAN , P., M
ELUCCI , D., T
ORSI , G. & Z
ATTONI , A. 2000 Standardless method for quantitativeparticle-size distribution studies by gravitational field-flow fractionation. application to silica particles.
Chromatographia (1-2), 87–94.S AFFMAN , P. G. 1965 The lift on a small sphere in a slow shear flow.
J. Fluid Mech. , 385–400.S CHMIESCHEK , S., B
ELYAEV , A. V., H
ARTING , J. & V
INOGRADOVA , O. I. 2012 Tensorial slip of super-hydrophobic channels.
Phys. Rev. E , 016324.S CHONBERG , J. A. & H
INCH , E. J. 1989 Inertial migration of a sphere in Poiseuille flow.
J. Fluid Mech. , 517–524.S
EGR ´ E , G. & S ILBERBERG , A. 1962 Behaviour of macroscopic rigid spheres in Poiseuille flow. Part 2.Experimental results and interpretation.
J. Fluid Mech. , 136–157.V ASSEUR , P. & C OX , R. G. 1976 The lateral migration of a spherical particle in two-dimensional shearflows. J. Fluid Mech. , 385–413.V INOGRADOVA , O. I. 1996 Hydrodynamic interaction of curved bodies allowing slip on their surfaces.
Langmuir , 5963–5968.V INOGRADOVA , O. I. 1999 Slippage of water over hydrophobic surfaces.
Int. J. Miner. Proc. , 31–60.V INOGRADOVA , O. I. & B
ELYAEV , A. V. 2011 Wetting, roughness and flow boundary conditions.
J. Phys.:Condens. Matter , 184104.W AKIYA , S., D
ARABANER , C.L. & M
ASON , S.G. 1967 Particle motions in sheared suspensions XXI:Interactions of rigid spheres (theoretical).
Rheol. Acta (3), 264–273.Y AHIAOUI , S. & F
EUILLEBOIS , F. 2010 Lift on a sphere moving near a wall in a parabolic flow.
J. FluidMech. , 447–474.Z
HANG , J., Y AN , S., A LICI , G., N
GUYEN , N.-T., D I C ARLO , D. & L I , W. 2014 Real-time control ofinertial focusing in microfluidics using dielectrophoresis (dep). RSC Adv. (107), 62076–62085. nertial focusing of finite-size particles in microchannels Z HANG , J., Y AN , S., Y UAN , D., A
LICI , G., N
GUYEN , N.-T., W
ARKIANI , M. E. & L I , W. 2016 Funda-mentals and applications of inertial microfluidics: a review. Lab. Chip16