2. Equations of motion

We model flow through an infinitely long square channel of side length $\ell$ . A 3D Poiseuille flow $\bar{\boldsymbol{u}}^{\prime }$ flowing in the $z^{\prime }$ -direction is disturbed by a rigid sphere of radius $a$ (figure 1 a). Here we use primes to denote dimensional variables. We denote the fluid viscosity by ${\it\mu}$ , fluid density by ${\it\rho}$ and the centre-line velocity of the background flow by $U_{m}$ . The particle is located at $(x_{0}^{\prime },y_{0}^{\prime },0)$ and is allowed to translate in the $z^{\prime }$ -direction with velocity $\boldsymbol{U}_{\!p}^{\prime }=U_{p}^{\prime }\boldsymbol{e}_{\boldsymbol{z}^{\prime }}$ , and rotate with angular velocity ${\it\bf\Omega}_{p}^{\prime }$ , until it is drag-free and torque-free. The objective of this paper is to calculate the lift forces acting on the particle in the $x^{\prime }$ - and $y^{\prime }$ -directions.

There are three important dimensionless parameters: (i) the dimensionless ratio of particle radius to channel diameter ${\it\alpha}=a/\ell$ , (ii) the channel Reynolds number $\mathit{Re}=U_{m}\ell /{\it\nu}$ , and (iii) the particle Reynolds number $\mathit{Re}_{p}=U_{m}a^{2}/\ell {\it\nu}$ . Here we write ${\it\nu}={\it\mu}/{\it\rho}$ for the kinematic viscosity. In common with previous theory (Cox & Brenner Reference Cox and Brenner1968; Ho & Leal Reference Ho and Leal1974; Schonberg & Hinch Reference Schonberg and Hinch1989), we will perform dual perturbation expansions in $\mathit{Re}_{p}$ and ${\it\alpha}$ , assuming that both quantities are asymptotically small. In inertial microfluidic experiments (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009), particle diameters may be comparable with the channel dimensions. We will show that our expansions converge even at the moderate values of ${\it\alpha}$ accessed in these experiments.

The background flow, $\bar{\boldsymbol{u}}^{\prime }$ , is square channel Poiseuille flow (Papanastasiou, Georgiou & Alexandrou Reference Papanastasiou, Georgiou and Alexandrou1999), and takes the form $\bar{\boldsymbol{u}}^{\prime }=\bar{u}^{\prime }(x^{\prime },y^{\prime })\boldsymbol{e}_{z^{\prime }}$ , where $\bar{u}^{\prime }$ is defined by

(2.1) $$\begin{eqnarray}\displaystyle \bar{u}^{\prime }(x^{\prime },y^{\prime }) & = & \displaystyle U_{m}\left[\vphantom{\displaystyle \frac{-4\ell ^{2}(-1)^{n}\cosh \left(\displaystyle \frac{(2n+1){\rm\pi}ax^{\prime }}{\ell }\right)}{(2n+1)^{3}{\rm\pi}^{3}a^{2}\cosh \left(\displaystyle \frac{(2n+1){\rm\pi}}{2}\right)}}-\frac{1}{2}\left(y^{^{\prime }2}-\left(\frac{\ell }{2a}\right)^{2}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\mathop{\sum }_{n=0}^{\infty }\frac{-4\ell ^{2}(-1)^{n}\cosh \left(\displaystyle \frac{(2n+1){\rm\pi}ax^{\prime }}{\ell }\right)}{(2n+1)^{3}{\rm\pi}^{3}a^{2}\cosh \left(\displaystyle \frac{(2n+1){\rm\pi}}{2}\right)}\cos \left(\frac{(2n+1){\rm\pi}ay^{\prime }}{\ell }\right)\right]\!.\end{eqnarray}$$

The velocity $\bar{\boldsymbol{u}}^{\prime }$ and pressure $\bar{p}^{\prime }$ solve the Stokes equations with boundary condition $\bar{\boldsymbol{u}}^{\prime }=\mathbf{0}$ on the channel walls. We will also need the Taylor series expansion for $\bar{u}^{\prime }$ around the centre of the particle:

(2.2) $$\begin{eqnarray}\displaystyle \bar{u}^{\prime }(x^{\prime },y^{\prime }) & = & \displaystyle {\it\beta}^{\prime }+{\it\gamma}_{x}^{\prime }(x^{\prime }-x_{0}^{\prime })+{\it\gamma}_{y}^{\prime }(y^{\prime }-y_{0}^{\prime })+{\it\delta}_{xx}^{\prime }(x^{\prime }-x_{0}^{\prime })^{2}+{\it\delta}_{xy}^{\prime }(x^{\prime }-x_{0}^{\prime })(y^{\prime }-y_{0}^{\prime })\nonumber\\ \displaystyle & & \displaystyle +\,{\it\delta}_{yy}^{\prime }(y^{\prime }-y_{0}^{\prime })^{2}+O(r^{^{\prime }3}).\end{eqnarray}$$

To illustrate the reference frame of the equations we will use later, we first list the dimensionless equations of motion and boundary conditions for the velocity and pressure fields $\boldsymbol{u}^{\prime \prime }$ and $p^{\prime \prime }$ expressed in particle-fixed coordinates. We non-dimensionalize these equations by scaling velocities by $U_{m}a/\ell$ , lengths by $a$ and pressures by ${\it\mu}U_{m}/\ell$ :

(2.3a ) $$\begin{eqnarray}\displaystyle {\rm\nabla}^{2}\boldsymbol{u}^{\prime \prime }-\boldsymbol{{\rm\nabla}}p^{\prime \prime } & = & \displaystyle \mathit{Re}_{p}[(\boldsymbol{u}^{\prime \prime }+\boldsymbol{U}_{\!p})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\prime \prime }],\end{eqnarray}$$

(2.3b ) $$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\prime \prime } & = & \displaystyle 0,\end{eqnarray}$$

(2.3c ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{\prime \prime } & = & \displaystyle {\it\bf\Omega}_{p}\times \boldsymbol{r}^{\prime \prime }\quad \text{on }r^{\prime \prime }=1,\end{eqnarray}$$

(2.3d ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{\prime \prime } & = & \displaystyle -\boldsymbol{U}_{\!p}\quad \text{on the walls},\end{eqnarray}$$

(2.3e ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{\prime \prime } & = & \displaystyle \bar{\boldsymbol{u}}-\boldsymbol{U}_{\!p}\quad \text{as }z^{\prime \prime }\rightarrow \pm \infty .\end{eqnarray}$$

Now we introduce the disturbance velocity and pressure fields $\boldsymbol{u}=\boldsymbol{u}^{\prime \prime }-\bar{\boldsymbol{u}}+\boldsymbol{U}_{\!p}$ and $p=p^{\prime \prime }-\bar{p}$ , in which the background flow $\bar{\boldsymbol{u}}-\boldsymbol{U}_{\!p}$ (as measured in this reference frame) is subtracted from $\boldsymbol{u}^{\prime \prime }$ . For reference, the fluid velocity in the laboratory frame is given by $\boldsymbol{v}=\boldsymbol{u}+\bar{\boldsymbol{u}}$ . We then obtain the equations of motion and boundary conditions that will be used throughout this paper:

(2.4a ) $$\begin{eqnarray}\displaystyle {\rm\nabla}^{2}\boldsymbol{u}-\boldsymbol{{\rm\nabla}}p & = & \displaystyle \mathit{Re}_{p}(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{u}}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}),\end{eqnarray}$$

(2.4b ) $$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u} & = & \displaystyle 0,\end{eqnarray}$$

(2.4c ) $$\begin{eqnarray}\displaystyle \boldsymbol{u} & = & \displaystyle {\it\bf\Omega}_{p}\times \boldsymbol{r}-\bar{\boldsymbol{u}}+\boldsymbol{U}_{\!p}\quad \text{on }r=1,\end{eqnarray}$$

(2.4d ) $$\begin{eqnarray}\displaystyle \boldsymbol{u} & = & \displaystyle \mathbf{0}\quad \text{on the walls},\end{eqnarray}$$

(2.4e ) $$\begin{eqnarray}\displaystyle \boldsymbol{u} & = & \displaystyle \mathbf{0}\quad \text{as }z\rightarrow \pm \infty .\end{eqnarray}$$

We call the variables that appear in (2.4) the inner variables. Appendix A summarizes the notations used for dimensionless and dimensional variables.

