Effect of inertial lift on a spherical particle suspended in flow through a curved duct

TL;DR

This paper models the inertial lift forces on spherical particles in curved duct flow at low Reynolds numbers, predicting stable equilibrium positions and how they vary with particle size, bend radius, and duct shape.

Contribution

It extends existing inertial lift models to curved ducts, incorporating effects of duct shape and curvature on particle migration.

Findings

Stable equilibrium positions depend on particle size and duct curvature.

Focusing position scales with a dimensionless parameter involving particle and duct dimensions.

Cross-sectional shape influences particle migration dynamics.

Abstract

We develop a model of the forces on a spherical particle suspended in flow through a curved duct under the assumption that the particle Reynolds number is small. This extends an asymptotic model of inertial lift force previously developed to study inertial migration in straight ducts. Of particular interest is the existence and location of stable equilibria within the cross-sectional plane towards which particles migrates. The Navier-Stokes equations determine the hydrodynamic forces acting on a particle. A leading order model of the forces within the cross-sectional plane is obtained through the use of a rotating coordinate system and a perturbation expansion in the particle Reynolds number of the disturbance flow. We predict the behaviour of neutrally buoyant particles at low flow rates and examine the variation in focusing position with respect to particle size and bend radius, independent of the flow rate. In this regime, the lateral focusing position of particles approximately collapses with respect to a dimensionless parameter dependent on three length scales, specifically the particle radius, duct height, and duct bend radius. Additionally, a trapezoidal shaped cross-section is considered in order to demonstrate how changes in the cross-section design influence the dynamics of particles.

Tables

Table 1. Table 1 : Summary of nomenclature used throughout the paper.

symbol	description
$a$	radius of the (spherical) particle
$𝒟,\partial 𝒟$	duct interior and its boundary (in lab frame)
$𝑭$	force on a particle (in lab frame)
${𝑭}_p^{\prime }$	cross-sectional forces on particle from our model
$ℱ,\partial ℱ$	fluid domain and its boundary when a particle is present (in lab frame)
$G$	pressure gradient along duct centre-line which drives the fluid flow
$H$	height of the duct cross-section
$\ell$	characteristic dimension of duct cross-section (taken to be $\min \{H,W\}$)
$I_p$	moment of inertia of the (spherical) particle
$𝕀$	identity tensor
$K$	perturbation parameter of the background fluid flow (taken to be $ϵ{\text{Rey}}^2/4$)
$m_p$	mass of the (spherical) particle
$𝒏$	normal vector on a prescribed surface pointing outwards relative to the fluid
$p,\overline{p}$	fluid pressure with and without particle present respectively (in lab frame)
$q$	disturbance fluid pressure
$r$	radial coordinate with respect to the $x$–$y$ plane (in the lab frame)
$r_p$	radial coordinate of particle position (frame independent)
$R$	bend radius of the duct (measured to centre of cross-section)
$\text{Rey}$	Reynolds number (taken to be $\left(\rho /\mu \right)U_m\ell$)
${\text{Rey}}_p$	particle Reynolds number (taken to be ${\left(a/\ell \right)}^2\text{Rey}$)
$t$	time variable
$𝑻$	torque on a particle (in lab frame)
$𝒖,\overline{𝒖}$	fluid velocity with and without particle present respectively (in lab frame)
${\overline{𝒖}}_a,{\overline{𝒖}}_s$	axial and secondary components of the background fluid velocity
${𝒖}_p$	particle velocity (in lab frame)
$U_m$	characteristic velocity of background flow (taken as the max of ${\overline{𝒖}}_a$)
$U_p$	axial component of the particle velocity
$𝒗$	disturbance fluid velocity
$W$	width of duct cross-section
$𝒙$	spatial coordinates (in lab frame)
${𝒙}_p$	particle position (in lab frame)
$x_p,y_p,z_p$	Cartesian coordinates of particle position (in lab frame)
$ϵ$	half the characteristic duct length divided by the bend radius, that is $\ell /2R$
$\theta$	angular coordinate with respect to the $x$–$y$ plane (in lab frame)
${\theta }_p$	angular coordinate of particle position (in lab frame)
$\Theta$	angular velocity of the rotating reference frame ($z$ component)
$\kappa$	dimensionless scale of the secondary flow drag (relative to the inertial lift force)
$\mu$	fluid viscosity
$\rho$	fluid density
${\rho }_p$	particle density (taken to be uniform)
${𝛀}_p$	particle spin (in lab frame)

Equations

$$\mathbf{x}(r,\theta ,z)=(R+r)\cos (\theta ){\mathbf{e}}_x+(R+r)\sin (\theta ){\mathbf{e}}_y+z{\mathbf{e}}_z.$$

$$-\nabla \bar{p}+\mu {\nabla }^2\bar{\mathbf{u}}$$

$$\nabla \cdot \bar{\mathbf{u}}$$

$$\bar{\mathbf{u}}$$

$$\bar{p}(\mathbf{x}(r,\theta ,z))=-GR\theta +\rho \mathbf{x}\cdot \mathbf{g}+\tilde{p}(r,z),$$

$$\bar{\mathbf{u}}={\bar{u}}_{\theta }(r,z){\hat{\mathbf{e}}}_{\theta }+{\bar{u}}_r(r,z){\hat{\mathbf{e}}}_r+{\bar{u}}_z(r,z){\mathbf{e}}_z,$$

$$-\nabla p+\mu {\nabla }^2\mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$\mathbf{u}$$

$$\mathbf{F}(t)$$

$$\mathbf{T}(t)$$

$${\mathbf{x}}^{\prime }(r^{\prime },{\theta }^{\prime },z^{\prime })=(R+r^{\prime })\cos ({\theta }^{\prime }){\mathbf{e}}_{x^{\prime }}+(R+r^{\prime })\sin ({\theta }^{\prime }){\mathbf{e}}_{y^{\prime }}+z^{\prime }{\mathbf{e}}_{z^{\prime }}.$$

$${\mathbf{e}}_{x^{\prime }}=\cos ({\theta }_p){\mathbf{e}}_x+\sin ({\theta }_p){\mathbf{e}}_y,{\mathbf{e}}_{y^{\prime }}=-\sin ({\theta }_p){\mathbf{e}}_x+\cos ({\theta }_p){\mathbf{e}}_y,{\mathbf{e}}_{z^{\prime }}={\mathbf{e}}_z.$$

$${\mathbf{e}}_x=\cos ({\theta }_p){\mathbf{e}}_{x^{\prime }}-\sin ({\theta }_p){\mathbf{e}}_{y^{\prime }},{\mathbf{e}}_y=\sin ({\theta }_p){\mathbf{e}}_{x^{\prime }}+\cos ({\theta }_p){\mathbf{e}}_{y^{\prime }},{\mathbf{e}}_z={\mathbf{e}}_{z^{\prime }}.$$

$$\Theta {\mathbf{e}}_z:=\frac{\partial {\theta }_p}{\partial t}{\mathbf{e}}_z=\frac{U_p}{R+r_p}{\mathbf{e}}_z.$$

