2 Theory

In this section we propose an analytical expression for the lift force on neutrally buoyant and slightly non-neutrally buoyant particles of radius $a$ , which translate parallel to a channel wall. Our expression is valid for $a/H\ll 1$ at any distance $z$ from the channel wall.

We consider a pressure-driven flow in a flat inclined microchannel of thickness $H$ . An inclination angle $\unicode[STIX]{x1D6FC}\geqslant 0$ is defined relative to the horizontal. The coordinate axis $x$ is parallel to the channel wall, and the normal to the wall coordinate is denoted by $z$ . The geometry is shown in figure 1. The undisturbed velocity profile in such a channel is given by

(2.1) $$\begin{eqnarray}U^{\prime }(z)=4U_{m}^{\prime }z(1-z/H)/H.\end{eqnarray}$$

Let us now introduce a dimensionless slip velocity $V_{s}=V_{s}^{\prime }/(aG_{m})$ , where $G_{m}=4U_{m}^{\prime }/H$ is the maximum shear rate at the channel wall. We can then rewrite equation (1.1) as

(2.2) $$\begin{eqnarray}F_{l}^{\prime }=\unicode[STIX]{x1D70C}a^{4}G_{m}^{2}c_{l},\end{eqnarray}$$

with the lift coefficient

(2.3) $$\begin{eqnarray}c_{l}=c_{l0}+c_{l1}V_{s}+c_{l2}V_{s}^{2},\end{eqnarray}$$

which depends on the slip velocity, $V_{s}$ , which in turn can be determined from the Stokes equations (a zero-order solution). Therefore, to construct a general expression for a lift force acting on finite-size particles in a 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 $\unicode[STIX]{x1D70C}_{p}$ equal to that of liquid, $\unicode[STIX]{x1D70C}$ . The expression for $V_{s}$ in a linear shear flow near a single wall has been derived before (Goldman, Cox & Brenner Reference Goldman, Cox and Brenner1967) and can be used to calculate the slip velocity in the near-wall region of our channel. The fits for $V_{s}$ are given in appendix A, equations (A 2)–(A 4). We first note that depending on $z/a$ one can distinguish between two different regimes of behaviour of $V_{s}$ . In the central part of the channel, i.e. when $z/a\gg 1$ , the slip contribution to the lift decays as $(a/z)^{3}$  (Wakiya, Darabaner & Mason Reference Wakiya, Darabaner and Mason1967), being always very small, but finite. In contrast, when the gap between the sphere and the wall is small, $z/a-1\ll 1$ , the slip velocity varies very rapidly with $z/a$  (Goldman et al. Reference Goldman, Cox and Brenner1967):

(2.4) $$\begin{eqnarray}V_{s}^{nb}=-1+\frac{0.7431}{0.6376-0.200\log (z/a-1)}.\end{eqnarray}$$

As a side note we should like to mention here that a logarithmic singularity in (2.4) implies that in the near-wall region the lift coefficient, equation (2.3), cannot be fitted by any power law $(a/z)^{n}$ as has been previously suggested (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Hood et al. Reference Hood, Lee and Roper2015; Liu et al. Reference Liu, Xue, Sun and Hu2016).

It follows from (2.4) that for an immobile particle in a contact with the wall, $z=a$ , the slip velocity is largest, $V_{s}=-1$ . In this limiting case the lift coefficient also takes its maximum value, $c_{l}^{KL}\simeq 9.257$  (Krishnan & Leighton Jr. Reference Krishnan and Leighton1995). Far from the wall, the slip velocity is much smaller and can be neglected, so that we can consider $c_{l}\simeq c_{l0}$ . Therefore, when $a\ll z\ll H$ , the value of $c_{l}$ in (2.3) is equal to $c_{l0}^{CV}|_{z/H\rightarrow 0}=55\unicode[STIX]{x03C0}/96\simeq 1.8$  (Cox & Hsu Reference Cox and Hsu1977), i.e. it becomes much smaller than for a particle at the wall. This illustrates that $c_{l}$ varies significantly in the vicinity of the wall due to a finite slip.