We formulated (2.4) as a finite element model (FEM) with approximately $650\,000$ linear tetrahedral elements, and solved for $\boldsymbol{u}$ and $p$ using Comsol Multiphysics (COMSOL, Los Angeles) in a rectangular domain with dimensions $\ell /a\times \ell /a\times 5(\ell /a)$ , prescribing $\boldsymbol{u}$ at the inlet $z=-5(\ell /a)$ , and imposing neutral boundary conditions (vanishing stress) at the outlet $z=+5(\ell /a)$ . In the FEM, we vary $U_{p}$ and ${\it\bf\Omega}_{p}$ until there is no drag force or torque on the particle. The FEM Lagrange multipliers, which enforce the velocity boundary condition on the particle, are used to compute the lift force $F_{L}$ on the drag-free and force-free particle. Bramble (Reference Bramble1981) rigorously demonstrates the accuracy of flux calculations from Lagrange multipliers for a Poisson’s equation with Dirichlet boundary conditions. Additionally, we discuss accuracy tests of the FEM discretization for our problem in appendix B.

First, we consider the lift force for particles located on the line of symmetry $x_{0}=0$ . Fixing particle position $y_{0}$ , we found that curves of lift force $F_{L}$ against particle size $a$ collapsed for different Reynolds numbers. Particles in different positions have different apparent scalings for $F_{L}$ as a function of $a$ (figure 5 a,b). By assaying a large range of particle sizes ${\it\alpha}$ , we see that the empirical fit $F_{L}\sim {\it\rho}U^{2}a^{3}$ observed by Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009) is not asymptotic as $a\rightarrow 0$ . The data for smallest particle sizes ( ${\it\alpha}<0.07$ ) are consistent with a scaling law of $F_{L}\sim {\it\rho}U^{2}a^{4}/\ell ^{2}$ as predicted by Ho & Leal (Reference Ho and Leal1974) and Schonberg & Hinch (Reference Schonberg and Hinch1989), but extrapolation of the asymptotic force law to the moderate particle sizes used in real inertial microfluidic devices ( ${\it\alpha}\approx 0.1{-}0.3$ ) over-predicts the lift force by more than an order of magnitude.

3. Dominant balances in the equations of motion

The governing equation (2.4) is a balance between momentum flux and the pressure and viscous stresses. Testing the hypothesis that two of these three contributions might form a dominant balance within the equation, we plotted the resultants of the three fluxes as functions of distance from the particle. Specifically, we integrate the $\ell _{2}$ norm of each flux over spherical control surfaces centred at the particle. Let $S_{r}$ be the boundary of a sphere of radius $r$ centred at the origin, and define the $\ell _{2}$ norm by $\Vert \boldsymbol{u}\Vert _{2}=\sqrt{u^{2}+v^{2}+w^{2}}$ . Then the dimensionless viscous stress resultant acting on the sphere $S_{r}$ is defined by

(3.1) $$\begin{eqnarray}V(r)=\int _{S_{r}}\Vert \boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\Vert _{2}\,\text{d}s,\end{eqnarray}$$

and the dimensionless inertial term $I(r)$ stress resultant by

(3.2) $$\begin{eqnarray}I(r)=\mathit{Re}_{p}\int _{S_{r}}\Vert [(\bar{\boldsymbol{u}}-\boldsymbol{U}_{\!p})\boldsymbol{u}+\boldsymbol{u}(\bar{\boldsymbol{u}}-\boldsymbol{U}_{\!p})+\boldsymbol{u}\boldsymbol{u}]\boldsymbol{\cdot }\boldsymbol{n}\Vert _{2}\,\text{d}s.\end{eqnarray}$$

The integrand in $I(r)$ is chosen to have divergence equal to the right-hand side of (2.4), and we pick a form of the inertial flux that decays in $\ell _{2}$ norm as $r\rightarrow \infty$ .

Numerically evaluating these two terms as well as $(1/r)\int _{S_{r}}\Vert p\boldsymbol{n}\Vert _{2}\,\text{d}s$ , we find that, contrary to the predictions of Ho & Leal (Reference Ho and Leal1974) and Schonberg & Hinch (Reference Schonberg and Hinch1989), at moderate channel Reynolds numbers, the viscous and pressure stress resultants are numerically larger than the momentum flux. In particular, there is no region in which $V(r)$ and $I(r)$ are co-dominant at $\mathit{Re}=10$ (figure 2). Indeed, even at higher Reynolds numbers ( $\mathit{Re}=50$ , 80) for which inertial stresses are numerically larger than viscous stresses, inertial stresses can be collapsed onto a single curve (see insets of figure 2 a,b) by rescaling with $\mathit{Re}$ . This scaling suggests that the underlying dynamics, even at moderate values of $\mathit{Re}$ , are inherited from the small- $\mathit{Re}$ dominant balance of pressure and viscous stresses. Dominance of viscous stresses over inertial stresses is surprising because, as Ho & Leal (Reference Ho and Leal1974) noticed, the resulting dominant balance equations are not self-consistent for isolated particles in unbounded fluid flow.

Figure 2. The dominant balance of the NSE for (a) a particle near the channel centre ( $y_{0}=0.15/{\it\alpha}$ and ${\it\alpha}=0.11$ ) and (b) a particle near the channel walls ( $y_{0}=0.35/{\it\alpha}$ and ${\it\alpha}=0.06$ ). The viscous stresses $V(r)$ for various Reynolds numbers are plotted as thin (red) lines, and the inertial stresses $I(r)$ are plotted as thick (black) lines. Reynolds numbers are indicated by line style: $\mathit{Re}=10$ , solid line; $\mathit{Re}=50$ , dashed line; and $\mathit{Re}=80$ , dotted line. The insets show that the inertial stresses $I(r)$ collapse when scaled by $\mathit{Re}$ , suggesting that the high-Reynolds-number dynamics are determined by the low-Reynolds-number dynamics.

We will now present a first-order estimate of the size of the domain in which inertial stresses may be expected to be dominant. The slowest-decaying component of the disturbance flow associated with a force-free particle on the plane of symmetry ( $x_{0}=0$ ) is given by the stresslet flow (Batchelor Reference Batchelor1967; Kim & Karrila Reference Kim and Karrila2005):

(3.3) $$\begin{eqnarray}\boldsymbol{u}_{\mathit{stresslet}}=\frac{5{\it\gamma}_{y}(y-y_{0})z\boldsymbol{r}}{2r^{5}}=O\left(\frac{1}{r^{2}}\right).\end{eqnarray}$$

Recall that ${\it\gamma}_{y}$ is the strain rate, defined in (2.2). For this flow field, the viscous stress term in (2.4) decays with distance like

(3.4) $$\begin{eqnarray}V(r)\sim O(\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{\mathit{stresslet}})\sim O\left(\frac{1}{r^{3}}\right),\end{eqnarray}$$

whereas the inertial stresses vary with distance like

(3.5) $$\begin{eqnarray}I(r)\sim O(\mathit{Re}_{p}(\boldsymbol{u}_{\mathit{stresslet}})(\bar{\boldsymbol{u}}-\boldsymbol{U}_{\!p}))\sim O\left(\frac{\mathit{Re}_{p}}{r}\right).\end{eqnarray}$$

We define the cross-over radius, $r_{\ast }$ , to be the distance at which the viscous and inertial stresses are comparable,

(3.6) $$\begin{eqnarray}r_{\ast }=O\left(\frac{1}{\mathit{Re}_{p}^{1/2}}\right).\end{eqnarray}$$

In order to compare the cross-over radius to the width of the channel, we consider when ${\it\alpha}r_{\ast }=O(1/\mathit{Re}^{1/2})$ is equal to one. To ensure that viscous stresses dominate over inertial stresses over the channel cross-section (i.e.  ${\it\alpha}r_{\ast }\gg 1$ ), Ho & Leal (Reference Ho and Leal1974) restrict to cases where $\mathit{Re}\ll 1$ . The asymptotic analysis of Schonberg & Hinch (Reference Schonberg and Hinch1989) allows that $\mathit{Re}=O(1)$ , but at the cost of needing to separately model and match the flows at $O(1)$ distances from the particle where viscous stresses are dominant, and at $O(1/\mathit{Re}^{1/2})$ distances where inertial and viscous stresses must both be included in the dominant balance.

