2.1. Mathematical model

Harding et al. (Reference Harding, Stokes and Bertozzi2019) developed an asymptotic model to calculate the forces acting on a neutrally buoyant spherical particle suspended in fluid flowing through a constant curvature duct at sufficiently small flow rates. By balancing these leading-order forces with viscous Stokes’ drag they constructed a first-order quasi-steady model for the trajectory of the particle. In this work, we adapt their model and apply it to spiral ducts with slowly varying curvature. By slowly varying curvature we mean that the background fluid flow through any cross-section of the spiral duct does not differ significantly from that through a circular duct with identical curvature. We expect the quasi-steady approximation of inertial migration to remain valid for such spiral ducts when approximating the background flow accordingly. The quality of this approximation will be demonstrated in § 2.3. Below, we briefly review the model of Harding et al. (Reference Harding, Stokes and Bertozzi2019); further details may be obtained from that paper.

As shown in figure 1, consider a neutrally buoyant spherical particle of radius $a$ suspended in fluid flowing through a spiral duct with a uniform rectangular cross-section. This specific example illustrates an Archimedean spiral for which the bend radius, $R(\theta )$, varies with the azimuthal angle, $\theta$. The rectangular cross-section has width $W$ and height $H$. We note that our model neglects the effects of any inlets and outlets, i.e. we analyse the dynamics of a particle assuming that the flow is fully developed throughout the duct. In a circular duct (with constant bend radius $R(\theta )=R$), and in the absence of the particle, a constant pressure gradient along the duct drives a steady incompressible fluid flow (Dean Reference Dean1927; Dean & Hurst Reference Dean and Hurst1959) which has been studied extensively for ducts with a rectangular cross-section (Winters Reference Winters1987; Yamamoto et al. Reference Yamamoto, Wu, Hyakutake and Yanase2004; Harding & Stokes Reference Harding and Stokes2018). In addition to the axial component of the fluid flow, the curvature induces a secondary flow within the cross-section consisting of a counter-rotating vortex pair. The presence of a particle in the duct disturbs this background flow and one can define disturbance pressure and velocity fields by taking the difference between the fields in the presence and absence of the particle. The resulting dimensionless dynamical equations for the disturbance pressure field $q$ and disturbance fluid velocity field $\boldsymbol {v}$ in a rotating frame (moving with the particle) are given by

(2.1)\begin{align} \left.\begin{gathered} - \boldsymbol{\nabla} q + \nabla^2\boldsymbol{v}=Re_p ((\boldsymbol{v}+\bar{\boldsymbol{u}}-\varTheta (\boldsymbol{e}_{z}\times \boldsymbol{x}))\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} \bar{\boldsymbol{u}} +\varTheta(\boldsymbol{e}_z\times\boldsymbol{v})) \quad \text{on}\ \boldsymbol{x} \in \mathcal{F}, \\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} = 0 \quad \text{on}\ \boldsymbol{x} \in \mathcal{F},\\ \boldsymbol{v}=\boldsymbol{0} \quad \text{on}\ \boldsymbol{x} \in \partial \mathcal{D},\\ \boldsymbol{v}={-}\bar{\boldsymbol{u}}+\varTheta(\boldsymbol{e}_z\times\boldsymbol{x})+ \boldsymbol{\varOmega}_p\times (\boldsymbol{x}-\boldsymbol{x}_p) \quad \text{on}\ \boldsymbol{x} \in \partial \mathcal{F} \setminus \partial \mathcal{D}. \end{gathered}\right\} \end{align}

Here, $\bar {\boldsymbol {u}}$ is the background fluid flow velocity in the absence of the particle, $\boldsymbol {\varOmega }_p$ is the particle spin and $\boldsymbol {e}_{z}$ is the unit vector in the vertical $z$ direction; all quantities are in a reference frame which is rotating about the $z$ axis at a rate $\varTheta :=\partial \theta _p/\partial t$ which is assumed to be (nearly) constant. The domain $\mathcal {D}$ denotes the interior of the duct, $\partial \mathcal {D}$ denotes the boundaries of the duct, $\mathcal {F}:=\{\boldsymbol {x}\in \mathcal {D}:|\boldsymbol {x}-\boldsymbol {x}_p|\geq 1\}$ is the dimensionless fluid domain in the presence of the particle and $\partial \mathcal {F} \setminus \partial \mathcal {D}=\{\boldsymbol {x}:|\boldsymbol {x}-\boldsymbol {x}_p|=1\}$ is the surface (boundary) of the particle. The dimensionless parameter $Re_p$ in (2.1) is the particle Reynolds number defined as $Re_p=Re\,(a/l)^2$, where $Re=\rho U_m l/\mu$ is the flow Reynolds number, $\rho$ is the particle/fluid density, $\mu$ is the fluid's dynamic viscosity, $l=\min \{W,H\}$ is the characteristic length scale associated with the duct cross-section and $U_m$ is a characteristic velocity scale chosen to be the maximum axial velocity of the background fluid flow. The variables in (2.1) have been scaled with physical parameters as in Harding et al. (Reference Harding, Stokes and Bertozzi2019): the particle radius $a$ for $\boldsymbol {x}$, $l/U_m$ for time $t$, $U_m a/l$ for the velocities $\boldsymbol {v}$ and $\bar {\boldsymbol {u}}$, $U_m/l$ for the particle spin $\boldsymbol {\varOmega }_p$ and the angular velocity $\varTheta$, and $\mu U_m/l$ for the disturbance pressure $q$.