We now remark that the Stokeslet contribution (the second and the third terms in (2.3)) is finite for $z\sim a$ only and vanishes in the central part of the channel. Within the close proximity to the wall we may neglect the corrections to the slip and the lift of order $a/H$ due to parabolic flow (Pasol, Sellier & Feuillebois Reference Pasol, Sellier and Feuillebois2006; Yahiaoui & Feuillebois Reference Yahiaoui and Feuillebois2010) and due to the second wall. Therefore, in this region one can use the results by Cherukat & McLaughlin (Reference Cherukat and McLaughlin1994) for the lift coefficients $c_{li}^{CM}$ . The stresslet contribution to the lift (first term in (2.3)) is finite for any $z$ . Close to the wall, the effect of particle size for this term is negligible as the coefficient $c_{l0}^{CM}(z/a)$ is nearly constant (Cherukat & McLaughlin Reference Cherukat and McLaughlin1994). So we may describe the stresslet contribution by the coefficient $c_{l0}^{VC}$ obtained by Vasseur & Cox (Reference Vasseur and Cox1976). This enables us to construct the following formula for the lift coefficient:

(2.5) $$\begin{eqnarray}c_{l}=c_{l0}^{VC}(z/H)+\unicode[STIX]{x1D6FE}c_{l1}^{CM}(z/a)V_{s}+c_{l2}^{CM}(z/a)V_{s}^{2},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}=G(z)/G_{m}=1-2z/H\leqslant 1$ is a dimensionless local shear rate at the particle position. The fitting expressions for three lift coefficients are summarized in appendix A. We, therefore, use (A 6) to calculate $c_{l0}^{VC}$ , equation (A 9) to calculate $c_{l1}^{CM}$ and equation (A 10) for $c_{l2}^{CM}$ . Note that in the second term of (2.5) we have introduced a correction factor $\unicode[STIX]{x1D6FE}$ , which takes into account the variation of $G$ in the second term of (1.1) and ensures the lift remains zero at the channel centreline.

We recall, that equation (2.5) is asymptotically valid for any $z$ when $a/H\ll 1$ and $Re\ll 1$ . However, one can argue that it should be accurate enough at moderate Reynolds numbers. Indeed, the contribution of undisturbed flow to inertial terms in the Navier–Stokes equations remains relatively small when $Re\leqslant 20$ . With this reason constructed for $Re\ll 1$ regular-perturbation methods (Ho & Leal Reference Ho and Leal1974; Vasseur & Cox Reference Vasseur and Cox1976; Cherukat & McLaughlin Reference Cherukat and McLaughlin1994) have successfully predicted the lift force on a point-like neutrally buoyant particle at a moderate $Re$ . For larger $Re$ , when a contribution of inertial terms becomes significant, the equilibrium positions should be shifted towards the wall with the increase in $Re$  (Schonberg & Hinch Reference Schonberg and Hinch1989; Asmolov Reference Asmolov1999).

We now turn to non-neutrally buoyant particles, for which the density is different from that of the liquid, so that they experience an external gravity force, $F_{g}^{\prime }$ , which in dimensionless form can be expressed as

(2.6) $$\begin{eqnarray}F_{g}={\displaystyle \frac{F_{g}^{\prime }}{\unicode[STIX]{x1D70C}a^{4}G_{m}^{2}}}={\displaystyle \frac{4\unicode[STIX]{x03C0}g}{3aG_{m}^{2}}}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C},\end{eqnarray}$$

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C})/\unicode[STIX]{x1D70C}$ . The gravity influences both the particle migration and equilibrium position when $F_{g}=O(1)$ . It also induces an additional slip velocity which is of the order of the Stokes settling velocity,