However, the predicted cross-over radius falls short of the numerical cross-over radius (figure 2, table 2). There are two explanations for the dominance of viscous stresses over inertial in these experimental geometries. First, the above estimates do not consider the coefficients in the stresslet, merely the order of magnitude of the terms. Second, although the stresslet describes the flow disturbance for a force-free particle in an unbounded fluid, the leading-order flow is considerably altered by the presence of the channel walls. Below, we demonstrate that both explanations contribute to the dominance of viscous stresses throughout the channel cross-section, pushing the cross-over radius $r^{\ast }$ out beyond the channel walls (table 2).

Table 2. The cross-over radius ${\it\alpha}r_{\ast }$ at which $I(r)\geqslant V(r)$ computed for $\mathit{Re}=10$ , 50, 80 using the following methods: (i) Ho & Leal’s (Reference Ho and Leal1974) calculation using the stresslet, (ii) our calculation using the stresslet, (iii) our calculation using the stresslet and first wall correction, and (iv) our calculation using the numerical solution to the full NSE.

3.1. Role of the stresslet constants

We compute $I(r)$ and $V(r)$ numerically for the stresslet flow field (i.e. substitute $\boldsymbol{u}=\boldsymbol{u}_{\mathit{stresslet}}$ in (3.1) and (3.2)). We examine two representative cases: a medium-sized particle near the channel centre ( $y_{0}=0.15/{\it\alpha}$ , ${\it\alpha}=0.11$ ) (figure 3 a), and a small particle near the channel wall ( $y=0.35/{\it\alpha}$ , ${\it\alpha}=0.06$ ) (figure 3 b). For $\mathit{Re}=10$ , in both cases the inertia is significantly smaller than the viscous stress throughout the channel. At larger values of $\mathit{Re}$ , $I(r)$ eventually exceeds $V(r)$ , but the cross-over radius $r_{\ast }$ is much larger than simple order-of-magnitude estimates would suggest (table 2).

Figure 3. The dominant balance that arises from the stresslet approximation of the flow for (a) a particle near the channel centre ( $y_{0}=0.15/{\it\alpha}$ and ${\it\alpha}=0.11$ ) and (b) a particle near the channel walls ( $y_{0}=0.35/{\it\alpha}$ and ${\it\alpha}=0.06$ ). The viscous stresses $V(r)$ for various Reynolds numbers are plotted as thin (red) lines, and the inertial stresses $I(r)$ are plotted as thick (black) lines. Reynolds numbers are indicated by line style: $\mathit{Re}=10$ , solid line; $\mathit{Re}=50$ , dashed line; and $\mathit{Re}=80$ , dotted line.

3.2. Role of wall effects

To estimate how wall modifications of the disturbance flow affect the dominant balances in (2.4), we numerically computed the first wall correction. That is, we substitute $\boldsymbol{u}=\boldsymbol{u}_{\mathit{stresslet}}+\boldsymbol{u}_{\mathit{image}}$ into (3.1) and (3.2), where $\boldsymbol{u}_{\mathit{image}}$ is a solution of Stokes’ equations with boundary condition $\boldsymbol{u}_{\mathit{image}}=-\boldsymbol{u}_{\mathit{stresslet}}$ on the channel walls. We examine the same two representative cases as in § 3.1: $(y_{0}=0.15/{\it\alpha},{\it\alpha}=0.11)$ and $(y_{0}=0.35/{\it\alpha},{\it\alpha}=0.06)$ (figure 4 a,b). For $\mathit{Re}=10$ , in both cases the inertia is significantly smaller than the viscous stress throughout the channel. At larger values of $\mathit{Re}$ , $I(r)$ eventually exceeds $V(r)$ , but the cross-over radius $r_{\ast }$ is larger than that predicted from the stresslet coefficients (table 2).

Figure 4. The dominant balance that arises from the stresslet and first wall correction of the flow for (a) a particle near the channel centre ( $y_{0}=0.15/{\it\alpha}$ and ${\it\alpha}=0.11$ ) and (b) a particle near the channel walls ( $y_{0}=0.35/{\it\alpha}$ and ${\it\alpha}=0.06$ ). The viscous stresses $V(r)$ for various Reynolds numbers are plotted as thin (red) lines, and the inertial stresses $I(r)$ are plotted as thick (black) lines. Reynolds numbers are indicated by line style: $\mathit{Re}=10$ , solid line; $\mathit{Re}=50$ , dashed line; and $\mathit{Re}=80$ , dotted line.

We can rationalize the larger values of the cross-over radius ${\it\alpha}r_{\ast }$ by considering the boundary conditions on the channel walls. Because the velocity field $\boldsymbol{u}$ vanishes on the channel walls, the inertial stresses vanish there. Therefore $I(r)$ is suppressed at larger radii. We see less suppression of $V(r)$ , presumably because viscous stresses do not need to vanish on the channel walls. Suppression of $I(r)$ increases the cross-over radius at which inertial stresses must be considered in the dominant balance.

4. A series expansion for the inertial lift force

Our careful evaluation of the stresslet prefactors and wall contributions shows that viscous stresses are dominant over inertial stresses over much of the fluid-filled domain, including at much greater distances from the particle than previous estimates have suggested. We therefore develop an asymptotic theory, based on Cox & Brenner (Reference Cox and Brenner1968) and Ho & Leal (Reference Ho and Leal1974), in which the flow field $\boldsymbol{u}$ , pressure $p$ , particle velocity $\boldsymbol{U}_{\!p}$ and rotation ${\it\bf\Omega}_{p}$ are expanded in powers of $\mathit{Re}_{p}$ , with inertia completely neglected in the leading-order equations:

(4.1a,b ) $$\begin{eqnarray}\boldsymbol{u}=\boldsymbol{u}^{(0)}+\mathit{Re}_{p}\boldsymbol{u}^{(1)}+\cdots \!,\quad p=p^{(0)}+\mathit{Re}_{p}p^{(1)}+\cdots \!,\quad \text{etc.}\end{eqnarray}$$

Notice that this is an expansion in the particle Reynolds number $\mathit{Re}_{p}$ and not the channel Reynolds number $\mathit{Re}$ . Although in experiments the channel Reynolds number is typically large, the expansion is formally valid provided that ${\it\alpha}^{2}$ is small enough that $\mathit{Re}_{p}={\it\alpha}^{2}\mathit{Re}\lesssim 1$ . In fact, when we compare our theory with numerical simulations in § 4.5, we find that the perturbative series gives a good approximation to the lift force even for $\mathit{Re}_{p}=7$ (figure 1 b).

First we compute the first two terms in the perturbative series $\boldsymbol{u}^{(0)}+\mathit{Re}_{p}\boldsymbol{u}^{(1)}$ numerically, showing that retaining these two terms gives the lift force quantitatively accurately over the entire dynamical range of experiments.

Series-expanding (2.4) and collecting like terms in $\mathit{Re}_{p}$ we arrive at equations for ( $\boldsymbol{u}^{(0)}$ , $p^{(0)}$ ), the first-order velocity and pressure:

(4.2a ) $$\begin{eqnarray}\displaystyle & {\rm\nabla}^{2}\boldsymbol{u}^{(0)}-\boldsymbol{{\rm\nabla}}p^{(0)}=\mathbf{0},\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{(0)}=0, & \displaystyle\end{eqnarray}$$

(4.2b ) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}^{(0)}=\boldsymbol{U}_{\!p}^{(0)}+{\it\bf\Omega}_{p}^{(0)}\times \boldsymbol{r}-\bar{\boldsymbol{u}}\quad \text{on }r=1, & \displaystyle\end{eqnarray}$$

(4.2c ) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}^{(0)}=\mathbf{0}\quad \text{on channel walls and as }z\rightarrow \pm \infty . & \displaystyle\end{eqnarray}$$

$(\boldsymbol{u}^{(1)},p^{(1)})$

(4.3a ) $$\begin{eqnarray}\displaystyle & {\rm\nabla}^{2}\boldsymbol{u}^{(1)}-\boldsymbol{{\rm\nabla}}p^{(1)}=(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(0)}+\boldsymbol{u}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{u}}+\boldsymbol{u}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(0)}),\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{(1)}=0, & \displaystyle\end{eqnarray}$$

(4.3b ) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}^{(1)}=\boldsymbol{U}_{\!p}^{(1)}+{\it\bf\Omega}_{p}^{(1)}\times \boldsymbol{r}\quad \text{on }r=1, & \displaystyle\end{eqnarray}$$

(4.3c ) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}^{(1)}=\mathbf{0}\quad \text{on channel walls and as }z\rightarrow \pm \infty . & \displaystyle\end{eqnarray}$$