$$\frac{\partial {\mathbf{x}}_p^{\prime }}{\partial t}=\frac{\partial x_p^{\prime }}{\partial t}{\mathbf{e}}_{x^{\prime }}+\frac{\partial z_p^{\prime }}{\partial t}{\mathbf{e}}_{z^{\prime }}=\frac{\partial r_p}{\partial t}{\mathbf{e}}_{x^{\prime }}+\frac{\partial z_p}{\partial t}{\mathbf{e}}_{z^{\prime }}=0.$$

$${\mathbf{\Omega }}_p^{\prime }:={\mathbf{\Omega }}_p-\Theta {\mathbf{e}}_{z^{\prime }}.$$

$$\mathbf{u}(\mathbf{x},t)={\mathbf{u}}^{\prime }({\mathbf{x}}^{\prime },t)+\Theta ({\mathbf{e}}_{z^{\prime }}\times {\mathbf{x}}^{\prime }).$$

$$\frac{d\mathbf{u}}{dt}=\frac{D{\mathbf{u}}^{\prime }}{Dt}+2\Theta ({\mathbf{e}}_{z^{\prime }}\times {\mathbf{u}}^{\prime })+{\Theta }^2({\mathbf{e}}_{z^{\prime }}\times ({\mathbf{e}}_{z^{\prime }}\times {\mathbf{x}}^{\prime })),$$

$$-{\nabla }^{\prime }{p}^{\prime }+\mu {\nabla }^{\prime 2}{\mathbf{u}}^{\prime }$$

$$\displaystyle\qquad\qquad+\Theta^{2}(\boldsymbol{e}_{z^{\prime}}\times(\boldsymbol{e}_{z^{\prime}}\times\boldsymbol{x}^{\prime}))+g\boldsymbol{e}_{z}\bigg{)}$$

$${\nabla }^{\prime }\cdot {\mathbf{u}}^{\prime }$$

$${\mathbf{u}}^{\prime }$$

$$\mathbf{F}$$

$$\mathbf{T}$$

$${\mathbf{F}}^{\prime }$$

$$+{\int }_{|{\mathbf{x}}^{\prime }-{\mathbf{x}}_p^{\prime }|=a}(-\mathbf{n})\cdot (-p^{\prime }\mathbb{I}+\mu ({\nabla }^{\prime }{\mathbf{u}}^{\prime }+{{\nabla }^{\prime }{\mathbf{u}}^{\prime }}^{⊺}))d{S}^{\prime },$$

$${\mathbf{T}}^{\prime }$$

$$+{\int }_{|{\mathbf{x}}^{\prime }-{\mathbf{x}}_p^{\prime }|=a}({\mathbf{x}}^{\prime }-{\mathbf{x}}_p^{\prime })\times \left((-\mathbf{n})\cdot (-p^{\prime }\mathbb{I}+\mu ({\nabla }^{\prime }{\mathbf{u}}^{\prime }+{{\nabla }^{\prime }{\mathbf{u}}^{\prime }}^{⊺}))\right)d{S}^{\prime }.$$

$${\bar{\mathbf{u}}}^{\prime }({\mathbf{x}}^{\prime })=\bar{\mathbf{u}}({\mathbf{x}}^{\prime })-\Theta {\mathbf{e}}_{z^{\prime }}\times {\mathbf{x}}^{\prime },{\bar{p}}^{\prime }({\mathbf{x}}^{\prime })=\bar{p}(\mathbf{x}).$$

$${\mathbf{v}}^{\prime }={\mathbf{u}}^{\prime }-{\bar{\mathbf{u}}}^{\prime },q^{\prime }=p^{\prime }-{\bar{p}}^{\prime },$$

$$-{\nabla }^{\prime }{q}^{\prime }+\mu {\nabla }^{\prime 2}{\mathbf{v}}^{\prime }$$

$$\displaystyle\qquad\qquad+\boldsymbol{v}^{\prime}\cdot\nabla^{\prime}\bar{\boldsymbol{u}}+\Theta(\boldsymbol{e}_{z^{\prime}}\times\boldsymbol{v}^{\prime})\big{)}$$

$${\nabla }^{\prime }\cdot {\mathbf{v}}^{\prime }$$

$${\mathbf{v}}^{\prime }$$

$${\mathbf{F}}^{\prime }$$

$$+{\int }_{|{\mathbf{x}}^{\prime }-{\mathbf{x}}_p^{\prime }|=a}(-\mathbf{n})\cdot (-q^{\prime }\mathbb{I}+\mu ({\nabla }^{\prime }{\mathbf{v}}^{\prime }+{{\nabla }^{\prime }{\mathbf{v}}^{\prime }}^{⊺}))d{S}^{\prime },$$

Full text

1 Introduction

Inertial lift force influences the motion of particles suspended in fluid flow through a duct. The effect was first demonstrated in the classical experiment of Segre & Silberberg (1961) where particles suspended in flow through a (straight) cylindrical pipe were observed to migrate towards an annulus approximately $0.6$ times the radius of the pipe. This sparked many studies of the simplified scenario of a spherical particle suspended in Poiseuille flow between two plane parallel walls (e.g. Ho & Leal (1974), Schonberg & Hinch (1989) and Asmolov (1999)). Subsequently, particle migration in a circular pipe was investigated at much larger Reynolds numbers ($\Rey$) with studies showing that the radius of the focusing annulus grows with increasing $\Rey$ (Matas et al., 2004, 2009). In recent years inertial lift has been used in the field of microfluidics for the separation of cells by size and has been extensively studied experimentally (Di Carlo, 2009; Martel & Toner, 2014; Warkiani et al., 2014; Geislinger & Franke, 2014). Much of the behaviour, especially in complex geometries, is still only understood at an empirical level.

In the case of fully enclosed non-circular ducts the inertial lift force is generally much more difficult to approach from an analytical perspective. Several experimental studies of straight ducts with square or rectangular cross-sections have shown that particles typically tend to migrate towards one of four (stable) equilibrium positions located a finite distance from the centre of each wall (Abbas et al., 2014; Amini et al., 2014; Di Carlo et al., 2009; Ciftlik et al., 2013). Recent work by Hood et al. (2015) extended the analysis of Schonberg & Hinch (1989) from the case of Poiseuille flow between two walls to Poiseuille flow through square ducts. They develop a model of the inertial lift force at small particle Reynolds number which is calculated via the Lorentz reciprocal theorem and utilises the finite element method (FEM) to approximate terms in the expansion that are analytically intractable. Their results agree well with experimental data and demonstrate that much of the behaviour can be characterised by effects that occur for small particle Reynolds number. Their approach was also applied to rectangular ducts where it is shown that the relative size of the attraction zone for equilibria near the side walls becomes increasingly small as the aspect ratio increases, and even disappears entirely for large enough particles (Hood, 2016). Other studies have considered the behaviour at higher Reynolds numbers, in particular Nakagawa et al. (2015) used an immersed boundary method to approximate the full Navier–Stokes equations and demonstrate that particles may additionally focus near the corners at high flow rates in agreement with experiments of Miura et al. (2014). These studies also clarify that particles are less focused at high flow rates suggesting that lower flow rates are more useful for cell separation and sorting.