The disturbed fluid flow exerts force and torque on the spherical particle. The dimensionless force (scaled by $\rho U_m^2 a^4/l^2$) and torque (scaled by $\rho U_m^2 a^5/l^2$) acting on the particle in the rotating frame are given by

(2.2)\begin{align} \boldsymbol{F}&={-}\frac{4{\rm \pi}}{3}\varTheta^2 (\boldsymbol{e}_z\times (\boldsymbol{e}_z\times \boldsymbol{x}_p))+ \int_{|\boldsymbol{x}-\boldsymbol{x}_p|<1} \bar{\boldsymbol{u}}\boldsymbol{\cdot} \boldsymbol{\nabla}\bar{\boldsymbol{u}}\,\text{d}V \nonumber\\ &\quad +\frac{1}{Re_p} \int_{|\boldsymbol{x}-\boldsymbol{x}_p|=1} (-\boldsymbol{n})\boldsymbol{\cdot} ({-}q\boldsymbol{I}+\boldsymbol{\nabla}\boldsymbol{v}+ \boldsymbol{\nabla}\boldsymbol{v}^{\rm T})\,\text{d}S, \end{align}

(2.3)\begin{align} \boldsymbol{T}&={-}\frac{8{\rm \pi}}{15}\varTheta (\boldsymbol{e}_z\times\boldsymbol{\varOmega}_p)+ \int_{|\boldsymbol{x}-\boldsymbol{x}_p|<1} (\boldsymbol{x}-\boldsymbol{x}_p)\times (\bar{\boldsymbol{u}}\boldsymbol{\cdot} \boldsymbol{\nabla}\bar{\boldsymbol{u}})\,\text{d}V \nonumber\\ &\quad +\frac{1}{Re_p} \int_{|\boldsymbol{x}-\boldsymbol{x}_p|=1} (\boldsymbol{x}-\boldsymbol{x}_p)\times((-\boldsymbol{n}) \boldsymbol{\cdot}({-}q\boldsymbol{I}+\boldsymbol{\nabla}\boldsymbol{v}+ \boldsymbol{\nabla}\boldsymbol{v}^{\rm T}))\,\text{d}S. \end{align}

Taking the particle Reynolds number as a small parameter and performing a perturbation expansion of the disturbance flow in powers of $Re_p$ we get

(2.4)\begin{equation} \left.\begin{gathered} \boldsymbol{v} = \boldsymbol{v}_0+Re_p \boldsymbol{v}_1+{O}(Re^2_p), \\ q = q_0+Re_p q_1+{O}(Re^2_p). \end{gathered}\right\} \end{equation}

Substituting into (2.1) gives us a zeroth-order system for $q_0,\boldsymbol {v}_0$, and a first-order system for $q_1,\boldsymbol {v}_1$. The zeroth-order system captures all the non-zero boundary conditions and the first-order system captures the most significant inertial contribution to the complete equations. Substituting the perturbation expansions in (2.2) results in the following equation for the leading-order cross-sectional force on the quasi-static particle:

(2.5)\begin{equation} \boldsymbol{F}_s=(\boldsymbol{e}_r\boldsymbol{\cdot}(Re^{{-}1}_p\boldsymbol{F}_{{-}1,s}+ \boldsymbol{F}_0))\boldsymbol{e}_r+(\boldsymbol{e}_z \boldsymbol{\cdot}(Re^{{-}1}_p \boldsymbol{F}_{{-}1,s}+\boldsymbol{F}_0))\boldsymbol{e}_z. \end{equation}

Here, $\boldsymbol {F}_{-1}$ is the leading-order force consisting of drag from the background flow, with the subscript ‘s’ denoting the secondary (cross-sectional) component, while $\boldsymbol {F}_0$ is the leading-order contribution of the inertial lift force. Then, using Stokes’ law, we balance the cross-sectional force in (2.5) with viscous Stokes’ drag and thus construct a first-order quasi-steady model for the trajectory of the particle giving the following dynamical equations of motion:

(2.6a–c)\begin{equation} \frac{\text{d} r_p}{\text{d} t}={-}Re_p\frac{F_{s,r}}{C_r},\quad \frac{\text{d} z_p}{\text{d} t}={-}Re_p\frac{F_{s,z}}{C_z}\quad \text{and}\quad \frac{\text{d} \theta_p}{\text{d} t}=\frac{\bar{u}_a}{R/a+r_p}, \end{equation}

where $\bar {u}_a$ is the axial component of the background fluid flow, $F_{s,r}=\boldsymbol {F}_s\boldsymbol {\cdot }\boldsymbol {e}_r$ and $F_{s,z}=\boldsymbol {F}_s\boldsymbol {\cdot }\boldsymbol {e}_z$ are the radial and the vertical components of the cross-sectional force, respectively, with corresponding drag coefficients $C_r$ and $C_z$ that vary with the particle's position in the cross-section. Thus, this model estimates the particle velocities within a circular duct with each of the forces implicitly depending on the radius of curvature (which coincides with the bend radius $R$) in addition to the particle radius $a$. We note that our model does not include particle–particle interactions. Hence, the results presented in this manuscript are valid in the regime of low particle density where particle–particle interactions can be safely neglected.