For both cases, we only need to solve the Stokes equations with a known body force term. In (4.2), the body force term is equal to $\mathbf{0}$ ; in (4.3) the body force term is equal to the inertia of flow $\boldsymbol{u}^{(0)}$ .

In fact, we can apply Lorentz’s reciprocal theorem (Leal Reference Leal1980) to calculate the lift force associated with $\boldsymbol{u}^{(1)}$ without needing to directly solve (4.3). We define the test fluid flow $(\hat{\boldsymbol{u}},\hat{p})$ representing Stokes flow around a sphere moving with unit velocity in the $y$ -direction: viz. satisfying (4.2) with the velocity condition on the sphere replaced by $\hat{\boldsymbol{u}}=\boldsymbol{e}_{y}$ . If ${\bf\sigma}^{(1)}$ and $\hat{{\bf\sigma}}$ are the viscous stress tensors associated with the flow fields $(\boldsymbol{u}^{(1)},p^{(1)})$ and $(\hat{\boldsymbol{u}},\hat{p})$ respectively, ${\bf\sigma}^{(1)}=\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(1)}+(\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(1)})^{\text{T}}-p^{(1)}\mathbf{1}$ , etc., and $\hat{\unicode[STIX]{x1D65A}}$ and $\unicode[STIX]{x1D65A}^{(1)}$ the respective rate-of-strain tensors, $\unicode[STIX]{x1D65A}^{(1)}=[\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(1)}+(\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(1)})^{\text{T}}]/2$ , etc., then by the divergence theorem, the following relation is valid for any volume $V$ enclosed by a surface  $S$ :

(4.4) $$\begin{eqnarray}\int _{S}(\boldsymbol{n}\boldsymbol{\cdot }\hat{{\bf\sigma}}\boldsymbol{\cdot }\boldsymbol{u}^{(1)}-\boldsymbol{n}\boldsymbol{\cdot }{\bf\sigma}^{(1)}\boldsymbol{\cdot }\hat{\boldsymbol{u}})\,\text{d}s=\int _{V}[\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\hat{{\bf\sigma}}\boldsymbol{\cdot }\boldsymbol{u}^{(1)})-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\bf\sigma}^{(1)}\boldsymbol{\cdot }\hat{\boldsymbol{u}})]\text{d}v.\end{eqnarray}$$

By setting $V$ equal to the fluid-filled domain and substituting boundary conditions from (2.4), we deduce that

(4.5) $$\begin{eqnarray}\displaystyle & & \displaystyle \boldsymbol{U}_{\!p}^{(1)}\boldsymbol{\cdot }\int _{S}(\hat{{\bf\sigma}}\boldsymbol{\cdot }\boldsymbol{n})\,\text{d}s+\int _{S}({\it\bf\Omega}_{p}^{(1)}\times \boldsymbol{r})\boldsymbol{\cdot }\hat{{\bf\sigma}}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}s-\boldsymbol{e}_{\boldsymbol{y}}\boldsymbol{\cdot }\int _{S}{\bf\sigma}^{(1)}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}s\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{V}[(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\sigma}^{(1)})\boldsymbol{\cdot }\hat{\boldsymbol{u}}+{\bf\sigma}^{(1)}\,\boldsymbol{ :}\,\hat{\unicode[STIX]{x1D65A}}-(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{{\bf\sigma}})\boldsymbol{\cdot }\boldsymbol{u}^{(1)}-\hat{{\bf\sigma}}\,\boldsymbol{ :}\,\unicode[STIX]{x1D65A}^{(1)}]\,\text{d}v.\end{eqnarray}$$

On the left-hand side of the equation, the first term is zero by symmetry. Similarly, the integrand of the second term can be rearranged to

(4.6) $$\begin{eqnarray}({\it\bf\Omega}_{p}^{(1)}\times \boldsymbol{r})\boldsymbol{\cdot }\hat{{\bf\sigma}}\boldsymbol{\cdot }\boldsymbol{n}={\it\bf\Omega}_{p}^{(1)}\boldsymbol{\cdot }(\boldsymbol{r}\times \hat{{\bf\sigma}}\boldsymbol{\cdot }\boldsymbol{n}),\end{eqnarray}$$

which also integrates to zero. On the right-hand side of (4.5), the third term is zero by definition (since $\hat{\boldsymbol{u}}$ solves the Stokes equations). Furthermore, we can rearrange the second and fourth terms to

(4.7) $$\begin{eqnarray}{\bf\sigma}^{(1)}\,\boldsymbol{ :}\,\hat{\unicode[STIX]{x1D65A}}-\hat{{\bf\sigma}}\,\boldsymbol{ :}\,\unicode[STIX]{x1D65A}^{(1)}=2\unicode[STIX]{x1D65A}^{(1)}\,\boldsymbol{ :}\,\hat{\unicode[STIX]{x1D65A}}-p^{(1)}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\hat{\boldsymbol{u}}-2\hat{\unicode[STIX]{x1D65A}}\,\boldsymbol{ :}\,\unicode[STIX]{x1D65A}^{(1)}+\hat{p}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{(1)}=0,\end{eqnarray}$$

since both flows are incompressible. So, on the right-hand side of (4.5), only the first term of the volume integral remains. Using the definitions of ${\bf\sigma}^{(1)}$ and $\hat{{\bf\sigma}}$ , we obtain the following formula, which we refer to as the reciprocal theorem:

(4.8) $$\begin{eqnarray}\boldsymbol{e}_{\boldsymbol{y}}\boldsymbol{\cdot }\boldsymbol{F}_{\boldsymbol{L}}=\int _{V}\hat{\boldsymbol{u}}\boldsymbol{\cdot }(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(0)}+\boldsymbol{u}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{u}}+\boldsymbol{u}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(0)})\,\text{d}v.\end{eqnarray}$$

We have now reduced our calculation of the lift force to that of solving two homogeneous Stokes equations and performing a volume integral. Numerically, we let $V$ be the truncated numerical domain modelled by our FEM. Next we solve numerically for $\boldsymbol{u}^{(0)}$ from (4.2) and $\hat{\boldsymbol{u}}$ . Again, we choose $U_{p}^{(0)}$ and ${\rm\Omega}_{p}^{(0)}$ so that the particle travels force-free and torque-free. We compute the lift force using the reciprocal theorem in (4.8) for particles at two different channel positions (figure 5 a,b). We see close quantitative agreement between the lift force computed from the full NSEs and the lift force computed from the reciprocal theorem using the two-term expansion in $\mathit{Re}_{p}$ . The comparison is accurate even when, as for $y_{0}=0.15/{\it\alpha}$ , there is no simple scaling law for the dependence of $F_{L}$ upon $a$ (figure 5 a). In the next section, we develop a model that nevertheless allows analytic evaluation of the lift force.

Figure 5. Numerical computation of the scaled lift force $F_{L}/{\it\rho}U_{m}^{2}\ell ^{2}$ using the NSEs in (2.4) as a function of particle size ${\it\alpha}$ for various channel Reynolds numbers: $\mathit{Re}=10$ (triangles), $\mathit{Re}=50$ (circles) and $\mathit{Re}=80$ (crosses). The dashed (black) line represents a scaling law with exponent 4, i.e.  $F_{L}\sim {\it\rho}U_{m}^{2}{\it\alpha}^{2}a^{2}$ as in Ho & Leal (Reference Ho and Leal1974), while the dotted (blue) line represents a scaling law with exponent 3, i.e.  $F_{L}\sim {\it\rho}U_{m}^{2}{\it\alpha}a^{2}$ , which is the line of best fit computed in Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009). The solid line represents the regular perturbation expansion computed numerically using the reciprocal theorem in (4.8). We compare all of these force predictions at two locations in the channel: (a) a particle near the channel centre ( $y_{0}=0.15/{\it\alpha}$ and ${\it\alpha}=0.11$ ); and (b) a particle near the channel walls ( $y_{0}=0.35/{\it\alpha}$ and ${\it\alpha}=0.06$ ).

4.1. Approximation of $\boldsymbol{u}^{(0)}$ and $\hat{\boldsymbol{u}}$ by method of images

