2.1 Governing equations in the lab frame of reference

Figure 1. Curved duct with rectangular cross-section containing a spherical particle located at $\boldsymbol{x}_{p}=\boldsymbol{x}(\unicode[STIX]{x1D703}_{p},r_{p},z_{p})$ . The enlarged view of the cross-section around the particle illustrates the origin of the local $r,z$ coordinates at the centre of the duct. The bend radius $R$ is with respect to the centreline of the duct and is of modest size for illustration.

Consider a duct which is curved, having a constant bend radius and uniform cross-sectional shape, an example of which is depicted in figure 1. Let $\boldsymbol{x}=(x,y,z)$ denote a Cartesian coordinate system in the lab reference frame with the duct positioned so that the $z$ axis is normal to the plane of the bend. Let the bend radius $R$ be the distance from the origin to the centre of the duct cross-section at any axial/angular position. We introduce $r,\unicode[STIX]{x1D703}$ such that $r,z$ are coordinates that describe the location within the duct cross-section relative to its centre and $\unicode[STIX]{x1D703}$ is the angle around the bend measured from the $x$ axis (see figure 1). Note that we do not consider the flow near the inlet/outlet. With this description $\boldsymbol{x}$ is parametrised as

(2.1) $$\begin{eqnarray}\boldsymbol{x}(r,\unicode[STIX]{x1D703},z)=(R+r)\cos (\unicode[STIX]{x1D703})\boldsymbol{e}_{x}+(R+r)\sin (\unicode[STIX]{x1D703})\boldsymbol{e}_{y}+z\boldsymbol{e}_{z}.\end{eqnarray}$$

For example, the point on the centreline $\boldsymbol{x}=R\boldsymbol{e}_{x}$ is equivalently described by $(r,\unicode[STIX]{x1D703},z)=(0,0,0)$ . For a duct having a rectangular cross-section with width $W$ and height $H$ , the walls of the duct are taken to be located at $r=\pm W/2$ and $z=\pm H/2$ . For the non-rectangular ducts considered later we take $H$ to be the height at $r=0$ and the origin of the cross-section to be such that for $r=0$ the top and bottom walls are located at $z=\pm H/2$ , respectively. Let $\ell =H$ be a characteristic minimum width of the duct cross-section, noting that the ducts considered in this paper all have $H\leqslant W$ .

Consider a steady fluid flow, through the duct, satisfying the incompressible Navier–Stokes equations. Letting $\bar{\boldsymbol{u}}$ denote the fluid velocity and $\bar{p}$ denote the pressure, we have

(2.2a ) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D735}\bar{p}+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\bar{\boldsymbol{u}} & = & \displaystyle \unicode[STIX]{x1D70C}(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{\boldsymbol{u}}-\boldsymbol{g})\quad \text{on }\boldsymbol{x}\in {\mathcal{D}},\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\boldsymbol{u}} & = & \displaystyle 0\quad \hphantom{(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{\boldsymbol{u}}-\boldsymbol{g})}\text{on }\boldsymbol{x}\in {\mathcal{D}},\end{eqnarray}$$

(2.2c ) $$\begin{eqnarray}\displaystyle \bar{\boldsymbol{u}} & = & \displaystyle \mathbf{0}\quad \hphantom{(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{\boldsymbol{u}}-\boldsymbol{g})}\text{on }\boldsymbol{x}\in \unicode[STIX]{x2202}{\mathcal{D}},\end{eqnarray}$$

$\unicode[STIX]{x1D70C}$

$\unicode[STIX]{x1D707}$

$\boldsymbol{g}$

${\mathcal{D}}$

$\unicode[STIX]{x2202}{\mathcal{D}}$

$\bar{p}$

(2.3) $$\begin{eqnarray}\bar{p}(\boldsymbol{x}(r,\unicode[STIX]{x1D703},z))=-GR\unicode[STIX]{x1D703}+\unicode[STIX]{x1D70C}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{g}+\tilde{p}(r,z),\end{eqnarray}$$

where $G$ denotes the drop in pressure per unit arclength along the centreline of the duct. Notice that we have assumed gravity does not influence the flow velocity by introducing the hydrostatic pressure $\unicode[STIX]{x1D70C}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{g}$ to counteract $\unicode[STIX]{x1D70C}\boldsymbol{g}$ in (2.2a ). For simplicity our study is restricted to the specific case where gravity acts in the $z$ direction – that is, $\boldsymbol{g}=-g\boldsymbol{e}_{z}$ .

The resulting flow, often referred to as Dean flow, has been well-studied for ducts having circular or rectangular cross-sections (Dean & Hurst Reference Dean and Hurst1959; Winters Reference Winters1987; Yanase, Goto & Yamamoto Reference Yanase, Goto and Yamamoto1989; Yamamoto et al. Reference Yamamoto, Wu, Hyakutake and Yanase2004). The pressure gradient drives an axial flow from which the inertia of the fluid around the bend then leads to the development of a secondary flow consisting of two counter-rotating vortices within the duct cross-section. (At high Dean numbers there can exist multiple solutions, some having several vortices; however, such flow conditions are far beyond those of practical interest in the context of microfluidics). Sufficiently far from the inlet/outlet region, the flow velocity components are independent of the angular coordinate (up to their direction) and can be expressed as

(2.4) $$\begin{eqnarray}\bar{\boldsymbol{u}}=\bar{u}_{\unicode[STIX]{x1D703}}(r,z)\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}+\bar{u}_{r}(r,z)\hat{\boldsymbol{e}}_{r}+\bar{u}_{z}(r,z)\boldsymbol{e}_{z},\end{eqnarray}$$

where $\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}:=-\text{sin}(\unicode[STIX]{x1D703})\boldsymbol{e}_{x}+\cos (\unicode[STIX]{x1D703})\boldsymbol{e}_{y}$ and $\hat{\boldsymbol{e}}_{r}:=\cos (\unicode[STIX]{x1D703})\boldsymbol{e}_{x}+\sin (\unicode[STIX]{x1D703})\boldsymbol{e}_{y}$ . The $\bar{u}_{\unicode[STIX]{x1D703}}$ component is the axial flow velocity, whereas the $\bar{u}_{r},\bar{u}_{z}$ components describe the secondary flow. Throughout this paper $\bar{\boldsymbol{u}}$ and $\bar{p}$ are referred to as the background velocity and pressure, respectively. In the context of studying the dynamics of a particle suspended in flow through a curved duct $\bar{\boldsymbol{u}},\bar{p}$ will be treated as being known/prescribed. Details of the specific approximation of $\bar{\boldsymbol{u}}$ used to obtain our results is deferred to § 3, while the treatment throughout the remainder of this section is quite general.

Consider a single spherical particle within the fluid flow through the duct. The axial velocity of the fluid is faster than its lateral velocity and, consequently, the same is expected for the velocity of a particle suspended in the flow. Therefore, we consider the particle to be travelling along a path traced out by the axial component of the background flow (i.e. a ‘streamline’ of $\bar{u}_{\unicode[STIX]{x1D703}}\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}}$ ) and that its motion is perturbed in the plane of the cross-section by drag force from the secondary fluid motion, inertia due to the curvature of the path, and inertial lift force. By separating the lateral dynamics from the (faster) axial dynamics we develop a framework in which particle focusing is more readily analysed.

Let $a$ denote the radius of the spherical particle, $\unicode[STIX]{x1D70C}_{p}$ denote its density (which we assume to be uniform), and $\boldsymbol{x}_{p}(t)=(x_{p}(t),y_{p}(t),z_{p}(t))$ denote the location of its centre at time $t$ . It is convenient to parametrise the particle position in cylindrical coordinates as $\boldsymbol{x}_{p}(t)=\boldsymbol{x}(r_{p}(t),\unicode[STIX]{x1D703}_{p}(t),z_{p}(t))$ (see figure 1). The particle follows an axial flow streamline (from which it will be perturbed) and thus $\unicode[STIX]{x1D6E9}:=\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{p}/\unicode[STIX]{x2202}t$ is constant, $\unicode[STIX]{x2202}r_{p}/\unicode[STIX]{x2202}t=\unicode[STIX]{x2202}z_{p}/\unicode[STIX]{x2202}t=0$ , and $\unicode[STIX]{x2202}\unicode[STIX]{x1D734}_{p}/\unicode[STIX]{x2202}t=\mathbf{0}$ where $\unicode[STIX]{x1D734}_{p}$ denotes the spin of the particle (i.e. the angular velocity with respect to its centre). The assumption $\unicode[STIX]{x2202}r_{p}/\unicode[STIX]{x2202}t=\unicode[STIX]{x2202}z_{p}/\unicode[STIX]{x2202}t=0$ may seem counter-intuitive (i.e. we assume that $r_{p},z_{p}$ is fixed only to calculate non-zero forces in these directions such that they cannot remain fixed). However, this is consistent with the approaches in previous studies of inertial lift forces in uni-directional flows (see, for example, Schonberg & Hinch (Reference Schonberg and Hinch1989) and Hood et al. (Reference Hood, Lee and Roper2015)). Furthermore, it can be shown that any additional effects from particle motion in the $r,z$ directions is much smaller than the inertial lift force and can therefore be neglected (given the other assumptions used here). Without loss of generality we can additionally assume that $\unicode[STIX]{x1D703}_{p}=0$ at time $t=0$ such that $\unicode[STIX]{x1D703}_{p}(t)=U_{p}t/(R+r_{p})$ (or equivalently $\unicode[STIX]{x1D6E9}=U_{p}/(R+r_{p})$ ), where $U_{p}$ here denotes the (linear) velocity of the particle.

The presence of the particle alters the fluid flow and we introduce $\boldsymbol{u},p$ to denote the fluid velocity and pressure in the presence of the particle. The fluid motion is again assumed to be governed by the incompressible Navier–Stokes equations; that is, in general

(2.5a ) $$\begin{eqnarray}\hspace{16.99998pt}-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}=\unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}+g\boldsymbol{e}_{z}\right)\qquad \qquad \hspace{5.0pt}\text{on }\boldsymbol{x}\in {\mathcal{F}},\end{eqnarray}$$

(2.5b,c ) $$\begin{eqnarray}\hphantom{-\unicode[STIX]{x1D735}p+\qquad \qquad }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0\quad \text{on }\boldsymbol{x}\in {\mathcal{F}},\qquad \qquad \hspace{15.00002pt}\boldsymbol{u}=\mathbf{0}\quad \text{on }\boldsymbol{x}\in \unicode[STIX]{x2202}{\mathcal{D}},\end{eqnarray}$$

(2.5d ) $$\begin{eqnarray}\hspace{30.00005pt}\boldsymbol{u}=\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{p}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D734}_{p}\times (\boldsymbol{x}-\boldsymbol{x}_{p})\hspace{60.00009pt}\text{on }\boldsymbol{x}\in \unicode[STIX]{x2202}{\mathcal{F}}\backslash \unicode[STIX]{x2202}{\mathcal{D}},\hspace{-40.00006pt}\end{eqnarray}$$

where ${\mathcal{F}}:=\{\boldsymbol{x}\in {\mathcal{D}}:|\boldsymbol{x}-\boldsymbol{x}_{p}|\geqslant a\}$ is the fluid domain given the presence of the spherical particle and $\unicode[STIX]{x2202}{\mathcal{F}}\backslash \unicode[STIX]{x2202}{\mathcal{D}}=\{\boldsymbol{x}:|\boldsymbol{x}-\boldsymbol{x}_{p}|=a\}$ is the surface of the particle. Note that the moving particle results in unsteady fluid flow in the lab frame of reference. In particular, ${\mathcal{F}}$ depends on time because $\boldsymbol{x}_{p}$ is non-stationary. Thus the $\unicode[STIX]{x2202}\boldsymbol{u}/\unicode[STIX]{x2202}t$ in (2.5a ) cannot be discarded (as it was for the background flow in (2.2a )). Note that the boundary condition (2.5d ) has a constant $\unicode[STIX]{x1D734}_{p}$ that denotes the spin (or angular velocity) of the particle (as observed in the lab frame). Far from the particle we expect $\boldsymbol{u},p$ to converge to the background flow $\bar{\boldsymbol{u}},\bar{p}$ .

The particle motion is driven solely by the force and torque exerted on the particle by the surrounding fluid in addition to the gravitational body force. The former are obtained by integrating the fluid stress over the particle surface so that the total force and torque on the particle is

(2.6a ) $$\begin{eqnarray}\displaystyle \boldsymbol{F}(t) & = & \displaystyle -m_{p}g\boldsymbol{e}_{z}+\int _{|\boldsymbol{x}-\boldsymbol{x}_{p}|=a}(-\boldsymbol{n})\boldsymbol{\cdot }(-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\intercal }))\,\text{d}S,\end{eqnarray}$$

(2.6b ) $$\begin{eqnarray}\displaystyle \boldsymbol{T}(t) & = & \displaystyle \int _{|\boldsymbol{x}-\boldsymbol{x}_{p}|=a}(\boldsymbol{x}-\boldsymbol{x}_{p})\times \left((-\boldsymbol{n})\boldsymbol{\cdot }(-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\intercal }))\right)\,\text{d}S,\end{eqnarray}$$

$m_{p}:=(4/3)\unicode[STIX]{x03C0}a^{3}\unicode[STIX]{x1D70C}_{p}$

$\boldsymbol{n}$

${\mathcal{F}}$

$-\boldsymbol{n}$

2.2 Equations of motion in a rotating frame of reference