We are able to extend this quasi-steady model to spiral ducts with slowly varying curvature by considering $R$ as a function of $\theta$ to describe the local radius of curvature (which must be distinguished from any notion of bend radius in non-circular duct geometries, see § 2.3 for more details). The radius of curvature is the appropriate quantity to determine local fluid behaviour as it is a well-defined geometrical construct which applies to any curve, whereas the notion of bend radius relies on having a well-defined origin (or other point of reference) which does not apply to generalised curves. Thus, while it is convenient to describe circular and Archimedean spiral ducts via a ‘bend radius’, we need to consider the radius of curvature when modelling the fluid.

The circular duct model is solved numerically using a finite element method to compute the forces acting on the particle. We refer the reader to Harding et al. (Reference Harding, Stokes and Bertozzi2019) for additional details on the numerical framework. Once the relevant forces have been pre-computed at numerous sample points over the cross-section (and for numerous system parameter values), smooth interpolants of $C_r$, $C_z$, $F_{s,r}$, $F_{s,z}$ are constructed and the particle dynamics is simulated using the MATLAB solver ode45.

For the results presented in this manuscript, we use the following non-dimensional scales for the system variables: dimensionless radius of curvature (or bend radius, as appropriate) $\tilde {R}=2R/H$, dimensionless particle size $\tilde {a}=2a/H$, dimensionless time $\tilde {t}=U_m t/H$ and cross-sectional co-ordinates $\tilde {r}=2r/H$ and $\tilde {z}=2z/H$.

2.2. Spiral duct geometry

Our illustrations thus far have utilised an Archimedean spiral, since these are readily recognised and understood. However, to extend the circular duct model to spiral ducts, it is convenient to instead consider an Archimedean-like spiral duct whose centreline is described by the parametric equation

(2.7)\begin{equation} \boldsymbol{\tilde{r}}_I(\varphi)=[\tilde{x}(\varphi),\tilde{y}(\varphi)]= \tilde{R}(\varphi)[\cos(\varphi),\sin(\varphi)]+\beta[-\sin(\varphi),\cos(\varphi)-1], \end{equation}

where

(2.8)$$\begin{gather} \tilde{R}(\varphi)=\tilde{R}_{start}+\beta\varphi, \end{gather}$$

(2.9)$$\begin{gather}\beta = \frac{\tilde{R}_{end}-\tilde{R}_{start}}{2{\rm \pi} N_{turns}}, \end{gather}$$

with $\tilde {R}_{start},\tilde {R}_{end}$ denoting the radius of curvature at the start and end of the spiral duct, respectively, and $N_{turns}$ denoting the number of turns. This curve corresponds to the involute of a circle (having radius $\beta$) and, although visually similar, differs in several ways from a traditional Archimedean spiral described by

(2.10)\begin{equation} \boldsymbol{\tilde{r}}_A(\theta)=[\tilde{x}(\theta),\tilde{y}(\theta)]= \tilde{R}(\theta)[\cos(\theta),\sin(\theta)], \end{equation}

with $\tilde {R}(\theta )$ as in (2.8).

For the traditional Archimedean spiral of (2.10), the radial distance from the origin, often referred to as the bend radius, varies linearly with respect to $\theta$. Accordingly, $\tilde {R}_{start},\tilde {R}_{end},\tilde {R}(\theta )$ should be interpreted as bend radii in the context of the Archimedean spiral. However, it is important to note that the bend radius of an Archimedean spiral does not equal its radius of curvature, and it is the latter which is fundamental to describing/modelling the local fluid flow. In contrast, for the involute spiral of (2.7) it is the radius of curvature that varies linearly with respect to $\varphi$ and $\tilde {R}_{start},\tilde {R}_{end},\tilde {R}(\varphi )$ can be interpreted as radii of curvature at the appropriate locations. The trade-off is that $\varphi$ no longer describes the polar angle.

There are a number of additional differences between the two curves worth noting. The involute spiral has a closed-form arclength parametrisation whereas the Archimedean spiral does not. The principal unit normal of the involute spiral is $\boldsymbol {N}_I(\varphi )=[-\cos (\varphi ),-\sin (\phi )]$ and is directed from $\tilde {\boldsymbol {r}}_I(\varphi )$ to the point $\beta [-\sin (\varphi ),\cos (\varphi )-1]$. The principal unit normal for the Archimedean spiral is more complex and, importantly, is not directed from $\boldsymbol {\tilde {r}}_A(\theta )$ towards the origin. These differences between the two curves are subtle when $\beta$ is small, but they become important when considering how to update the value of $\tilde {R}$ as a particle traverses a spiral duct having positive width. The simplest way one might describe the sidewalls in both spirals is to replace $\tilde {R}({\cdot })$ with $\tilde {R}({\cdot })\pm \tilde {W}/2$ in each of (2.7) and (2.10). For the involute spiral this leads to a consistent width of $\tilde {W}$ as measured with respect to the principal unit normal of $\tilde {\boldsymbol {r}}_I(\varphi )$, but the same cannot be said for the Archimedean spiral. Moreover, when considering how to update $\tilde {R}$ as a particle traverses the involute spiral we have

(2.11)\begin{equation} \frac{{\rm d}}{{\rm d}\tilde{t}}\tilde{R}(\varphi)=\beta\frac{{\rm d}\varphi}{{\rm d}\tilde{t}}= \beta\frac{{\rm d}\varphi}{{\rm d} s_I}\frac{{\rm d} s_I}{{\rm d}\tilde{t}}, \end{equation}