In the previous section we showed that a single regular perturbation in $\mathit{Re}_{p}$ of Stokes equations agrees excellently with the numerically computed lift force. We calculated the terms in this perturbation series numerically; but to rationally design inertial microfluidic devices, we need an asymptotic theory for how the lift force and the inertial focusing points depend on the size of the particle and its position within the channel. We derive this theory from asymptotic expansion of $\boldsymbol{u}^{(0)}$ and $\hat{\boldsymbol{u}}$ in powers of ${\it\alpha}$ , the dimensionless particle size. We follow Ho & Leal (Reference Ho and Leal1974) and use the method of reflections to generate expansions in powers of ${\it\alpha}$ for the Stokes flow fields appearing in (4.8) (Happel & Brenner Reference Happel and Brenner1982),

(4.9) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{(0)}=\boldsymbol{u}_{1}^{(0)}+\boldsymbol{u}_{2}^{(0)}+\boldsymbol{u}_{3}^{(0)}+\boldsymbol{u}_{4}^{(0)}+\cdots \!, & & \displaystyle\end{eqnarray}$$

with similar expansions for $p$ , $\hat{\boldsymbol{u}}$ and $\hat{p}$ . Here, $\boldsymbol{u}_{1}^{(0)}$ is the Stokes solution for a particle in unbounded flow, $\boldsymbol{u}_{2}^{(0)}$ is the Stokes solution with boundary condition $\boldsymbol{u}_{2}^{(0)}=-\boldsymbol{u}_{1}^{(0)}$ on the channel walls, and $\boldsymbol{u}_{3}^{(0)}$ is the Stokes solution with boundary condition $\boldsymbol{u}_{3}^{(0)}=-\boldsymbol{u}_{2}^{(0)}$ on the particle surface, etc. Odd terms impose the global boundary conditions on the particle, whereas even terms impose the global boundary conditions on the channel walls. We will show below that the terms in this series constitute a power series in  ${\it\alpha}$ .

Since the odd terms in the expansion, $\boldsymbol{u}_{2i-1}^{(0)}$ , are prescribed on the sphere’s surface, they can be calculated using Lamb’s method for solving the flow external to a sphere (Lamb Reference Lamb1945; Happel & Brenner Reference Happel and Brenner1982). This method expands the velocity field as a sum of multipoles located at the sphere centre, namely,

(4.10) $$\begin{eqnarray}\boldsymbol{u}_{2i-1}^{(0)}=\mathop{\sum }_{n=0}^{\infty }\frac{1}{r^{n+1}}\boldsymbol{f}_{n}^{i}\left(\frac{x-x_{0}}{r},\frac{y-y_{0}}{r},\frac{z}{r}\right)\!,\end{eqnarray}$$

where each term $\boldsymbol{f}_{n}^{i}/r^{n+1}$ is a combination of the stokeslet $n$ -pole and the source $(n-1)$ -pole. We can similarly expand the odd terms of $\hat{\boldsymbol{u}}$ as:

(4.11) $$\begin{eqnarray}\hat{\boldsymbol{u}}_{2i-1}=\mathop{\sum }_{n=0}^{\infty }\frac{1}{r^{n+1}}\boldsymbol{g}_{n}^{i}\left(\frac{x-x_{0}}{r},\frac{y-y_{0}}{r},\frac{z}{r}\right)\!.\end{eqnarray}$$

The full analytic forms for the $\boldsymbol{f}_{n}^{1}$ and $\boldsymbol{g}_{n}^{1}$ are listed in appendix C. From the analytic form of $\boldsymbol{u}_{1}^{(0)}$ , we can find $\boldsymbol{u}_{2}^{(0)}$ by solving the associated Stokes problem numerically. Given $-\boldsymbol{u}_{2}^{(0)}$ on the particle surface, we can appeal to Lamb’s solution to find $\boldsymbol{u}_{3}^{(0)}$ , and so on. The same sequence of reflections can be used to expand the reference velocity  $\hat{\boldsymbol{u}}$ .

4.2. Approximation to the reciprocal theorem integral

Given the Stokes velocities $\boldsymbol{u}^{(0)}$ and $\hat{\boldsymbol{u}}$ , we can compute the inertial lift force $F_{L}$ up to terms of $O(\mathit{Re}_{p})$ using the reciprocal theorem (4.8). As in Ho & Leal (Reference Ho and Leal1974), it is advantageous to divide the fluid-filled domain $V$ into two subdomains, $V_{1}$ and $V_{2}$ , where

(4.12a,b ) $$\begin{eqnarray}V_{1}=\{\boldsymbol{r}\in V:r\leqslant {\it\xi}\}\quad \text{and}\quad V_{2}=\{\boldsymbol{r}\in V:r\geqslant {\it\xi}\}.\end{eqnarray}$$

The intermediate radius ${\it\xi}$ is any parameter satisfying $1\ll {\it\xi}\ll 1/{\it\alpha}$ . Call the corresponding integrals the inner integral and the outer integral, and identify their contributions to the lift force as $F_{L_{1}}$ and $F_{L_{2}}$ , respectively (  $F_{L}=F_{L_{1}}+F_{L_{2}}$ ). The division of the integral into inner and outer regions allows one to incorporate varying length scales ( ${\it\alpha}$ for the inner region and $\ell$ for the outer region) into our model. Note that, distinct from Schonberg & Hinch (Reference Schonberg and Hinch1989), inertia remains subdominant even in the outer region $V_{2}$ . In the next two sections, we will separately consider the contributions from the inner and outer integrals.

4.3. The inner integral

For the inner integral, we continue to scale lengths by $a$ , so that $1\leqslant r\leqslant {\it\xi}\ll {\it\alpha}^{-1}$ . The inner integral can be expressed as the following expansion in ${\it\alpha}$ :

(4.13) $$\begin{eqnarray}F_{L_{1}}={\it\rho}U_{m}^{2}a^{2}(h_{4}{\it\alpha}^{2}+h_{5}{\it\alpha}^{3}+\cdots ).\end{eqnarray}$$

In order to calculate the terms $h_{4}$ and $h_{5}$ , we sort the terms of the Stokes velocities by leading order in ${\it\alpha}$ . The terms contributing at $O({\it\alpha}^{2})$ in the inner region are

(4.14a−c ) $$\begin{eqnarray}\boldsymbol{u}_{1}^{(0)}\sim {\it\alpha}\left(\frac{1}{r^{2}}\boldsymbol{f}_{1}^{1}+\frac{1}{r^{4}}\boldsymbol{f}_{3}^{1}\right),\quad \hat{\boldsymbol{u}}_{1}\sim \frac{1}{r}\boldsymbol{g}_{0}^{1}+\frac{1}{r^{3}}\boldsymbol{g}_{2}^{1},\quad \bar{\boldsymbol{u}}\sim {\it\gamma}{\it\alpha}r.\end{eqnarray}$$

All of these terms are known analytically (see appendix C), and it can be shown that their contribution to the inner integral evaluates to zero, i.e.  $h_{4}=0$ .

At $O({\it\alpha}^{3})$ the velocity terms contributing to calculation of $h_{5}$ are

(4.15) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{u}_{1}^{(0)}\sim {\it\alpha}\left(\frac{1}{r^{2}}\boldsymbol{f}_{1}^{1}+\frac{1}{r^{4}}\boldsymbol{f}_{3}^{1}\right)+{\it\alpha}^{2}\left(\frac{1}{r^{3}}\boldsymbol{f}_{2}^{1}+\frac{1}{r^{5}}\boldsymbol{f}_{4}^{1}\right),\quad \bar{\boldsymbol{u}}\sim {\it\gamma}{\it\alpha}r+{\it\delta}{\it\alpha}^{2}r^{2},\\ \displaystyle \hat{\boldsymbol{u}}_{1}\sim \frac{1}{r}\boldsymbol{g}_{0}^{1}+\frac{1}{r^{3}}\boldsymbol{g}_{2}^{1},\quad \hat{\boldsymbol{u}}_{2}\sim {\it\alpha}\mathscr{S}\left[\frac{1}{r}\boldsymbol{g}_{0}^{1}\right]_{0},\quad \hat{\boldsymbol{u}}_{3}\sim {\it\alpha}\left(\frac{1}{r}\boldsymbol{g}_{0}^{3}+\frac{1}{r^{3}}\boldsymbol{g}_{2}^{3}\right)\!,\end{array}\!\right\}\end{eqnarray}$$