Although the fluid motion is not steady in the lab frame, we view it as steady relative to the motion of the particle centre (assumed to be moving with constant angular velocity $\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{p}/\unicode[STIX]{x2202}t=\unicode[STIX]{x1D6E9}$ and fixed $r_{p},z_{p}$ ). Therefore, in the particle reference frame (noting that the particle may still have a constant spin) we can eliminate the $\unicode[STIX]{x2202}\boldsymbol{u}/\unicode[STIX]{x2202}t$ term.

Consider a reference frame rotating about the $z$ axis at the (constant) rate $\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{p}/\unicode[STIX]{x2202}t:=\unicode[STIX]{x1D6E9}$ . We introduce the coordinate system $\boldsymbol{x}^{\prime }=(x^{\prime },y^{\prime },z^{\prime })$ , which we parametrise by $(r^{\prime },\unicode[STIX]{x1D703}^{\prime },z^{\prime })$ such that

(2.7) $$\begin{eqnarray}\boldsymbol{x}^{\prime }(r^{\prime },\unicode[STIX]{x1D703}^{\prime },z^{\prime })=(R+r^{\prime })\cos (\unicode[STIX]{x1D703}^{\prime })\boldsymbol{e}_{x^{\prime }}+(R+r^{\prime })\sin (\unicode[STIX]{x1D703}^{\prime })\boldsymbol{e}_{y^{\prime }}+z^{\prime }\boldsymbol{e}_{z^{\prime }}.\end{eqnarray}$$

The unit vectors $\boldsymbol{e}_{x^{\prime }},\boldsymbol{e}_{y^{\prime }},\boldsymbol{e}_{z^{\prime }}$ in the rotating frame are related to $\boldsymbol{e}_{x},\boldsymbol{e}_{y},\boldsymbol{e}_{z}$ in the lab reference frame by

(2.8a-c ) $$\begin{eqnarray}\boldsymbol{e}_{x^{\prime }}=\cos (\unicode[STIX]{x1D703}_{p})\boldsymbol{e}_{x}+\sin (\unicode[STIX]{x1D703}_{p})\boldsymbol{e}_{y},\quad \boldsymbol{e}_{y^{\prime }}=-\sin (\unicode[STIX]{x1D703}_{p})\boldsymbol{e}_{x}+\cos (\unicode[STIX]{x1D703}_{p})\boldsymbol{e}_{y},\quad \boldsymbol{e}_{z^{\prime }}=\boldsymbol{e}_{z}.\end{eqnarray}$$

or conversely

(2.9a-c ) $$\begin{eqnarray}\boldsymbol{e}_{x}=\cos (\unicode[STIX]{x1D703}_{p})\boldsymbol{e}_{x^{\prime }}-\sin (\unicode[STIX]{x1D703}_{p})\boldsymbol{e}_{y^{\prime }},\quad \boldsymbol{e}_{y}=\sin (\unicode[STIX]{x1D703}_{p})\boldsymbol{e}_{x^{\prime }}+\cos (\unicode[STIX]{x1D703}_{p})\boldsymbol{e}_{y^{\prime }},\quad \boldsymbol{e}_{z}=\boldsymbol{e}_{z^{\prime }}.\end{eqnarray}$$

The angular velocity of the rotating reference frame (relative to the lab frame) is

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D6E9}\boldsymbol{e}_{z}:=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{p}}{\unicode[STIX]{x2202}t}\boldsymbol{e}_{z}=\frac{U_{p}}{R+r_{p}}\boldsymbol{e}_{z}.\end{eqnarray}$$

Let $\boldsymbol{x}_{p}^{\prime }=(x_{p}^{\prime },y_{p}^{\prime },z_{p}^{\prime })$ be the centre of the particle in the rotating frame of reference. It is always the case that $x_{p}^{\prime }=R+r_{p}$ , $y_{p}^{\prime }=0$ and $z_{p}^{\prime }=z_{p}$ . Consequently, the particle velocity in the rotating reference frame is

(2.11) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{p}^{\prime }}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}x_{p}^{\prime }}{\unicode[STIX]{x2202}t}\boldsymbol{e}_{x^{\prime }}+\frac{\unicode[STIX]{x2202}z_{p}^{\prime }}{\unicode[STIX]{x2202}t}\boldsymbol{e}_{z^{\prime }}=\frac{\unicode[STIX]{x2202}r_{p}}{\unicode[STIX]{x2202}t}\boldsymbol{e}_{x^{\prime }}+\frac{\unicode[STIX]{x2202}z_{p}}{\unicode[STIX]{x2202}t}\boldsymbol{e}_{z^{\prime }}=0.\end{eqnarray}$$

The particle spin as viewed from the rotating reference frame is reduced by the solid-body rotation of the coordinate system and is therefore denoted

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D734}_{p}^{\prime }:=\unicode[STIX]{x1D734}_{p}-\unicode[STIX]{x1D6E9}\boldsymbol{e}_{z^{\prime }}.\end{eqnarray}$$

Let $\boldsymbol{u}^{\prime }(\boldsymbol{x}^{\prime },t)$ denote the velocity of a fluid parcel in the rotating frame, which is related to the velocity $\boldsymbol{u}(\boldsymbol{x},t)$ in the lab frame via

(2.13) $$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{u}^{\prime }(\boldsymbol{x}^{\prime },t)+\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime }).\end{eqnarray}$$

Furthermore, the accelerations of a fluid parcel in the two frames are related via

(2.14) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{u}}{\text{d}t}=\frac{\text{D}\boldsymbol{u}^{\prime }}{\text{D}t}+2\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{u}^{\prime })+\unicode[STIX]{x1D6E9}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime })),\end{eqnarray}$$

where $\text{D}\boldsymbol{u}^{\prime }/\text{D}t:=\unicode[STIX]{x2202}\boldsymbol{u}^{\prime }/\unicode[STIX]{x2202}t+\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\boldsymbol{u}^{\prime }$ denotes the material derivative. Letting $p^{\prime }$ denote the fluid pressure in the rotating frame (which is such that $p^{\prime }(\boldsymbol{x}^{\prime },t)=p(\boldsymbol{x},t)$ and thus $\unicode[STIX]{x1D735}^{\prime }p^{\prime }=\unicode[STIX]{x1D735}p$ ), it follows that the equations of motion for the fluid in the rotating reference frame are

(2.15a ) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D735}^{\prime }p^{\prime }+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{\prime \,2}\boldsymbol{u}^{\prime } & = & \displaystyle \unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}^{\prime }}{\unicode[STIX]{x2202}t}+\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\boldsymbol{u}^{\prime }+2\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{u}^{\prime })\right.\qquad \nonumber\\ \displaystyle & & \displaystyle \left.+\,\unicode[STIX]{x1D6E9}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime }))+g\boldsymbol{e}_{z}\right)\quad \text{on }\boldsymbol{x}^{\prime }\in {\mathcal{F}}^{\prime },\qquad\end{eqnarray}$$

(2.15b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}^{\prime } & = & \displaystyle 0\hspace{125.00018pt}\text{on }\boldsymbol{x}^{\prime }\in {\mathcal{F}}^{\prime },\end{eqnarray}$$

(2.15c ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{\prime } & = & \displaystyle -\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime })\hspace{77.00008pt}\text{on }\boldsymbol{x}^{\prime }\in \unicode[STIX]{x2202}{\mathcal{D}},\qquad\end{eqnarray}$$

(2.15d ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{\prime } & = & \displaystyle \unicode[STIX]{x1D734}_{p}^{\prime }\times (\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime })\hspace{70.0001pt}\text{on }\boldsymbol{x}^{\prime }\in \unicode[STIX]{x2202}{\mathcal{F}}^{\prime }\backslash \unicode[STIX]{x2202}{\mathcal{D}},\qquad\end{eqnarray}$$

${\mathcal{F}}^{\prime }:=\{\boldsymbol{x}^{\prime }\in {\mathcal{D}}:|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|\geqslant a\}$

$-\unicode[STIX]{x1D6E9}\boldsymbol{e}_{z}$

${\mathcal{F}}^{\prime }$

$\unicode[STIX]{x1D6E9},\unicode[STIX]{x1D734}_{p}^{\prime },g$

$\unicode[STIX]{x2202}\boldsymbol{u}^{\prime }/\unicode[STIX]{x2202}t$

$\boldsymbol{e}_{z}$

$\unicode[STIX]{x1D703}_{p}(t)$

The force and torque on the particle need to be expressed in terms of the rotating reference frame as well. Let $\boldsymbol{F}^{\prime },\boldsymbol{T}^{\prime }$ denote the force and torque on the particle in the rotating reference frame, respectively. Then one has

(2.16a ) $$\begin{eqnarray}\displaystyle \boldsymbol{F} & = & \displaystyle \boldsymbol{F}^{\prime }+m_{p}\unicode[STIX]{x1D6E9}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}_{p}^{\prime })),\end{eqnarray}$$

(2.16b ) $$\begin{eqnarray}\displaystyle \boldsymbol{T} & = & \displaystyle \boldsymbol{T}^{\prime }+I_{p}\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \unicode[STIX]{x1D734}_{p}^{\prime }),\end{eqnarray}$$

$m_{p}:=(4/3)\unicode[STIX]{x03C0}a^{3}\unicode[STIX]{x1D70C}_{p}$

$I_{p}:=(2/5)a^{2}m_{p}$

$\unicode[STIX]{x1D6E9}$

$p^{\prime },\boldsymbol{u}^{\prime }$

(2.17a ) $$\begin{eqnarray}\displaystyle \boldsymbol{F}^{\prime } & = & \displaystyle -m_{p}g\boldsymbol{e}_{z^{\prime }}-m_{p}\unicode[STIX]{x1D6E9}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}_{p}^{\prime }))\nonumber\\ \displaystyle & & \displaystyle +\,\int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|=a}(-\boldsymbol{n})\boldsymbol{\cdot }(-p^{\prime }\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}^{\prime }\boldsymbol{u}^{\prime }+{\unicode[STIX]{x1D735}^{\prime }\boldsymbol{u}^{\prime }}^{\intercal }))\,\text{d}S^{\prime },\end{eqnarray}$$

(2.17b ) $$\begin{eqnarray}\displaystyle \boldsymbol{T}^{\prime } & = & \displaystyle -I_{p}\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \unicode[STIX]{x1D734}_{p}^{\prime })\nonumber\\ \displaystyle & & \displaystyle +\,\int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|=a}(\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime })\times ((-\boldsymbol{n})\boldsymbol{\cdot }(-p^{\prime }\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}^{\prime }\boldsymbol{u}^{\prime }+\unicode[STIX]{x1D735}^{\prime }\boldsymbol{u}^{\prime \intercal })))\,\text{d}S^{\prime }.\end{eqnarray}$$

$-m_{p}\unicode[STIX]{x1D6E9}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}_{p}^{\prime }))$

We also need the background flow in the rotating reference frame. Since the background flow is steady in time and independent of angular position then, letting $\bar{p}^{\prime },\bar{\boldsymbol{u}}^{\prime }$ denote pressure and velocity of the background fluid flow, respectively, they are related to those in the stationary/lab frame via

(2.18a,b ) $$\begin{eqnarray}\bar{\boldsymbol{u}}^{\prime }(\boldsymbol{x}^{\prime })=\bar{\boldsymbol{u}}(\boldsymbol{x}^{\prime })-\unicode[STIX]{x1D6E9}\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime },\quad \bar{p}^{\prime }(\boldsymbol{x}^{\prime })=\bar{p}(\boldsymbol{x}).\end{eqnarray}$$

We are now equipped to consider the disturbance flow in the rotating reference frame – that is, the difference between flow $\boldsymbol{u}^{\prime },p^{\prime }$ in the presence of the particle and the flow $\bar{\boldsymbol{u}}^{\prime },\bar{p}^{\prime }$ in the absence of the particle.

2.3 Disturbance flow in the rotating frame

In the rotating reference frame we introduce the disturbance velocity and pressure

(2.19a,b ) $$\begin{eqnarray}\boldsymbol{v}^{\prime }=\boldsymbol{u}^{\prime }-\bar{\boldsymbol{u}}^{\prime },\quad q^{\prime }=p^{\prime }-\bar{p}^{\prime },\end{eqnarray}$$

respectively. Substituting $\boldsymbol{u}^{\prime }=\boldsymbol{v}^{\prime }+\bar{\boldsymbol{u}}^{\prime }=\boldsymbol{v}^{\prime }+\bar{\boldsymbol{u}}-\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime })$ and $p^{\prime }=q^{\prime }+\bar{p}^{\prime }$ into (2.15) one obtains

(2.20a ) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D735}^{\prime }q^{\prime }+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{\prime \,2}\boldsymbol{v}^{\prime } & = & \displaystyle \unicode[STIX]{x1D70C}((\boldsymbol{v}^{\prime }+\bar{\boldsymbol{u}}-\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime }))\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\boldsymbol{v}^{\prime }\hphantom{-\unicode[STIX]{x1D735}^{\prime }q^{\prime }+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{\prime \,2}\boldsymbol{v}^{\prime }=}\qquad \qquad \nonumber\\ \displaystyle & & \displaystyle +\,\boldsymbol{v}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}+\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{v}^{\prime })\!)\hspace{48.00009pt}\text{on }\boldsymbol{x}^{\prime }\in {\mathcal{F}}^{\prime },\qquad \qquad\end{eqnarray}$$

(2.20b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime } & = & \displaystyle 0\hspace{145.00021pt}\text{on }\boldsymbol{x}^{\prime }\in {\mathcal{F}}^{\prime },\qquad \qquad\end{eqnarray}$$