(2.7) $$\begin{eqnarray}V^{St}={\displaystyle \frac{F_{g}^{\prime }}{6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a^{2}G_{m}}}={\displaystyle \frac{Re_{p}F_{g}}{6\unicode[STIX]{x03C0}}},\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ is the dynamic viscosity. The effect of this velocity on the lift is comparable to $F_{l}^{nb}$ when $V^{St}=O(1)$ , i.e. at large gravity, $F_{g}\sim 6\unicode[STIX]{x03C0}Re_{p}^{-1}\gg 1$ , and is very important for vertical or nearly vertical channels. For horizontal channels, the slip velocity is equal to that of a neutrally buoyant sphere since $F_{x}=0$ . Equation (2.5) can also be applied in this case since the slip velocity remains small far from walls. Equilibrium positions of particles, $z_{eq}$ , can then be deduced from the balance between the lift and the gravity,

(2.8) $$\begin{eqnarray}c_{l}(z_{eq})=F_{g}.\end{eqnarray}$$

Equation (2.8) may have two, one or no stable equilibrium points depending on $F_{g}$ , and the sensitivity of the equilibrium positions to the value of $a$ or $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ is defined by the value $\unicode[STIX]{x2202}c_{l}/\unicode[STIX]{x2202}z$ . Thus, when the 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}^{\prime }$ .

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, Succi & Vergassola Reference Benzi, Succi and Vergassola1992; Kunert, Harting & Vinogradova Reference Kunert, Harting and Vinogradova2010; Dubov et al. Reference Dubov, Schmieschek, Asmolov, Harting and Vinogradova2014) which has been successfully employed earlier to simulate a motion of particles in the channel flow. We use a simulation cell confined by two impermeable no-slip walls located at $z=0$ and $z=79\unicode[STIX]{x1D6FF}$ , so that in all simulations $H=79\unicode[STIX]{x1D6FF}$ , and two periodic boundaries with $N_{x}=N_{y}=256\unicode[STIX]{x1D6FF}$ , where $\unicode[STIX]{x1D6FF}$ is the lattice spacing. Spherical particles of radii $a=4\unicode[STIX]{x1D6FF}-12\unicode[STIX]{x1D6FF}$ are implemented as moving no-slip boundaries (Ladd & Verberg Reference Ladd and Verberg2001; Janoschek, Toschi & Harting Reference Janoschek, Toschi and Harting2010; Harting et al. Reference Harting, Frijters, Ramaioli, Robinson, Wolf and Luding2014), where the chosen radii are sufficient to keep discretization effects of the order of a few per cent (Janoschek Reference Janoschek2013). A Poiseuille flow is generated by applying a body force, which is equivalent to a pressure gradient $-\unicode[STIX]{x1D735}p$ . We use a three-dimensional, 19 velocity, single relaxation time implementation of the lattice Boltzmann method, where the relaxation time $\unicode[STIX]{x1D70F}$ is kept to 1 throughout this paper. Different flow rates are obtained by changing the fluid forcing. We use two channel Reynolds numbers, $Re=11.3$ and $22.6$ . To simulate the migration in an inclined channel we apply the gravity force directed at an angle $\unicode[STIX]{x1D6FC}$ relative to the $z$ -axis at the centre of the particle. In our simulations the values of dimensionless $F_{g}$ vary from $0$ (neutrally buoyant particle) to $13.91$ .

In our computer experiments we determine the lift by using two different strategies. In the first method we extract the lift from the migration velocity. We measure the $x$ - and $z$ -components of the particle velocity to find the dimensionless slip, $V_{s}=(V_{x}^{\prime }-U^{\prime }(z))/(aG_{m})$ , and migration velocities, $V_{m}=V_{z}^{\prime }/(aG_{m})$ . To suppress the fluctuations arising from the discretization artefacts we average the velocities over approximately 4000 time steps. The error does not exceed $3\,\%$ for the particles with $a=4$ and rapidly decreases with $a$ . The lift force can then be found from these calculations, by assuming that the particle motion is quasi-stationary. The lift is balanced by the $z$ -component of the drag, $F_{l}^{\prime }=-F_{dz}^{\prime }$ . Following Dubov et al. (Reference Dubov, Schmieschek, Asmolov, Harting and Vinogradova2014) we use an expression

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle F_{dz}^{\prime }\approx -6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}aV_{m}^{\prime }\,f_{z}(z/H,a/H), & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle f_{z}=1+{\displaystyle \frac{a}{z-a}}+{\displaystyle \frac{a}{H-a-z}}, & \displaystyle\end{eqnarray}$$