where we define $\boldsymbol{v}\equiv \mathscr{S}[\boldsymbol{u}]$ as the image of the function $\boldsymbol{u}$ , and we define $\boldsymbol{v}_{0}\equiv \mathscr{S}[\boldsymbol{u}]_{0}$ as the velocity $\boldsymbol{v}$ evaluated at the particle centre. That is, $\boldsymbol{v}$ solves the Stokes equations with boundary condition $\boldsymbol{v}=-\boldsymbol{u}$ on the channel walls, and $\boldsymbol{v}_{0}=\boldsymbol{v}(x_{0},y_{0},0)$ . We determine $\mathscr{S}[(1/r)\boldsymbol{g}_{0}^{1}]$ numerically, by discretizing Stokes equations as an FEM, with quadratic elements for the velocity field and linear elements for the pressure field, and solving the FEM in Comsol Multiphysics.

The $O({\it\alpha}^{3})$ contribution to the inner integral is

(4.16) $$\begin{eqnarray}\displaystyle h_{5} & = & \displaystyle \int _{\mathbb{R}^{3}}(\hat{\boldsymbol{u}}_{1}+\hat{\boldsymbol{u}}_{2}+\hat{\boldsymbol{u}}_{3})\boldsymbol{\cdot }(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{1}^{(0)}+\boldsymbol{u}_{1}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{u}}+\boldsymbol{u}_{1}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{1}^{(0)})\text{d}v\nonumber\\ \displaystyle & = & \displaystyle \int _{\mathbb{R}^{3}}\hat{\boldsymbol{u}}_{1}\boldsymbol{\cdot }(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{1}^{(0)}+\boldsymbol{u}_{1}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{u}}+\boldsymbol{u}_{1}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{1}^{(0)})\,\text{d}v\end{eqnarray}$$

where we have made use of the fact that the contributions to the integral from $\hat{\boldsymbol{u}}_{2}$ and $\hat{\boldsymbol{u}}_{3}$ evaluate to zero. Since all of the terms in the integrand are $O(r^{3})$ as $r\rightarrow \infty$ , we can take ${\it\xi}\rightarrow \infty$ ; viz., replace integration over $V_{1}$ by integration over $\mathbb{R}^{3}$ . In doing so, we pick up an error that is $O(1/{\it\xi})$ . We neglect this contribution, since ${\it\xi}\gg 1$ ; in fact the error terms can be shown to cancel with corresponding contributions from the outer integral if expansions are continued to higher-order powers of ${\it\alpha}$ . Evaluating the final integral, we obtain

(4.17) $$\begin{eqnarray}F_{L_{1}}=\frac{{\it\rho}U_{m}^{2}h_{5}a^{5}}{\ell ^{3}}+O(a^{6}),\end{eqnarray}$$

where

(4.18) $$\begin{eqnarray}h_{5}=-\frac{26\,171{\rm\pi}{\it\gamma}_{y}^{2}}{277\,200}-\frac{53{\rm\pi}{\it\gamma}_{y}{\it\delta}_{xx}}{1728}-\frac{283{\rm\pi}{\it\gamma}_{y}{\it\delta}_{yy}}{3150}\end{eqnarray}$$

is $O(1)$ , and depends only on the location of the particle. Recall that the constants ${\it\gamma}_{y}$ , ${\it\delta}_{xx}$ and ${\it\delta}_{yy}$ were defined in the expansion of $\bar{u}$ in (2.2), and depend on the particle position.

4.4. The outer integral

For the outer integral, we will consider alternative dimensionless variables, by using the rescaled distance $\boldsymbol{R}={\it\alpha}\boldsymbol{r}$ . This corresponds to using $\ell$ to non-dimensionalize lengths, rather than $a$ . We call these variables the outer variables, and we will denote them with upper-case roman letters. A detailed comparison of the dimensionless variables is given in appendix A.

In the outer region $V_{2}$ , we must express our functions in terms of $\boldsymbol{R}$ and rearrange our functions by order of magnitude in ${\it\alpha}$ . These expansions are listed in full in appendix D. In the outer region, the reciprocal theorem integral takes the dimensional form

(4.19) $$\begin{eqnarray}F_{L_{2}}={\it\rho}U_{m}^{2}\ell ^{2}\int _{V_{C}}\hat{\boldsymbol{U}}\boldsymbol{\cdot }(\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}^{(0)}+\boldsymbol{U}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{U}}+\boldsymbol{U}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}^{(0)})\,\text{d}v,\end{eqnarray}$$

where we have expanded our domain of integration from $V_{2}=\{\boldsymbol{R}\in V:R\geqslant {\it\xi}\}$ to the entire empty channel $V_{C}$ . This expansion of the domain is justified since the contribution from the region that we add to the integral $\{\boldsymbol{R}:0\leqslant R\leqslant {\it\alpha}{\it\xi}\}$ is $O({\it\alpha}^{4}{\it\xi})$ , and ${\it\xi}\ll 1/{\it\alpha}$ . In fact, this residue (which would show up in the $O({\it\alpha}^{3})$ inner integral) is exactly zero.

As we did for the inner integral, we can write the outer integral as an expansion in ${\it\alpha}$ :

(4.20) $$\begin{eqnarray}F_{L_{2}}={\it\rho}U_{m}^{2}\ell ^{2}(k_{4}{\it\alpha}^{4}+k_{5}{\it\alpha}^{5}+\cdots ).\end{eqnarray}$$

The velocity terms that contribute to $k_{4}$ are the following:

(4.21) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{U}_{1}^{(0)}\sim {\it\alpha}^{3}\frac{1}{R^{2}}\boldsymbol{f}_{1}^{1},\quad \boldsymbol{U}_{2}^{(0)}\sim {\it\alpha}^{3}\mathscr{S}\left[\frac{1}{R^{2}}\boldsymbol{f}_{1}^{1}\right],\\ \displaystyle \hat{\boldsymbol{U}}_{1}\sim {\it\alpha}\frac{1}{R}\boldsymbol{g}_{0}^{1},\quad \hat{\boldsymbol{U}}_{2}\sim {\it\alpha}\mathscr{S}\left[\frac{1}{R}\boldsymbol{g}_{0}^{1}\right],\quad \bar{\boldsymbol{U}}\sim {\it\gamma}R+{\it\delta}R^{2}+\cdots \,.\end{array}\!\right\}\end{eqnarray}$$

Again, we define $\boldsymbol{V}=\mathscr{S}[\boldsymbol{U}]$ as the image of the function $\boldsymbol{U}$ , and we compute $\mathscr{S}[(1/R^{2})\boldsymbol{f}_{1}^{1}]$ and $\mathscr{S}[(1/R)\boldsymbol{g}_{0}^{1}]$ numerically. Furthermore, we can approximate the term $\boldsymbol{f}_{1}^{1}$ by the stresslet terms, since the rotlet terms have coefficients that are order $O({\it\alpha}^{2})$ higher than the coefficients of the stresslet terms. The $O(a^{4})$ contribution to the reciprocal theorem integral takes the following form:

(4.22) $$\begin{eqnarray}k_{4}=\int _{V_{C}}(\hat{\boldsymbol{U}}_{1}+\hat{\boldsymbol{U}}_{2})\boldsymbol{\cdot }[\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(\boldsymbol{U}_{1}^{(0)}+\boldsymbol{U}_{2}^{(0)})+(\boldsymbol{U}_{1}^{(0)}+\boldsymbol{U}_{2}^{(0)})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{U}}]\,\text{d}v.\end{eqnarray}$$

We run into a problem numerically evaluating the integral in (4.22) when considering only the first terms in the series expansions, $\boldsymbol{U}_{1}^{(0)}$ and $\hat{\boldsymbol{U}}_{1}$ . The problem arises because $\boldsymbol{U}_{1}^{(0)}$ and $\hat{\boldsymbol{U}}_{1}$ have singularities of the form

(4.23a,b ) $$\begin{eqnarray}\boldsymbol{U}_{1}^{(0)}\approx -\frac{5{\it\gamma}_{y}}{2}\frac{(Y-Y_{0})Z\boldsymbol{R}}{R^{5}},\quad \hat{\boldsymbol{U}}_{1}\approx \frac{3}{4}\left(\boldsymbol{e}_{Y}+\frac{(Y-Y_{0})\boldsymbol{R}}{R^{2}}\right)\frac{1}{R},\end{eqnarray}$$

which are respectively the stresslet and stokeslet components of the two velocity fields. When the singularities are integrated against the shear term of $\bar{\boldsymbol{U}}$ , that is $\bar{\boldsymbol{U}}_{{\it\gamma}}\approx {\it\gamma}_{y}Y\boldsymbol{e}_{\boldsymbol{Z}}$ , the result is an integral that is undefined near $R=0$ :