(2.20c ) $$\begin{eqnarray}\displaystyle \boldsymbol{v}^{\prime } & = & \displaystyle \mathbf{0}\hspace{145.00021pt}\text{on }\boldsymbol{x}^{\prime }\in \unicode[STIX]{x2202}{\mathcal{D}},\qquad \qquad\end{eqnarray}$$

(2.20d ) $$\begin{eqnarray}\displaystyle \boldsymbol{v}^{\prime } & = & \displaystyle -\bar{\boldsymbol{u}}+\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime })+\unicode[STIX]{x1D734}_{p}^{\prime }\times (\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime })\quad \text{on }\boldsymbol{x}^{\prime }\in \unicode[STIX]{x2202}{\mathcal{F}}^{\prime }\backslash \unicode[STIX]{x2202}{\mathcal{D}},\qquad \qquad\end{eqnarray}$$

1. (i) the background flow $\bar{\boldsymbol{u}},\bar{p}$ satisfies (2.2),

2. (ii) for any $\boldsymbol{w}$ one has $\boldsymbol{w}\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime })=\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{w}$ ,

3. (iii) using (2.4) it can be shown $(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime })\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}=\boldsymbol{e}_{z^{\prime }}\times \bar{\boldsymbol{u}}$ ,

4. (iv) for $\boldsymbol{x}^{\prime }\in \unicode[STIX]{x2202}{\mathcal{D}}$ one has $\bar{\boldsymbol{u}}=\mathbf{0}$ .

Similarly the force and torque in (2.17) can be decomposed as

(2.21a ) $$\begin{eqnarray}\displaystyle \boldsymbol{F}^{\prime } & = & \displaystyle -m_{p}g\boldsymbol{e}_{z}-m_{p}\unicode[STIX]{x1D6E9}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}_{p}^{\prime }))+\int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|=a}(-\boldsymbol{n})\boldsymbol{\cdot }(-\bar{p}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}+{\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}}^{\intercal }))\,\text{d}S^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|=a}(-\boldsymbol{n})\boldsymbol{\cdot }(-q^{\prime }\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}^{\prime }\boldsymbol{v}^{\prime }+{\unicode[STIX]{x1D735}^{\prime }\boldsymbol{v}^{\prime }}^{\intercal }))\,\text{d}S^{\prime },\end{eqnarray}$$

(2.21b ) $$\begin{eqnarray}\displaystyle \boldsymbol{T}^{\prime } & = & \displaystyle -I_{p}\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \unicode[STIX]{x1D734}_{p}^{\prime })+\int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|=a}(\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime })\times ((-\boldsymbol{n})\boldsymbol{\cdot }(-\bar{p}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}+{\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}}^{\intercal })))\,\text{d}S^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|=a}(\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime })\times ((-\boldsymbol{n})\boldsymbol{\cdot }(-q^{\prime }\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}^{\prime }\boldsymbol{v}^{\prime }+{\unicode[STIX]{x1D735}^{\prime }\boldsymbol{v}^{\prime }}^{\intercal })))\,\text{d}S^{\prime },\end{eqnarray}$$

$\unicode[STIX]{x1D735}^{\prime }(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime })+{\unicode[STIX]{x1D735}^{\prime }(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}^{\prime })}^{\intercal }=0$

$q^{\prime },\boldsymbol{v}^{\prime }$

$\bar{p},\bar{\boldsymbol{u}}$

$-m_{p}g\boldsymbol{e}_{z}-m_{p}\unicode[STIX]{x1D6E9}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}_{p}^{\prime }))$

$-I_{p}\unicode[STIX]{x1D6E9}(\boldsymbol{e}_{z^{\prime }}\times \unicode[STIX]{x1D734}_{p}^{\prime })$

$\bar{\boldsymbol{u}},\bar{p}$

(2.22) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|=a}(-\boldsymbol{n})\boldsymbol{\cdot }(-\bar{p}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}+{\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}}^{\intercal }))\,\text{d}S^{\prime }\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D70C}\int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|\leqslant a}\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}+g\boldsymbol{e}_{z}\,\text{d}V^{\prime }\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{4}{3}\unicode[STIX]{x03C0}a^{3}\unicode[STIX]{x1D70C}g\boldsymbol{e}_{z}+\unicode[STIX]{x1D70C}\int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|\leqslant a}\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}\,\text{d}V^{\prime },\end{eqnarray}$$

recalling that $\boldsymbol{n}$ points in to the centre of the particle. Note that the most significant component of $\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}$ is $-\bar{u}_{\unicode[STIX]{x1D703}}^{2}\boldsymbol{e}_{r}/(R+r)$ . Thus the contribution from the background flow $\bar{p},\bar{\boldsymbol{u}}$ to $\boldsymbol{F}^{\prime }$ is essentially just a buoyancy term plus a centrifugal force to counter the (rotational) inertia of the fluid that would otherwise take up the volume occupied by the particle. For a neutrally buoyant particle (i.e. $\unicode[STIX]{x1D70C}_{p}=\unicode[STIX]{x1D70C}$ ) travelling with a velocity close to that of the axial velocity of the background flow we would then expect $m_{p}\unicode[STIX]{x1D6E9}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{x}_{p}^{\prime }))$ to be close to the inertia of the displaced fluid so that the difference between the two is negligible. If the particle has significantly lower or higher density than the fluid then the net effect of the two terms is a force directed towards the inside or outside wall, respectively. The contribution of $q,\boldsymbol{v}^{\prime }$ to $\boldsymbol{F}^{\prime }$ in (2.21a ) consists of both the viscous forces on the particle from the surrounding fluid (i.e. essentially fluid drag) and an inertial contribution (which will be the inertial lift force).

Similarly, for the background contribution to the torque, it may be shown that

(2.23) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|=a}(\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime })\times ((-\boldsymbol{n})\boldsymbol{\cdot }(-\bar{p}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}+{\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}}}^{\intercal })))\,\text{d}S^{\prime }\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D70C}\int _{|\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime }|\leqslant a}(\boldsymbol{x}^{\prime }-\boldsymbol{x}_{p}^{\prime })\times (\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\bar{\boldsymbol{u}})\,\text{d}V^{\prime }.\end{eqnarray}$$

At first glance this may appear like a misuse of the divergence theorem; however, it indeed holds by first re-arranging the order of the cross and dot products so that the integrand has the form $-\boldsymbol{n}\boldsymbol{\cdot }(\ast )$ , after which the divergence theorem can be applied and further re-arrangement and use of (2.2a ) leads to (2.23). In any case this term is sufficiently small that it can be neglected for the purpose of estimating the inertial lift force to leading order, as will be demonstrated upon non-dimensionalising our equations.

2.4 Non-dimensionalisation of the disturbance flow

Let us now non-dimensionalise the variables and equations of motion (working in the rotating reference frame). Spatial variables are non-dimensionalised using the radius of the particle $a$ and velocity variables are non-dimensionalised via the characteristic velocity $U:=U_{m}a/\ell$ , where $U_{m}$ denotes a characteristic velocity of the background flow $\bar{\boldsymbol{u}}$ , which we take to be the value of its axial component at $r=\unicode[STIX]{x1D703}=z=0$ . This is essentially the same scale as in Hood et al. (Reference Hood, Lee and Roper2015). The non-dimensionalisation of the remaining variables follows from these via the standard approach for flows in which viscous stresses are dominant. Thus we introduce dimensionless variables denoted by  $\hat{\cdot }$ ,

(2.24) $$\begin{eqnarray}\left.\begin{array}{@{}rclrclrcl@{}}\boldsymbol{x}^{\prime }\ & =\ & a\hat{\boldsymbol{x}}^{\prime },\quad & \boldsymbol{v}^{\prime }\ & =\ & (U_{m}a/\ell )\hat{\boldsymbol{v}}^{\prime },\quad & q^{\prime }\ & =\ & (\unicode[STIX]{x1D707}U_{m}/\ell )\hat{q}^{\prime },\\ \boldsymbol{x}_{p}^{\prime }\ & =\ & a\hat{\boldsymbol{x}}_{p}^{\prime },\quad & \bar{\boldsymbol{u}}\ & =\ & (U_{m}a/\ell )\hat{\bar{\boldsymbol{u}}},\quad & \bar{p}\ & =\ & (\unicode[STIX]{x1D707}U_{m}/\ell )\hat{\bar{p}},\\ t\ & =\ & (\ell /U_{m})\hat{t},\quad & \unicode[STIX]{x1D734}_{p}^{\prime }\ & =\ & (U_{m}/\ell )\hat{\unicode[STIX]{x1D734}}_{p}^{\prime },\quad & \unicode[STIX]{x1D6E9}\ & =\ & (U_{m}/\ell )\hat{\unicode[STIX]{x1D6E9}},\\ \unicode[STIX]{x1D735}^{\prime }\ & =\ & a^{-1}\hat{\unicode[STIX]{x1D735}}^{\prime },\quad & g\ & =\ & (U_{m}^{2}a/\ell ^{2}){\hat{g}}.\quad & \ & \ & \end{array}\right\}\end{eqnarray}$$

With these scalings, the dimensionless form of (2.20) is

(2.25a ) $$\begin{eqnarray}\displaystyle -\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{q}^{\prime }+\hat{\unicode[STIX]{x1D6FB}}^{\prime \,2}\hat{\boldsymbol{v}}^{\prime } & = & \displaystyle Re_{p} ((\hat{\boldsymbol{v}}^{\prime }+\hat{\bar{\boldsymbol{u}}}-\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}^{\prime }))\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\boldsymbol{v}}^{\prime }\qquad \quad \nonumber\\ \displaystyle & & \displaystyle +\,\hat{\boldsymbol{v}}^{\prime }\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}+\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{v}}^{\prime })\!)\hspace{47.50006pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \hat{{\mathcal{F}}}^{\prime },\qquad \quad\end{eqnarray}$$

(2.25b ) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{\cdot }\hat{\boldsymbol{v}}^{\prime } & = & \displaystyle 0\hspace{145.00021pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \hat{{\mathcal{F}}}^{\prime },\qquad \quad\end{eqnarray}$$

(2.25c ) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{v}}^{\prime } & = & \displaystyle \mathbf{0}\hspace{145.00021pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{D}}},\qquad \quad\end{eqnarray}$$

(2.25d ) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{v}}^{\prime } & = & \displaystyle -\hat{\bar{\boldsymbol{u}}}+\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}^{\prime })+\hat{\unicode[STIX]{x1D734}}_{p}^{\prime }\times (\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime })\quad \text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{F}}}^{\prime }\backslash \unicode[STIX]{x2202}\hat{{\mathcal{D}}},\qquad \quad\end{eqnarray}$$

$\hat{{\mathcal{F}}}^{\prime }:=\{\hat{\boldsymbol{x}}^{\prime }:a\hat{\boldsymbol{x}}^{\prime }\in {\mathcal{F}}^{\prime }\}$

$\hat{{\mathcal{D}}}^{\prime }:=\{\hat{\boldsymbol{x}}^{\prime }:a\hat{\boldsymbol{x}}^{\prime }\in {\mathcal{D}}^{\prime }\}$

(2.26) $$\begin{eqnarray}Re_{p}:=\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}Ua=\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}U_{m}\frac{a^{2}}{\ell },\end{eqnarray}$$

is the particle Reynolds number. Writing the duct Reynolds number, $Re:=(\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D707})U_{m}\ell$ , we have $Re_{p}=Re(a/\ell )^{2}$ , so that the particle Reynolds number can be small even if the duct Reynolds number is not (note that $a<\ell /2$ is necessary for the particle to fit in the duct and $a\lesssim \ell /10$ is typical).

The forces on the particle in the rotating frame must also be non-dimensionalised. Noting that the dimensional force $\boldsymbol{F}^{\prime }$ and torque $\boldsymbol{T}^{\prime }$ are

(2.27a,b ) $$\begin{eqnarray}\boldsymbol{F}^{\prime }=m_{p}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{x}_{p}^{\prime }}{\unicode[STIX]{x2202}t^{2}}\quad \text{and}\quad \boldsymbol{T}^{\prime }=I_{p}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D734}_{p}^{\prime }}{\unicode[STIX]{x2202}t},\end{eqnarray}$$

we introduce the dimensionless quantities $\hat{\boldsymbol{F}}^{\prime },\hat{\boldsymbol{T}}^{\prime }$ defined by

(2.28a,b ) $$\begin{eqnarray}\hat{\boldsymbol{F}}^{\prime }:=\frac{\ell ^{2}}{\unicode[STIX]{x1D70C}U_{m}^{2}a^{4}}\boldsymbol{F}^{\prime }=\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}}\frac{4\unicode[STIX]{x03C0}}{3}\frac{\unicode[STIX]{x2202}^{2}\hat{\boldsymbol{x}}_{p}^{\prime }}{\unicode[STIX]{x2202}\hat{t}^{2}}\quad \text{and}\quad \hat{\boldsymbol{T}}^{\prime }:=\frac{\ell ^{2}}{\unicode[STIX]{x1D70C}U_{m}^{2}a^{5}}\boldsymbol{T}^{\prime }=\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}}\frac{8\unicode[STIX]{x03C0}}{15}\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D734}}_{p}^{\prime }}{\unicode[STIX]{x2202}\hat{t}}.\end{eqnarray}$$