where $s_I(\varphi )$ is the arclength along the involute spiral. The term $\textrm {d} s_I/\textrm {d}\tilde {t}$ is simply the velocity of the particle down the main axis (i.e. aligned with the tangent vector $\hat {\boldsymbol {T}}_I(\varphi )=[-\sin (\varphi ),\cos (\varphi )]$). For $\textrm {d}\varphi /\textrm {d} s_I$ we note that $\textrm {d} s_I/\textrm {d}\varphi =\tilde {R}(\varphi )+\tilde {r}_p$, where the addition of $\tilde {r}_p$ accounts for the radial coordinate of the particle within the cross-section. For the Archimedean spiral, we similarly have

(2.12)\begin{equation} \frac{{\rm d}}{{\rm d}\tilde{t}}\tilde{R}(\theta)= \beta\frac{{\rm d}\theta}{{\rm d}\tilde{t}}=\beta\frac{{\rm d}\theta}{{\rm d} s_A}\frac{{\rm d} s_A}{{\rm d}\tilde{t}}, \end{equation}

where $s_A(\varphi )$ is the arclength along the Archimedean spiral. The term $\textrm {d} s_A/\textrm {d}\tilde {t}$ is again the particle velocity down the main axis; however, one must recognise that the tangent vector is given by the more complex expression

(2.13)\begin{equation} \hat{\boldsymbol{T}}_A(\theta)=\frac{\tilde{R}(\theta)[-\sin(\theta),\cos(\theta)]+ \beta[\cos(\theta),\sin(\theta)]}{\sqrt{\tilde{R}(\theta)^2+\beta^2}}. \end{equation}

Moreover, for $\textrm {d}\theta /\textrm {d} s_A$ we note that $\textrm {d} s_A/\textrm {d}\theta =\sqrt {\tilde {R}(\theta )^2+\beta ^2}+\tilde {r}_p$. Upon accounting for these subtleties, one must still determine the radius of curvature of the Archimedean spiral to ensure the correct migration forces are sampled from our circular duct model. All together, these differences make the involute spiral much better suited to extending our curved duct model of inertial particle migration. We will use $\tilde {\boldsymbol {r}}(\varphi )=\tilde {\boldsymbol {r}}_I(\varphi )$ and $\tilde {R}=\tilde {R}(\varphi )$ going forwards, where $\tilde {R},\tilde {R}_{start},\tilde {R}_{end}$ should always be interpreted as radii of curvature (with the exception of appendix where we compare results with an Archimedean spiral).

When designing a spiral duct, upper and lower bounds on the number of turns, ${N}_{turns}$, are determined by both geometrical and physical/practical constraints. From a purely geometrical point of view, the maximum number of turns, $N_{\max }$, is constrained by the fact that the spiral duct, having non-zero width along the spiral curve, should not intersect itself. This means that, after one turn of the spiral, $\tilde {R}(\varphi )$ should change at least by the width of the cross-section, resulting in the constraint $2{\rm \pi} \beta >\tilde {W}$, or equivalently $N_{\max }<|\tilde {R}_{end}-\tilde {R}_{start}|/\tilde {W}$, where $\tilde {W}=W/H$ is the aspect ratio of the rectangular cross-section. In addition to this geometrical constraint, there may be additional practical experimental/design constraints that further limit the maximum number of turns. For example, the required pressure difference to drive fluid flow through the duct increases with the length of the duct and there may be an upper limit on the pressure difference to preserve the structural integrity of a microfluidic device. The minimum number of turns $N_{min}$ is constrained by two factors: (i) the maximum change in radius of curvature that can be considered reasonable in our slowly varying curvature model, and (ii) the time/length required for the particles (initially randomly distributed in the cross-section) to focus (sufficiently close) to their particle attractors. Both of these considerations are also influenced by the flow rate at which the device is expected to operate at.

2.3. Numerical test for the slowly varying curvature approximation

Steady flow through an Archimedean spiral duct having a rectangular cross-section with large aspect ratio was considered in Harding & Stokes (Reference Harding and Stokes2018). (Non-rectangular cross-sections were also considered but in the context of this work only the rectangular cross-sections are of interest.) An approximation was developed via a regular perturbation expansion with respect to both the curvature parameter and the height-to-width ratio. It was found that the leading-order solution depended on the local curvature of the spiral but not its rate of change. First-order corrections with respect to both perturbation parameters were also considered, but these too did not introduce any non-trivial dependence on the rate of change of the spiral curvature. This provides strong evidence that spiral duct flow does not differ significantly from circular duct flow (when using the appropriate curvature parameter at any given point in the spiral). More specifically, it suggests that non-trivial dependence of the flow with respect to our $\beta$ parameter only occurs at order $O(\epsilon ^2)$ or smaller (relative to the magnitude of each leading-order component).

To test this further we conducted a numerical simulation of fluid flow through a relatively short spiral duct having a somewhat large value of $\beta$. Specifically, we considered a set-up with $\tilde {R}_{start}=100$, $\tilde {R}_{end}=20$, $N_{turns}=1/2$, $W=4$ and $H=2$, i.e. such that $|\beta |=80/{\rm \pi} \approx 25.5$. We generated a structured mesh of this domain consisting of roughly 780 000 tetrahedra and then solved the Navier–Stokes equations using the finite element method with Taylor–Hood elements (utilising the FEniCS computing platform Logg, Mardal & Wells Reference Logg, Mardal and Wells2012; Alnaes et al. Reference Alnaes, Blechta, Hake, Johansson, Kehlet, Logg, Richardson, Ring, Rognes and Wells2015). At different cross-sections we then compared the velocity field with that which is obtained by solving the equations of steady flow through a circular duct (Harding Reference Harding2019, Reference Harding2022). At the cross-section of the spiral duct for which the radius of curvature is $50$, and with a flow Reynolds number of $100$, we measured a relative $L_2$ difference in the axial, radial and vertical velocity components of the circular and spiral duct flows as $3.3\,\%$, $5.5\,\%$ and $3.1\,\%$, respectively. The error decreased for smaller Reynolds numbers, for example with $Re=50$ we observe relative differences of $0.53\,\%$, $0.99\,\%$ and $2.0\,\%$, respectively.