This paper extends the work of Hood et al. (2015) to the case of curved ducts. The focusing behaviour of curved and spiral ducts has been largely confined to experimental studies which have shown it to be an effective mechanism to separate particles/cells by size (Martel & Toner, 2012, 2013; Ramachandraiah et al., 2014; Russom et al., 2009; Nivedita et al., 2017). There have been some studies utilising direct numerical simulation, see for example Liu et al. (2016) and Ookawara et al. (2006), in which the inertial lift effect is often added as an external force on the particle via an over-simplified model. Curved ducts are well-known to develop a secondary flow consisting of two counter-rotating vortices within the duct cross-section commonly referred to as Dean flow in reference to the early studies of Dean (1927) and Dean & Hurst (1959). The experimental literature demonstrates that size based separation is caused by an interplay between the inertial lift force and the drag forces coming from the secondary fluid flow. Here we carefully analyse the Navier–Stokes equations to develop a model of the forces, acting on a particle, within the cross-sectional plane. The model is then applied to the special case of low flow rate to gain some insight into how the inertial lift force and secondary flow drag interact and lead to focusing positions that differ with respect to particle size. Motivated by the experiments of Guan et al. (2013) and Warkiani et al. (2014) we also consider a curved duct having a trapezoidal cross-section to better understand how the shape of a duct influences the separation of different size particles.

Microfluidic devices used for particle separation often have a spiral design in order to accommodate the required focusing length relative to the bend radius. Although the bend radius in such devices is continuously changing, it has been demonstrated that the slowly changing curvature in such devices has negligible impact on the fluid flow in the sense that it is reasonable to treat the bend radius as being locally constant at each point within the spiral (Harding & Stokes, 2018). Therefore, ducts with constant bend radius are considered in this study and, by smoothly interpolating results obtained for several different bend radii (using the same cross-section), we are able to approximate the particle dynamics at any location within a spiral duct.

A leading order model of the relevant cross-sectional forces is developed in Section 2. There are several key steps in the case of a curved duct. The first is the introduction of a rotating reference frame in which the fluid and particle motion is approximately steady. This allows us to neglect time dependence of the cross-sectional forces on the particle, but also introduces inertia and Coriolis forces into the Navier–Stokes equations and force model. The second step is the introduction of disturbance flow variables followed by a non-dimensionalisation of the problem and a perturbation expansion with respect to the particle Reynolds number. This is similar in spirit to the approach of (Hood et al., 2015) for a straight duct. The Lorentz reciprocal theorem is used to re-express first order forces, which include the inertial lift force, as a volume integral depending on the leading order approximation of the disturbance flow. Third, we break up the background flow velocity into its axial and secondary components and use this to further expand the force and torque experienced by the particle into distinct parts. Lastly, a model for the trajectory of a particle is developed based on its terminal velocity resulting from the forces driving its motion within the plane of the duct cross-section.

In Section 3 we apply the general model from the preceding section to the specific case of particle behaviour in the limit of low flow rate. This is an interesting case to consider for two reasons. Firstly it allows us to consider a simplified approximation of the background fluid flow and eliminate dependence on flow velocity so that we may focus on the effect of bend radius and particle size on the focusing behaviour. In particular, by utilising an expansion of the background flow developed in (Harding, 2018a) we demonstrate that it is sufficient to take a leading order approximation of the axial and secondary components of the flow. It then becomes clear that the interaction between the secondary flow drag and inertial lift force converges towards a limiting behaviour (for a fixed particle size and duct dimensions) that is different from that obtained for straight ducts. Second, it allows us to validate the model by demonstrating that the focusing behaviour observed in straight ducts is recovered for large bend radii. This provides new insights into how the focusing behaviour in curved ducts differs from that in straight ducts.

In Section 4 we present and discuss the focusing behaviour resulting from our model for ducts having square, rectangular and trapezoidal cross-sections. Curved square ducts are not be particularly good at focusing or separating particles but illustrate some interesting dynamics which we comment on briefly. On the other hand, rectangular and trapezoidal ducts are quite effective at focusing particles and the focusing position depends on the size of the particle and bend radius of the duct. Furthermore, the behaviour approximately collapses to a single curve when plotted with respect to a dimensionless parameter which describes the relative size of the inertial lift and secondary drag forces. These results provide insights into recently published experimental results.

2 Governing equations and perturbation expansion

Here we derive the equations for our leading order force model.

2.1 Governing equations in the lab frame of reference

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,\theta$ such that $r,z$ are coordinates that describe the location within the duct cross-section relative to its centre and $\theta$ 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

[TABLE]

For example, the point on the centre-line $\boldsymbol{x}=R\boldsymbol{e}_{x}$ is equivalently described by $(r,\theta,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 average height 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\leq 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

[TABLE]

where $\rho$ and $\mu$ denote the fluid density and viscosity respectively (both assumed to be constant), $\boldsymbol{g}$ is the (constant) acceleration due to gravity, $\mathcal{D}$ denotes the interior of the duct and $\partial\mathcal{D}$ denotes its boundary/surface (i.e. the walls of the duct). The fluid motion is driven by a steady pressure gradient such that $\bar{p}$ can be decomposed into

[TABLE]

where $G$ denotes the drop in pressure per unit arc-length along the centreline of the duct. Notice that we have assumed gravity does not influence the flow velocity by introducing the hydrostatic pressure $\rho\boldsymbol{x}\cdot\boldsymbol{g}$ to counteract $\rho\boldsymbol{g}$ in (1a). 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, 1959; Winters, 1987; Yanase et al., 1989; Yamamoto et al., 2004). 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

[TABLE]

where $\hat{\boldsymbol{e}}_{\theta}:=-\sin(\theta)\boldsymbol{e}_{x}+\cos(\theta)\boldsymbol{e}_{y}$ and $\hat{\boldsymbol{e}}_{r}:=\cos(\theta)\boldsymbol{e}_{x}+\sin(\theta)\boldsymbol{e}_{y}$. The $\bar{u}_{\theta}$ 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 Section 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. There are two ways one could analyse the trajectory of a particle within the flow. If the secondary component of the background flow is relatively strong (i.e. such that inertial lift force will be small compared to fluid drag from the secondary flow) then one might view the particle as behaving like a passive tracer (with respect to the background flow velocity $\bar{\boldsymbol{u}}$) which is then perturbed by additional drag (due to the nearby walls) and inertial lift forces. Alternatively, if the secondary component of the background flow is moderate or weak (i.e. such that inertial lift force is similar in magnitude to drag from the secondary fluid motion) one could consider the particle to be travelling along a streamline of the axial component of the background flow (i.e. streamlines of $\bar{u}_{\theta}\hat{\boldsymbol{e}}_{\theta}$) and that its motion is perturbed by the drag forces from the secondary fluid motion, increased drag from the walls, and inertial lift forces. We consider applications in which the drag forces from the secondary fluid motion is not significantly larger than the inertial lift force and thus take the latter point of view throughout this paper.