Noting that $\unicode[STIX]{x1D70C}U_{m}^{2}a^{4}/\ell ^{2}=Re_{p}\unicode[STIX]{x1D707}U_{m}a^{2}/\ell$ , one has from (2.21a )

(2.29) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{F}}^{\prime } & = & \displaystyle -\frac{\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}}\frac{4\unicode[STIX]{x03C0}}{3}{\hat{g}}\boldsymbol{e}_{z}-\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}}\frac{4\unicode[STIX]{x03C0}}{3}\hat{\unicode[STIX]{x1D6E9}}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}_{p}^{\prime }))+\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|<1}\hat{\bar{\boldsymbol{u}}}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}\,\text{d}\hat{V}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{Re_{p}}\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(-\boldsymbol{n})\boldsymbol{\cdot }(-\hat{q}^{\prime }\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\boldsymbol{v}}^{\prime }+\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\boldsymbol{v}}^{\prime \intercal })\,\text{d}{\hat{S}}^{\prime },\end{eqnarray}$$

and similarly from (2.21b )

(2.30) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{T}}^{\prime } & = & \displaystyle -\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}}\frac{8\unicode[STIX]{x03C0}}{15}\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \hat{\unicode[STIX]{x1D734}}_{p}^{\prime })+\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|<1}(\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime })\times (\hat{\bar{\boldsymbol{u}}}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}})\,\text{d}\hat{V}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{Re_{p}}\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime })\times ((-\boldsymbol{n})\boldsymbol{\cdot }(-\hat{q}^{\prime }\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\boldsymbol{v}}^{\prime }+{\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\boldsymbol{v}}^{\prime }}^{\intercal }))\,\text{d}{\hat{S}}^{\prime }.\end{eqnarray}$$

At this stage, given $a,\hat{\boldsymbol{x}}_{p}^{\prime },\hat{\unicode[STIX]{x1D6E9}},\hat{\unicode[STIX]{x1D734}}_{p}^{\prime }$ and an approximation of $\hat{\bar{\boldsymbol{u}}}$ (for a desired duct shape/size and flow rate), one could solve the nonlinear problem (2.25) and then determine the resulting force and torque on the particle via (2.29) and (2.30), respectively. In the absence of a secondary component of the background flow the inertial lift force only becomes the dominant force once the particle has reached ‘terminal’ velocity and spin in relation to its main axial motion. As such, it is typical to determine $\hat{\unicode[STIX]{x1D6E9}},\hat{\unicode[STIX]{x1D734}}_{p}^{\prime }$ such that $\hat{\boldsymbol{T}}^{\prime }$ and the axial component of $\hat{\boldsymbol{F}}^{\prime }$ (i.e. $\hat{\boldsymbol{F}}^{\prime }\boldsymbol{\cdot }\boldsymbol{e}_{y^{\prime }}$ ) are negligible through an iterative process. In the simpler case of flow through a straight duct, the remaining components of $\hat{\boldsymbol{F}}^{\prime }$ provide the inertial lift force which leads to drift/migration within the cross-section. However, for curved ducts the secondary component of the background flow results in $\hat{\boldsymbol{F}}^{\prime }$ being influenced by the drag force from the secondary flow motion. A location in the cross-section where the net force from both effects is zero is an equilibrium to/from which a particle may be attracted/repelled, depending on the eigenvalues of the Jacobian of the net force (as a function of the particle position within the cross-section – that is, $\hat{\boldsymbol{x}}_{p}^{\prime }$ ).

In the case of a straight duct, the inertial lift force can be approximated from a perturbation expansion of the disturbance flow for small $Re_{p}$ (Hood et al. Reference Hood, Lee and Roper2015). Through this expansion, the nonlinear equation (2.25) can be broken up into a sequence of linear Stokes problems that are more easily solved. Furthermore, the resulting decomposition informs the way in which inertial lift influences the particle motion, particularly in our case where the problem is further complicated by the secondary fluid motion. The remainder of this section describes how the approach of Hood et al. (Reference Hood, Lee and Roper2015) is adopted to the curved duct geometry and, additionally, how separating the background flow into axial and secondary components allows us to separate and identify the terms that influence particle motion within the cross-section.

2.5 Perturbation expansion of the disturbance flow and forces on the particle

We consider a perturbation expansion of the disturbance flow in powers of the particle Reynolds number:

(2.31a,b ) $$\begin{eqnarray}\hat{\boldsymbol{v}}^{\prime }=\boldsymbol{v}_{0}+Re_{p}\boldsymbol{v}_{1}+O(Re_{p}^{2}),\quad \hat{q}^{\prime }=q_{0}+Re_{p}q_{1}+O(Re_{p}^{2}).\end{eqnarray}$$

Note that carets are dropped from the perturbation variables since these will always be taken to be dimensionless quantities. Substituting (2.31) into (2.25) and matching powers of $Re_{p}$ , one obtains the zeroth-order equations

(2.32a ) $$\begin{eqnarray}\displaystyle -\hat{\unicode[STIX]{x1D735}}^{\prime }q_{0}+\hat{\unicode[STIX]{x1D6FB}}^{\prime \,2}\boldsymbol{v}_{0} & = & \displaystyle \mathbf{0},\hspace{136.00026pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \hat{{\mathcal{F}}}^{\prime },\end{eqnarray}$$

(2.32b ) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}_{0} & = & \displaystyle 0,\hspace{136.00026pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \hat{{\mathcal{F}}}^{\prime },\end{eqnarray}$$

(2.32c ) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{0} & = & \displaystyle \mathbf{0},\hspace{136.00026pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{D}}}^{\prime },\end{eqnarray}$$

(2.32d ) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{0} & = & \displaystyle \hat{\unicode[STIX]{x1D734}}_{p}^{\prime }\times (\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime })+\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}^{\prime })-\hat{\bar{\boldsymbol{u}}},\quad \text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{F}}}^{\prime }\backslash \unicode[STIX]{x2202}\hat{{\mathcal{D}}}^{\prime },\qquad \quad\end{eqnarray}$$

(2.33a ) $$\begin{eqnarray}\displaystyle -\hat{\unicode[STIX]{x1D735}}^{\prime }q_{1}+\hat{\unicode[STIX]{x1D6FB}}^{\prime \,2}\boldsymbol{v}_{1} & = & \displaystyle \hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{v}_{0})+\boldsymbol{v}_{0}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}\qquad \nonumber\\ \displaystyle & & \displaystyle +\,(\boldsymbol{v}_{0}+\hat{\bar{\boldsymbol{u}}}-\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}^{\prime }))\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0},\quad \text{on }\hat{\boldsymbol{x}}^{\prime }\in \hat{{\mathcal{F}}}^{\prime },\qquad\end{eqnarray}$$

(2.33b ) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}_{1} & = & \displaystyle 0,\hspace{130.0002pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \hat{{\mathcal{F}}}^{\prime },\qquad\end{eqnarray}$$

(2.33c ) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{1} & = & \displaystyle \mathbf{0},\hspace{130.0002pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{D}}}^{\prime },\qquad\end{eqnarray}$$

(2.33d ) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{1} & = & \displaystyle \mathbf{0},\hspace{130.0002pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{F}}}^{\prime }\backslash \unicode[STIX]{x2202}\hat{{\mathcal{D}}}^{\prime }.\qquad\end{eqnarray}$$

Equations (2.32) and (2.33) are similar to those obtained for the straight duct (see Hood et al. (Reference Hood, Lee and Roper2015)) but with a few key differences. First, the right-hand side of (2.33a ) has additional terms introduced in the change to a rotating coordinate system. Second, the background velocity $\hat{\bar{\boldsymbol{u}}}$ is no longer a simple Poiseuille flow through the duct. Last, the boundary condition (2.25d ) is captured entirely in (2.32d ) so that $\boldsymbol{v}_{1}=\mathbf{0}$ in (2.33d ).

The perturbation expansion (2.31) can also be substituted into (2.29) and (2.30) to obtain the expansion

(2.34a ) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{F}}^{\prime } & = & \displaystyle Re_{p}^{-1}\boldsymbol{F}_{-1}+\boldsymbol{F}_{0}+O(Re_{p}),\end{eqnarray}$$

(2.34b ) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{T}}^{\prime } & = & \displaystyle Re_{p}^{-1}\boldsymbol{T}_{-1}+\boldsymbol{T}_{0}+O(Re_{p}),\end{eqnarray}$$

(2.35a ) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{-1} & := & \displaystyle \int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(-\boldsymbol{n})\boldsymbol{\cdot }\left(-q_{0}\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0}+{\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0}}^{\intercal }\right)\,\text{d}{\hat{S}}^{\prime },\end{eqnarray}$$

(2.35b ) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{0} & := & \displaystyle -\frac{\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}}\frac{4\unicode[STIX]{x03C0}}{3}{\hat{g}}\boldsymbol{e}_{z^{\prime }}-\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}}\frac{4\unicode[STIX]{x03C0}}{3}\hat{\unicode[STIX]{x1D6E9}}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}_{p}^{\prime }))+\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|<1}\hat{\bar{\boldsymbol{u}}}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}\,d\hat{V}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(-\boldsymbol{n})\boldsymbol{\cdot }\left(-q_{1}\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{1}+{\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{1}}^{\intercal }\right)\,\text{d}{\hat{S}}^{\prime },\end{eqnarray}$$

(2.36a ) $$\begin{eqnarray}\displaystyle \boldsymbol{T}_{-1} & := & \displaystyle \int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime })\times ((-\boldsymbol{n})\boldsymbol{\cdot }(-q_{0}\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0}+{\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0}}^{\intercal }))\,\text{d}{\hat{S}}^{\prime },\end{eqnarray}$$

(2.36b ) $$\begin{eqnarray}\displaystyle \boldsymbol{T}_{0} & := & \displaystyle -\frac{\unicode[STIX]{x1D70C}_{p}}{\unicode[STIX]{x1D70C}}\frac{8\unicode[STIX]{x03C0}}{15}\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \hat{\unicode[STIX]{x1D734}}_{p}^{\prime })+\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|<1}(\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime })\times (\hat{\bar{\boldsymbol{u}}}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}})\,\text{d}\hat{V}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime })\times ((-\boldsymbol{n})\boldsymbol{\cdot }(-q_{1}\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{1}+{\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{1}}^{\intercal }))\,\text{d}{\hat{S}}^{\prime }.\end{eqnarray}$$

In the case of a straight duct Hood et al. (Reference Hood, Lee and Roper2015) showed that the integral over the stress from $q_{1},\boldsymbol{v}_{1}$ can be computed without explicitly solving for these terms by utilising the Lorentz reciprocal theorem. The same result holds in the case of a curved duct; specifically, it is straightforward to show

(2.37) $$\begin{eqnarray}\displaystyle & & \displaystyle \boldsymbol{e}_{\ast }\boldsymbol{\cdot }\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(-\boldsymbol{n})\boldsymbol{\cdot }(-q_{1}\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{1}+{\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{1}}^{\intercal })\,\text{d}{\hat{S}}^{\prime }\nonumber\\ \displaystyle & & \displaystyle \quad =-\int _{\hat{{\mathcal{F}}}^{\prime }}\hat{\boldsymbol{u}}_{\ast }\boldsymbol{\cdot }(\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{v}_{0})+\boldsymbol{v}_{0}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}+(\boldsymbol{v}_{0}+\hat{\bar{\boldsymbol{u}}}-\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}^{\prime }))\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0})\,\text{d}\hat{V}^{\prime },\qquad\end{eqnarray}$$

for $\ast =x^{\prime },y^{\prime },z^{\prime }$ , where each velocity $\hat{\boldsymbol{u}}_{\ast }$ , along with a corresponding pressure $\hat{p}_{\ast }$ , solves

(2.38a,b ) $$\begin{eqnarray}-\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{p}_{\ast }+\hat{\unicode[STIX]{x1D6FB}}^{\prime \,2}\hat{\boldsymbol{u}}_{\ast }=\mathbf{0}\quad \text{on }\hat{\boldsymbol{x}}^{\prime }\in \hat{{\mathcal{F}}}^{\prime },\quad \hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{\cdot }\hat{\boldsymbol{u}}_{\ast }=0\quad \text{on }\hat{\boldsymbol{x}}^{\prime }\in \hat{{\mathcal{F}}}^{\prime },\end{eqnarray}$$

(2.38c,d ) $$\begin{eqnarray}\hspace{83.00015pt}\hat{\boldsymbol{u}}_{\ast }=\mathbf{0}\quad \text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{D}}}^{\prime },\hspace{20.00003pt}\hat{\boldsymbol{u}}_{\ast }=\boldsymbol{e}_{\ast }\hspace{8.00003pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{F}}}^{\prime }\backslash \unicode[STIX]{x2202}\hat{{\mathcal{D}}}^{\prime }.\end{eqnarray}$$

Note that we only need to consider the $x^{\prime },z^{\prime }$ components since the $y^{\prime }$ component, which perturbs the particle motion in the $y^{\prime }$ direction at $O(1)$ , has no effect on the cross-section forces at the same order of magnitude. This approach cuts down on computation time (since each $\hat{p}_{\ast },\hat{\boldsymbol{u}}_{\ast }$ is computed to estimate drag coefficients, see § 2.7), and also provides opportunities to examine how different components of the flow influence $\boldsymbol{F}_{0}$ in more detail (as discussed in § 2.6). Additionally, note that $\boldsymbol{T}_{0}$ can be ignored since perturbations of the particle spin have no influence on the (linear) forcing to this order of approximation. Given the boundary condition (2.32d ), $\boldsymbol{F}_{-1},\boldsymbol{T}_{-1}$ include drag and torque components that come from the secondary fluid motion. In order to separate the effect of the axial and secondary flow components and better understand how each influences the particle motion it is convenient to break up $\hat{\bar{\boldsymbol{u}}},q_{0},\boldsymbol{v}_{0}$ and (2.32) into two distinct parts. This will be done in the following section.