In a second computation we generated a structured mesh of a spiral duct with ${\tilde {R}_{start}=50}$, $\tilde {R}_{end}=300$, $N_{turns}=1$ and other parameters as above, such that, now, $|\beta |=125/{\rm \pi} \approx 40.8$. The same number of tetrahedra were used in this case. At the cross-section for which the radius of curvature is $200$, and again with a flow Reynolds number of $100$, we measured a relative $L_2$ difference in the three velocity components as $0.25\,\%$, $2.6\,\%$ and $0.46\,\%$. In this instance the difference did not decrease significantly when decreasing the Reynolds number, which suggests that the dominant error is due to the mesh resolution. A similar observation was made for other cross-sections (which are not too close or far from the ends where there are inlet/outlet effects). A precise analysis of the difference between circular and spiral duct flows is the subject of ongoing work, but together these results indicate that $\beta =O(10)$ can be considered small enough that the spiral geometry is sufficiently slowly varying that the fluid flow through any given cross-section can be well approximated by circular duct flow with appropriately chosen curvature parameter, for Reynolds numbers up to $O(100)$.

4.1. Spiral ducts with a single point attractor

We start by analysing the dynamics of particles in spiral ducts where there is a single particle attractor. In ducts with rectangular cross-sections, this regime arises at smaller radii of curvature, where secondary drag forces dominate inertial lift forces and one obtains either: (i) a single stable node near the inner wall of the curve duct (e.g. see $\tilde {R}=3000$ in figure 2c), or (ii) a vertically symmetric pair of stable nodes/spirals (e.g. see $\tilde {R}=1000$ in figure 2c).

A typical example of particle focusing dynamics in an in-spiral duct (i.e. the particle travels in the direction of decreasing radius of curvature) consisting of a vertically symmetric pair of stable spirals at the same horizontal $\tilde {r}$ location, is shown in figure 3. Here, we compare the particle focusing dynamics for spirals having the same starting ($\tilde {R}_{start}$) and ending ($\tilde {R}_{end}$) radii of curvature but a different number of turns ($N_{turns}$). Note that, even in the absence of bifurcations, the particle equilibria are not static throughout the spiral ducts; they move in response to the changing radius of curvature of each cross-section. We observe that, for the spiral having fewer turns, the combination of insufficient duct length and continuous movement of the point attractors inhibits complete particle focusing and a spread in the particle distribution is observed at the end of the spiral duct (see figure 3a). Alternatively, with more turns, one observes complete focusing of particles at the end of the spiral duct (see figure 3b). For cross-sections G through to H, the focused particle cluster tracks the moving attractors closely and there is a small ‘lag’ between the centre of the focused clusters and the centre of the attractors. Hence, in a typical spiral duct with a single horizontal attractor and sufficient turns, the particle focusing is similar to circular ducts whose radius of curvature coincides with the local radius of curvature at the end of the spiral; in both cases the particles focus in close proximity to the single attractor whose location is determined by the radius of curvature at the end of the spiral. Nevertheless, spiral duct geometries have an advantage over circular ducts since one can increase the duct length by having multiple non-intersecting turns (which is not practically possible in circular ducts, although circular ducts might be wrapped around in the form of a helix to obtain multiple turns while avoiding self-intersection, in which case the effects of torsion may become important and our present model may not be applicable). This becomes particularly useful at smaller radii of curvature where one turn of a circular duct is not sufficient for particle focusing (e.g. see figure 10a).

Figure 3.

Particle focusing dynamics in an in-spiral duct with a $2\times 1$ rectangular cross-section that has a vertically symmetric pair of point attractors. The system parameters are $R_{start}=300$, $R_{end}=50$, $Re=50$ and $\tilde {a}=0.10$. Snapshots of the cross-section are shown at (A,E) $\tilde {R}=300$, (B,F) $200$, (C,G) $100$ and (D,H) $50$, for column (a) $N_{turns}=1$ and column (b) $N_{turns}=4$, respectively. The coloured circles denote the type of particle equilibria with cyan for stable spirals (point attractors) and yellow for a saddle point. The grey circles denote the particle positions while the grey curves denote their trajectories. If the centre of a particle lies on the dashed rectangle, it will touch at least one wall of the duct.

4.2. Spiral ducts with a single limit-cycle pair

Next, we consider an in-spiral duct with square cross-section having a single one-dimensional attractor in the horizontal direction, as shown in figure 4. Specifically, we have a vertically symmetric pair of stable limit cycles (black curves) which persists throughout the entire length of the spiral. However, the decreasing radius of curvature throughout the in-spiral duct leads to a continuous change in the size of the limit cycle and the dynamics along it. The limit cycle is large in extent near the inlet of the spiral and shrinks progressively along the spiral, while the period to traverse the limit cycle simultaneously decreases. We make an interesting observation regarding the phase of the limit cycle occupied by particles based on their initial location in the cross-section. We find that most particles that do not start near the unstable spiral equilibria or along the horizontal centreline of the duct, aggregate around the same phase of the limit cycle (see Movie $1$). One of the key factors for this particle aggregation is that initially, the majority of particles are ‘squeezed’ along the horizontal centreline where their dynamics slow down. This results in aggregation of particles as they approach the stable limit cycles. This phenomenon is robust with respect to the number of turns, as can be seen in figure 4(a,b) respectively (compare Movies $1$ and $2$), which allows the phase occupied by the majority of particles within the final cross-section to be manipulated by modifying the number of turns accordingly.