where the second and the third terms are corrections to the Stokes drag due to hydrodynamic interactions with two channel walls. In what follows

(3.3) $$\begin{eqnarray}c_{l}=6\unicode[STIX]{x03C0}V_{m}\,f_{z}Re_{p}^{-1}.\end{eqnarray}$$

The second method to calculate the lift (and to check the validity of the first approach) uses the balance of the lift and the gravity forces described by (2.8). By varying the gravity force $F_{g}$ one can, therefore, comprehend the whole range of equilibrium positions within the channel to obtain $c_{l}(z)$ . The advantage of such an approach is that it does not require the particle motion to be quasi-stationary. However, the disadvantage of this method is that the convergence to 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.

4 Results and discussion

In this section, we present the lattice Boltzmann simulation results and compare them with theoretical predictions.

Figure 2.

(a) Dimensionless migration velocity computed as a function of $z/H$ for particles of $a=4\unicode[STIX]{x1D6FF}$ (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 point-like particles. Solid curve shows a polynomial fit to simulation data. (b) Dimensionless slip velocities computed for the same particles (symbols). Solid curve plots 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.1 Neutrally buoyant particles

We start with neutrally buoyant particles and first calculate their migration $V_{m}^{nb}$ and the slip $V_{s}^{nb}$ velocities as a function of $z/H$ . Figure 2 plots simulation data obtained for particles of radius $a=4\unicode[STIX]{x1D6FF}$ . Here we show only a half of the channel since the curves are antisymmetric with respect to the channel axis $z=H/2$ . These results demonstrate that the migration velocity differs significantly from the velocity $c_{l0}Re_{p}/(6\unicode[STIX]{x03C0})$ , where $Re_{p}=a^{2}G_{m}/\unicode[STIX]{x1D708}$ , predicted theoretically for point-like particles (Vasseur & Cox Reference Vasseur and Cox1976). We also see that the equilibrium position, $V_{m}^{nb}=0$ , of finite-size particles is shifted towards the channel axis compared to that of point-like particles, which is obviously due to their interactions with the wall resulting in a finite slip velocity. Indeed, figure 2(b) demonstrates that the computed $V_{s}^{nb}$ grows rapidly near the wall, being close to the theoretical predictions for a linear shear flow near a single wall (Goldman et al. Reference Goldman, Cox and Brenner1967). Unlike theoretical predictions by Goldman et al. (Reference Goldman, Cox and Brenner1967), the computed slip velocity does not vanish in the central part of the channel. Its value is roughly twice larger than the Faxen correction $4U_{m}^{\prime }a^{2}/(3H^{2})$  (Happel & Brenner Reference Happel and Brenner1965). Note that a similar difference has been obtained in simulations of the migration of finite-size particles based on the force coupling method (Loisel et al. Reference Loisel, Abbas, Masbernat and Climent2015). These deviations from the Faxen corrections are likely also caused by hydrodynamic interactions of particles with the wall in a parabolic flow.

Figure 3.

Lift coefficient, $c_{l}$ , for neutrally buoyant particles of $a=4\unicode[STIX]{x1D6FF}$ (circles) and $8\unicode[STIX]{x1D6FF}$ (squares) obtained from the migration velocity at $Re=11.3$ (grey symbols) and $22.6$ (white symbols). Solid and dash-dotted curves show predictions of (2.5) for $a=4\unicode[STIX]{x1D6FF}$ and $8\unicode[STIX]{x1D6FF}$ , dashed curve plots predictions for point-like particles. Vertical dotted lines show $z=a$ . Black symbols show $c_{l}$ obtained for non-neutrally buoyant particles of $a=4\unicode[STIX]{x1D6FF}$ from the force balance at $Re=22.6$ . The inset plots $c_{l}$ in the central part of the channel.

Figure 4.

Equilibrium positions for neutrally buoyant finite-size (grey circles) and point-like (white circle) particles. Dashed curve shows predictions of (2.5).

Figure 3 shows $c_{l}$ for particles of $a=4\unicode[STIX]{x1D6FF}$ and $8\unicode[STIX]{x1D6FF}$ . The lift coefficient has been obtained from the migration velocity and from the force balance as specified above, and simulations have been made for two moderate Reynolds numbers, $Re=11.3$ and $22.6$ . As we discussed above, if $Re\leqslant 20$ a potential dependence of $c_{l}$ on $Re$ could be ruled out a priori, and this is indeed confirmed by our simulations. Therefore, below we provide a detailed comparison of our simulation data with asymptotic solutions obtained for $Re\ll 1$ , which should be applicable for finite moderate $Re$ . Figure 3 also includes theoretical predictions by Vasseur & Cox (Reference Vasseur and Cox1976) and curves calculated with (2.5). One can see that simulation results show a strong discrepancy with the point-particle approximation, especially in the near-wall region, where hydrodynamic interactions are significant. This discrepancy increases with the size of particles. We can, however, conclude that predictions of our equation (2.5) are generally in good agreement with simulation results. Thus, for smaller particles, of $a=4\unicode[STIX]{x1D6FF}$ , equation (2.5) perfectly fits the simulation data in the near-wall region, where the theory for point-like particles fails. Simulation results slightly deviate from predictions of (2.5) near the equilibrium positions and in the central part of the channel. For bigger particles, of $a=8\unicode[STIX]{x1D6FF}$ , these deviations are more pronounced. We emphasize, however, that they are still much smaller than from the point-particle theory by Vasseur & Cox (Reference Vasseur and Cox1976).

To examine a significance of the particle size in more detail, we plot in figure 4(a) the computed equilibrium position, $z_{eq}/H$ , as a function of $a/H$ . We recall that the lift $c_{l}^{nb}(z)$ is antisymmetric with respect to the midplane of the channel axis, so that neutrally buoyant particles have a second equilibrium position at $H-z_{eq}$ . In a point-particle approximation $z_{eq}/H\simeq 0.19$  (Vasseur & Cox Reference Vasseur and Cox1976). We see that for finite-size particles $z_{eq}/H$ is always larger, and increases with the particle size. Note that the increase in $z_{eq}/H$ is nearly linear when $a/H\leqslant 0.1$ . Also included in figure 4 are predictions of (2.5). One can conclude that the theory correctly predicts the trend observed in simulations, but slightly deviates from the simulation data. A possible explanation for this discrepancy could be effects of parabolic flow (which are of the order of $O(a/H)$ ) on the slip velocity and the stresslet (see Yahiaoui & Feuillebois Reference Yahiaoui and Feuillebois2010; Hood et al. Reference Hood, Lee and Roper2015), which are neglected in our theory.

4.2 Non-neutrally buoyant particles

We now turn to the particle migration under both inertial lift and gravity forces.

Figure 5.

Migration velocity of non-neutrally buoyant particles in a horizontal channel. Symbols show simulation data. Solid curve is calculation with (4.1) using data for neutrally buoyant particles.

4.2.1 Horizontal channel

Let us start with the investigation of migration of particles in the most relevant experimental case of a horizontal channel ( $\unicode[STIX]{x1D6FC}=0^{\circ }$ ).

We first fix a weak gravity force, $F_{g}=0.694$ , and compute the migration velocity of particles of radii $a=4\unicode[STIX]{x1D6FF}$ 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. The migration velocity can be calculated as

(4.1) $$\begin{eqnarray}V_{m}=V_{m}^{nb}-V^{St}/f_{z},\end{eqnarray}$$

where we use a fit for $V_{m}^{nb}$ computed for neutrally buoyant particles (see figure 2 a). The agreement between simulation data and calculations using (4.1) is excellent, which confirms that (2.5) remains valid in the case of slightly non-neutrally buoyant particles. We remark that due to gravity $V_{m}$ is shifted downwards relative to $V_{m}^{nb}(z)$ , shown in figure 2. As a result, 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 the wall by varying $F_{g}$ . We have included these force balance results in figure 3 and can conclude that they agree very well with data obtained by using another computational method and for neutrally buoyant particles. This suggests again that above results could be used at moderate Reynolds numbers, $Re\leqslant 20$ , since in this case the lift coefficient does not depend on $Re$ .

Figure 6.

(a) Equilibrium positions of non-neutrally buoyant particles ( $a=4\unicode[STIX]{x1D6FF}$ ) in a horizontal channel; (b) trajectories of the same particles released at $z_{0}=1.125a$ computed at different $F_{g}=0.268$ (squares), $0.804$ (diamonds), $1.340$ (triangles), $2.681$ (circles), $9.383$ (diamonds).

Figure 6(a) shows $z_{eq}/a$ computed at several $F_{g}$ . It can be seen that when the gravity force is getting larger, the equilibrium positions decrease rapidly. This trend can be used to separate particles even when $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ is very small. To illustrate this we now fix $Re=11.3$ , inject particles of $a=4\unicode[STIX]{x1D6FF}$ close to the bottom of the channel, $z_{0}=1.125a$ 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\unicode[STIX]{x1D708}$ . The data show 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 any $Re\leqslant 20$ (see appendix B). Therefore, even in the case of very small $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ , one can always tune the value of $Re$ to induce the required separation difference in $F_{g}$ . For example, we have to separate particles of $a=2~\unicode[STIX]{x1D707}$ m and different $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ in a channel of $H=40~\unicode[STIX]{x1D707}$ m. If we chose $Re=0.3$ , the separation length $L=50xG_{m}a\unicode[STIX]{x1D708}$ of figure 6(b) will be approximately $3.3$  cm. By evaluating $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ with (2.6), we can immediately conclude that trajectories plotted in figure 6(b) from top to bottom correspond to $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=0.007$ , $0.022$ , $0.037$ and $0.073$ , which are indeed extremely small.

4.2.2 Inclined channel

When $F_{g}$ is large enough, it can also influence the slip velocity, and therefore, change the lift itself. This effect is especially important for vertical channels. Note that due to the linearity of the Stokes equations, which govern a disturbance flow at small particle Reynolds numbers, we can decouple the contributions of the particle–wall interaction and of the gravity force into the slip velocity:

(4.2) $$\begin{eqnarray}V_{s}=V_{s}^{nb}+\unicode[STIX]{x0394}V_{s}\sin \unicode[STIX]{x1D6FC},\end{eqnarray}$$

where $\unicode[STIX]{x0394}V_{s}=V^{St}/f_{x}$ is the gravity-induced slip velocity for a vertical channel ( $\unicode[STIX]{x1D6FC}=90^{\circ }$ ) 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=4\unicode[STIX]{x1D6FF}$ in a vertical channel computed by using several values of $F_{g}$ are shown in figures 7(a) and 8(a). Note that the slip velocity, $V_{s}$ , grows with $F_{g}$ since the Stokes velocity, $V^{St}$ , is linearly proportional to $F_{g}$ (see (2.7)). We now use simulation data presented in figures 2(a) and 7(a) to compute $\unicode[STIX]{x0394}V_{s}$ , and then $\unicode[STIX]{x0394}V_{s}/F_{g}$ . The results for $\unicode[STIX]{x0394}V_{s}/F_{g}$ are shown in figure 7(b), and we see that all data collapse into a single curve, which confirms the validity of (4.2). Figure 7(b) also shows that $\unicode[STIX]{x0394}V_{s}/F_{g}$ is nearly constant in the central region of the channel, being smaller than $V^{St}$ , but the deviations from $V^{St}$ grow when particles approach the wall. These results again illustrate that hydrodynamic interactions with the walls significantly affect the motion of particles in the channel.

Figure 7.

Slip velocities (a) and $\unicode[STIX]{x0394}V_{s}/F_{g}$ (b) computed for non-neutrally buoyant particles of $a=4\unicode[STIX]{x1D6FF}$ in a vertical channel. The data sets correspond to $F_{g}=3.475$ (circles), $6.956$ (triangles), $10.44$ (squares) and $13.91$ (diamonds). Dashed line shows $V^{St}/F_{g}$ , vertical dotted lines plot $z=a$ .

Figure 8.

Migration velocities (a) and $\unicode[STIX]{x0394}V_{m}/F_{g}$ (b) computed for non-neutrally buoyant particles of $a=4\unicode[STIX]{x1D6FF}$ in a vertical channel. The data sets correspond to $F_{g}=3.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.

We recall that the variation of the slip velocity caused by gravity is small for slightly non-neutrally buoyant particles (see figure 7), so that (2.5) can be linearized with respect to $\unicode[STIX]{x0394}V_{s}$ :

(4.3) $$\begin{eqnarray}c_{l}\simeq c_{l}^{nb}+\unicode[STIX]{x0394}V_{s}\frac{\unicode[STIX]{x2202}c_{l}(V_{s}^{nb})}{\unicode[STIX]{x2202}V_{s}},\end{eqnarray}$$

where $c_{l}^{nb}=c_{l}(V_{s}^{nb})$ is the lift coefficient for neutrally buoyant particles. By using (3.3) we can then calculate the migration velocity

(4.4) $$\begin{eqnarray}V_{m}=V_{m}^{nb}+\unicode[STIX]{x0394}V_{m}=V_{m}^{nb}+\unicode[STIX]{x0394}V_{s}\frac{\unicode[STIX]{x2202}c_{l}(V_{s}^{nb})}{\unicode[STIX]{x2202}V_{s}}\frac{Re_{p}}{6\unicode[STIX]{x03C0}f_{z}}.\end{eqnarray}$$

The computed migration velocity is shown in figure 8(a). We see that it decreases with $F_{g}$ , and the equilibrium position shifts towards the wall, since $\unicode[STIX]{x0394}V_{s}/F_{g}$ is positive while $\unicode[STIX]{x2202}c_{l}/\unicode[STIX]{x2202}V_{s}$ is negative.

We can now evaluate $\unicode[STIX]{x0394}V_{m}/F_{g}$ by using the simulation data presented in figures 2 and 8(a), and these 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).