2.6 Separation of axial and secondary flow effects

Let

(2.39) $$\begin{eqnarray}\hat{\bar{\boldsymbol{u}}}=\hat{\bar{\boldsymbol{u}}}_{a}+\hat{\bar{\boldsymbol{u}}}_{s},\end{eqnarray}$$

where

(2.40) $$\begin{eqnarray}\displaystyle \hat{\bar{\boldsymbol{u}}}_{a} & := & \displaystyle \hat{\bar{u}}_{\unicode[STIX]{x1D703}}(r^{\prime },z^{\prime })\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}^{\prime }}=\frac{1}{U_{m}a/\ell }\bar{u}_{\unicode[STIX]{x1D703}}(r^{\prime },z^{\prime })\hat{\boldsymbol{e}}_{\unicode[STIX]{x1D703}^{\prime }},\end{eqnarray}$$

(2.41) $$\begin{eqnarray}\displaystyle \hat{\bar{\boldsymbol{u}}}_{s} & := & \displaystyle \hat{\bar{u}}_{r}(r^{\prime },z^{\prime })\hat{\boldsymbol{e}}_{r^{\prime }}+\hat{\bar{u}}_{z}(r^{\prime },z^{\prime })\hat{\boldsymbol{e}}_{z^{\prime }}=\frac{1}{U_{m}a/\ell }(\bar{u}_{r}(r^{\prime },z^{\prime })\hat{\boldsymbol{e}}_{r^{\prime }}+\bar{u}_{z}(r^{\prime },z^{\prime })\hat{\boldsymbol{e}}_{z^{\prime }}),\end{eqnarray}$$

denote the axial and secondary components, respectively. Similarly, we separate $q_{0},\boldsymbol{v}_{0}$ into two parts: specifically

(2.42a,b ) $$\begin{eqnarray}\boldsymbol{v}_{0}=\boldsymbol{v}_{0,a}+\boldsymbol{v}_{0,s},\quad q_{0}=q_{0,a}+q_{0,s},\end{eqnarray}$$

such that the pairs $q_{0,a},\boldsymbol{v}_{0,a}$ and $q_{0,s},\boldsymbol{v}_{0,s}$ each solve (2.32) except with the boundary condition (2.32d ) replaced by

(2.43a ) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{0,a} & = & \displaystyle \hat{\unicode[STIX]{x1D734}}_{p}^{\prime }\times (\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime })+\hat{\unicode[STIX]{x1D6E9}}(\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}^{\prime })-\hat{\bar{\boldsymbol{u}}}_{a},\quad \text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{F}}}^{\prime }\backslash \unicode[STIX]{x2202}\hat{{\mathcal{D}}}^{\prime },\end{eqnarray}$$

(2.43b ) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{0,s} & = & \displaystyle -\hat{\bar{\boldsymbol{u}}}_{s},\hspace{129.00012pt}\text{on }\hat{\boldsymbol{x}}^{\prime }\in \unicode[STIX]{x2202}\hat{{\mathcal{F}}}^{\prime }\backslash \unicode[STIX]{x2202}\hat{{\mathcal{D}}}^{\prime },\end{eqnarray}$$

$\boldsymbol{F}_{-1}=\boldsymbol{F}_{-1,a}+\boldsymbol{F}_{-1,s}$

$\boldsymbol{T}_{-1}=\boldsymbol{T}_{-1,a}+\boldsymbol{T}_{-1,s}$

$\ast =a,s$

(2.44a ) $$\begin{eqnarray}\displaystyle \boldsymbol{F}_{-1,\ast } & = & \displaystyle \int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(-\boldsymbol{n})\boldsymbol{\cdot }(-q_{0,\ast }\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0,\ast }+{\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0,\ast }}^{\intercal })\,\text{d}{\hat{S}}^{\prime },\end{eqnarray}$$

(2.44b ) $$\begin{eqnarray}\displaystyle \boldsymbol{T}_{-1,\ast } & = & \displaystyle \int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime })\times ((-\boldsymbol{n})\boldsymbol{\cdot }(-q_{0,\ast }\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0,\ast }+{\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0,\ast }}^{\intercal }))\,\text{d}{\hat{S}}^{\prime }.\end{eqnarray}$$

This separation of axial and secondary flow effects yields a clear methodology for solving the leading-order disturbance equations. One first solves for $q_{0,a},\boldsymbol{v}_{0,a}$ and in doing so determines $\hat{\unicode[STIX]{x1D6E9}},\hat{\unicode[STIX]{x1D734}}_{p}^{\prime }$ such that $\boldsymbol{F}_{-1,a}=\boldsymbol{T}_{-1,a}=\mathbf{0}$ . The $q_{0,s},\boldsymbol{v}_{0,s}$ components are then solved independently, from which $\boldsymbol{F}_{-1,s},\boldsymbol{T}_{-1,s}$ are calculated to provide the drag and torque on the particle from the secondary fluid motion. The drag term $\boldsymbol{F}_{-1,s}$ can then be added to $\boldsymbol{F}_{0}$ to obtain the net perturbing forces on the particle (up to $O(1)$ ). While this drag from the secondary flow motion features as an $O(Re_{p}^{-1})$ factor it is typical that the magnitude of $\hat{\bar{\boldsymbol{u}}}_{s}$ is sufficiently small that the resulting force is comparable to $\boldsymbol{F}_{0}$ (since if this was not the case it is unlikely that focusing of particles would be observed since they would instead follow the streamlines of the secondary flow).

The $\boldsymbol{F}_{0}$ component of the force can also be expanded based on the separation of axial and secondary background flow components. The volume integral over $\hat{\bar{\boldsymbol{u}}}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}$ in (2.35b ) can be expanded as

(2.45) $$\begin{eqnarray}\displaystyle \int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|<1}\hat{\bar{\boldsymbol{u}}}_{a}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}_{a}\,\text{d}\hat{V}^{\prime }+\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|<1}\hat{\bar{\boldsymbol{u}}}_{s}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}_{a}+\hat{\bar{\boldsymbol{u}}}_{a}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}_{s}\,\text{d}\hat{V}^{\prime }+\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|<1}\hat{\bar{\boldsymbol{u}}}_{s}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}_{s}\,\text{d}\hat{V}^{\prime }, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where the $y^{\prime },z^{\prime }$ components of the first integral are exactly zero (and the $x^{\prime }$ component is comparable to $(4\unicode[STIX]{x03C0}/3)\hat{\unicode[STIX]{x1D6E9}}^{2}(\boldsymbol{e}_{z^{\prime }}\times (\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}_{p}^{\prime }))$ for a small neutrally buoyant particle), the second integral consists of only a small $y^{\prime }$ component and the third integral provides a small force in the $x^{\prime },z^{\prime }$ directions. Utilising (2.37), the $\ast =x^{\prime },z^{\prime }$ components of the surface integral over the fluid stress from $q_{1},\boldsymbol{v}_{1}$ in (2.35b ) can each be expanded as

(2.46a ) $$\begin{eqnarray}\displaystyle & & \displaystyle -\int _{\hat{{\mathcal{F}}}^{\prime }}\hat{\boldsymbol{u}}_{\ast }\boldsymbol{\cdot }(\hat{\unicode[STIX]{x1D6E9}}\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{v}_{0,a}+\boldsymbol{v}_{0,a}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}_{a}+(\boldsymbol{v}_{0,a}+\hat{\bar{\boldsymbol{u}}}_{a}-\hat{\unicode[STIX]{x1D6E9}}\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}^{\prime })\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0,a})\,\text{d}\hat{V}^{\prime }\qquad \quad\end{eqnarray}$$

(2.46b ) $$\begin{eqnarray}\displaystyle & & \displaystyle -\int _{\hat{{\mathcal{F}}}^{\prime }}\hat{\boldsymbol{u}}_{\ast }\boldsymbol{\cdot }(\boldsymbol{v}_{0,s}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}_{s}+(\boldsymbol{v}_{0,s}+\hat{\bar{\boldsymbol{u}}}_{s})\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0,s})\,\text{d}\hat{V}^{\prime }\end{eqnarray}$$

(2.46c ) $$\begin{eqnarray}\displaystyle & & \displaystyle -\int _{\hat{{\mathcal{F}}}^{\prime }}\hat{\boldsymbol{u}}_{\ast }\boldsymbol{\cdot }(\!\hat{\unicode[STIX]{x1D6E9}}\boldsymbol{e}_{z^{\prime }}\times \boldsymbol{v}_{0,s}-(\hat{\unicode[STIX]{x1D6E9}}\boldsymbol{e}_{z^{\prime }}\times \hat{\boldsymbol{x}}^{\prime })\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0,s}+\boldsymbol{v}_{0,a}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}_{s}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\boldsymbol{v}_{0,s}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\bar{\boldsymbol{u}}}_{a}+(\boldsymbol{v}_{0,a}+\hat{\bar{\boldsymbol{u}}}_{a})\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0,s}+(\boldsymbol{v}_{0,s}+\hat{\bar{\boldsymbol{u}}}_{s})\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}}^{\prime }\boldsymbol{v}_{0,a}\! )\,\text{d}\hat{V}^{\prime }.\end{eqnarray}$$

$x^{\prime },z^{\prime }$

$y^{\prime }$

$y^{\prime }$

$\boldsymbol{F}_{0}$

In summary, the force on the particle within the cross-section up to $O(1)$ is given by

(2.47) $$\begin{eqnarray}\boldsymbol{F}_{p}^{\prime }:=(\boldsymbol{e}_{x^{\prime }}\boldsymbol{\cdot }(Re_{p}^{-1}\boldsymbol{F}_{-1,s}+\boldsymbol{F}_{0}))\boldsymbol{e}_{x^{\prime }}+(\boldsymbol{e}_{z^{\prime }}\boldsymbol{\cdot }(Re_{p}^{-1}\boldsymbol{F}_{-1,s}+\boldsymbol{F}_{0}))\boldsymbol{e}_{z^{\prime }}.\end{eqnarray}$$

Given a fixed particle size and duct size/shape, then $\boldsymbol{F}_{p}^{\prime }$ can be considered as a function of the particle position within the cross-section (i.e. $\hat{x}_{p,r^{\prime }},\hat{x}_{p,z^{\prime }}$ , recalling that at each point each of $\hat{\unicode[STIX]{x1D6E9}},\hat{\unicode[STIX]{x1D734}}_{p}^{\prime }$ are chosen such that $\boldsymbol{F}_{-1,a}=\boldsymbol{T}_{-1,a}=\mathbf{0}$ ).

2.7 A first-order model of particle trajectories

Here we consider a simple model for the trajectory of a particle within the cross-section based on the terminal velocity; that is, we take the particle velocity in the $r^{\prime },z^{\prime }$ directions to be such that the drag force is equal and opposite to the perturbing force $\boldsymbol{F}_{p}^{\prime }(\hat{x}_{p,r^{\prime }},\hat{x}_{p,z^{\prime }})$ . In order to use this model we must first determine the drag coefficients for a particle with non-zero velocity in the cross-sectional plane.

Suppose the spherical particle is moving with a velocity $\hat{u} _{p,r^{\prime }}^{\prime },\hat{u} _{p,z^{\prime }}^{\prime }$ in the $r^{\prime },z^{\prime }$ directions, respectively, then the cross-sectional drag force on the particle (non-dimensionalised with respect to the force scale $\unicode[STIX]{x1D70C}U_{m}^{2}a^{4}/\ell ^{2}$ ) can be expressed as $Re_{p}^{-1}(C_{r^{\prime }}\hat{u} _{p,r^{\prime }}\hat{\boldsymbol{e}}_{r^{\prime }}+C_{z^{\prime }}\hat{u} _{p,z^{\prime }}\hat{\boldsymbol{e}}_{z^{\prime }})$ , where $C_{r^{\prime }},C_{z^{\prime }}$ are (dimensionless) drag coefficients. For an asymptotically small particle away from the walls, the (dimensionless) drag coefficients could be estimated via Stokes drag law as simply $6\unicode[STIX]{x03C0}$ . However, in the context of a microfluidic duct, the finite particle size and proximity of the duct walls have an effect such that the true drag coefficients are noticeably larger than this. A better estimate of the drag coefficients is obtained by explicitly computing the drag force on the particle. Note that the $\hat{p}_{\ast },\hat{\boldsymbol{u}}_{\ast }$ from § 2.5 which each solve (2.38) describe the Stokes flow resulting from a particle moving through a stationary fluid in the duct with velocity $\boldsymbol{e}_{\ast }$ . Thus, we can re-use these solutions to compute the drag coefficients. Specifically, for each $C_{\ast }$ we need only integrate the resulting fluid stress over the surface of the particle, that is

(2.48) $$\begin{eqnarray}C_{\ast }=\hat{\boldsymbol{e}}_{\ast }\boldsymbol{\cdot }\int _{|\hat{\boldsymbol{x}}^{\prime }-\hat{\boldsymbol{x}}_{p}^{\prime }|=1}(-\boldsymbol{n})\boldsymbol{\cdot }(-\hat{p}_{\ast }\unicode[STIX]{x1D644}+\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\boldsymbol{u}}_{\ast }+{\hat{\unicode[STIX]{x1D735}}^{\prime }\hat{\boldsymbol{u}}_{\ast }}^{\intercal })\,\text{d}{\hat{S}}^{\prime }.\end{eqnarray}$$