Figure 4.

Particle focusing dynamics in a spiral duct with square cross-section having a vertically symmetric pair of limit-cycle attractors. (a,b) Particle dynamics for in-spiral ducts having (a) $N_{turns}=4$ (see supplementary Movie 1 available at https://doi.org/10.1017/jfm.2024.487) and (b) $N_{turns}=5$ (see Movie 2). Cross-sectional images at the (A,C) start and (B,D) end of the spiral show particle equilibria (dark coloured circles with purple for unstable spirals and yellow for saddle points) along with particle positions (light coloured circles) whose colour is based on the phase angle $\phi$ occupied by the particles along the limit cycle (black curves) at the end of the spiral. If the centre of a particle lies on the dashed square, it will touch at least one wall of the duct. The system parameters are fixed to $R_{start}=1250$, $R_{end}=500$, $\tilde {a}=0.05$ and $Re=50$.

4.3. Spiral ducts with multiple particle attractors

Multiple point attractors with different horizontal locations can coexist in circular ducts with rectangular cross-sections. When this occurs, the cross-sectional domain can be partitioned based on the basin of attraction for each point attractor. This allows one to readily identify where a particle will converge based on its initial location. However, when spiral ducts contain multiple point attractors throughout (with no bifurcations), both the point attractors and the basins of attraction move in response to the continuously changing radius of curvature. This can lead to some particles finding themselves in a different basin of attraction to the one they started in and this thereby induces non-trivial focusing effects.

Figure 5(a,b) compares a particle's trajectory in two spiral duct geometries having different numbers of turns. The duct has a $2\times 1$ rectangular cross-section and the spiral is such that multiple point attractors persist throughout. Two horizontally distinct stable nodes are present, one on the horizontal centreline near the inner wall, and a vertically symmetric pair of stable nodes a bit further away from the inner wall. Two saddle points between the three stable node attractors separate their respective basins of attraction. In figure 5(a) we show the trajectory of a particle starting out near the top inner (left) corner of the cross-section, specifically within the basin of the upper stable node (with respect to the starting radius of curvature) but near the boundary of that basin. As the particle is carried through the spiral duct, within the cross-section it is initially attracted towards the upper saddle point within a neighbourhood of its stable manifold. By the time the particle reaches its shortest distance to the saddle point, it finds itself in the basin of the left-most stable node due to the upward motion of the saddle point (see yellow arrows). Conversely, with more turns of the spiral as shown in figure 5(b), a particle starting from the same initial position remains in the basin of attraction of the upper-most stable node throughout the spiral. Thus, this particle–fluid system exhibits rate-induced tipping (R-tipping) where, depending on the rate at which the radius of curvature changes (e.g. via manipulating $N_{turns}$), particles starting near a basin boundary can approach a stable node with a basin of attraction different from the one they started in. Similarly, rate-induced tipping occurs in a spiral duct with a tall $1\times 2$ rectangular cross-section, as shown in figure 5(c,d).

Figure 5.

Rate-induced tipping (R-tipping) for small particles of (a,b) radius $\tilde {a}=0.05$ in an in-spiral duct with $\tilde {R}_{start}=4500$, $\tilde {R}_{end}=3100$, $Re=50$ and a $2\times 1$ rectangular cross-section. Particles which start near a saddle point that separates the basins of attraction of adjacent stable node attractors can either fall in the basin of attraction of the left-most stable node for (a) $N_{turns}=1$, or one of the stable node pair to the right of this for (b) $N_{turns}=2$. Similar behaviour is observed in (c,d) an in-spiral duct with $\tilde {R}_{start}=1000$, $\tilde {R}_{end}=500$, $Re=50$ and a $1\times 2$ rectangular cross-section for larger particles of radius $\tilde {a}=0.15$ with (c) $N_{turns}=0.5$ and (d) $N_{turns}=1$. The cross-sections of the ducts show particle equilibria (coloured circles) and their motion as the radius of curvature changes through the spiral (coloured arrows) with unstable nodes in red, stable nodes (point attractors) in green and saddle points in yellow. The particle location at the end of each spiral duct (grey filled circle), its motion (grey arrows) and trajectory (black curve) are also shown. If the centre of a particle lies on the dashed rectangle it will touch at least one wall of the duct.

We note that, since the saddle-point equilibria, which separate the basins of attraction, do not move significantly between the starting and ending cross-sections of the spirals shown in figure 5, the R-tipping phenomenon only affects a small fraction of particles with initial positions near the basin boundaries. Hence, we do not expect R-tipping by itself to have major consequences for inertial particle focusing and separation. However, if R-tipping occurs in combination with bifurcations in particle equilibria, then significant changes in focusing behaviour can occur, as will be discussed in § 5.2.

5.1. Spiral ducts with bifurcations at a single radius of curvature

We first investigate the particle dynamics in spiral ducts where bifurcations take place at a single critical radius of curvature, $\tilde {R}_{cr}$. Figure 6(a) shows the focusing dynamics in an in-spiral duct with a square cross-section containing small particles ($\tilde {a}=0.05$). At the starting radius of curvature, there are two stable node attractors in the cross-section located on the horizontal centreline near the inner and outer walls. Randomly distributed particles in the cross-section initially approach these two point attractors along the slow manifold. However, at $\tilde {R}_{cr}\approx 3000$, the stable node near the outer wall takes part in a subcritical pitchfork bifurcation which ultimately makes this location unstable. All of the particles which initially migrated towards the stable node near the outer wall, which is replaced with a saddle after the bifurcation occurs, now migrate towards the stable node located near the inner wall (along heteroclinic orbits which form a slow manifold). Hence, we observe bifurcation-induced tipping (also known as B-tipping) where a focused particle cluster at an attractor tips due to a bifurcation that makes this attractor unstable. A second example of B-tipping due to saddle-node bifurcations in an in-spiral duct with square cross-section for larger particles ($\tilde {a}=0.15$) is depicted in figure 6(b). In both of these examples, since the bifurcations occur on a slow manifold, the particle migration is slow after the B-tipping takes place.

Figure 6.

Bifurcation-induced tipping (B-tipping) at a flow Reynolds number of $Re=50$. (a) An in-spiral duct with $\tilde {R}_{start}=4000$, $\tilde {R}_{end}=1600$, $N_{turns}=5$ and a square cross-section, in which particles of radius $\tilde {a}=0.05$ undergo B-tipping following a subcritical pitchfork bifurcation at $\tilde {R}_{cr}\approx 3000$; snapshots (A–E) are shown at $\tilde {R}=4000, 3500, 3170, 2000, 1600$. (b) An in-spiral duct with $\tilde {R}_{start}=1500$, $\tilde {R}_{end}=500$, $N_{turns}=3$ and a square cross-section, in which particles of radius $\tilde {a}=0.15$ undergo B-tipping following a pair of saddle-node bifurcations at $\tilde {R}_{cr}\approx 900$; snapshots (F–J) are shown at $\tilde {R}=1500, 1000, 800, 600, 500$. (c) An in-spiral duct with $\tilde {R}_{start}=120$, $\tilde {R}_{end}=50$, $N_{turns}=5$ and a rectangular $1\times 2$ cross-section, in which particles of radius $\tilde {a}=0.15$ undergo B-tipping following a saddle-node infinite period (SNIPER) bifurcation at $\tilde {R}_{cr}\approx 100$; snapshots (K–O) are shown at $\tilde {R}=120, 100, 80, 60, 50$. The cross-sectional images show the particle equilibria as coloured circles with unstable nodes in red, stable nodes (point attractor) in green, saddle points in yellow and unstable spirals in purple. Stable limit cycles are shown as black curves. Particle locations (grey circles), their motion (black arrows) and trajectories (grey curves) are also shown. If the centre of a particle lies on the dashed square/rectangle it will touch at least one wall of the duct. Red cross on spirals denotes the location of the bifurcation.

Bifurcation-induced tipping can also cause a cluster of particles which first focus to a point attractor, to ultimately migrate towards a one-dimensional limit-cycle attractor. Figure 6(c) illustrates this for particles of size $\tilde {a}=0.15$ in an in-spiral duct with $1\times 2$ rectangular cross-section. At the starting radius of curvature, the system has only one attractor, a point attractor on the horizontal centreline near the inner wall to which all particles are initially attracted. However, this stable node vanishes at $\tilde {R}_{cr}\approx 100$ and a symmetric pair of stable limit cycles emerge in what appears to be a saddle-node infinite period (SNIPER) bifurcation (Strogatz Reference Strogatz2015). Although the bifurcation takes place at a relatively short distance from the spiral inlet, the effects of the initial point attractor persist even after the bifurcation due to a saddle-node ghost (Strogatz Reference Strogatz2015), leading to slow migration and further clustering of particles near the inner wall along the horizontal centreline (see panel L in figure 6c). We note that, since the radius of curvature value at the inlet of the spirals in figure 4 is just after the SNIPER bifurcation, the dynamical origin of particle aggregation in figure 4 is also due to similar effects. After the bifurcation takes place, as shown in panel M, the pre-focused cluster of particles migrates towards the saddle point near the outer wall along its stable manifold (horizontal direction). Then, the pre-focused cluster of particles become more separated along the unstable manifold (vertical direction) of the saddle point as they move away from the saddle point (see panel N). After leaving the saddle point along its unstable manifold, the particles re-cluster within a narrow range of phases along one of the pair of vertically symmetric limit cycles, similar to the behaviour observed in figure 4.

5.2. Spiral ducts with bifurcations at multiple radii of curvature

Due to rich bifurcations taking place in particle equilibria as a function of the radius of curvature, spiral ducts can cover a range of radii of curvature where more than one bifurcation takes place. This gives rise to intricate tipping phenomena that may involve a combination of bifurcation-induced, rate-induced and/or phase-induced effects.

Figure 7 shows an example of complex tipping phenomena arising from a combination of bifurcation-induced and rate-induced effects. We consider two out-spiral ducts with matching radius of curvature at each end but a different number of turns. The particle attractors associated with the starting radius of curvature consists of a vertically symmetric pair of stable spiral attractors. These point attractors undergo a sequence of bifurcations (a saddle node and a supercritical pitchfork, as shown in the inset of figure 2a) along each spiral duct near $\tilde {R}_{cr}\approx 2850$ leading to the formation of a stable node attractor and an unstable node on the horizontal centreline of the duct cross-section; see panels C and H of figure 7. Further along the spiral, a vertically symmetric pair of saddle-node bifurcations take place off the horizontal centreline near $\tilde {R}_{cr}\approx 3050$ leading to the birth of two additional stable node attractors. For the spiral duct with more turns (figure 7a), the particles initially focus to the stable spiral attractor pair and then closely track the moving attractors. As these point attractors undergo bifurcations, the closely following particles focus to the stable node attractor on the horizontal centreline (see also Movie 3). In the spiral duct with less turns (figure 7b), the particles are unable to closely track the moving stable spiral attractor pair. The lagging particles further slow down near the location of forthcoming off-centred saddle-node bifurcations. The critical slowing prior to a saddle-node bifurcation is due to the presence of a saddle-node ghost and follows a square-root scaling law (Strogatz Reference Strogatz2015). After the saddle-node bifurcation takes place, the majority of the slowed particles are attracted towards the pair of newly formed off-centred stable node attractors (see also Movie 4). Independently, the tipping phenomena arising in each case are due to bifurcation-induced effects. However, if we compare the two spiral ducts whose only difference is the number of turns, we see a rate-induced effect where, depending on the rate at which the radius of curvature changes, the particles ultimately focus to different point attractors of the system. Thus, the tipping observed in figure 7 is a combination of bifurcation-induced and rate-induced effects.

Figure 7.

Complex tipping phenomena with both bifurcation-induced and rate-induced effects for a flow Reynolds number of $Re=65$, in an out-spiral duct with $\tilde {R}_{start}=1500$, $\tilde {R}_{end}=4000$ and a rectangular $2\times 1$ cross-section, and particles of radius $\tilde {a}=0.05$. Panels show (a) $N_{turns}=5$ (see Movie 3); all particles focus to the stable node on the horizontal centreline near the inner wall and (b) $N_{turns}=1.5$ (see Movie 4); particles focus to all three stable nodes with majority focusing near the pair of off-centred stable nodes. The difference is due to the multiple bifurcations in particle equilibria near $\tilde {R}_{cr}\approx 2850$ and $3050$ together with rate-induced effects from the different number of turns. Snapshots of the cross-section are shown at (A,F) $\tilde {R}=1500$, (B,G) $2500$, (C,H) $3000$, (D,I) $3500$ and (E,J) $4000$. The cross-sectional images show the particle equilibria as coloured circles with unstable nodes in red, stable nodes/spirals (point attractors) in green/cyan and saddle points in yellow. The particle locations (grey circles) and their motion (black arrows) are also shown. If the centre of a particle lies on the dashed rectangle it will touch at least one wall of the duct. Red crosses on spirals denote the location of the bifurcations.

Another example of a complex tipping phenomenon that involves limit-cycle attractors is shown in figure 8 for an out-spiral duct with a $1\times 2$ rectangular cross-section containing particles of size $\tilde {a}=0.15$. A vertically symmetric pair of limit cycles exists at the starting radius of curvature. As the spiral is traversed, the limit-cycle attractor pair undergoes a SNIPER bifurcation near $\tilde {R}_{cr}\approx 100$, resulting in the emergence of a stable node attractor on the horizontal centreline near the inner wall. Further along the spiral near $\tilde {R}_{cr}\approx 400$, a subcritical pitchfork bifurcation takes place on the horizontal centreline near the outer wall, resulting in the emergence of another stable node attractor. Thus, the pair of limit-cycle attractors at the inlet is replaced with two stable node point attractors at the outlet of the spiral duct. Randomly distributed particles in the $1\times 2$ cross-section of this spiral initially focus on the limit cycles. However, since the period of this limit cycle diverges as the SNIPER bifurcation is approached (due to the square-root singularity arising from the saddle-node ghost Strogatz Reference Strogatz2015), many particles accumulate on phases of the limit cycles near the inner wall while few particles are on the phases of the limit cycles located near the outer wall (see panel B of figure 8). After the SNIPER bifurcation, the accumulated particles near the inner wall focus to the newly emerged stable node attractor and more particles approach this stable point along the slow manifold. However, after the subcritical pitchfork bifurcation near the outer wall, a new stable node attractor emerges there too. Hence, some of the particles that were initially on phases of the limit cycle near the outer wall, fall in the basin of attraction of the stable node attractor near the outer wall and focus to this point attractor. Hence, we observe a complex tipping behaviour where in addition to the bifurcations, the phase of particles on each limit cycle determines the stable node attractor they will focus towards.

Figure 8.

Complex tipping phenomena with bifurcation-induced and phase-induced effects for a flow Reynolds number of $Re=25$ in an out-spiral duct with $\tilde {R}_{start}=50$, $\tilde {R}_{end}=750$, a rectangular $1\times 2$ cross-section, $N_{turns}=3$ and particles of radius $\tilde {a}=0.15$. These are governed by bifurcations of particle equilibria at $\tilde {R}_{cr}\approx 100$ and $400$ as well as the phase of the particles on the limit-cycle attractors. Snapshots of the cross-section are shown at (A) $\tilde {R}=50$, (B) $100$, (C) $150$, (D) $500$ and (E) $750$. The cross-sectional images show the particle equilibria as coloured circles with unstable nodes in red, stable nodes (point attractors) in green, saddle points in yellow and unstable spirals in purple. Stable limit cycles are shown as black curves. The particle locations (grey circles), their motion (black arrows) and trajectories (grey curves) are also shown. If the centre of a particle lies on the dashed rectangle it will touch at least one wall of the duct. Red crosses on spiral denote the location of the bifurcations.