Let $a$ denote the radius of the spherical particle, $\rho_{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),\theta_{p}(t),z_{p}(t))$, see Figure 1. The particle follows an axial flow streamline (from which it will be perturbed) and thus $\Theta:=\partial\theta_{p}/\partial t$ is constant, $\partial r_{p}/\partial t=\partial z_{p}/\partial t=0$, and $\partial\boldsymbol{\Omega}_{p}/\partial t=\boldsymbol{0}$ where $\boldsymbol{\Omega}_{p}$ denotes the spin of the particle (i.e. the angular velocity with respect to its centre). The assumption $\partial r_{p}/\partial t=\partial z_{p}/\partial 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 (1989) and Hood et al. (2015). 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 $\theta_{p}=0$ at time $t=0$ such that $\theta_{p}(t)=U_{p}t/(R+r_{p})$ (or equivalently $\Theta=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

[TABLE]

where $\mathcal{F}:=\{\boldsymbol{x}\in\mathcal{D}:|\boldsymbol{x}-\boldsymbol{x}_{p}|\geq a\}$ is the fluid domain given the presence of the spherical particle and $\partial\mathcal{F}\backslash\partial\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 $\partial\boldsymbol{u}/\partial t$ in (3a) cannot be discarded (as it was for the background flow in (1a)). Note that the boundary condition (3c) has a constant $\boldsymbol{\Omega}_{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 stresses over the particle surface:

[TABLE]

where $m_{p}:=(4/3)\pi a^{3}\rho_{p}$ is the mass of the particle. Note that we take the normal $\boldsymbol{n}$ with respect to the fluid domain $\mathcal{F}$ and thus on the surface of the spherical particle $-\boldsymbol{n}$ is the normal vector pointing outwards from the centre of the particle.

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 $\partial\theta_{p}/\partial t=\Theta$ 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 $\partial\boldsymbol{u}/\partial t$ term.

Consider a reference frame rotating about the $z$ axis at the (constant) rate $\partial\theta_{p}/\partial t=:\Theta$. We introduce the coordinate system $\boldsymbol{x}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})$ which we parametrise by $(r^{\prime},\theta^{\prime},z^{\prime})$ such that

[TABLE]

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

[TABLE]

or conversely

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where $D\boldsymbol{u}^{\prime}/Dt:=\partial\boldsymbol{u}^{\prime}/\partial t+\boldsymbol{u}^{\prime}\cdot\nabla^{\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 $\nabla^{\prime}p^{\prime}=\nabla p$), it follows that the equations of motion for the fluid in the rotating reference frame are

[TABLE]