Similar to $\boldsymbol{F}_{p}^{\prime }$ , given a fixed particle size and duct size/shape the drag coefficients $C_{r^{\prime }},C_{z^{\prime }}$ depend primarily on the particle position within the cross-section.

Therefore, with the addition of these drag terms, the cross-sectional force on the particle is

(2.49) $$\begin{eqnarray}(Re_{p}^{-1}C_{r^{\prime }}\hat{u} _{p,r^{\prime }}^{\prime }+\hat{F}_{p,r^{\prime }}^{\prime })\hat{\boldsymbol{e}}_{r^{\prime }}+(Re_{p}^{-1}C_{z^{\prime }}\hat{u} _{p,z^{\prime }}^{\prime }+\hat{F}_{p,z^{\prime }}^{\prime })\hat{\boldsymbol{e}}_{z^{\prime }}.\end{eqnarray}$$

Thus, equating this to zero and writing this out in terms of $\hat{x}_{p,r^{\prime }}^{\prime },\hat{x}_{p,z^{\prime }}^{\prime }$ , we obtain the first-order system of ordinary differential equations

(2.50a,b ) $$\begin{eqnarray}\frac{\text{d}\hat{x}_{p,r^{\prime }}^{\prime }}{\text{d}t^{\prime }}=-Re_{p}\frac{\hat{F}_{p,r^{\prime }}^{\prime }(\hat{x}_{p,r^{\prime }}^{\prime },\hat{x}_{p,z^{\prime }}^{\prime })}{C_{r^{\prime }}(\hat{x}_{p,r^{\prime }}^{\prime },\hat{x}_{p,z^{\prime }}^{\prime })},\quad \frac{\text{d}\hat{x}_{p,z^{\prime }}^{\prime }}{\text{d}t^{\prime }}=-Re_{p}\frac{\hat{F}_{p,z^{\prime }}^{\prime }(\hat{x}_{p,r^{\prime }}^{\prime },\hat{x}_{p,z^{\prime }}^{\prime })}{C_{z^{\prime }}(\hat{x}_{p,r^{\prime }}^{\prime },\hat{x}_{p,z^{\prime }}^{\prime })}.\end{eqnarray}$$

The solution of these coupled equations approximates the particle trajectory within the cross-section. Interpreting $\hat{\unicode[STIX]{x1D6E9}}$ as a function of $\hat{x}_{p,r^{\prime }}^{\prime },\hat{x}_{p,z^{\prime }}^{\prime }$ (recalling that $\hat{\unicode[STIX]{x1D6E9}},\hat{\unicode[STIX]{x1D734}}_{p}^{\prime }$ are chosen at each $(\hat{x}_{p,r^{\prime }}^{\prime },\hat{x}_{p,z^{\prime }}^{\prime })$ such that $\boldsymbol{F}_{-1,a}=\boldsymbol{T}_{-1,a}=\mathbf{0}$ ), the angular position of the particle can also be estimated via $\text{d}\unicode[STIX]{x1D703}_{p}/\text{d}t^{\prime }=\hat{\unicode[STIX]{x1D6E9}}(\hat{x}_{p,r^{\prime }}^{\prime },\hat{x}_{p,z^{\prime }}^{\prime })$ .

4.1 Ducts having a square cross-section

Figure 3. The force $\boldsymbol{F}_{p}^{\prime }$ on neutrally buoyant particles of radius $a$ as shown at different locations in a curved duct with square cross-section ( $W=H=2$ ) and bend radius $R=160$ ((a) $a=0.05$ , (b) $a=0.10$ , (c) $a=0.15$ , (d) $a=0.20$ ). The left wall is on the inside of the bend. For comparison, results for a straight square duct (obtained in the limit $R\rightarrow \infty$ ) are also included ((e) $a=0.05$ , (f) $a=0.10$ , (g) $a=0.15$ , (h) $a=0.20$ ). The colour background shows the magnitude of $\boldsymbol{F}_{p}^{\prime }$ . Black and white contours are the zero-level set curves of the horizontal and vertical components, respectively, whereas the arrows indicate the sign of each component in the area bounded by the respective contour. The dashed red line shows where the centre of the particle lies when its surface touches the wall.

Figure 4. Approximate trajectories of neutrally buoyant particles within a curved duct having a square cross-section ( $W=H=2$ ) with bend radius $R=160$ ((a)  $a=0.05$ , (b) $a=0.10$ , (c) $a=0.15$ , (d) $a=0.20$ ). The left wall is on the inside of the bend. For comparison, results for a straight square duct (obtained in the limit $R\rightarrow \infty$ ) are also included ((e) $a=0.05$ , (f) $a=0.10$ , (g) $a=0.15$ , (h) $a=0.20$ ). Trajectories from several starting positions have been superimposed. Green, orange and red markers show the location of stable, (unstable) saddle and unstable equilibria, respectively (with marker size indicative of particle size). Stable orbits are coloured green (where they exist). The dashed red line shows where the centre of the particle lies when its surface touches the wall.

For a duct having a square cross-section with $W=H=2$ the cross-sectional force and drag on particles is computed for the particle sizes $a\in \{0.20,0.15,0.10,0.05\}$ and bend radii $R\in \{640,320,160,80,40\}$ (or equivalently $\unicode[STIX]{x1D716}^{-1}\in \{640,320,160,80,40\}$ ). Figures 3 and 4 show the force and trajectory plots, respectively, for the specific bend radius $R=160$ . Additionally, to help illustrate the effect of duct curvature, results for a square duct (obtained as $R\rightarrow \infty$ ) are also included.

We start by discussing the results in the case of a curved duct. Results for the smallest particle ( $a=0.05$ ) are depicted in figures 3(a) and 4(a). Here $\boldsymbol{F}_{p}^{\prime }$ is dominated by the secondary flow drag since $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R=200$ is reasonably large. Three equilibria can be identified in the force plot: one near the right (outer) wall and two which are near the centre, horizontally and vertically symmetric about the $r$ axis. The equilibrium near the outer wall is an unstable saddle (it attracts horizontally but repels vertically). The remaining two are difficult to infer visually from the force plot but are also unstable (both eigenvalues of the Jacobian have positive real part, albeit quite small). The trajectory plot makes it clear that particles are strongly affected by the vortices of the secondary flow. It is notable the particles starting too close to the walls migrate inwards a little before becoming trapped in the vortex motion.

Results for the second smallest particle ( $a=0.10$ and $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R=25$ ) are depicted in figures 3(b) and 4(b). Three equilibria can again be identified from the force plot, each of which is unstable, similar to those in figures 3(a) and 4(a). Notice that the two symmetric equilibria are now closer to the inside wall. The trajectory plot again shows that the vortex motion of the secondary flow has a strong influence on the trajectories; however, in this case the particles appear to migrate onto one of two (symmetrical) stable orbits relatively quickly. This suggests the inertial lift force and secondary flow drag are similar in magnitude.

Results for the second largest particle ( $a=0.15$ and $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R\approx 7.41$ ) are depicted in figures 3(c) and 4(c). Three equilibria can be identified from the force plot, each centred vertically within the cross-section: one near the centre of the outside wall, one slightly left of centre and another near the centre of the inside wall. The equilibrium near the inside wall is stable while the remaining two are unstable (the one nearest to the outside being a saddle point). The trajectory plot illustrates that inertial lift force is the dominant effect with complete focusing onto the stable equilibrium. Particles initially migrate onto a ‘slow manifold’ along which they more slowly migrate towards the unique stable equilibrium.

Results for the largest particle ( $a=0.20$ and $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R=3.125$ ) are depicted in figures 3(d) and 4(d) and are similar to those in figures 3(c) and 4(c). The trajectory plot shows a more direct convergence onto the slow manifold, given the weaker effect of the secondary flow drag in this case. Migration along the slow manifold is quite slow, as evidenced by the close proximity of the black and white contours in the force plot. Note that for a particle which is a little larger, and/or a larger bend radius, we expect these contours to cross over, leading to additional equilibria. In particular, a stable equilibria may form near the centre of the outside wall, although relatively few particles would be expected to migrate towards it.

For other bend radii the results are qualitatively similar to those described above given a similar value of the ratio $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R$ . We briefly comment on some general trends. For very large $\unicode[STIX]{x1D705}$ , particles effectively follow the vortex motion of the background flow (apart from some initial migration away from the walls). In the limit $R\rightarrow \infty$ (or equivalently $\unicode[STIX]{x1D705}\rightarrow 0$ given fixed $a,\ell$ ) we recover the known behaviour for a straight duct (specifically, the four stable equilibria near the centre of each wall as in Hood et al. (Reference Hood, Lee and Roper2015)). It is notable that $R\gg 10^{3}$ is required before particle focusing similar to that in a straight duct is observed, because even a very small influence from secondary flow motion is enough to modify the stability of some equilibria.

We take a moment to compare with the results for a straight square duct. Figure 3(e–h) shows the inertial lift force (noting there is no rotational inertia or secondary flow drag). There is clear symmetry with respect to both axes, as one would expect, and only subtle changes with respect to particle size. Comparing with figure 3(a–d) it is clear that adding curvature to the duct has the most significant effect on the horizontal component of the force, as evident in the changes to the corresponding zero-level contour (black). Only when $\unicode[STIX]{x1D705}$ is relatively large do we start to notice significant changes to the zero-level contours of the vertical force component (in white). Figure 4(e–h) shows the trajectory plots for the straight square duct. Without curvature it is clear the behaviour does not change significantly for different size particles. The stability of the equilibria can be understood in terms of the relative positioning of the two zero-level contours in figure 3(e–h). Note that as the $a$ increases the quasi-circular parts of the two zero-level contours change in aspect ratio. It is expected this trend continues such that the two zero-level curves cross/switch, and once this occurs the stability of equilibria around the curves will change such that the corner points become stable and the side points become saddles. This supports previous work that demonstrates larger particles focus towards the corners of a square duct (Prohm & Stark Reference Prohm and Stark2014). A similar switching of stability may also occur in curved ducts when $\unicode[STIX]{x1D705}$ is sufficiently small (albeit for particles larger than considered here).

Figure 5. Approximate trajectories of neutrally buoyant particles with radius (a,b) $a=0.10$ and (c,d) $a=0.15$ within a curved duct with square cross-section ( $W=H=2$ ) having different bend radii $R$ as shown in each case. Trajectories from several starting positions have been superimposed. Green, orange and red markers show the location of stable, (unstable) saddle and unstable equilibria, respectively (with marker size indicative of particle size). Stable orbits are also coloured green (where they exist). The left wall is on the inside of the bend. The dashed red line shows where the centre of the particle lies when its surface touches the wall.

Returning to the case of curved ducts, the transition from non-focusing to focusing behaviour is quite interesting. As $R$ increases (and $\unicode[STIX]{x1D705}$ decreases), the particles begin to focus towards trapping orbits, as is the case observed for $a=0.10$ and $R=160$ . For the specific case of $a=0.10$ , as $R$ continues to increase the size of the trapping orbits grows until they eventually meet and break up (around $\unicode[STIX]{x1D705}\in [18,19]$ ) leaving behind the slow manifold and a stable equilibrium, as depicted in figure 5(a,b). Interestingly this is somewhat different for the slightly larger particle with $a=0.15$ , in which the trapping orbits instead shrink while migrating towards what becomes the stable equilibrium (around $\unicode[STIX]{x1D705}\in [13,15]$ ), as depicted in figure 5(c,d). At this stage the significance of the trapping orbits and their structure are unclear and warrant further investigation.

Notice that the larger particles appear to focus near the inside wall while smaller ones become trapped in orbits a finite distance away from the inside wall. This provides a potential mechanism for separating larger particles from smaller ones in a fluid sample that contains a low concentration of both. For example, by splitting the flow at the outlet into two streams which separate the inside $1/3$ and outside $2/3$ (measured laterally) one would expect the majority of larger particles to be collected from the inside stream (with a reduced portion of smaller particles). However, the orbits onto which smaller particles are trapped are perhaps too close to the focusing location of the larger particles in a square duct for this to provide robust and efficient separation in practice. Rectangular ducts of small aspect ratio are more common in experiments that demonstrate separation, so we move onto examine two such cross-sections.

4.2 Ducts having a rectangular cross-section

Two rectangular cross-sections are considered: specifically, with $H=2$ and $W\in \{4,8\}$ . For each we computed $\boldsymbol{F}_{p}^{\prime }$ for particles of size $a\in \{0.20,0.15,0.10,0.05\}$ and bend radii $R\in \{640,320,160,80\}$ (or equivalently $\unicode[STIX]{x1D716}^{-1}\in \{640,320,160,80\}$ ). To illustrate the effect of duct curvature the straight duct case ( $R\rightarrow \infty$ , or equivalently $\unicode[STIX]{x1D716}\rightarrow 0$ ) is also computed, with results provided in appendix C. The force and trajectory plots for $W=4$ and $R=160$ are shown in figures 6 and 7, respectively. Results for the smallest particle ( $a=0.05$ and $\ell ^{4}/4a^{3}R=200$ ) are depicted in figures 6(a) and 7(a). As in the case of the square duct, the effect of secondary flow drag is the dominant cross-sectional force. The force plot is similar to that in the square duct case for $a=0.05$ (only stretched widthwise) and, in particular, three equilibria can be identified. The equilibrium near the right wall is again an unstable saddle. However, unlike the square duct case, the remaining (symmetric) equilibria pair is stable in this case. This is evident in the trajectory plot, where particles are observed to slowly converge towards the equilibria despite the relatively strong vortex motion. Note that the location of the stable pair is slightly left of the centre. Thus, despite a strong influence from the secondary flow drag, the particles eventually becomes focused.