Finally, we briefly discuss the case of an arbitrary inclination angle $\unicode[STIX]{x1D6FC}$ , where the $z$ -component of the force can be written as

(4.5) $$\begin{eqnarray}F_{z}=c_{l}(V_{s})+F_{g}\cos \unicode[STIX]{x1D6FC}.\end{eqnarray}$$

By using (4.2), (4.4) and (4.5), we can express the migration velocity as

(4.6) $$\begin{eqnarray}V_{m}=V_{m}^{nb}+\unicode[STIX]{x0394}V_{m}\sin \unicode[STIX]{x1D6FC}+F_{g}\cos \unicode[STIX]{x1D6FC}\frac{Re_{p}}{6\unicode[STIX]{x03C0}f_{z}},\end{eqnarray}$$

where $\unicode[STIX]{x0394}V_{m}$ is evaluated for a vertical channel (see figure 8 b). The equilibrium positions can be found by using the condition $V_{m}=0$ , where $V_{m}$ is calculated with (4.6). The results of these calculations, made at a fixed $F_{g}=3.475$ and different $\unicode[STIX]{x1D6FC}$ , are plotted in figure 9 together with direct simulation data, and one can see that they practically coincide. Our results show that in a vertical channel two stable equilibrium positions coexist. They are symmetric relative to the midplane and 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 Reference Vasseur and Cox1976; Asmolov Reference Asmolov1999). If we slightly reduce $\unicode[STIX]{x1D6FC}$ both stable equilibrium positions become shifted towards the lower wall due to gravity, as seen in figure 9. These two positions coexist only for $\unicode[STIX]{x1D6FC}\geqslant 85.7^{\circ }$ . On reducing $\unicode[STIX]{x1D6FC}$ further the upper equilibrium position disappears, and only one, a lower, equilibrium position remains. This obviously indicates that the inertial lift cannot balance gravity anymore. We note that this remaining single equilibrium position becomes insensitive to the inclination angle when $\unicode[STIX]{x1D6FC}\leqslant 60^{\circ }$ .