(4.24) $$\begin{eqnarray}\int _{R<{\it\epsilon}}\hat{\boldsymbol{U}}_{1}\boldsymbol{\cdot }[\bar{\boldsymbol{U}}_{{\it\gamma}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}_{1}^{(0)}+\boldsymbol{U}_{1}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{U}}_{{\it\gamma}}]\,\text{d}v.\end{eqnarray}$$

However, converting to spherical coordinates, we find that the angular dependence forces the integral in (4.24) to be zero:

(4.25) $$\begin{eqnarray}\int _{0}^{{\rm\pi}}\!\!\int _{0}^{2{\rm\pi}}\!\!\int _{0}^{{\it\epsilon}}\!\left(\frac{15{\it\gamma}_{y}^{2}(1+2\cos 2{\it\theta})\sin ^{4}{\it\theta}\sin ^{3}{\it\phi}}{4R}\right)\text{d}R\,\text{d}{\it\phi}\,\text{d}{\it\theta}=0.\end{eqnarray}$$

This angular behaviour is difficult to capture numerically, especially if the mesh is not symmetric. Instead, we propose a regularization of the outer integral, where we integrate the problematic terms analytically in a small region near $R=0$ . Now considering the full expansion of $\bar{u}$ , we derive the following analytic form for the integral in the region near the origin:

(4.26) $$\begin{eqnarray}\int _{R<{\it\epsilon}}\hat{\boldsymbol{U}}_{1}\boldsymbol{\cdot }[\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}_{1}^{(0)}+\boldsymbol{U}_{1}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{U}}]\,\text{d}v=-{\rm\pi}{\it\gamma}_{y}({\it\delta}_{xx}+3{\it\delta}_{yy}){\it\epsilon}.\end{eqnarray}$$

Recall that the constants ${\it\gamma}_{y}$ , ${\it\delta}_{xx}$ and ${\it\delta}_{yy}$ were defined in the expansion of $\bar{\boldsymbol{u}}$ in (2.2). Using this analytic expression, we split up the rest of the reciprocal theorem integral (4.22) into the following parts:

(4.27) $$\begin{eqnarray}\displaystyle k_{4} & = & \displaystyle \int _{V_{C}}\hat{\boldsymbol{U}}_{2}\boldsymbol{\cdot }[\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(\boldsymbol{U}_{1}^{(0)}+\boldsymbol{U}_{2}^{(0)})+(\boldsymbol{U}_{1}^{(0)}+\boldsymbol{U}_{2}^{(0)})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{U}}]\,\text{d}v\nonumber\\ \displaystyle & & \displaystyle +\,\int _{V_{C}}\hat{\boldsymbol{U}}_{1}\boldsymbol{\cdot }[\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}_{2}^{(0)}+\boldsymbol{U}_{2}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{U}}]\,\text{d}v\nonumber\\ \displaystyle & & \displaystyle +\,\int _{\{\boldsymbol{r}\in V_{C}\,:\,R\geqslant {\it\epsilon}\}}\hat{\boldsymbol{U}}_{1}\boldsymbol{\cdot }[\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}_{1}^{(0)}+\boldsymbol{U}_{1}^{(0)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\bar{\boldsymbol{U}}]\,\text{d}v\nonumber\\ \displaystyle & & \displaystyle -\,{\rm\pi}{\it\gamma}_{y}({\it\delta}_{xx}+3{\it\delta}_{yy}){\it\epsilon}.\end{eqnarray}$$

The first three lines in (4.27) are evaluated numerically using the FEM. Evaluating the integral in (4.27), we arrive at the scaling law

(4.28) $$\begin{eqnarray}F_{L_{2}}=\frac{{\it\rho}U_{m}^{2}k_{4}a^{4}}{\ell ^{2}}+O(a^{5}),\end{eqnarray}$$

where $k_{4}=O(1)$ is a constant that depends on the location of the particle in the channel, and is computed numerically.

Similarly, the $O({\it\alpha}^{5})$ correction to the outer integral comes from terms

(4.29) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{U}_{1}^{(0)}\sim {\it\alpha}^{3}\frac{1}{R^{2}}\boldsymbol{f}_{1}^{1},\quad \boldsymbol{U}_{2}^{(0)}\sim {\it\alpha}^{3}\mathscr{S}\left[\displaystyle \frac{1}{R^{2}}\boldsymbol{f}_{1}^{1}\right],\\ \displaystyle \hat{\boldsymbol{U}}_{3}\sim {\it\alpha}^{2}\frac{1}{R}\boldsymbol{g}_{0}^{3},\quad \hat{\boldsymbol{U}}_{4}\sim {\it\alpha}^{2}\mathscr{S}\left[\displaystyle \frac{1}{R}\boldsymbol{g}_{0}^{3}\right],\quad \bar{\boldsymbol{U}}\sim {\it\gamma}R+{\it\delta}R^{2}+\cdots \,.\end{array}\right\}\end{eqnarray}$$

Again, we must regularize the outer integral, since $\hat{\boldsymbol{U}}_{3}$ also has a stokeslet singularity. We use the same regularization as before, replacing $\hat{\boldsymbol{U}}_{1}$ and $\hat{\boldsymbol{U}}_{2}$ with $\hat{\boldsymbol{U}}_{3}$ and $\hat{\boldsymbol{U}}_{4}$ , respectively.

Finally, combining terms at $O(a^{4})$ and $O(a^{5})$ , we obtain

(4.30) $$\begin{eqnarray}f_{L_{2}}=\frac{{\it\rho}U_{m}^{2}k_{4}a^{4}}{\ell ^{2}}+\frac{{\it\rho}U_{m}^{2}k_{5}a^{5}}{\ell ^{3}}+O(a^{6}),\end{eqnarray}$$

where $k_{5}=O(1)$ is a constant that depends on the location of the particle in the channel. We have now calculated the $V_{2}$ contribution to the reciprocal theorem integral up to order $O(a^{5})$ .

4.5. Results

In the last section, we described our method of computing the correction to the scaling law made by Ho & Leal (Reference Ho and Leal1974). Combining the inner and outer integrals, the result is a new approximation of the form

(4.31) $$\begin{eqnarray}F_{L}=\frac{{\it\rho}U_{m}^{2}c_{4}a^{4}}{\ell ^{2}}+\frac{{\it\rho}U_{m}^{2}c_{5}a^{5}}{\ell ^{3}}+O(a^{6}),\end{eqnarray}$$

where $c_{4}=k_{4}$ from (4.27), and $c_{5}=h_{5}+k_{5}$ from (4.18) and (4.30). The prefactors $c_{4}$ and $c_{5}$ are $O(1)$ in ${\it\alpha}$ , and depend only on the location of the particle in the channel. The extended series agrees well with numerical data for particle sizes up to ${\it\alpha}=0.2{-}0.3$ (figure 6). This calculation could in principle be extended by computing the contributions from higher-order terms. Completing the series (Hinch Reference Hinch1991), i.e. approximating

(4.32) $$\begin{eqnarray}F_{L}\approx \frac{{\it\rho}U_{m}^{2}c_{4}a^{4}}{\ell ^{2}\left(\displaystyle 1-\frac{c_{5}a}{c_{4}\ell }\right)},\end{eqnarray}$$

produces a modest increase in the accuracy of the asymptotic approximation (figure 6).

Figure 6. Numerical computation of the lift force $F_{L}$ using the NSEs in (2.4) and plotted as a function of particle radius $a$ for channel Reynolds numbers $\mathit{Re}=10$ (triangles), $\mathit{Re}=50$ (circles) and $\mathit{Re}=80$ (crosses). The dashed (blue) line represents a scaling law of particle radius to the fourth power, $F_{L}={\it\rho}U_{m}^{2}c_{4}{\it\alpha}^{2}a^{2}$ , the solid (black) line represents the sum of the fourth and fifth power terms in (4.31), and the dotted (green) line represents the completion of series in (4.32), with (a) particle displacement $y_{0}=0.15/{\it\alpha}$ and (b) particle displacement $y_{0}=0.4/{\it\alpha}$ .