where $\mathcal{F}^{\prime}:=\{\boldsymbol{x}^{\prime}\in\mathcal{D}:|\boldsymbol{x}^{\prime}-\boldsymbol{x}_{p}^{\prime}|\geq a\}$ denotes the fluid domain with respect to the rotating reference frame. The boundary condition (6c) arises from the no slip condition on the walls, rotating with angular velocity $-\Theta\boldsymbol{e}_{z}$. Notice that, expressed in this frame, the fluid flow can be considered to be steady since $\mathcal{F}^{\prime}$ is stationary and each of $\Theta,\boldsymbol{\Omega}_{p}^{\prime},g$ are constant. Consequently $\partial\boldsymbol{u}^{\prime}/\partial t$ is zero and can be dropped from (6a). Observe that it is important for gravity to acts in the $\boldsymbol{e}_{z}$ direction since otherwise the relative direction of gravity in the rotating frame would be dependent on $\theta_{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

[TABLE]

where $m_{p}:=(4/3)\pi a^{3}\rho_{p}$ and $I_{p}:=(2/5)a^{2}m_{p}$ are the mass and moment of inertia of the particle respectively. Note there is no Coriolis like force since the particle velocity is taken to be zero in the rotating reference frame, nor is there an Euler force since $\Theta$ is assumed to be constant. Thus, upon re-expressing (4) in terms of $p^{\prime},\boldsymbol{u}^{\prime}$, and then combining with (7) one obtains

[TABLE]

Note that $-m_{p}\Theta^{2}(\boldsymbol{e}_{z^{\prime}}\times(\boldsymbol{e}_{z^{\prime}}\times\boldsymbol{x}_{p}^{\prime}))$ provides the centripetal force experienced by the particle due to its rotational motion as it travels through the curved duct.

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

[TABLE]

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

[TABLE]

respectively. Substituting $\boldsymbol{u}^{\prime}=\boldsymbol{v}^{\prime}+\bar{\boldsymbol{u}}^{\prime}=\boldsymbol{v}^{\prime}+\bar{\boldsymbol{u}}-\Theta(\boldsymbol{e}_{z^{\prime}}\times\boldsymbol{x}^{\prime})$ and $p^{\prime}=q^{\prime}+\bar{p}^{\prime}$ into (6) one obtains

[TABLE]

where we have additionally used the fact that

• (i)

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

• (ii)

for any $\boldsymbol{w}$ one has $\boldsymbol{w}\cdot\nabla^{\prime}(\boldsymbol{e}_{z^{\prime}}\times\boldsymbol{x}^{\prime})=\boldsymbol{e}_{z^{\prime}}\times\boldsymbol{w}$,

• (iii)

using (2) it can be shown $(\boldsymbol{e}_{z^{\prime}}\times\boldsymbol{x}^{\prime})\cdot\nabla^{\prime}\bar{\boldsymbol{u}}=\boldsymbol{e}_{z^{\prime}}\times\bar{\boldsymbol{u}}$,

• (iv)

for $\boldsymbol{x}^{\prime}\in\partial\mathcal{D}$ one has $\bar{\boldsymbol{u}}=\boldsymbol{0}$.

Similarly the force and torque can be expressed in terms of the disturbance flow as

[TABLE]

respectively (noting that $\nabla^{\prime}(\boldsymbol{e}_{z^{\prime}}\times\boldsymbol{x}^{\prime})+{\nabla^{\prime}(\boldsymbol{e}_{z^{\prime}}\times\boldsymbol{x}^{\prime})}^{\intercal}=0$). Because $\bar{\boldsymbol{u}}^{\prime},\bar{p}^{\prime}$ are well-defined on the interior of the particle then one may use the divergence theorem and (1a) to obtain

[TABLE]

recalling that $\boldsymbol{n}$ points in to the centre of the particle. Thus the contribution from the background flow to the force is essentially just the inertia of the fluid that would otherwise take up the volume occupied by the particle. For a neutrally buoyant particle (i.e. $\rho_{p}=\rho$) travelling with a velocity close to that of the axial velocity of the background flow we would then expect $m_{p}\Theta^{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 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. Similarly, for the background contribution to the torque, it may be shown that

[TABLE]

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}\cdot(\ast)$ after which the divergence theorem can be applied and then further re-arrangement and use of (1a) leads to (13). 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=\theta=z=0$. This is essentially the same scale as in Hood et al. (2015). The non-dimensionalisation of the remaining variables follow from these via the standard approach for flows in which viscous stresses are dominant. Thus we introduce dimensionless variables denoted by $\hat{\cdot}$,

[TABLE]

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

[TABLE]

where $\hat{\mathcal{F}}^{\prime}:=\{\hat{\boldsymbol{x}}^{\prime}:a\hat{\boldsymbol{x}}^{\prime}\in\mathcal{F}^{\prime}\}$ and $\hat{\mathcal{D}}^{\prime}:=\{\hat{\boldsymbol{x}}^{\prime}:a\hat{\boldsymbol{x}}^{\prime}\in\mathcal{D}^{\prime}\}$ denote the rescaled domains and

[TABLE]

is the particle Reynolds number. Writing the duct Reynolds number, $\Rey:=(\rho/\mu)U_{m}\ell$, we have $\Rey_{p}=\Rey(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

[TABLE]

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

[TABLE]

Noting that $\rho U_{m}^{2}a^{4}/\ell^{2}=\Rey_{p}\mu U_{m}a^{2}/\ell$, one has from (11a)

[TABLE]

and similarly from (11b)

[TABLE]

At this stage, given $a,\hat{\boldsymbol{x}}_{p}^{\prime},\hat{\Theta},\hat{\boldsymbol{\Omega}}_{p}^{\prime}$ and an approximation of $\hat{\bar{\boldsymbol{u}}}$ (for a desired duct shape/size and flow rate), one could solve the non-linear problem (14) and then determine the resulting force and torque on the particle via (16) and (17) 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{\Theta},\hat{\boldsymbol{\Omega}}_{p}^{\prime}$ such that $\hat{\boldsymbol{T}}^{\prime}$ and the axial component of $\hat{\boldsymbol{F}}^{\prime}$ (i.e. $\hat{\boldsymbol{F}}^{\prime}\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, i.e. $\hat{\boldsymbol{x}}_{p}^{\prime}$).

In the case of a straight duct, the the inertial lift force can be approximated from a perturbation expansion of the disturbance flow for small $\Rey_{p}$ (Hood et al., 2015). Through this expansion, the non-linear equation (14) 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. (2015) 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:

[TABLE]

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

[TABLE]

and the first order equations

[TABLE]

Equations (19) and (20) are similar to those obtained for the straight duct (see Hood et al. (2015)) but with a few key differences. Firstly, the right side of (20a) has additional terms introduced in the change to a rotating coordinate system. Secondly, the background velocity $\hat{\bar{\boldsymbol{u}}}$ is no longer a simple Poiseuille flow through the duct. Lastly, the boundary condition (14d) is captured entirely in (19d) so that $\boldsymbol{v}_{1}=\boldsymbol{0}$ in (20d).

The perturbation expansion (18) can also be substituted into (16) and (17) to obtain the expansion

[TABLE]

where

[TABLE]

and similarly

[TABLE]

In the case of a straight duct Hood et al. (2015) 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

[TABLE]

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

[TABLE]

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 Section 2.7), and also provides opportunities to examine how different components of the flow influence $\boldsymbol{F}_{0}$ in more detail (as discussed in Section 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 (19d), $\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 (19) into two distinct parts. This will be done in the following section.

2.6 Separation of axial and secondary flow effects

Let

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

respectively. The force and torque is similarly separated as $\boldsymbol{F}_{-1}=\boldsymbol{F}_{-1,a}+\boldsymbol{F}_{-1,s}$ and $\boldsymbol{T}_{-1}=\boldsymbol{T}_{-1,a}+\boldsymbol{T}_{-1,s}$. Specifically, for each $\ast=a,s$ one has

[TABLE]

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 determine $\hat{\Theta},\hat{\boldsymbol{\Omega}}_{p}^{\prime}$ such that $\boldsymbol{F}_{-1,a}=\boldsymbol{T}_{-1,a}=\boldsymbol{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(\Rey_{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}}}\cdot\hat{\nabla}^{\prime}\hat{\bar{\boldsymbol{u}}}$ in (22b) can be expanded as

[TABLE]

where $y^{\prime},z^{\prime}$ components of the first integral are exactly zero (and the $x^{\prime}$ component is comparable to $(4\pi/3)\hat{\Theta}^{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 (24), the $\ast=x^{\prime},z^{\prime}$ components of the surface integral over the fluid stress from $q_{1},\boldsymbol{v}_{1}$ in (22b) can each be expanded as

[TABLE]

Here the first and second integrals provide a force in the $x^{\prime},z^{\prime}$ directions and the third integral provides a small force in the $y^{\prime}$ directions. Note that, in the context of understanding the particle behaviour within the cross-section, the $y^{\prime}$ components of $\boldsymbol{F}_{0}$ can essentially be ignored since the small perturbation to the axial motion has no further influence on the cross-section forcing to this order of magnitude. As such, one can discard the second and third component of (31) and (32) respectively.

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

[TABLE]

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{\Theta},\hat{\boldsymbol{\Omega}}_{p}^{\prime}$ are chosen such that such that $\boldsymbol{F}_{-1,a}=\boldsymbol{T}_{-1,a}=\boldsymbol{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 $\rho U_{m}^{2}a^{4}/\ell^{2}$) can be expressed as $\Rey_{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\pi$. 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 Section 2.5 which each solve (25) 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

[TABLE]

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

[TABLE]

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

[TABLE]

The solution of these coupled equations approximates the particle trajectory within the cross-section. Interpreting $\hat{\Theta}$ as a function of $\hat{x}_{p,r^{\prime}}^{\prime},\hat{x}_{p,z^{\prime}}^{\prime}$ (recalling that $\hat{\Theta},\hat{\boldsymbol{\Omega}}_{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{0}$), the angular position of the particle can also be estimated via $d\theta_{p}/dt^{\prime}=\hat{\Theta}(\hat{x}_{p,r^{\prime}}^{\prime},\hat{x}_{p,z^{\prime}}^{\prime})$.

3 Estimation of particle behaviour at low flow rates

Experiments involving curved and spiral ducts have been done with $U_{m},\Rey$ as high as $O(1),O(100)$ respectively (in order to achieve reasonable throughput/flow-rate). That said, the low flow rate case it is particularly interesting for several reasons. First, the behaviour at low flow rate in a curved duct differs from that in a straight duct, i.e. the behaviour as $U_{m}\rightarrow 0$ (or equivalently $\Rey\rightarrow 0$) differs from that as $R\rightarrow\infty$. As such, understanding what happens at low flow rates is a good first step towards understanding the difference between straight and curved ducts. Second, it allows us to reduce parameter space since the focusing behaviour becomes independent of flow rate; we are able to thoroughly examine how focusing behaviour is influenced by the bend radius and particle size in this regime. Third, we can additionally consider the limit $R\rightarrow\infty$ to check that the focusing positions in a straight duct are recovered in order to validate the model. Lastly, the assumption is not as restrictive as it may first seem and we expect the predictions may be applicable up to at least $\Rey=O(10)$.

The scaling employed to study the inertial lift force differs from the characteristic scales in the background flow; we return to the physical background flow velocity $\bar{\boldsymbol{u}}$ to examine the scale of its components. Recalling that $U_{m}$ denotes the maximum of $\bar{u}_{\theta}(r,z)$, and letting $\epsilon:=\ell/(2R)$, then it can be shown that if $K:=\epsilon\Rey^{2}/4$ is sufficiently small then

[TABLE]

This scaling is discussed Harding (2018a) in which the governing equations (1) are explicitly non-dimensionalised (a brief account is also provided in Appendix B for completeness). Furthermore, the dimensionless velocity components can be expressed in terms of a perturbation expansion in $K$, specifically

[TABLE]

This series leads to a linearisation of the governing equations in which each $\bar{u}_{\theta,i}$ and pair $\bar{u}_{r,i},\bar{u}_{z,i}$ (which are combined into a stream function $\Phi_{i}$) can be solved in an alternating fashion.

It has been demonstrated that the terms in the expansion (37) decay in a manner such that the expansion generally converges whenever $K\lessapprox 200$ (Harding, 2018a). Another observation that can be made from the results of Harding (2018a) is that the $i=0$ terms provide a good approximation for each velocity component (e.g. within a few percent) whenever $\Rey^{2}\lessapprox 10/\epsilon$. Since $\epsilon\ll 10^{-1}$ is typical in spiral microfluidic ducts one can generally expect the $i=0$ terms to provide a reasonable approximation up to at least $\Rey=O(10)$. Furthermore, although (36) implies the secondary flow scales with $U_{m}^{2}$, and therefore vanishes quickly in the limit of small flow rate, it is important to note that the inertial lift force also scales with $U_{m}^{2}$. A consequence of this is the $i=0$ terms in the expansion of the background flow are sufficient to approximate the cross-sectional force on the particle when the flow rate is small. To make this clear we further examine the expansion of $\boldsymbol{F}_{p}^{\prime}$ into its different parts depending on the axial and secondary components of the background flow.

Returning to the non-dimensional setting from Section 2 we have

[TABLE]

Given it is now clear that $\hat{\bar{\boldsymbol{u}}}_{s}\propto\ell^{2}\Rey/4aR$ it follows that $\boldsymbol{v}_{0,s}\propto\ell^{2}\Rey/4aR$, and consequently

[TABLE]

Thus the contribution of the secondary flow drag to $\boldsymbol{F}_{p}^{\prime}$ scales with the dimensionless constant $\kappa:=\ell^{4}/4a^{3}R$ which depends only on physical length scales. At the same time it can be seen that the $\boldsymbol{F}_{0}$ contribution does not change in magnitude with respect to physical lengths. In particular, note that $\boldsymbol{F}_{0}$ can be approximated as

[TABLE]

since the additional terms in (31) and (32) involving $\hat{\bar{\boldsymbol{u}}}_{s}$ vanish much faster than the terms involving $\hat{\bar{\boldsymbol{u}}}_{a}$, which we have retained. Hence, the relative magnitude of the secondary flow drag $\Rey_{p}^{-1}\boldsymbol{F}_{-1,s}$ relative to $\boldsymbol{F}_{0}$ is $\kappa=\ell^{4}/4a^{3}R$ is consistent with a simplified study of the forces on the particle in the limit of large $R$ and small flow rate (Harding, 2018b). Furthermore, it follows that we need only a leading order approximation of $\hat{\bar{\boldsymbol{u}}}_{s}$ to estimate $\boldsymbol{F}_{-1,s}$ and a leading order approximation of $\hat{\bar{\boldsymbol{u}}}_{a}$ to estimate $\boldsymbol{F}_{0}$. That is, it suffices to take

[TABLE]

The scaling of $\Rey_{p}^{-1}\boldsymbol{F}_{-1,s}$ with $\kappa$ makes it clear that, as $\epsilon=\ell/(2R)\rightarrow 0$, the relative effect of the secondary flow drag becomes insignificant leaving only $\boldsymbol{F}_{0}$. Furthermore, $\boldsymbol{F}_{0}$ approaches the inertial lift force obtained in a straight duct since

[TABLE]

where $\hat{\bar{\boldsymbol{u}}}_{a}$ converges to a classical Poiseuille flow profile and $\boldsymbol{v}_{0,a}$ converges to the appropriate disturbance flow. On the other hand, with fixed physical dimensions (and thus fixed $\epsilon$) the relative size of the secondary drag force in relation to the other forces making up $\boldsymbol{F}_{p}^{\prime}$ remains constant. This implies that in the limit of small flow rate the interaction between the inertial lift force and secondary flow drag force converge to something which is different than that in the case of a straight duct. In Section 4 we examine this particular limit in detail for several cross-sections, particle sizes and bend radii.

4 Computational results and discussion

A finite element method is used to solve each Stokes problem, implemented using the FEniCS software (Logg et al., 2012; Alnæs et al., 2015). The computational domain consists of only a portion of the curved duct, in a neighbourhood of the particle, discretised by a mesh consisting of $O(10^{6})$ tetrahedra refined such that those elements near the particle have much smaller scale than those far away (since the disturbance varies most near the particle). Continuous Taylor–Hood elements are used, specifically a quadratic and linear polynomial basis for the velocity and pressure respectively. The Dirichlet boundary conditions on the duct walls and particle surface are enforced explicitly (i.e. at the linear algebraic level) and natural (stress-free) boundary conditions are applied at the ends of the (truncated) curved duct. The approximation of background flow terms $\hat{\bar{\boldsymbol{u}}}_{a},\hat{\bar{\boldsymbol{u}}}_{s}$ as in (38) are precomputed using the Rayleigh–Ritz method described in Harding (2018a) and subsequently interpolated in the finite element space. With $\boldsymbol{F}_{p}^{\prime}$ and drag coefficients $C_{r^{\prime}},C_{z^{\prime}}$ estimated at a large number of (equidistant) points in the cross-section we then construct a bivariate cubic spline to interpolate the data. The smooth interpolants are then used for subsequent analysis, for example the estimation of particle trajectories via the solution of (35). The behaviour of particles suspended in flow through curved ducts having square, rectangular and trapezoidal shaped cross-sections is examined. We use the approximation for low flow rate developed in Section 3 and particles are assumed to be neutrally buoyant (i.e. $\rho_{p}=\rho$).

4.1 Ducts having a square cross-section

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 $\epsilon^{-1}\in\{640,320,160,80,40\}$). Figures 2 and 3 show the force and trajectory plots respectively for the specific bend radius $R=160$. Results for the smallest particle ($a=0.05$) are depicted in Figures 2(a) and 3(a). Here $\boldsymbol{F}_{p}^{\prime}$ is dominated by the secondary flow drag since $\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 $\ell^{4}/4a^{3}R=25$) are depicted in Figures 2(b) and 3(b). Three equilibria can again be identified from the force plot, each of which is unstable, similar to those in (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 $\ell^{4}/4a^{3}R\approx 7.41$) are depicted in Figures 2(c) and 3(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 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 $\ell^{4}/4a^{3}R=3.125$) are depicted in Figures 2(d) and 3(d) and are similar to those in Figures 2(c) and 3(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 $\kappa=\ell^{4}/4a^{3}R$. We briefly comment on some general trends. For very large $\kappa$ 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 $\kappa\rightarrow 0$ given fixed $a,\ell$) we recover the results for a straight duct, specifically the four stable equilibria near the centre of each wall as in Hood et al. (2015). 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.

The transition from non-focusing to focusing behaviour is quite interesting. As $R$ increases (and $\kappa$ 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 $\kappa\in[18,19]$) leaving behind the slow manifold and a stable equilibrium as depicted in Figures 4(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 $\kappa\in[13,14]$) as depicted in Figures 4(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 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 $\epsilon^{-1}\in\{640,320,160,80\}$). The force and trajectory plots for $W=4$ and $R=160$ are shown in Figures 5 and 6 respectively. Results for the smallest particle ($a=0.05$ and $\ell^{4}/4a^{3}R=200$) are depicted in Figures 5(a) and 6(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 width-wise) 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 (the green marker is hidden by the spiral trajectories). Thus, despite a strong influence from the secondary flow drag, the particles eventually becomes focused.

Results for the second smallest particle ($a=0.10$ and $\ell^{4}/4a^{3}R=25$) are depicted in Figures 5(b) and 6(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 5(c) and 6(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 5(d) and 6(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.

The results for different values of $R$ are similar to those described above given a similar value of $\kappa=\ell^{4}/4a^{3}R$ and so we remark on some general observations and trends. In the limit $R\rightarrow\infty$ (equivalently $\kappa=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. For the smaller size particles one additionally finds stable equilibria near the centre of the side walls, however very few particles are likely to be found here in practice because the attractive region is much smaller (Hood, 2016). 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 $\kappa$ 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 centre 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 $\kappa$ 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 7 and 8 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 7 and 8 with Figures 5(a,c) and 6(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.

To study how the focusing position of different size particles changes with respect to the bend radius, we interpolate $\boldsymbol{F}_{p}^{\prime}$ 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 $\epsilon^{-1}=2R/\ell$. Despite the wide range of $R$, each of $\bar{u}_{\theta,0},\bar{u}_{r,0},\bar{u}_{z,0}$ (used to approximate $\hat{\bar{\boldsymbol{u}}}_{a},\hat{\bar{\boldsymbol{u}}}_{s}$ as in (38)) 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 $\kappa=\ell^{4}/4a^{3}R$; changes in the remaining components of $\boldsymbol{F}_{p}^{\prime}$ are generally subtle. For a much finer sampling of $R$ we then analyse the force on the particle $\boldsymbol{F}_{p}^{\prime}$ and identify/track the horizontal location of stable equilibria. The results of this analysis is shown in Figure 9 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 $\epsilon^{-1}$ shown, suggesting that some care should be taken in choosing the bend radius depending on the size of particles to be separated.

For the duct with $W/H=2$ and with $\epsilon^{-1}\lessapprox 80$ the ordering going from the inside wall to the centre is from largest to smallest particle. As $\epsilon^{-1}$ increases there are many changes to the ordering as largest particles begin to migrate towards the centre whereas smaller ones migrate towards the wall. Note that for $\epsilon^{-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 $\epsilon^{-1}$. This is one of the extra equilibria that is expected to occur as the behaviour approaches that of a straight duct. For $\epsilon^{-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 $\epsilon^{-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 $\epsilon^{-1}$ until for $\epsilon^{-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 don’t appear to get as close to the centre over the range of $\epsilon^{-1}$ shown.

Figure 10 shows the same data plotted against the ratio $\kappa=\ell^{4}/4a^{3}R$ for both rectangular ducts. Recall that small $\kappa$ indicates that inertial lift forces are dominant and large $\kappa$ indicates that secondary flow drag is dominant. The focusing horizontal focusing location approximately collapses onto a single curve, particularly for $\kappa\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 $\kappa$ 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 $\kappa$ 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.

Experiments using a spiral duct with two different rectangular cross-sections were performed by Wu et al. (2012). Both ducts had a width of $W=500\mu$m but with different heights $H=90\mu$m and $H=120\mu$m. Beads with diameter $2a=10\mu$m and $2a=6\mu$m were suspended in flow through the ducts at flow rates of $1$mL$/$min and $2$mL$/$min 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 $\kappa$ (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 10(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. 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 sized 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 $\kappa=\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 $\kappa$. 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.

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. (2013) and Warkiani et al. (2014) 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 side walls). We consider a cross-section width $W=8$ and heights $H_{\text{left}}=1.5$ and $H_{\text{right}}=2.5$ at the left (inside) and right (outside) walls respectively (see figures below). 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\}$. Figures 11 and 12 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 11(a) and 12(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 (the green marker is somewhat hidden by the trajectories). 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 11(b) and 12(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 converge 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 11(c) and 12(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 slightly left of centre.

We again consider how the focusing location changes with respect to $R$ (or equivalently $\epsilon^{-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 $\epsilon^{-1}$ is shown in Figure 13(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 13 (a) there are three main regions that can be identified over the $\epsilon^{-1}$ considered. For $\epsilon^{-1}\lessapprox 80$ the smaller two particle sizes are focused near $r=2$, the largest 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\epsilon^{-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.

Looking carefully around the discontinuities in Figure 13(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 $\epsilon^{-1}\approx 400$ and the second smallest particle ($a=0.10$, green triangles) at $\epsilon^{-1}\approx 80$. The force and trajectory plots for the latter case ($a=0.10$ and $\epsilon^{-1}=80$) is shown in Figure 14. Comparing Figure 14(a) with Figure 11(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 $\epsilon^{-1}$ (or equivalently $R$) the relative location of the black and white contours in the mid-section shift further such that only the right most intersections remain thus leaving only the two stable equilibria on the right. The trajectories in Figure 14(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 13 (b) the focusing data is plotted against the ratio $\kappa=\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 $\kappa$. Additionally, for $\kappa<10$ and $\kappa>100$ the behaviour approximately collapses onto a single curve. However, for intermediate $\kappa$ 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 $\kappa$ 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 $\kappa$ is a little larger than where the minimum is achieved. Furthermore, this rapid change occurs at larger $\kappa$ 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 $\kappa$ (and is only briefly left of centre). Furthermore, its focusing behaviour appears to be reasonably smooth with respect to $\kappa$ 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 sized 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\mu$m and $2a=6\mu$m were performed by Wu et al. (2012). Both ducts had a width of $W=500\mu$m but with the different heights of $H_{\text{left}}=70\mu$m to $H_{\text{right}}=100\mu$m from the inner to outer wall in the smaller cross-section and $H_{\text{left}}=90\mu$m to $H_{\text{right}}=120\mu$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 $\kappa$ (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 13(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 $\kappa$ 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.

5 Conclusions

This paper develops a general model for the forces that govern the motion of a spherical particle suspended in flow through a curved duct. A key component in extending the approach of Hood et al. (2015) from straight ducts to curved ducts is the use of a rotating coordinate system as a frame of reference in which the flow is approximately steady. Additionally, an expansion of the background flow into axial and secondary components identifies how different components of the force on a particle are affected by the background flow. We performed further analysis on the special case of low flow rate and neutrally buoyant particles. We found that the the secondary flow drag scales with $\kappa=\ell^{4}/4a^{3}R$ relative to the inertial lift force. We computed the position and stability of equilibria for several different cross-sectional geometries. Further, a simple first order model of particle trajectories allows us to plot approximate trajectories and identify a slow manifold in many cases.

An analysis of the location of stable equilibria in rectangular and trapezoidal cross-sections demonstrates that $\kappa$ plays an important role in the general focusing behaviour of particles. We observed that a stable equilibria pair exists over a large range of bend radii and particle diameter. While several changes in the lateral ordering of focused particles are observed we find that the focusing behaviour approximately collapse onto a single curve when plotted against $\kappa$, particularly for rectangular cross-sections. The results suggest that this is the mechanism for the size based particle/cell separation observed in the experimental literature, providing a dimensionless parameter that may aid in the design and operation of microfluidic ducts under appropriate conditions.

While our study assumes the bend radius to be constant the results can be applied to provide some insight into the focusing behaviour within spiral microfluidic ducts. In particular, the final location of a particle in a spiral device can be estimated by looking at the focusing behaviour in a curved duct with (constant) bend radius matching that near the outlet of the spiral. Furthermore, our results that illustrate the change in particle focusing position with respect to $\epsilon^{-1}=2R/\ell$ illustrate where particles will focus towards as the bend radius changes throughout a spiral device. Lastly we note that the results could be reasonably applied to dilute suspensions of particles/cells provided the concentration is small enough that interactions between particles/cells can be neglected and the particles/cells are sufficiently rigid and sphere like in shape.

There are several ways in which this work could be extended. This study considered a single trapezoidal shaped duct and an interesting extension is to examine how the focusing dynamics evolve as the slope of the top wall is slowly increased starting from a rectangular cross-section. Further, one could investigate what happens when the trapezoidal cross-section is tallest at the inside wall as in the experiment of Sofela et al. (2018). Exploration of other cross-sectional geometries may also be interesting, for example a symmetric trapezoidal shape which we expect to eliminate the staggering of equilibria observed with the asymmetric trapezoidal shape. One could even go further and consider cross-sections in which the top and bottom walls have a more general polynomial shape.

It would be interesting to explore how a change in density of the particle influences the existence and stability of equilibria. In particular, if the density of the particle differs enough from that of the fluid that gravity becomes relevant then this will break the symmetry of the force on the particles in rectangular cross-sections. Additionally, there will be a noticeable change in the centrifugal force for a non-neutrally buoyant particle which will influence the horizontal location of equilibria. The focusing behaviour at higher flow rates is another important aspect to consider in order to achieve a reasonable throughput while maintaining a high separation efficiency. Finally, the simple first order trajectory model might be extended to a more complete second order model which takes into account axial acceleration/deceleration of the particle as it migrates within the cross-section. In particular it is reasonable to expect changes to the axial velocity to have an appreciable effect on the inertial lift force, particularly when $\kappa$ is large.

Acknowledgements

This research was supported by the Australian Research Council via a Discovery Project (DP160102021) and a Future Fellowship (FT160100108) to YMS, and by the Simons Foundation via a Math + X grant (510776) to ALB. Funding from the University of Adelaide via an establishment grant to YMS is also gratefully acknowledged. The results were computed using supercomputing resources provided by the Phoenix HPC service at the University of Adelaide.

Appendix A Nomenclature

Table 1 summarises the nomenclature used throughout the paper. We note that a large number of variables listed are also introduced in the context of a rotating reference frame which is always denoted with a prime (e.g. as in $\ast^{\prime}$). For brevity these variants are excluded from the table. Dimensionless forms of variables are denoted throughout the paper with a hat/caret (i.e. as in $\hat{\ast}$) and are omitted from the table for brevity. Perturbation expansions are also introduced throughout the paper via a subscript, for example $\hat{\boldsymbol{v}}^{\prime}=\boldsymbol{v}_{0}+\Rey_{p}\boldsymbol{v}_{1}+O(\Rey_{p}^{2})$, and are also omitted from the table. We additionally note that the caret and prime are dropped from perturbation variables for ease of readability since these are always taken to refer to dimensionless quantities in the rotating reference frame. Lastly, some terms are decomposed into separate parts dependent on the axial and secondary components of the background flow typically denoted as $\ast_{a},\ast_{s}$. With the exception of $\bar{\boldsymbol{u}}_{a},\bar{\boldsymbol{u}}_{s}$ these too are omitted from the table.

Appendix B Estimation of the background flow

We provide a brief account of the derivation of the equations governing the background flow from Harding (2018a) which we utilise in Sections 3 and 4. In the notation of Section 2.1 the Navier–Stokes equations for steady flow through a curved duct in cylindrical coordinates are

[TABLE]

Solving these equations is most easily done via a non-dimensionalisation that reflects the scaling of the background flow rather than the flow near a particle. In particular, the spatial coordinates are non-dimensionalised as $(r,z)=(\ell/2)(\grave{r},\grave{z})$ and the fluid velocity components are non-dimensionalised as $\bar{u}_{\theta}=U_{m}\grave{\bar{u}}_{\theta}$ and $(\bar{u}_{r},\bar{u}_{z})=(\epsilon\Rey/2)U_{m}(\grave{\bar{u}}_{r},\grave{\bar{u}}_{z})$. Observe that $R+r=R(1+\epsilon\grave{r})$ where $\epsilon=\ell/2R$ and for convenience we define $\grave{R}:=(1+\epsilon\grave{r})$ such that $R+r=R\grave{R}$. We additionally introduce the (dimensionless) stream function $\Phi$ such that $\partial\Phi/\partial\grave{z}=-\grave{R}\grave{\bar{u}}_{r}$ and $\partial\Phi/\partial\grave{r}=-\grave{R}\grave{\bar{u}}_{z}$. The equations (39) can then be reduced to

[TABLE]

where $K=\epsilon\Rey^{2}/4$ and $\Delta:=\partial^{2}/\partial\grave{r}^{2}+\partial^{2}/\partial\grave{z}^{2}$. For small $K$ we can introduce a perturbation expansion of $\grave{\bar{u}}_{\theta}$ and $\Phi$ with respect to $K$, specifically

[TABLE]

For $K\ll 1$ it is sufficient to estimate the cross-sectional forces on a particle using the approximation $\grave{\bar{u}}_{\theta}\approx\grave{\bar{u}}_{\theta,0}$ and $\Phi\approx\Phi_{0}$ where

[TABLE]

A Rayleigh–Ritz method for approximating these components, which is used in the numerics in Section 4, is described in Harding (2018a).