Figure 9.

Equilibrium positions $z_{eq}/H$ for $a=4\unicode[STIX]{x1D6FF}$ and $F_{g}=3.475$ . Circles show simulation data. Solid curve plots results obtained using $V_{m}=0$ , where $V_{m}$ is calculated with (4.6).

5 Concluding remarks

In this paper we have studied the inertial migration of finite-size particles in a plane channel flow at moderate Reynolds numbers, $Re\leqslant 20$ . We have shown that the slip velocity, $V_{s}$ , which is finite even for neutrally buoyant particles, contributes to the lift and determines the equilibrium positions in the channel. We have proposed an expression for the lift which generalizes theories, originally applied for some cases of limited guidance, to finite-size particles in a channel flow. When the size of the particle is zero, our formula recovers the known expression of the point-particle approximation (Vasseur & Cox Reference Vasseur and Cox1976). For particles close to the walls we recover earlier predictions for finite-size particles in a linear shear flow (Cherukat & McLaughlin Reference Cherukat and McLaughlin1994). Our theoretical model, which is probably the simplest realistic model for lift in a channel that one might contemplate, provides considerable insight into inertial migration of finite-size particles in microchannels. In particular, it provides a simple explanation of a significant increase in the lift near walls. It also allows one to predict the number of equilibrium positions and determine their locations in various situations.