Figure 6. The force $\boldsymbol{F}_{p}^{\prime }$ on neutrally buoyant particles at different locations in a curved duct with rectangular cross-section ( $W/2=H=2$ ) and bend radius $R=160$ . The colour background shows the magnitude of $\boldsymbol{F}_{p}^{\prime }$ . Black and white contours are the zero level set curves of the horizontal and vertical components, respectively, whereas the arrows indicate the sign of each component in the area bounded by the respective contour. The left wall is on the inside of the bend. The dashed red line shows where the centre of the particle lies when its surface touches the wall.

Figure 7. Approximate trajectories of neutrally buoyant particles in a curved duct with rectangular cross-section ( $W/2=H=2$ ) and bend radius $R=160$ . Trajectories from several starting positions have been superimposed. Green, orange and red markers show the location of stable, (unstable) saddle and unstable equilibria, respectively (with marker size indicative of particle size). The left wall is on the inside of the bend. The dashed red line shows where the centre of the particle lies when its surface touches the wall.

Results for the second smallest particle ( $a=0.10$ and $\ell ^{4}/4a^{3}R=25$ ) are depicted in figures 6(b) and 7(b). Again, the force plot could be interpreted as a stretched version of the square duct case. The equilibrium near the right wall is an unstable saddle, whereas the remaining (symmetric) pair is stable. Notably the stable pair is much closer to the inside wall than in the case of the smallest particle. The trajectory plot demonstrates that although the secondary flow drag and inertial lift force are expected to have similar magnitude, the particle converges onto a slow manifold before completing a full orbit around an equilibrium, and then proceeds to migrate to the equilibrium. This too is quite different from the square duct case despite the seemingly similar structure of the force plots.

Results for the second largest particle ( $a=0.15$ and $\ell ^{4}/4a^{3}R\approx 7.41$ ) are depicted in figures 6(c) and 7(c). Three equilibria can again be identified, the unstable saddle near the outside wall and a stable pair. The stable pair is slightly closer to the centre than in the case of $a=0.10$ , but still closer to the inside wall than the case $a=0.05$ . The trajectory plot shows particles converge onto the slow manifold with minimal impact from the secondary flow motion before then migrating towards the stable equilibria.

Results for the largest particle ( $a=0.20$ and $\ell ^{4}/4a^{3}R=3.125$ ) are depicted in figures 6(d) and 7(d). Two new unstable equilibria can be identified in the force plot, located on the horizontal symmetry line on the left side. Additionally, the stable equilibria pair has shifted further towards the right. Inertial lift force is the dominant effect in this case and the trajectory plot shows a more direct convergence onto the slow manifold prior to convergence towards the stable equilibria.

While isolines of the total shear rate of $\bar{\boldsymbol{u}}$ have been found to coincide with the slow manifold in straight ducts (Lashgari et al. Reference Lashgari, Ardekani, Banerjee, Russom and Brandt2017), such a correlation does not exist in the case of curved ducts. Specifically, note that the slow manifold deforms and eventually breaks down with decreasing particle size (and all else fixed), which demonstrates that the focusing mechanism in curved ducts cannot be understood solely in terms of the total shear rate.

The results for different values of $R$ are similar to those described above, given a similar value of $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R$ , and so we remark on some general observations and trends. In the limit $R\rightarrow \infty$ (equivalently $\unicode[STIX]{x1D705}=0$ for fixed $a,\ell$ ) we recover the focusing behaviour of particles in straight ducts; specifically, particles will focus towards stable equilibria located a small distance from the centre of the two longer walls (see appendix C). For the smaller size particles one additionally finds stable equilibria near the centre of the sidewalls; however, very few particles are likely to be found here in practice because the attractive region is much smaller (Hood Reference Hood2016). Similar to the square duct case these generally only become evident for $R\gg 10^{3}$ , since even a small influence from the secondary flow is enough to modify/destroy these equilibria. On the other hand, when $\unicode[STIX]{x1D705}$ is large the effect of the secondary flow drag is dominant. Unlike the square duct case, we find that the effect of the inertial lift force in the rectangular duct is enough that the centres of the vortices are stable equilibria. However, it would be reasonable to expect this to break down for smaller particles (that is with $a<0.05$ ), as the relative influence of the secondary flow motion becomes more significant. In the two limits discussed we note the stable equilibria pair is close to the centre horizontally. For intermediate values of $\unicode[STIX]{x1D705}$ on the other hand we observe the focusing position to be much closer to the inside wall. The further left a particle can migrate varies with size, and is between $0.3H$ and $0.4H$ for the radii considered herein. Also noteworthy is that our results suggest that stable focusing positions can never occur in the right half of the cross-section. In contrast, some experiments have shown focusing can occur in this region; although we believe this is likely to be an effect of higher flow rates, and hence is not captured by the low flow rate approximation. In particular, at higher flow rates the location of the maximum of the axial flow and the centre of the secondary flow vortices are pushed towards the outer wall by inertia, and it is reasonable to expect particle behaviour to be perturbed accordingly.

Force and trajectory plots for select particle sizes within the higher aspect ratio rectangular cross-section with $W=8$ , $H=2$ and $R=160$ are shown in figures 8 and 9, respectively. Generally speaking, the behaviour in the wider duct is qualitatively similar to that observed in the case with $W=4$ . This is evident in comparing figures 8 and 9 with figures 6(a,c) and 7(a,c), respectively. The main advantage of the wider duct is a greater separation distance for particles focused towards the inside wall from those that focus towards the centre. Variations in focusing position with respect to particle size and bend radius are critical in determining how effective a microfluidic duct is able to separate particles. To better compare the differences for the two aspect ratios the focusing positions will be examined in detail.

Figure 8. The force $\boldsymbol{F}_{p}^{\prime }$ on neutrally buoyant particles at different locations in a curved duct with rectangular cross-section ( $W/4=H=2$ ) and bend radius $R=160$ . The colour background shows the magnitude of $\boldsymbol{F}_{p}^{\prime }$ . Black and white contours are the zero level set curves of the horizontal and vertical components, respectively, whereas the arrows indicate the sign of each component in the area bounded by the respective contour. The left wall is on the inside of the bend. The dashed red line shows where the centre of the particle lies when its surface touches the wall.

Figure 9. Approximate trajectories of neutrally buoyant particles in a curved duct with rectangular cross-section ( $W/4=H=2$ ) and bend radius $R=160$ . Trajectories from several different starting positions have been superimposed. Green and orange markers show the location of stable and (unstable) saddle equilibria, respectively (with marker size indicative of particle size). The left wall is on the inside of the bend. The dashed red line shows where the centre of the particle lies when its surface touches the wall.

To study how the focusing position of different size particles changes with respect to the bend radius, we interpolate the components that make up $\boldsymbol{F}_{p}^{\prime }$ in (2.47) between $R\in \{640,320,160,80\}$ (and extrapolate for $R$ outside this range). Since $\ell =2$ is fixed, changes in $R$ can be interpreted as changes to the dimensionless ratio $\unicode[STIX]{x1D716}^{-1}=2R/\ell$ . Despite the wide range of $R$ , each of $\bar{u}_{\unicode[STIX]{x1D703},0},\bar{u}_{r,0},\bar{u}_{z,0}$ (used to approximate $\hat{\bar{\boldsymbol{u}}}_{a},\hat{\bar{\boldsymbol{u}}}_{s}$ as in (3.3)) change only modestly. The most significant effect is a change in scale of the secondary flow drag relative to the inertial lift force, which is captured separately by the factor $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R$ ; changes in the remaining components of $\boldsymbol{F}_{p}^{\prime }$ are generally small and can therefore be interpolated/extrapolated with a reasonable degree of confidence. For a much finer sampling of $R$ we then reconstruct and analyse the total force on the particle $\boldsymbol{F}_{p}^{\prime }$ and identify/track the horizontal location of stable equilibria. The results of this analysis are shown in figure 10 for both rectangular ducts. The vertical axis shows the horizontal location of the focusing position with the inside (left) wall at the bottom and the centre at the top. The focusing position of particles undergoes a transition from being near the centre to being a finite distance away from the inside wall. Importantly, this differs significantly with particle size, thereby suggesting a mechanism for size-based separation of particles. However, note that the relative ordering of particles changes many times over the range of $\unicode[STIX]{x1D716}^{-1}$ shown, suggesting that some care should be taken in choosing the bend radius, depending on the size of particles to be separated.

Figure 10. Horizontal location of the (stable) focusing positions of particles versus $\unicode[STIX]{x1D716}^{-1}=2R/\ell$ for the two rectangular ducts. The particle sizes are $a=0.05$ (blue, circles), $a=0.10$ (green, triangles), $a=0.15$ (red, squares) and $a=0.20$ (cyan, diamonds). Note the vertical axis includes only the left half of the duct ( $0$ being the centre and $-2,-4$ being the inside wall in (a,b), respectively).

For the duct with $W/H=2$ and with $\unicode[STIX]{x1D716}^{-1}\lessapprox 80$ the ordering going from the inside wall to the centre is from largest to smallest particle. As $\unicode[STIX]{x1D716}^{-1}$ increases there are many changes to the ordering as larger particles begin to migrate towards the centre whereas smaller ones migrate towards the wall. Note that for $\unicode[STIX]{x1D716}^{-1}\gtrapprox 500$ an additional stable equilibrium appears near the inside wall for the particle with radius $a=0.10$ , and varies very little with respect to $\unicode[STIX]{x1D716}^{-1}$ . This is one of the extra equilibria that is expected to occur as the behaviour approaches that of a straight duct (see appendix C). For $\unicode[STIX]{x1D716}^{-1}\gtrapprox 640$ the ordering of the dominant equilibria is from smallest to largest particle. Similar trends are observed for the duct with $W/H=4$ . In particular, the ordering of the focusing location for $\unicode[STIX]{x1D716}^{-1}\lessapprox 80$ is the same, with the exception of the two larger particles being switched. There are then several changes to the ordering for increasing $\unicode[STIX]{x1D716}^{-1}$ until for $\unicode[STIX]{x1D716}^{-1}\gtrapprox 640$ the ordering is from smallest to largest particle (noting we again observe the emergence of an additional stable equilibrium for $a=0.15$ near the inside wall). One notable difference for the two different aspect ratios is that particles generally do not appear to get as close to the centre over the range of $\unicode[STIX]{x1D716}^{-1}$ shown.

Figure 11 shows the same data plotted against the ratio $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R$ for both rectangular ducts. Recall that small $\unicode[STIX]{x1D705}$ indicates that inertial lift forces are dominant and large $\unicode[STIX]{x1D705}$ indicates that secondary flow drag is dominant. The focusing horizontal focusing location approximately collapses onto a single curve, particularly for $\unicode[STIX]{x1D705}\ll 10$ , in each case. Note that for the two smaller particle sizes ( $a\in \{0.05,0.10\}$ ) an additional stable equilibrium is observed when $\unicode[STIX]{x1D705}$ is small (these are the additional ‘tails’ to the main curve). These are again expected, since $R$ is comparatively large such that the focusing behaviour is approaching that of a straight duct, and are not part of the scaling theory. Going from small to large $\unicode[STIX]{x1D705}$ we can see that particles initially focus close to the centre of the duct, then move towards the left wall before reaching a minimum and then gradually moving back towards the centre again. Observe that the minimum achieved is increasing with respect to particle size (that is, larger particles do not get as close to the inside wall as smaller ones). This effect is amplified in the case of the largest particle in the wider of the two ducts.

Figure 11. Horizontal location of the (stable) focusing positions of particles versus $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R$ for the two rectangular ducts. The particle sizes are $a=0.05$ (blue, circles), $a=0.10$ (green, triangles), $a=0.15$ (red, squares) and $a=0.20$ (cyan, diamonds). Note the vertical axis includes only the left half of the duct.

Experiments using a spiral duct with two different rectangular cross-sections were performed by Wu et al. (Reference Wu, Guan, Hou, Bhagat and Han2012). Both ducts had a width of $W=500~\unicode[STIX]{x03BC}\text{m}$ , but with different heights $H=90~\unicode[STIX]{x03BC}\text{m}$ and $H=120~\unicode[STIX]{x03BC}\text{m}$ . The bend radius is smallest near the outlet and is approximately $R=4.5~\text{mm}$ . Beads with diameter $2a=10~\unicode[STIX]{x03BC}\text{m}$ and $2a=6~\unicode[STIX]{x03BC}\text{m}$ were suspended in flow through the ducts at flow rates of $1~\text{mL}~\text{min}^{-1}$ and $2~\text{mL}~\text{min}^{-1}$ for the smaller and larger cross-sections, respectively. In both cases the larger particles focused towards the inside wall, whereas the smaller ones were less focused in a region centred slightly towards the inside wall from the centre. The corresponding $\unicode[STIX]{x1D705}$ (near the outlet where $R$ is smallest) for the larger particle is approximately $29$ and $92$ for the smaller and larger cross-section, respectively, while for the smaller particle it is approximately $135$ and $427$ . Looking at the results in figure 11(b), noting that the particle sizes $a=0.05,0.10$ are closest to those of the experiment (given $H=2$ in our computations), we see that the larger particle is indeed expected to focus nearer the inside wall in each case. While the case $\unicode[STIX]{x1D705}=427$ is outside the range of figure 11(b), the trend for increasing $\unicode[STIX]{x1D705}$ is reasonably clear (i.e. particles focus closer to the centre). Additionally, since $2R/\ell =100$ and $2a/\ell \approx 0.067$ , the behaviour can also be estimated from our results for a particle with radius $a=0.05$ in figure 10(b). Thus our model shows qualitative agreement with these experiments.