By including two terms in our asymptotic expansion, we can describe how the particle equilibrium position depends on its size – a key prediction for rationally designing devices that use inertial lift forces to fractionate particles, or to transfer them between fluid streams (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Hur et al. Reference Hur, Tse and Di Carlo2010; Mach et al. Reference Mach, Kim, Arshi, Hur and Di Carlo2011; Chung et al. Reference Chung, Gossett and Di Carlo2013; Sollier et al. Reference Sollier, Go, Che, Gossett, O’Byrne, Weaver, Kummer, Rettig, Goldman, Nickols, McCloskey, Kulkarni and Di Carlo2014) (figure 8). We compare our asymptotic calculation predictions directly with experiments of Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009), finding good agreement in focusing positions up to $a=0.3$ (figure 8 b).

5. Three-dimensional asymptotic expansion

Previous asymptotic studies have considered inertial migration in 2D flows (Ho & Leal Reference Ho and Leal1974; Schonberg & Hinch Reference Schonberg and Hinch1989; Hogg Reference Hogg1994; Asmolov Reference Asmolov1999). At sufficiently small values of $a$ , there is qualitative agreement between the 2D theories and our theory, but only when the particle is located on a symmetry plane, e.g.  $x_{0}=0$ or $y_{0}=0$ . However, real inertial microfluidic devices focus in $x$ - and $y$ -directions, taking initially uniformly dispersed particles to four focusing positions. Our asymptotic approach allows us to compute the focusing forces for particles placed at arbitrary positions in the channel.

The calculation is very similar to the one outlined in § 4; we only need to add similar terms driven by the shear in the $x$ -direction, and allow for a reciprocal velocity $\hat{\boldsymbol{u}}$ associated with moving the particle in this direction. The full Lamb’s solution for $\boldsymbol{u}_{1}^{(0)}$ has additional terms from the shear in the $x$ -direction (i.e. the terms with coefficients ${\it\gamma}_{x}$ ), shown in appendix C. The only additional components of $\boldsymbol{u}_{1}^{(0)}$ that contribute to the 3D calculation are the stresslet and source quadrupole.

Figure 7. Numerical computation of the lift force $F_{L}$ from (2.4) as a function of particle radius for channel Reynolds numbers $\mathit{Re}=10$ (triangles), $\mathit{Re}=50$ (circles) and $\mathit{Re}=80$ (crosses). The dashed (blue) line represents a scaling law of particle radius to the fourth power, $F_{L}={\it\rho}U_{m}^{2}c_{4}{\it\alpha}^{2}a^{2}$ , the solid (black) line represents the fifth power correction term in (5.4), and the dotted (green) line represents the completion of series in (4.32), with particle displacement $x_{0}=0.2/{\it\alpha}$ and $y_{0}=0.15/{\it\alpha}$ for lift force (a) in the $x$ -direction and (b) in the $y$ -direction.

The inner integral in 3D evaluates to

(5.1) $$\begin{eqnarray}F_{L_{1}}^{(3D)}=\frac{{\it\rho}U_{m}^{2}h_{5}^{(3D)}a^{5}}{\ell ^{3}}+O(a^{6}),\end{eqnarray}$$

where

(5.2) $$\begin{eqnarray}\displaystyle h_{5}^{(3D)}=\frac{4381{\rm\pi}{\it\gamma}_{x}{\it\gamma}_{y}}{554\,400}-\frac{26\,171{\rm\pi}{\it\gamma}_{y}^{2}}{277\,200}+\frac{527{\rm\pi}{\it\psi}_{y}{\it\gamma}_{x}{\it\gamma}_{y}}{116\,424}-\frac{53{\rm\pi}{\it\gamma}_{y}{\it\delta}_{xx}}{1728}+\frac{19{\rm\pi}{\it\gamma}_{x}{\it\delta}_{yy}}{3150}-\frac{283{\rm\pi}{\it\gamma}_{y}{\it\delta}_{yy}}{3150}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

We define ${\it\psi}_{y}$ to be the value of the $y$ -component of the image of the stokeslet evaluated at the location of the particle:

(5.3) $$\begin{eqnarray}{\it\psi}_{y}=\left.\left[\mathscr{S}\left[\frac{1}{r}\boldsymbol{g}_{0}^{1}\right]\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{y}}\right]\right|_{(x,y,z)=(x_{0},y_{0},0)},\end{eqnarray}$$

where $\boldsymbol{g}_{0}^{1}$ is the stokeslet and the leading term of $\hat{\boldsymbol{u}}_{1}$ defined in (C 8) in appendix C. The outer integral remains the same; however, $\boldsymbol{u}_{1}^{(0)}$ and $\boldsymbol{u}_{3}^{(0)}$ each now include a stresslet contribution associated with shear in the $x$ -direction. Computing this integral gives a scaling law of the form

(5.4) $$\begin{eqnarray}F_{L}^{(3D)}=\frac{{\it\rho}U_{m}^{2}c_{4}^{(3D)}a^{4}}{\ell ^{2}}+\frac{{\it\rho}U_{m}^{2}c_{5}^{(3D)}a^{5}}{\ell ^{3}}+O(a^{6}).\end{eqnarray}$$

It remains true that, for particles located arbitrarily in the square channel, the lift force scales like $a^{4}$ in the asymptotic limit ${\it\alpha}\rightarrow 0$ . Additionally, our $O(a^{5})$ correction to the scaling law remains accurate for moderately large ${\it\alpha}$ , shown in figure 7(a,b) for the forces in the $x$ - and $y$ -direction, respectively. We provide the calculated values of the 3D lift force in a square channel in a Matlab code in the supplementary data available at http://dx.doi.org/10.1017/jfm.2014.739. In particular, we find that lift forces vanish only at eight symmetrically placed points around the channel, with four points being stable and four unstable, in good agreement with experimental observations (figure 8 a).

Figure 8. (a) Lift force calculated using (5.4) for locations in the lower right quadrant of the channel for a particle of radius ${\it\alpha}=0.11$ and $\mathit{Re}=80$ . The solid (black) circles mark stable equilibrium points, while the open (white) circles mark unstable equilibrium points. (b) Trajectories of particles calculated using (5.5) for particle size ${\it\alpha}=0.11$ and $\mathit{Re}=80$ . The solid (black) circles mark stable equilibrium points, while the open (white) circles mark unstable equilibrium points.

We can compute particle streamlines using the lift force prediction, and confirm that there are four stable focusing positions in the channel (figure 8 b). Particles are advected using a forward Euler time stepping scheme. We find the particle velocity by equating the $O(a^{5})$ lift force (5.4) with the $O(a)$ drag force (Happel & Brenner Reference Happel and Brenner1982). That is, $v_{L}$ , the $y$ -component of velocity $v$ , satisfies the equation

(5.5) $$\begin{eqnarray}6{\rm\pi}{\it\mu}a(v_{L}+{\it\psi}_{y})=\left.\left[\frac{{\it\rho}U_{m}^{2}c_{4}^{(3D)}a^{4}}{\ell ^{2}}+\frac{{\it\rho}U_{m}^{2}c_{5}^{(3D)}a^{5}}{\ell ^{3}}\right]\right|_{(x_{0},y_{0},0)},\end{eqnarray}$$

where ${\it\psi}_{y}$ is the image velocity of the stokeslet defined in (5.3). The velocity $u_{L}$ , the $x$ -component of velocity $u$ , is computed in the same way by substituting ${\it\psi}_{x}$ from the $x$ -stokeslet for ${\it\psi}_{y}$ .

In addition, the distance of the focusing positions from the channel centre-line can be predicted by solving the implicit equation $F_{L}^{(3D)}=0$ . Recall that the lift force coefficients depend on the location of the particle, i.e.  $c_{4}^{(3D)}=c_{4}^{(3D)}(x_{0},y_{0})$ and $c_{5}^{(3D)}=c_{5}^{(3D)}(x_{0},y_{0})$ . Since the lift force formula has both $O(a^{4})$ and $O(a^{5})$ terms, the focusing position will have a functional dependence on the particle size $a$ . This prediction of the focusing position compares well with experimental data by Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009), especially for particle sizes up to ${\it\alpha}\leqslant 0.3$ (figure 9 b).

Figure 9. (a) Inertial focusing position as a function of particle size predicted by our theory. The markers are data collected by Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Toner2009), the dashed (blue) line is the theory predicted by the first term of $O(a^{4})$ in (4.31), and the solid (red) line is the theory predicted by (4.31). (b) Schematic diagram plotting the outlines of particles at their predicted focusing position along the positive $x$ -axis. The particle sizes range between ${\it\alpha}=0.03$ and ${\it\alpha}=0.29$ .