To check the validity of our theory, we have employed lattice Boltzmann simulations. Generally, the simulation results have fully confirmed the theory, and have shown that many of our theoretical results have validity beyond the initial restrictions of our model. Thus, it has been confirmed that the predictions of our theory do not depend on Reynolds number when $Re\leqslant 20$ , that equilibrium positions of heavy particles in a horizontal channel can be accurately determined by using data for the neutrally buoyant case, and more.

Several of our theoretical predictions could be tested by experiment. In particular, we have shown that particles with a very small density contrast should follow divergent trajectories, so that channel flows with low Reynolds numbers $Re\sim 1$ can be used to separate such particles. We stress that our theory should correctly predict the lift in near-wall regions also in pipes or square channels, and we expect that for this geometry it could be accurate even at $Re\geqslant 20$ since the length scale of the disturbance flow would be the distance to the wall rather than the channel width. For this reason it would be possible to neglect the effects of other distant walls and parabolic flow on the lift. Note, however, that these effects should be taken into account in the central part of the channel.

Our model and computational approach can be extended to more complex situations, which include, for example, hydrophobic walls or particles allowing hydrodynamic slip at their surfaces (Vinogradova Reference Vinogradova1999; Neto et al. Reference Neto, Evans, Bonaccurso, Butt and Craig2005). In this case the hydrodynamic interaction in the near-wall region changes significantly (Davis, Kezirian & Brenner Reference Davis, Kezirian and Brenner1994; Vinogradova Reference Vinogradova1996), so that we expect that the lift force can also be dramatically modified. It would also be interesting to consider the case of an anisotropic superhydrophobic wall, which could induce secondary flows transverse to the direction of applied pressure gradient (Feuillebois, Bazant & Vinogradova Reference Feuillebois, Bazant and Vinogradova2010; Vinogradova & Belyaev Reference Vinogradova and Belyaev2011; Schmieschek et al. Reference Schmieschek, Belyaev, Harting and Vinogradova2012). It has been recently shown (Pimponi et al. Reference Pimponi, Chinappi, Gualtieri and Casciola2014; Asmolov et al. Reference Asmolov, Dubov, Nizkaya, Kuehne and Vinogradova2015) that particles translating in a superhydrophobic channel can be laterally displaced due to such a transverse flow. The use of this effect in combination with the inertial migration should be a fruitful direction, which could allow us to separate particles of different size or density contrast not only by their vertical but also by their transverse positions.