We summarise our results here for the rectangular ducts. With the exception of an additional focusing position for small $a$ and large $R$ , there is generally a unique horizontal focusing location of a particle suspended in flow through a curved rectangular duct. This is significantly different from the square duct case and explains why rectangular ducts are particularly useful for microfluidic applications requiring a focused stream of particles. Furthermore, the horizontal focusing position is sensitive to both the particle size and bend radius, providing a practical mechanism for size-based particle separation. However, our results demonstrate that the ordering of focused particles changes several times with bend radius, even in the simplified scenario of low flow rates. This suggests that much care should be taken with device design (specifically, the choice of bend radius, depending upon the size of particles to be separated). We have shown that the horizontal focusing position approximately collapses onto a single curve when plotted against $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R$ . Additionally, although there are some small differences for the two different aspect ratio ducts considered, the general focusing behaviour is qualitatively similar with respect to $\unicode[STIX]{x1D705}$ . This provides a tangible variable that can be used to guide the design of devices with rectangular cross-section operating at suitably low flow rates.

Figure 12. The force $\boldsymbol{F}_{p}^{\prime }$ on neutrally buoyant particles in a curved duct with trapezoidal cross-section ( $W=8$ , $H_{left}=3/2$ and $H_{right}=5/2$ ) and bend radius $R=160$ . The colour background shows the magnitude of $\boldsymbol{F}_{p}^{\prime }$ . Black and white contours are the zero level set curves of the horizontal and vertical components, respectively, whereas the arrows indicate the sign of each component in the area bounded by the respective contour. The left wall is on the inside of the bend. The dashed red line shows where the centre of the particle lies when its surface touches the wall.

Figure 13. Approximate trajectories of neutrally buoyant particles in a curved duct with trapezoidal cross-section ( $W=8$ , $H_{left}=3/2$ and $H_{right}=5/2$ ) and bend radius $R=160$ . Trajectories from several starting positions are superimposed. Green, orange and red markers show the location of stable, (unstable) saddle and unstable equilibria, respectively (with marker size indicative of particle size). The left wall is on the inside of the bend. The dashed red line shows where the centre of the particle lies when its surface touches the wall.

4.3 Ducts having a trapezoidal cross-section

We now consider the behaviour of particles in a curved duct with trapezoidal cross-section inspired by the experiments of Guan et al. (Reference Guan, Wu, Bhagat, Li, Chen, Chao, Ong and Han2013) and Warkiani et al. (Reference Warkiani, Guan, Luan, Lee, Bhagat, Kant Chaudhuri, Tan, Lim, Lee and Chen2014), in which the bottom wall remains perpendicular to the flow direction and parallel to the bend plane, whilst the top wall is sloped (with straight/vertical sidewalls). We consider a cross-section width $W=8$ and heights $H_{left}=1.5$ and $H_{right}=2.5$ at the left (inside) and right (outside) walls, respectively (see figure 12). This new shape has a significant effect on the background flow compared to the rectangular ducts; in particular, the location of the maximum of the axial component and the centre of the vortices of the secondary component are each shifted towards the outside wall. Furthermore, the asymmetry of the cross-section means that the axial and secondary components no longer have even and odd symmetry, respectively, with respect to $z$ . These changes have a significant effect on the focusing behaviour of particles within the duct.

The perturbing force $\boldsymbol{F}_{p}^{\prime }$ was computed for each particle size $a\in \{0.20,0.15,0.10,0.05\}$ and bend radius $R\in \{640,320,160,80\}$ . To illustrate the effect of duct curvature the straight duct case ( $R\rightarrow \infty$ , or equivalently $\unicode[STIX]{x1D716}\rightarrow 0$ ) is also computed with results provided in appendix C. Figures 12 and 13 show force and trajectory plots, respectively, in the specific case of $R=160$ for the particle sizes $a\in \{0.05,0.10,0.20\}$ . Results for the smallest particle ( $a=0.05$ and $\ell ^{4}/4a^{3}R=200$ ) are shown in figures 12(a) and 13(a). It is immediately evident from the zero level set curves that the asymmetric shape of the cross-section has a significant impact on the location of equilibria in comparison to the rectangular duct case. Three equilibria can be identified in the force plot: one very close to the centre of the outer wall, which is clearly a saddle equilibrium. The trajectory plot shows that, despite the secondary flow drag being the dominant effect, the particles gradually spiral in towards one of the two remaining equilibria, which are stable. In contrast to the results for a rectangular duct, these stable equilibria are in the outside half of the cross-section and are also slightly staggered with respect to their horizontal position.

The results for the second smallest particle ( $a=0.10$ and $\ell ^{4}/4a^{3}R=25$ ) are shown in figures 12(b) and 13(b). Three equilibria can again be identified, a saddle remains near the outside wall, whereas the stable pair has shifted to the opposite side of the cross-section. The magnitude of the secondary flow drag and inertial lift force is similar and their interaction leads to this significant shift in the focusing location. The particle dynamics shows an initial convergence onto a ‘slow manifold’ along which a slower migration towards one of the two stable equilibria takes place. Note that the horizontal location of the stable equilibria is staggered even more so than in the case of the smaller particle. The force and trajectories of the largest particle ( $a=0.20$ and $\ell ^{4}/4a^{3}R=3.125$ ) are shown in figures 12(c) and 13(c). Five equilibria can be identified in this case. Three of these lie on the white contour running along the centre (vertically) and are all unstable (the outer two are saddles). The remaining two are stable and are now located slightly right of centre. The inertial lift is dominant in this case and the trajectory plot demonstrates that particles converge more directly onto the ‘slow manifold’ before migrating towards a stable equilibrium. The results for $a=0.15$ have been omitted as they are similar to the case $a=0.20$ , but with the stable pair located on the ‘slow manifold’ slightly left of centre.

Figure 14. Horizontal location of (stable) focusing positions of particles with respect to (a) $\unicode[STIX]{x1D716}^{-1}=2R/\ell$ and (b) $\unicode[STIX]{x1D705}\ell ^{4}/4a^{3}R$ for a trapezoidal duct with $W=8$ , $H_{left}=3/2$ and $H_{right}=5/2$ . The particle sizes are $a=0.05$ (blue, circles), $a=0.10$ (green, triangles), $a=0.15$ (red, squares) and $a=0.20$ (cyan, diamonds).

We again consider how the focusing location changes with respect to $R$ (or equivalently $\unicode[STIX]{x1D716}^{-1}=2R/\ell$ given fixed $\ell$ ). Using the same approach as in the case of the rectangular cross-sections, the $\boldsymbol{F}_{p}^{\prime }$ are interpolated between the given $R$ samples and the horizontal location of the stable equilibria are determined over a finer sampling of $R$ . The horizontal focusing position versus $\unicode[STIX]{x1D716}^{-1}$ is shown in figure 14(a). The staggering of the equilibria pair due to the asymmetry of the cross-section is shown as two different, albeit close, focusing positions for each $a$ and $R$ . In figure 14(a) there are three main regions that can be identified over the $\unicode[STIX]{x1D716}^{-1}$ considered. For $\unicode[STIX]{x1D716}^{-1}\lessapprox 80$ the two smaller particle sizes are focused near $r=2$ , the largest is near $r=0$ , and the second largest shows some significant staggering but is generally on the inside half of the cross-section. For $80\lessapprox \unicode[STIX]{x1D716}^{-1}\lessapprox 400$ the second smallest particle effectively jumps and now focuses nearest to the inside wall while the relative order of the remaining three is the same. For $R\gtrapprox 400$ the smallest particle has now jumped and focuses nearest to the inside wall. These apparent discontinuities in the focusing location with respect to bend radius are a significant difference compared to the rectangular duct case, where the curve is smooth.

Figure 15. The force $\boldsymbol{F}_{p}^{\prime }$ and trajectories of a neutrally buoyant particle having radius $a=0.10$ in a curved duct having trapezoidal cross-section ( $W=8$ , $H_{left}=3/2$ and $H_{right}=5/2$ ) and bend radius $R=80$ . Interpretation of the two plots is the same as in figures 12 and 13, respectively. Many stable equilibria can be observed in this case.

Looking carefully around the discontinuities in figure 14(a) there is a small range for which more than two stable focusing positions exist. Specifically, this is evident for the smallest particle ( $a=0.05$ , blue circles) at $\unicode[STIX]{x1D716}^{-1}\approx 400$ and the second smallest particle ( $a=0.10$ , green triangles) at $\unicode[STIX]{x1D716}^{-1}\approx 80$ . The force and trajectory plots for the latter case ( $a=0.10$ and $\unicode[STIX]{x1D716}^{-1}=80$ ) are shown in figure 15. Comparing figure 15(a) with figure 12(b) the white contours along the mid-section have shifted relative to the black contours such that they each now cross three times, resulting in two stable equilibria and an unstable saddle in between. With further decrease in $\unicode[STIX]{x1D716}^{-1}$ (or equivalently $R$ ) the relative locations of the black and white contours in the mid-section shift further such that only the rightmost intersections remain, thus leaving only the two stable equilibria on the right. The trajectories in figure 15(b) suggest that the majority of starting positions migrate towards the two stable equilibria on the right. However, in the context of a spiral device, in which $R$ is continuously changing, this may not be the case and further investigation is required.

In figure 14(b) the focusing data are plotted against the ratio $\unicode[STIX]{x1D705}=\ell ^{4}/4a^{3}R$ . In a broad sense the general trend is similar to the rectangular duct case (specifically the focusing positions move towards the inside wall and back again with increasing $\unicode[STIX]{x1D705}$ ). Additionally, for $\unicode[STIX]{x1D705}<10$ and $\unicode[STIX]{x1D705}>100$ the behaviour approximately collapses onto a single curve. However, for intermediate $\unicode[STIX]{x1D705}$ some significant differences in focusing position occur, depending on the particle size. The three smaller particle sizes achieve a stable focusing position within two units of the left wall for $\unicode[STIX]{x1D705}$ in this intermediate range. Each also exhibits a rapid change in focusing position from the inside half of the duct to the outside half when $\unicode[STIX]{x1D705}$ is a little larger than where the minimum is achieved. Furthermore, this rapid change occurs at larger $\unicode[STIX]{x1D705}$ as the particle size decreases. On the other hand, the largest particle ( $a=0.20$ ) does not get close to the inside wall for any $\unicode[STIX]{x1D705}$ (and is only briefly left of centre). Furthermore, its focusing behaviour appears to be reasonably smooth with respect to $\unicode[STIX]{x1D705}$ compared to that observed for the smaller particles.

It is interesting to consider how these different features of focusing behaviour can affect the efficiency of such devices for the size-based separation of particles. The staggering of the horizontal location of the stable equilibria in the top half and bottom half of the duct is somewhat undesirable. To some degree this is offset by the larger separation distances provided by the trapezoidal duct, resulting from stable equilibria being able to exist in the outer half of the cross-section. An obvious modification of the duct would be to make the bottom wall slanted (in the opposite direction) to achieve a symmetric trapezoidal shape. This should eliminate the horizontal staggering of equilibria pairs while maintaining a large separation distance, although it may be difficult to reliably produce such cross-sections. We also observe that the focusing location of smaller particles can change rapidly with respect to small changes in $R$ . While this results in a narrow band of design parameters in which the particle may be found at multiple locations, it could also be exploited in carefully designed microfluidic devices to achieve reasonably high separation efficiency of particles whose size may differ by a comparatively small amount.

Experiments in two different spiral ducts with trapezoidal cross-section using beads with diameter $2a=10~\unicode[STIX]{x03BC}\text{m}$ and $2a=6~\unicode[STIX]{x03BC}\text{m}$ were performed by Wu et al. (Reference Wu, Guan, Hou, Bhagat and Han2012). Both ducts had a width of $W=500~\unicode[STIX]{x03BC}\text{m}$ but with the different heights of $H_{left}=70~\unicode[STIX]{x03BC}\text{m}$ to $H_{right}=100~\unicode[STIX]{x03BC}\text{m}$ from the inner to outer wall in the smaller cross-section and $H_{left}=90~\unicode[STIX]{x03BC}\text{m}$ to $H_{right}=120~\unicode[STIX]{x03BC}\text{m}$ in the larger cross-section. In both cases, the larger particles focused towards the inside wall while the smaller ones were less focused in a region slightly towards the outside wall from the centre. The corresponding $\unicode[STIX]{x1D705}$ (near the outlet where $R$ is smallest) for the larger particle is approximately $23$ and $54$ for the smaller and larger cross-section, respectively, while for the smaller particle it is approximately $107$ and $250$ . Looking at the results in figure 14(b), noting that the particle sizes $a=0.05,0.10$ are closest to those of the experiment, we see that the larger particle is indeed expected to focus near the inside wall for these $\unicode[STIX]{x1D705}$ , whereas the smaller particle will focus towards the outer wall. Interestingly the staggering of the equilibria pair is not evident in the experimental results, although this is potentially because the staggering is less than the width of the fluorescent streaks in the experimental data.