Inertial migration of a neutrally buoyant spheroid in plane Poiseuille flow

TL;DR

This study analyzes the inertial migration behavior of neutrally buoyant spheroids in plane Poiseuille flow, revealing how shape and inertial effects influence equilibrium positions and lift velocities at low Reynolds numbers.

Contribution

It provides a Jeffery-averaged analysis of inertial migration for spheroids with arbitrary aspect ratios, extending classical sphere results to more complex particle shapes.

Findings

Equilibrium positions are similar to classical sphere results, starting near 0.6H/2 from the centerline.

Inertial effects stabilize specific Jeffery orbits depending on aspect ratio.

Lift profiles are shape-independent and match those of spheres at low Re_c.

Abstract

We study the cross-stream inertial migration of a torque-free neutrally buoyant spheroid, of an arbitrary aspect ratio $\kappa$, in wall-bounded plane Poiseuille flow for small particle Reynolds numbers\,($Re_p\ll1$) and confinement ratios\,($\lambda\ll1$), with the channel Reynolds number, $Re_c = Re_p/\lambda^2$, assumed to be arbitrary; here, $\lambda=L/H$ where $L$ is the semi-major axis of the spheroid and $H$ denotes the separation between the channel walls. In the Stokes limit\,($Re_p =0)$ and for $\lambda \ll 1$, a spheroid rotates along any of an infinite number of Jeffery orbits parameterized by an orbit constant $C$, while translating with a time dependent speed along a given ambient streamline. Weak inertial effects stabilize either the spinning\,($C=0$) or the tumbling orbit\,($C=\infty$), or both, depending on $\kappa$. The separation of the Jeffery-rotation and orbital drift time scales, from that associated with cross-stream migration, implies that the latter occurs due to a Jeffery-averaged lift velocity. Although the magnitude of this averaged lift velocity depends on $\kappa$ and $C$, the shape of the lift profiles are identical to those for a sphere, regardless of $Re_c$. In particular, the equilibrium positions for a spheroid remain identical to the classical Segre-Silberberg ones for a sphere, starting off at a distance of about $0.6(H/2)$ from the channel centerline for small $Re_c$, and migrating wallward with increasing $Re_c$. For spheroids with $\kappa \sim O(1)$, the Jeffery-averaged analysis is valid for $Re_p\ll1$; for extreme aspect ratio spheroids, the regime of validity becomes more restrictive being given by $Re_p\,\kappa/\ln \kappa \ll 1$ and $Re_p/\kappa^2 \ll 1$ for $\kappa \rightarrow \infty$\,(slender fibers) and $\kappa \rightarrow 0$\,(flat disks), respectively.

Figures

Equations

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

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

$$\mathbf{u}$$

$${\mathbf{u}}^{\infty }=(\alpha +\beta r_2+\gamma {r_2}^2)1_1-{\mathbf{U}}_p,$$

$$\nabla \cdot {\mathbf{\sigma }}^{\prime }={\nabla }^2{\mathbf{u}}^{\prime }-\nabla p^{\prime }$$

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

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

$$\nabla \cdot {\mathbf{\sigma }}^t$$

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

$${\mathbf{u}}^t$$

$$(\nabla \cdot {\mathbf{\sigma }}^{\prime })\cdot {\mathbf{u}}^t-(\nabla \cdot {\mathbf{\sigma }}^t)\cdot {\mathbf{u}}^{\prime }$$

$${\int }_{V^F}\nabla \cdot ({\mathbf{\sigma }}^{\prime }\cdot {\mathbf{u}}^t-{\mathbf{\sigma }}^t\cdot {\mathbf{u}}^{\prime })dV+{\int }_{V^F}({\mathbf{\sigma }}^t\cdot {\nabla \mathbf{u}}^{\prime }-{\mathbf{\sigma }}^{\prime }\cdot {\nabla \mathbf{u}}^t)dV$$

$$-{\int }_{S_p+S_w+S_{\infty }}({\mathbf{\sigma }}^{\prime }\cdot {\mathbf{u}}^t-{\mathbf{\sigma }}^t\cdot {\mathbf{u}}^{\prime })\cdot \mathbf{n}dS$$

$$-{\int }_{S_p}({\mathbf{\sigma }}^{\prime }\cdot {\mathbf{u}}^t-{\mathbf{\sigma }}^t\cdot {\mathbf{u}}^{\prime })\cdot \mathbf{n}dS$$

$$-{\mathbf{U}}_p^t\cdot {\int }_{S_p}{\mathbf{\sigma }}^{\prime }\cdot \mathbf{n}dS-{\mathbf{\Omega }}_p^t\cdot {\int }_{S_p}\mathbf{r}\wedge ({\mathbf{\sigma }}^{\prime }\cdot \mathbf{n})dS+{\mathbf{\Omega }}_p\cdot {\int }_{S_p}\mathbf{r}\wedge ({\mathbf{\sigma }}^t\cdot \mathbf{n})dS$$

$$+{\mathbf{U}}_p\cdot {\int }_{S_p}{\mathbf{\sigma }}^t\cdot \mathbf{n}dS-{\int }_{S_p}({\mathbf{\sigma }}^t\cdot \mathbf{n})\cdot (\alpha +\beta r_2+\gamma r_2^2)1_{1}dS=R{e}_p{\int }_{V^F}{\mathbf{u}}^t\cdot \mathbf{f}dV.$$

$$-R{e}_p{\mathbf{U}}_p^t\cdot \frac{d{\mathbf{U}}_p}{dt}-R{e}_p{\mathbf{\Omega }}_p^t\cdot \frac{d({\mathbf{I}}_p\cdot {\mathbf{\Omega }}_p)}{dt}+{\mathbf{U}}_p\cdot {\int }_{S_p}{\mathbf{\sigma }}^t\cdot \mathbf{n}dS$$

$$-{\int }_{S_p}({\mathbf{\sigma }}^t\cdot \mathbf{n})\cdot (\alpha +\beta r_2+\gamma r_2^2)1_{1}dS=R{e}_p{\int }_{V^F}{\mathbf{u}}^t\cdot {\mathbf{f}}_{s}dV.$$

$${\mathbf{f}}_s=\left(\frac{\partial {\mathbf{u}}_s}{\partial t}+{\mathbf{u}}_s\cdot {\nabla \mathbf{u}}_s+{\mathbf{u}}_s\cdot {\nabla \mathbf{u}}^{\infty }+{\mathbf{u}}^{\infty }\cdot {\nabla \mathbf{u}}_s\right).$$

$${\mathbf{U}}_{p0}\cdot {\int }_{S_p}{\mathbf{\sigma }}^t\cdot \mathbf{n}dS={\int }_{S_p}({\mathbf{\sigma }}^t\cdot \mathbf{n})\cdot (\alpha +\beta r_2+\gamma r_2^2)1_{1}dS.$$

$$-{\mathbf{U}}_p^t\cdot \frac{d{\mathbf{U}}_{p0}}{dt}-{\mathbf{\Omega }}_p^t\cdot \frac{d({\mathbf{I}}_p\cdot {\mathbf{\Omega }}_{p0})}{dt}+{\mathbf{U}}_{p1}\cdot {\int }_{S_p}{\mathbf{\sigma }}^t\cdot \mathbf{n}dS={\int }_{V^F}{\mathbf{u}}^t\cdot {\mathbf{f}}_{s}dV.$$

$${\mathbf{\sigma }}^t\cdot \mathbf{n}=$$

$${\mathbf{U}}_{p0}=\left[\alpha +\frac{\gamma }{3{\kappa }^2}[{\cos }^2{\varphi }_j+({\cos }^2{\theta }_j+{\kappa }^2{\sin }^2{\theta }_j){\sin }^2{\varphi }_j]\right]1_1,$$

$${\mathbf{\Omega }}_{p0}=-\frac{\beta ({\kappa }^2-1)}{4({\kappa }^2+1)}\cos {\varphi }_j\sin 2{\theta }_j{1}_1+\frac{\beta ({\kappa }^2-1)}{4({\kappa }^2+1)}\sin {\varphi }_j\sin 2{\theta }_j{1}_2$$

$$\displaystyle+\frac{\beta}{2}\big{[}-1+\frac{(\kappa^{2}-1)}{(\kappa^{2}+1)}\cos 2\phi_{j}\sin^{2}\theta_{j}\big{]}\bm{1}_{3}.$$

$${\mathbf{u}}^{\infty }=\left[(\beta r_2+\gamma r_2^2)-\frac{\gamma }{3{\kappa }^2}[{\cos }^2{\varphi }_j+({\cos }^2{\theta }_j+{\kappa }^2{\sin }^2{\theta }_j){\sin }^2{\varphi }_j]\right]1_1.$$

$$V_p=R{e}_p{\mathbf{U}}_{p1}\cdot 1_2$$

$$-R{e}_p{\mathbf{U}}_p^{t2}\cdot \frac{d{\mathbf{U}}_{p0}}{dt}-R{e}_p{\mathbf{\Omega }}_p^{t2}\cdot \frac{d({\mathbf{I}}_p\cdot {\mathbf{\Omega }}_{p0})}{dt}$$

$$⟨V_p⟩=-R{e}_p{\int }_{V^F}⟨{\mathbf{u}}^{t2}\cdot \left(\frac{\partial {\mathbf{u}}_s}{\partial t}+{\mathbf{u}}_s\cdot {\nabla \mathbf{u}}_s+{\mathbf{u}}_s\cdot {\nabla \mathbf{u}}^{\infty }+{\mathbf{u}}^{\infty }\cdot \nabla {\mathbf{u}}_s\right)⟩dV,$$

$$\frac{d{\varphi }_j}{dt}=$$

$$\frac{d{\theta }_j}{dt}=$$

Full text

1 Introduction

In an experimental study of pressure-driven pipe flow of a suspension of neutrally buoyant spheres, Segre & Silberberg (1962a, b) reported the migration of particles to an intermediate annular location. This was termed the tubular pinch effect to indicate that the initially uniform distribution of particles over the pipe cross-section is ‘pinched’ to a narrow annulus with increasing downstream distance. The authors performed experiments for $0.03\leq d_{p}/D\leq 0.15$ and for $Re$ upto 700, where $d_{p}$ and $D$ are the particle and pipe diameters, respectively, and $Re$ is the Reynolds number based on $D$ and the mean velocity. The aforementioned pinch occurred at $0.6\times$pipe radius for the smaller $Re$’s, and shifted towards the walls with increasing $Re$. In any unidirectional shearing flow, including pressure-driven flow in the pipe and channel geometries, reversibility of the Stokes equations prohibits cross-stream migration of a neutrally buoyant sphere, and the pinch effect above is a consequence of migration driven by inertial lift forces.

Since the original experiments of Segre and Silberberg, inertial migration in varying flow configurations has been investigated in a number of theoretical (Saffman (1965), Cox & Brenner (1968), Ho & Leal (1974), Vasseur & Cox (1976), Cox & Hsu (1977), Schonberg & Hinch (1989), Asmolov (1999), Matas et al. (2009)), numerical (Chun & Ladd (2006), Shao et al. (2008), Morita et al. (2017), Nakayama et al. (2019), Pan et al. (2021)) and experimental (Repetti & Leonard (1964), Jeffrey & Pearson (1965), Karnis et al. (1966), Tachibana (1973), Aoki et al. (1979), Matas et al. (2004), Masaeli et al. (2012), Nakayama et al. (2019)) studies. From the rheological perspective, the pinch effect in the original pipe-flow experiments is undesirable, since it manifests as an apparent non-Newtonian behaviour for a dilute suspension of spheres. For instance, in a capillary viscometer, depending on the suspension flow rate, the residence time of the particles may or may not be sufficient for inertial migration to be complete; incomplete migration leads to a non-linear dependence of the inferred viscosity on the shear rate(Segre & Silberberg, 1963). On the other hand, inertial migration has recently been used as a tool in microfluidics to separate particles based on size (Di Carlo et al., 2007). Due to the robust fault-tolerant physical effects employed, and high rates of operation (the latter being a natural consequence of being in the inertial regime), inertial microfluidic systems are expected to have a broad range of applications in continuous bioparticle separation, cell and particle manipulation, and filtration systems(Di Carlo, 2009).

A first theoretical explanation of the phenomenon was given by Ho & Leal (1974) who determined the lift force on a sphere in plane Poiseuille flow for $Re_{p},Re_{c}\ll 1$, where $Re_{p}$ and $Re_{c}$ are the particle and channel Reynolds numbers, respectively. A more accurate calculation was performed in Vasseur & Cox (1976) using the framework developed earlier in Cox & Brenner (1968), wherein the lift velocity was expressed in terms of a volume integral involving the Green’s function for creeping flow in the presence of a pair of plane boundaries (the channel walls). The resulting lift force profile had a pair of zero crossings, symmetrically located on either side of the channel centerline, and that corresponded to stable equilibrium locations. These were interpreted as the analog of the intermediate annulus observed in the experiments, implicitly pointing to the similarity of the physics governing inertial migration in the channel and pipe geometries. Later, Schonberg & Hinch (1989) calculated the lift velocity of a sphere in plane Poiseuille flow for $Re_{p}\ll 1$ and for $Re_{c}$ upto $150$, finding the equilibrium locations to move towards the respective walls with increasing $Re_{c}$, consistent with the original observations. Asmolov (1999) confirmed these findings, and further extended the calculation to $Re_{c}=3000$111This upper bound was regarded as reasonable based on the threshold ($Re_{c}\sim 11544$) for the Tolmein-Schlichting instability. It is now known, however, that the actual transition of plane Poiseuille flow to turbulence has a subcritical character, occurring at much lower $Re_{c}$’s of $O(2000)$..

There are numerous instances in microfluidic and other settings where the particles of interest are anisotropic, for example, cancer cells (Suresh, 2007), blood cells (Toner & Irimia, 2005) or polymeric microstructures (Chung et al., 2008). For anisotropic particles, the inertial lift force, and any equilibria that arise as a result of this force vanishing at specific locations, are expected to depend on particle shape. Motivated by this, Masaeli et al. (2012) conducted experiments on spheroids of aspect ratio $1\leq\kappa\leq 5$, suspended in pressure-driven flow through a rectangular duct, for $Re_{c}$ upto $80$ ($Re_{c}$ defined based on the smaller cross-sectional dimension). The study confirmed the existence of shape-sensitive equilibria, with large-aspect-ratio spheroids migrating to locations near the channel centerline and those with order unity aspect ratios migrating towards the duct walls, over a range of cross-sectional aspect ratios.

Motivated by the experiments of Masaeli et al. (2012), we take a first step towards analyzing the inertial migration of a freely rotating neutrally buoyant spheroid in plane Poiseuille flow. Specifically, we calculate the leading order time-averaged lift velocity for $Re_{p}\ll 1$ within the framework of a point-particle approximation; $Re_{c}$, while much larger than $Re_{p}$, is otherwise arbitrary. $\S$2 below presents the governing equations and boundary conditions in the context of the problem formulation. Next, in $\S$3, we examine the small-$Re_{c}$ limit where the time-averaged lift velocity is determined semi-analytically using a generalized reciprocal theorem formulation originally used by Ho & Leal (1974), and that is derived in $\S$3.1. Scaling arguments given in $\S$3.2 show that the dominant contribution to the lift velocity in this limit comes from scales of order the channel width $H$ which is much smaller than the inertial screening length of $O(HRe_{c}^{-1/2})$. Inertia therefore has a regular character, with the lift velocity being $O(Re_{p})$, and its calculation requiring knowledge of only the Stokesian disturbance fields in the confined domain. These Stokesian fields are calculated in $\S$3.3, and are then used to obtain the time-averaged lift velocity for spheroids in $\S$3.4. This is followed by a presentation of the results in $\S$3.5. In $\S$4, the time-averaged lift velocity is calculated numerically for $Re_{c}\gtrsim O(1)$, with inertia now acting as a singular perturbation. Following Schonberg & Hinch (1989), this calculation involves a partial Fourier transform of the linearized Navier-Stokes equations, leading to coupled ODEs (in the transverse coordinate) for the transformed pressure and normal velocity fields, and their numerical solution using a shooting method. A brief outline of the aforementioned calculation procedure is given in $\S$4.1, with the results obtained discussed in $\S$4.2. In both $\S$3 and $\S$4, a time averaging is necessary on account of the separation between spheroid rotation and inertial migration time scales, and at leading order in $Re_{p}$, is based on the Jeffery angular velocity. This, together with the fact that the dominant contributions to the lift velocity come from scales much larger than the particle size for all $Re_{c}$, lead to the equilibrium locations being independent of the spheroid aspect ratio, and identical to those for spheres; the magnitude of the lift force does depend on aspect ratio. In $\S$5, we summarize the main results, and briefly discuss their implications for shape-sorting in microfluidic settings.

2 Problem Formulation

Figure 1a shows a neutrally buoyant spheroid of aspect ratio $\kappa=L/b$ ($L$ and $b$ are the semi-major and minor axes) freely suspended in a wall-bounded plane Poiseuille flow at a distance $d$ from the lower wall; $\kappa<1$ and $>1$ for oblate and prolate spheroids, respectively. The non-dimensional equations governing the velocity field $\bm{u}$ and pressure field $p$ are given by:

[TABLE]

where $\bm{u}$ satisfies the following boundary conditions:

[TABLE]

In (1a-2c), all variables are non-dimensionalized using $L$ and the velocity scale $V_{c}=V_{\text{max}}L/H$, this being the order of the velocity change across the ends of the spheroid. The particle Reynolds number $Re_{p}=V_{c}L/\nu=V_{\text{max}}L^{2}/H\nu$, $\nu$ being the kinematic viscosity of the fluid. $\bm{\Omega}_{p}$ in (2a) is the angular velocity of the spheroid whose surface is denoted by $S_{p}$. The coordinate system chosen in writing the above equations translates with the spheroid velocity $\bm{U}_{p}$, with its origin at the spheroid center. Thus, $\bm{u}^{\infty}$ in (2b), the ambient plane Poiseuille flow in this frame, is given by,

[TABLE]

where $\alpha\bm{1}_{1}-\bm{U}_{p}$, with $\alpha=4\lambda^{-1}s(1-s)$, is the ambient slip velocity at the spheroid center, $\beta=4(1-2s)$ is the local shear rate that varies linearly across the channel, and $\gamma=-4\lambda$ is the constant curvature of the plane Poiseuille flow; $s=d/H$ here being the (non-dimensional) spheroid location. $\lambda=L/H$ is the confinement ratio assumed small, so the channel Reynolds number $Re_{c}=Re_{p}/\lambda^{2}\gg Re_{p}$. As will be argued below, at leading order in $Re_{p}$ and $\lambda$, $\bm{U}_{p}$ is along the flow direction.

We now define the disturbance fields $\bm{u}^{\prime}=\bm{u}-\bm{u}^{\infty}$, $p^{\prime}=p-p^{\infty}$ which satisfy the following governing equations:,

[TABLE]

where $\bm{u}^{\prime}$ satisfies:

[TABLE]

Unlike a sphere, one needs to include the unsteady terms in the inertial acceleration on account of spheroid rotation even in the Stokes limit ($Re_{p}=0$). We analyze the inertial migration problem defined above in the limit $\lambda,Re_{p}\ll 1$ with $Re_{c}=Re_{p}/\lambda^{2}$ arbitrary. It is also assumed that $s,\,(1-s)\gg\lambda$, implying that the analysis is restricted to the spheroid being at distances from either wall that are much larger than $O(L)$.

Before embarking on a detailed analysis for weak inertia, it is worth summarizing the nature of spheroid motion in plane Poiseuille flow in the Stokes limit. Since the local linear flow approximation for the plane Poiseuille profile is simple shear flow, the neutrally buoyant spheroid must rotate along Jeffery orbits(Jeffery, 1922), this rotation being characterized by the polar ($\theta_{j}$) and azimuthal ($\phi_{j}$) angles of the spheroid symmetry axis (see Figure 1b) being functions of time; these equations are given later in $\S$3 (see 24a,b). Jeffery rotation is best described in $(C,\tau$) coordinates, where $C\in[0,\infty)$ is the orbit constant, and $\tau$ is the orbit phase that changes at a constant albeit $\kappa$-dependent rate from [math] to $2\pi$ over a single period(Leal & Hinch, 1971; Dabade et al., 2016); the two limiting orbits correspond to the spinning ($C=0$) and tumbling ($C=\infty$) modes. The Jeffery period is independent of $C$, being given by $T_{\text{jeff}}=2\pi(\beta V_{\text{max}}/H)^{-1}(\kappa+\kappa^{-1})$, where $\beta$ accounts for the linearly varying shear rate. The curvature of the plane Poiseuille profile only affects the spheroid translation velocity that is now a function of its orientation, and thence, of time. A neutrally buoyant spheroid in an unbounded plane Poiseuille flow continues to move along a given ambient streamline, with a speed that ranges from a maximum, corresponding to the flow-aligned orientation ($\phi_{j}=0$ and $\frac{\pi}{2}$ for $\kappa>1$ and $<1$), to a minimum when the spheroid is aligned orthogonal to the flow direction ($\phi_{j}=\frac{\pi}{2}$ and [math] for $\kappa>1$ and $<1$)(Chwang, 1975). For very large (small) $\kappa$, the spheroid spends an increasing fraction of a Jeffery period in the flow (gradient-vorticity plane)-aligned orientation, and the translational motion thereby acquires an increasingly jerky character owing to the spheroid abruptly slowing down during the brief periods of mis-alignment.

In presence of walls, the leading order correction to the motion above arises from interaction of the spheroid with time dependent image stresslets induced by each wall, resulting in an $O(\lambda^{2})$ lateral velocity component even in the Stokes limit, with an additional $O(\lambda^{3})$ correction to the Jeffery angular velocity. Note that this image-stresslet interaction does not lead to transverse motion for a sphere, as may be seen from the stresslet orientation (along the local extensional axis) and the associated purely radial velocity field; as already mentioned in $\S$1, this must be so, independent of $\lambda$, due to reversibility constraints. Thus, for $\lambda\ll 1$, the center of mass of a neutrally buoyant spheroid in wall-bounded plane Poiseuille flow exhibits a small-amplitude oscillatory motion about an ambient streamline, the amplitude being $O(\lambda^{2})$. In other words, unlike a sphere, Stokesian reversibility does not preclude an instantaneous lift force for an anisotropic particle. Evidence for such oscillation-cum-tumbling spheroid trajectories, for finite $\lambda$, is available from earlier computations, with there being a transition from tumbling (rotation) to angular oscillations beyond a threshold $\lambda$ close to unity, corresponding to sufficiently narrow channels,(Sugihara-Seki, 1993, 1996; Staben et al., 2003, 2006). However, Stokesian reversibility still forbids a net migration of the spheroid in the transverse direction, and to allow for such a motion, one needs inertia. It is this net cross-stream migration, for small $Re_{p}$ and $\lambda$, that is analyzed in the following sections. We obtain the time-averaged motion of a spheroid in this limit, and since the wall and inertia-induced modification of the primary translational and rotational motion are weak, the time average corresponds to an average over a Jeffery period.

It is well-known that the Stokes equations do not provide a uniformly valid leading order approximation for small but finite $Re_{p}$, and that in general one requires a matched asymptotic expansion approach (Proudman & Pearson, 1957) involving an inner expansion in the neighbourhood of the particle, and an outer expansion at distances of order an inertial screening length, to calculate inertial corrections. The screening length for the present problem is $LRe_{p}^{-1/2}$, or equivalently, $HRe_{c}^{-1/2}$. For small $Re_{c}$, the screening length is larger than $H$, so the channel walls lie in the inner Stokesian region where fluid inertia may be treated as a regular perturbation. The inertial migration problem for a sphere in this limit was investigated by Ho & Leal (1974), and in a series of papers by Brenner, Cox and collaborators (Cox & Brenner (1968), Vasseur & Cox (1976), Cox & Hsu (1977)). We examine the analogous problem for a spheroid in this limit in $\S$3 below. When $Re_{c}\gtrsim O(1)$, the inertial screening length is of $O(H)$ or smaller, and the solution procedure involves the outer expansion, the leading order term of which satisfies the linearized Navier-Stokes equations. This limit was first examined for a sphere in Schonberg & Hinch (1989), and we examine the same for a spheroid in $\S$4.

3 The inertial lift velocity for $Re_{c}\ll 1$

Although inertial effects are a regular perturbation for small $Re_{c}$, following Ho & Leal (1974), we don’t calculate this perturbation explicitly, and instead use the generalized reciprocal theorem that relates the velocity and stress fields of the problem of interest with the velocity and stress fields of a simpler test problem whose solution is known. The problem of interest ($\bm{u}^{\prime},\bm{\sigma}^{\prime}$) corresponds to the motion of a torque-free neutrally buoyant spheroid in a wall-bounded plane Poiseuille flow, taking into account the inertial acceleration of the suspending fluid. Since the quantity of interest is the inertial lift velocity, the test problem ($\bm{u}^{t},\bm{\sigma}^{t}$) is taken to be a torque-free spheroid, in a quiescent ambient between parallel walls, acted on by an arbitrarily oriented unit force; the test spheroid has the same instantaneous orientation as the one in the actual problem. The governing equations and boundary conditions for the actual problem have already been given in (4a,b) and (5a-c). Those for the test problem are:

[TABLE]

with the boundary conditions:

[TABLE]

where $\bm{U}_{p}^{t}$ and $\bm{\Omega}_{p}^{t}$ are the spheroid translational and angular velocities.

3.1 The Generalized Reciprocal theorem

To derive the generalized reciprocal theorem identity, we contract (4a) with $\bm{u}^{t}$ and (6a) with $\bm{u}^{\prime}$, and subtract the resulting expressions(Ho & Leal, 1974) to obtain:

[TABLE]

where $\bm{f}=\left(\frac{\partial\bm{u}^{\prime}}{\partial t}+\bm{u}^{\prime}\cdot\bm{\nabla u}^{\prime}+\bm{u}^{\prime}\cdot\bm{\nabla u}^{\infty}+\bm{u}^{\infty}\cdot\bm{\nabla u}^{\prime}\right)$. Integrating (8) over the fluid volume ($V^{F}$) between the channel walls:

[TABLE]

Using (4b) and (6b), along with the definitions of $\bm{\sigma}^{\prime}$ and $\bm{\sigma}^{t}$, the second integral on the LHS in (9) can be shown to be identically zero. Applying the divergence theorem to the first integral yields,

[TABLE]

where $\bm{n}$ is the unit normal at all bounding surfaces pointing into the fluid domain. The set of bounding surfaces include the particle surface ($S_{p}$), the channel walls ($S_{w}$) and the surface at infinity ($S_{\infty}$); the latter can be thought of as the curved surface of a cylinder of radius $R$, with its axis along the gradient direction, in the limit $R\to\infty$. For $R\gg H$, the disturbance fields in the actual and test problems are exponentially small (see (96) and (37)), and therefore, the integral over $S_{\infty}$ in (10) is vanishingly small. The integral $\int_{S_{w}}(\bm{\sigma}^{\prime}\cdot\bm{u}^{t}-\bm{\sigma}^{t}\cdot\bm{u}^{\prime})\cdot\bm{n}\,dS\,\,$ can be shown to be identically zero on using the no-slip conditions (5c) and (7c) on the channel walls. Thus, (10) reduces to:

[TABLE]

Applying the no-slip conditions (5a) and (7a) to (11),

[TABLE]

The hydrodynamic force and torque experienced by the particle may be written as $\int_{S_{p}}\bm{\sigma}^{\prime}\cdot\bm{n}\,dS=Re_{p}\,\,d\bm{U}_{p}/dt$ and $\int_{S_{p}}\bm{r}\wedge(\bm{\sigma}^{\prime}\cdot\bm{n})\,dS=Re_{p}\,\,d(\bm{I}_{p}\cdot\bm{\Omega}_{p})/dt$, $\bm{I}_{p}$ being the spheroid moment of inertia tensor. Next, to $O(Re_{p})$, one may replace $\bm{u}^{\prime}$ by $\bm{u}_{s}$ in the volume integral in (12), $\bm{u}_{s}$ being the corresponding Stokesian approximation. The validity of such a replacement requires the resulting volume integral to be convergent, this being related to inertia acting as a regular perturbation. Detailed arguments in this regard are given in the next subsection. Further, noting that $\int_{S_{p}}\bm{r}\wedge(\bm{\sigma}^{t}\cdot\bm{n})\,dS=0$ owing to the spheroid being torque-free in the test problem, one obtains:

[TABLE]

where $\bm{f}_{s}$ denotes the approximation of $\bm{f}$ based on replacing $\bm{u}^{\prime}$ by $\bm{u}_{s}$, being given by:

[TABLE]

One may now use the small-$Re_{p}$ expansions $\bm{U}_{p}=\bm{U}_{p0}+Re_{p}\bm{U}_{p1}+O(Re_{p}^{2})$ and $\bm{\Omega}_{p}=\bm{\Omega}_{p0}+Re_{p}\bm{\Omega}_{p1}+O(Re_{p}^{2})$, of the spheroid translational and angular velocities, in (13), and obtain the following relations at successive orders in $Re_{p}$:

• •

$O(1)$:

[TABLE]

• •

$O(Re_{p})$:

[TABLE]

The $O(1)$ problem:

For the $O(1)$ problem defined by (15), we choose $\int_{S_{p}}\bm{\sigma}^{t}\bm{\cdot n}\,dS=\bm{1}_{1}$, corresponding to the test spheroid translating due to a unit flow-aligned force. In the absence of boundaries, the induced velocity field is known in terms of vector spheroidal harmonics(Dabade et al., 2015), and the surface force density to be used on the RHS is given by(Kushch & Sangani, 2003):

[TABLE]

where $\xi>\xi_{0}>1$ is the coordinate characterizing a family of confocal spheroids with $\xi_{0}$ representing the spheroid surface, $\bm{1}_{\xi}$ is the unit normal to the spheroid surface and $p^{t}$ is the pressure field. Evaluating the surface integral in (15), using (17), one obtains:

[TABLE]

which is Faxen’s law for the translation of a force-free spheroid truncated to second order, the term proportional to $\gamma$ arising from the curvature of the ambient flow; higher order terms in the expansion are zero for an ambient quadratic flow(Happel & Brenner, 2012). (18) pertains to spheroid of a given orientation, and a time trajectory requires knowing $\theta_{j}$ and $\phi_{j}$ as functions of time. For this purpose, one may use Faxen’s law relating the torque and angular velocity of a spheroid. This is again an infinite expansion in the general case, but with only the leading order term, proportional to the ambient velocity gradient, surviving for a quadratic flow. Thus, spheroid rotation remains identical to that in an ambient linear flow, and application of the torque-free constraint yields:

[TABLE]

The equations governing Jeffery orbits arise as components of (19) in the body-aligned coordinate system, and are given below in (23ab). With $\theta_{j}$ and $\phi_{j}$ defined in this manner, $\bm{U}_{p0}$ as defined in (18) exhibits the time dependence described earlier in $\S$2. Substituting $\bm{U}_{p0}$ in (3), the ambient flow in the reference frame chosen for the reciprocal theorem is given by:

[TABLE]

Note that the $O(\lambda^{2})$ center-of-mass oscillations mentioned in $\S$2 can only arise from $\bm{U}_{p0}$ having a component along the gradient direction, and this requires an expression for $\bm{\sigma}^{t}\cdot\bm{n}$ that includes the influence of the plane boundaries.

The $O(Re_{p})$ problem:

At $O(Re_{p})$, one chooses the test force as $\int_{S_{p}}\bm{\sigma}^{t}\bm{\cdot n}\,dS=-\bm{1}_{2}$, which leads to:

[TABLE]

for the instantaneous lift velocity of a neutrally buoyant spheroid in plane Poiseuille flow, which is a function of the changing spheroid orientation via $\bm{u}^{t2}$ and $\bm{u}_{s}$, both of them depend on $\bm{p}$. The additional superscript ‘$2$’ for the test problem quantities indicate the choice of the test force orientation above.

As already mentioned, for small confinement ratios, a spheroid in wall-bounded plane Poiseuille flow rotates along Jeffery orbits (Jeffery, 1922) in the Stokes limit. For small but finite $Re_{p}$, there is an additional $O(Re_{p})$ drift across orbits which stabilizes the tumbling mode ($C=\infty$) for prolate spheroids, the spinning mode ($C=0$) for oblate spheroids with $0.14<\kappa<1$ and, depending on the initial orientation, either the spinning or tumbling mode for oblate spheroids with $\kappa<0.14$ (Einarsson et al., 2015; Dabade et al., 2016). The time scale for rotation along a Jeffery orbit is $T_{\text{jeff}}\sim O(H/V_{\text{max}})$ for $\kappa\sim O(1)$. The inertia-driven orbital drift above occurs on longer time scales $t_{\text{drift}}\sim O(Re_{p}^{-1}H/V_{\text{max}})$. Scaling arguments in $\S$3.2 below show that the volume integral in (21) is the dominant contribution to the cross-stream migration, being $O(1)$ for $\lambda\ll 1$, so the time scale for inertial migration is $t_{\text{lift}}\sim O(Re_{p}^{-1}\lambda^{-1}H/V_{\text{max}})$. Since $t_{\text{lift}}/t_{\text{drift}}\sim\lambda^{-1}\gg 1$, for purposes of the migration calculation, one may assume that the spheroid has settled into its stable Jeffery orbit. Moreover, since $t_{\text{lift}}/T_{\text{jeff}}\sim Re_{p}^{-1}\lambda^{-1}\gg 1$, the leading order migration is the result of an orientation-averaged lift velocity, the average being over the orientations sampled in the stabilized Jeffery orbit. As indicated below, this average may be approximated based on the Jeffery angular velocity for small $Re_{p}$. Superposed on this orientation-averaged migration trajectory would be small-amplitude oscillations of $O(Re_{p}\lambda)$, corresponding to the fluctuations in the instantaneous lift velocity arising from the rapidly changing spheroid orientation, as may be established formally using the method of multiple scales.

Owing to the aforementioned time scale separation in the limit $\lambda,Re_{p}\ll 1$, it is of interest to determine the Jeffery-averaged rather than the instantaneous lift velocity. Averaging both sides of (21) over a Jeffery period, one obtains:

[TABLE]

where the averaging operation is defined as $\langle.\rangle=\frac{1}{T_{\text{jeff}}}\int_{0}^{T_{\text{jeff}}}(.)\,dt$, and pertains to the inertially stabilized Jeffery orbit (that is, $C$ is either [math] or $\infty$) for times longer than $O(Re_{p}^{-1}H/V_{\text{max}})$.

In writing down (22), we have neglected the terms on the RHS of (21) involving the translational and angular accelerations in the actual problem in the Stokesian limit. These depend on $\lambda$ at leading order, and are therefore asymptotically smaller than the Jeffery-averaged volume integral that, as mentioned above, is independent of $\lambda$ for $\lambda\ll 1$. For a spinning spheroid both acceleration terms are trivially zero because, independent of $\lambda$, a spinning spheroid (like a sphere) translates with a constant speed along an ambient streamline, while rotating at a uniform rate about the ambient vorticity direction. For a spheroid rotating in any other Jeffery orbit, both terms do lead to non-zero Jeffery-averaged contributions. The first of these terms is non-zero because of a correlation between the time-periodic variations of $\frac{d\bm{U}_{p0}}{dt}$ along the flow direction, and the analogous variation of the flow-directed component of $\bm{U}_{p}^{t2}$, the latter arising due to the periodically varying test spheroid orientation for a fixed force (along $-\bm{1}_{2}$). The acceleration in the actual problem is only $O(\lambda^{2})$, however, as may be seen from the Faxen’s translation law given earlier. The smallness of the second term is because $\bm{\Omega}_{p}^{t2}$, at leading order, is driven by the velocity gradient associated with the equal and oppositely directed image-Stokeslet, and is again $O(\lambda^{2})$; the angular acceleration in the actual problem, $\frac{d(\bm{I}_{p}\cdot\bm{\Omega}_{p0})}{dt}$, is $O(1)$. Thus, both acceleration terms in (21) are $O(\lambda^{2})$ smaller than the volume integral.

The averaging operation in (22) corresponds to a fixed $C$, and is naturally evaluated out in $(C,\tau)$ coordinates(Leal & Hinch, 1971; Dabade et al., 2016). Here, $C=\frac{\tan\theta_{j}(\kappa^{2}\sin^{2}\phi_{j}+\cos^{2}\phi_{j})^{1/2}}{\kappa}$ and $\tau=\tan^{-1}\left[\frac{1}{\kappa\tan\phi_{j}}\right]$, where $\theta_{j}$ and $\phi_{j}$ are the polar and azimuthal angles characterizing the spheroid orientation. The latter is denoted by the unit vector $\bm{p}$ in Figure 1b with $\bm{p}=\sin\theta_{j}\cos\phi_{j}\bm{1}_{1}+\sin\theta_{j}\sin\phi_{j}\bm{1}_{2}+\cos\theta_{j}\bm{1}_{3}$. The equations:

[TABLE]

describe rotation along Jeffery orbits, and are obtained as the individual Cartesian components of the Faxen angular velocity relation, (19), given above ($\bm{\Omega}_{p0}\cdot\bm{1}_{r_{b_{2}}}=\frac{d\theta_{j}}{dt}$, $\bm{\Omega}_{p0}\cdot\bm{1}_{r_{b_{1}}}=-\sin\theta_{j}\frac{d\phi_{j}}{dt}$). In terms of $C$ and $\tau$, (23a) and (23b) take the form $\frac{dC}{dt}=0$ and $\frac{d\tau}{dt}=\frac{\beta}{\kappa+\kappa^{-1}}$. The former must be the case by definition, while the latter may be used to transform the time-averaged integral in (22) into a $\tau$-averaged one. Inverting the definitions above, one may write $\theta_{j}$ and $\phi_{j}$ as:

[TABLE]

which are used below in evaluating $\langle V_{p}\rangle$.

For $Re_{p}$ fixed, the Jeffery-averaged description of inertial migration, given by (22), breaks down for sufficiently small or large $\kappa$. This is because the leading order Jeffery angular velocity becomes small in these limits, being $O(\kappa^{-2})$ for fibers ($\kappa\gg 1$) close to flow-alignment, and $O(\kappa^{2})$ for flat disks ($\kappa\ll 1$) close to alignment with the gradient-vorticity plane. As a result, the $O(Re_{p})$ inertial correction becomes comparable in magnitude to the Jeffery contribution, leading to a slow down and eventual arrest of rotation, first shown by Subramanian & Koch (2005) for the case of a slender fiber. Herein, we will nevertheless confine ourselves to analyzing the Jeffery-averaged approximation, only noting that it remains valid provided $Re_{p}\,\kappa/\ln\kappa\ll 1\,(Re_{p}/\kappa^{2}\ll 1)$ for $\kappa\gg 1\,(\kappa\ll 1)$, an increasingly restrictive assumption for extreme-aspect-ratio particles(Subramanian & Koch, 2005; Marath & Subramanian, 2017). The consequence of an inertia-induced slow down at leading order, for the said particles, will be analyzed in a later communication.

A final point worth mentioning is that, for thin oblate spheroids with $\kappa<0.14$, the aforementioned inertial drift time scale of $O(Re_{p}^{-1}H/V_{\text{max}})$ only applies in the absence of stochastic orientation fluctuations. In presence of such fluctuations, either of a thermal origin or otherwise (for instance, due to pair-hydrodynamic interactions or weak turbulence; see Subramanian & Marath (2022)), there is a barrier-hopping time associated with the eventual equilibration between the numbers of spinning and tumbling spheroids in a manner independent of the initial orientation distribution (Dabade et al., 2016; Marath et al., 2017; Marath & Subramanian, 2018). This time scale increases exponentially with decreasing amplitude of the fluctuations, becoming much longer than the nominal drift time scale, and likely comparable to the time scale for inertial migration. Under these conditions, the migration dynamics will have a probabilistic rather than deterministic character, being described by a kinetic equation for the probability density that is a function of both $C$ and the transverse channel coordinate ($s$) - we briefly revisit this issue when calculating the lift velocity for arbitrary $C$ later in this section.

3.2 Scaling analysis and the Point-particle formulation

The dominant contribution to the volume integral in (22), for $\lambda\ll 1$, can arise from either scales of $O(L)$ or those of $O(H)$, the inertial screening length being irrelevant in the small-$Re_{c}$ limit. In order to assess the relative importance of these contributions, we consider the intermediate asymptotic interval $1\ll r\ll\lambda^{-1}$ ($r$ is measured in units of $L$), in which case $\bm{u}^{t2}\sim 1/r$, $\bm{u}_{s}\sim\beta/r^{2}+\gamma/r^{3}$, corresponding to the Stokeslet scaling for the test velocity field, and the stresslet-cum-force-quadrupole scaling for the velocity field in the actual problem. Using these $r$-scalings along with $\bm{u}^{\infty}\sim\beta r+\gamma r^{2}$ for the ambient flow, and separating the estimates for the linear and non-linear parts of the integrand in (22), one obtains:

• •

$\bm{u}^{t2}\cdot\left(\dfrac{\partial\bm{u}_{s}}{\partial t}+\bm{u}_{s}\cdot\bm{\nabla u}^{\infty}+\,\,\bm{u}^{\infty}\cdot\bm{\nabla u}_{s}\right)\sim\dfrac{\beta^{2}}{r^{3}}+\dfrac{\beta\gamma}{r^{2}}+\dfrac{\gamma\beta}{r^{4}}+\dfrac{\gamma^{2}}{r^{3}}\,\,\,\,$ (linear in $\bm{u}_{s}$),

• •

$\bm{u}^{t2}\cdot(\bm{u}_{s}\cdot\bm{\nabla}\bm{u}_{s})\sim\dfrac{\beta^{2}}{r^{6}}+\dfrac{\beta\gamma}{r^{7}}+\dfrac{\gamma^{2}}{r^{8}}\,\,\,\,\,\,\,$ (non-linear in $\bm{u}_{s}$).

Note that, over the range of scales under consideration, wall effects only contribute at a smaller order in $\lambda$, and hence, use of the unbounded domain estimates above. Next, using $dV\sim O(r^{2}dr)$, one obtains the following estimates for the contributions of the different terms to the lift velocity integral:

[TABLE]

The algebraically growing contribution in (25a) will be dominated by scales of $O(H)$ ($r\sim\lambda^{-1}$: the outer region), while contributions in (25a,b) that decay with $r$ will be dominated by length scales of $O(L)$ ($r\sim O(1)$: the inner region). The algebraic growth with $r$ will be cut off for $r\gtrsim\lambda^{-1}$ due to the more rapid decay of the disturbance velocity fields induced by wall-induced screening (recall that this more rapid decay was used in neglecting the integral over $S_{\infty}$ in (10)). Likewise, the apparent divergence for $r\to 0$, for the algebraically decaying terms, will be cut off at $r\sim O(1)$ by the finite size of the spheroid. The terms proportional to $\ln r$ in (25a) imply the dominance of the matching interval $1\ll r\ll\lambda$, resulting in a leading order contribution proportional to $\ln\lambda^{-1}$, with logarithmically smaller contributions arising from the inner and outer regions. Use of these cutoffs leads to the estimates:

[TABLE]

for the contributions from terms linear and nonlinear in $\bm{u}_{s}$. Using the definitions of $\beta$ and $\gamma$ given below (3), and reverting to the dimensional form, imply a leading order lift velocity of the form $V_{\text{max}}^{2}L^{3}/(\nu H^{2})[\ln(H/a)+O(1)]$. Here, the contribution to the $O(1)$ term within brackets arises from the $\beta\gamma\lambda^{-1}$ term in (26a) pertaining to the outer region, and from the $\beta^{2}$ terms in both (26a) and (26b) pertaining, respectively, to the outer and inner regions; the logarithmically larger lift contribution arises from the $\beta^{2}\ln\lambda^{-1}$ term in (26a) that pertains to the matching region. However, the matching-region $\beta^{2}\ln\lambda^{-1}$ contribution and the inner-region $\beta^{2}$-contribution must both vanish by symmetry because they cannot involve boundaries, and therefore relate to the time-averaged lift on a neutrally buoyant spheroid in an unbounded simple shear flow. Since such a spheroid cannot exhibit a net lateral drift, only the outer-region $\beta^{2}$-contribution survives. This and the outer-region $\beta\gamma$-contribution both yield a dimensional lift velocity of $V_{\text{max}}^{2}L^{3}/(\nu H^{2})$; in dimensionless terms, the Jeffery-averaged volume integral in (22) is $O(1)$, with the scaled lift velocity being $O(Re_{p})$. More detailed arguments along these lines given in Anand & Subramanian (2022) show that the next order contribution to the volume integral for a neutrally buoyant sphere arises from the inner region, involves an additional factor of $\lambda$, but leads to a qualitative alteration of the lift velocity profiles for large $Re_{c}$. Owing to the acceleration terms in (21) being $O(\lambda^{2})$, the aforementioned inner-region contribution will also be relevant to a neutrally buoyant spheroid.

Owing to the dominant contribution arising from length scales of $O(H)$ implied by the above arguments, the spheroids in the actual and test problems can be replaced by the corresponding point singularities. Thus, $\bm{u}_{s}$ and $\bm{u}^{t2}$ may be approximated as being induced by a time dependent stresslet ($\bm{u}_{\text{str}}$) and a Stokeslet due to a point force directed along the negative gradient direction ($\bm{u}_{\text{St}}$), respectively. The Jeffery-averaged lift velocity given by (22), at leading order in $\lambda$, may therefore be written as:

[TABLE]

where, on account of the subdominant nature of the length scales of $O(L)$, the domain of integration has now been extended to include the particle volume $V^{P}$. Thus, (27) is the leading order point particle approximation for the lift velocity for $Re_{c}\ll 1$.

In (27), $\bm{u}_{\text{St}}$ is time independent since the unit force points in a fixed direction, despite the changing orientation of the test spheroid. Further, $\langle\partial\bm{u}_{\text{str}}/\partial t\rangle=0$ for rotation along Jeffery orbits. Also noting that the ambient flow is steady in the chosen non-rotating reference frame, (27) reduces to:

[TABLE]

Using $r_{2}\sim O(\lambda^{-1})$ in (20) on account of the outer-region dominance, the Faxen’s correction to spheroid translation turns out to be $O(\lambda^{2})$ smaller than the terms linear and quadratic in $r_{2}$, and one may therefore approximate the ambient flow in (28) as $\bm{u}^{\infty}\approx(\beta r_{2}+\gamma r_{2}^{2})\bm{1}_{1}$.

With the spheroid volume neglected, the volume integral in (28) is most easily evaluated by Fourier transforming the flow ($r_{1}$) and vorticity ($r_{3}$) coordinates, the partial Fourier transform being defined as:

[TABLE]

Applying the convolution theorem along these coordinates (Arfken & Weber, 1999), and substituting the above approximate form of $\bm{u}^{\infty}$ in (28), one obtains:

[TABLE]

in terms of the partial Fourier transforms of the Stokeslet and stresslet velocity fields; here, $\bm{k}_{\perp}\equiv(k_{1},k_{3})$ and $y_{2}=s\lambda^{-1}$ is the transverse distance of the stresslet from the lower wall. The expressions for $\hat{\bm{u}}_{\text{St}}$ is derived in Appendix A, while that for $\langle\hat{\bm{u}}_{\text{str}}\rangle$ is developed in the next subsection.

3.3 Solution for $(\langle\bm{u}_{\text{str}}\rangle$, $\langle p_{\text{str}}\rangle)$

The Jeffery-averaged disturbance fields ($\langle\bm{u}_{\text{str}}\rangle$, $\langle p_{\text{str}}\rangle$), that appear in the point-particle approximation of the inertial lift velocity given by (28), satisfy:

[TABLE]

with the boundary conditions:

[TABLE]

where the stresslet singularity on the RHS of (31a) approximates the torque-free neutrally buoyant spheroid, in an unbounded simple shear flow, on scales much larger than $O(L)$. The stresslet coefficient $\bm{S}$ is a function of the spheroid orientation $\bm{p}$, being given by(Marath & Subramanian, 2017):

[TABLE]

where $E_{ij}=\frac{\beta}{2}(\delta_{i1}\delta_{j2}+\delta_{i2}\delta_{j1})$ is the rate of strain tensor associated with the local simple shear. The simple shear flow may be resolved into an axisymmetric extension aligned with $\bm{p}$, longitudinal planar extensions in a pair of orthogonal planes containing $\bm{p}$, and a pair of transverse planar extensions in the plane perpendicular to $\bm{p}$(Subramanian & Koch, 2006; Dabade et al., 2016). The $A_{i}$’s in (33) are the $\kappa$-dependent stresslet amplitudes corresponding to the aforesaid component flows, there being only three of these (in contrast to the five component flows) owing to the axisymmetry of the spheroidal geometry(Kim & Karrila, 1991; Marath & Subramanian, 2017). For $\kappa>1$, these amplitudes are given by:

[TABLE]

The corresponding expressions for an oblate spheroid ($\kappa<1$), may be obtained by first substituting $\kappa=\xi_{0}/(\xi_{0}^{2}-1)^{1/2}$ in terms of the coordinate $\xi_{0}\,(>1)$ labeling the surface of the spheroid, and then using the transformation $d\to-\iota d,\xi_{0}\to\iota(\xi_{0}^{2}-1)^{1/2}$ in the dimensional form of the stresslet obtained from multiplying $\bm{S}$ above by $\mu L^{3}V_{\text{max}}/H$ with $L=d\xi_{0}$(Marath & Subramanian, 2017). Note that $\bm{p}$, and thence $\bm{S}$, is a function of time owing to rotation along Jeffery orbits. In the limit of a sphere, $A_{1}=A_{2}=A_{3}=-\frac{20\pi}{3}$, and $\bm{S}$ reduces to $-20\pi\bm{E}/3$, independent of time, corresponding to the stresslet induced by a freely rotating sphere that leads to the well known Einstein coefficient in the suspension viscosity.

While one may, in principle, solve for ($\langle\bm{u}_{\text{str}}\rangle$, $\langle p_{\text{str}}\rangle$) in a manner similar to that for the bounded-domain Stokeslet fields (see Appendix A) by defining the disturbance fields as the sum of an unbounded-domain contribution, and a contribution that accounts for the confinement induced by channel walls, one may also obtain them directly as a gradient of the bounded-domain Stokeslet fields, derived in Appendix A; the gradient is with respect to the location $\bm{y}$ of the Stokeslet (Swan & Brady, 2010). Therefore,

[TABLE]

where the bounded-domain Stokeslet field corresponding to a point force $\bm{F}$ is given by $\bm{J}\cdot\bm{F}$; the Fourier transform ($\hat{\bm{J}}$) of the tensor $\bm{J}$ is defined in Appendix A (see text below (96)). Since it is $\langle\hat{\bm{u}}_{\text{str}}\rangle$ that appears in the final expression for the point-particle approximation of the lift velocity, (30), we Fourier transform (37) to obtain:

[TABLE]

where the third-order tensor $\hat{N}_{ijm}$ is defined as:

[TABLE]

In (38), $\langle\bm{S}\rangle=\frac{1}{2\pi}\int_{0}^{2\pi}\bm{S}\,d\tau$ denotes the time-averaged stresslet corresponding to a neutrally buoyant spheroid rotating along a given Jeffery orbit in an unbounded simple shear flow (with unit velocity gradient). The expression for $\langle\bm{S}\rangle$ is given below in the next subsection.

3.4 The time-averaged lift velocity $\langle V_{p}\rangle$

Using (38) for the stresslet velocity field in (30), the Jeffery-averaged lift velocity takes the form:

[TABLE]

Note that the averaging operation ($\langle.\rangle$) in (40) is independent of the shear rate since it involves the scaled Jeffery angular velocity which is only a function of $\kappa$ and $\phi_{j}$ for the inertially stabilized orbits. Thus, $\langle\bm{S}\rangle$ must be linear in $\bm{E}$, with a prefactor that is a function of $C$ and $\kappa$. For simple shear flow, this implies $\langle S_{12}\rangle=\langle S_{21}\rangle$ are the only non-zero components of the Jeffery-averaged stresslet. The following expressions for these components may be obtained by starting from (33), with $\bm{p}$ expressed in terms of $C$ and $\tau$ using (24ab), followed by an integration over $\tau$:

[TABLE]

for prolate spheroids, and

[TABLE]

for oblate spheroids. For both (41) and (42)$,\lim_{\kappa\rightarrow 1}\langle S_{12}\rangle=-10\pi/3$ independent of $C$, which yields the Einstein coefficient ($5/2$) in the $O(\phi)$ contribution to the suspension viscosity, $\phi$ being the sphere volume fraction; the $C$-independence arises owing to the sphere orientation being a degenerate degree of freedom. At the other extreme, one has $\lim_{\kappa\to\infty}\langle S_{12}\rangle|_{C=\infty}=-2\pi/(3\kappa\ln\kappa)$ in (41). Here, the factor $1/\ln\kappa$ arises from viscous slender body theory(Subramanian & Koch, 2005), while the additional $\kappa^{-1}$-factor reflects the probability of occurrence of non-aligned orientations that contribute dominantly to the viscosity of a dilute non-interacting suspension of slender fibers (Leal & Hinch, 1971; Dabade et al., 2016).

For times much longer than $O(Re_{p}^{-1}H/V_{\text{max}})$, only the inertially stabilized orbits are relevant to cross-stream migration. As mentioned earlier, for prolate spheroids, this is the tumbling orbit ($C=\infty$), in which case, $\langle S_{12}\rangle$ in (41) reduces to:

[TABLE]

For oblate spheroids, one can have either the spinning ($C=0$) or tumbling ($C=\infty$) orbits depending on $\kappa$, and accordingly, $\langle S_{12}\rangle$ in (42) reduces to:

[TABLE]

(43)-(45) will be used to determine $\langle V_{p}\rangle$ below. In light of the above, one may write (38) in the form:

[TABLE]

with the Jeffery-averaged lift velocity given by:

[TABLE]

Since the dominant contributions to the integral in (47) come from scales of $O(H)$ (or, equivalently, $k_{\perp}\sim O(H^{-1})$), it is natural to transform to rescaled variables $\bm{k}_{\perp}=\bm{k}_{\perp}^{\prime\prime}\lambda,r_{2}=r_{2}^{\prime\prime}/\lambda$. The original disturbance fields may be written in terms of the rescaled ones as: $\hat{\bm{u}}_{\text{St}}=\hat{\bm{u}}_{\text{St}}^{{}^{\prime\prime}}/\lambda$, $\hat{\bm{N}}=\hat{\bm{N}}^{\prime\prime}$, and further, noting that $y_{2}^{\prime\prime}=s$ and $d\bm{k}_{\perp}\,dr_{2}\equiv\lambda d\bm{k}_{\perp}^{\prime\prime}dr_{2}^{\prime\prime}$, one obtains,

[TABLE]

where $\gamma^{\prime\prime}=\gamma\lambda^{-1}=-4$, and $\langle S_{12}\rangle$ is given by (43), or (44) and (45), depending on $\kappa$. The essential consequence of the rescaling above is to show that the volume integral in (48) is independent of $\lambda$. In fact, the integral is only a function of $s$ (the particle location in the channel; see Figure 1), with the additional dependence on $\kappa$ entirely contained in $\langle S_{12}\rangle$. Accounting for the scaling $V_{\text{max}}\lambda$ used for $V_{p}$, (48) shows that the point-particle framework leads to a lift velocity of $O(V_{\text{max}}\lambda Re_{p})$, which was used earlier in section 3.1 to estimate the inertial migration time scale ($t_{\text{lift}}$).

A consequence of the integral in (48) being independent of $\kappa$ is that, within the Jeffery-averaged framework used, the inertial lift velocity of a sphere, and the time-averaged lift velocity of a spheroid, only differ by a multiplicative function of $C$ and $\kappa$ in general, and only by a function of $\kappa$ if one assumes the spheroid to have settled onto the inertially stabilized Jeffery orbit

• this multiplicative function is the ratio of the Jeffery-averaged spheroid stresslet to the sphere stresslet, being equal to $-\frac{3\langle S_{12}\rangle}{10\pi}$. An immediate consequence of this proportionality relation is that the shapes of the lift velocity profiles for the two cases must be identical, with the zero-crossings in particular (that correspond to the Segre-Silberberg equilibria for a sphere) being identical. In other words, a change in $\kappa$ only affects the magnitude of the inertial lift, the equilibrium positions being unaffected. In light of the obvious importance of the above conclusion for shape-sorting applications, it is worth documenting the reasons that leads to the simple proportionality relation:

The dominant contribution from the linearized inertial terms on scales of $O(H)$. The latter outer-region dominance led to the point-particle approximation, and thence, to the time independence of the test velocity field in the reciprocal theorem volume integral. 2. 2.

The asymptotic smallness, for small confinement ratios, of the additional contributions arising from the translational and angular accelerations associated with the motion of the neutrally buoyant spheroid in the Stokes limit, 3. 3.

The approximation of the time averaging based on the $Re_{p}-$independent Jeffery angular velocity, which leads to the time-averaged spheroid stresslet being proportional to $\bm{E}$ (just as the sphere stresslet).

The Jeffery-averaged lift velocity in (48) may be written in the compact form:

[TABLE]

where the functions $F(s)$ and $G(s)$ are given by:

[TABLE]

and satisfy $F(s)=-F(1-s)$ and $G(s)=G(1-s)$, consistent with the antisymmetry of the lift velocity profile about the channel centerline.

In (49), the $\beta^{2}$ contribution characterizes the effect of the asymmetrically located walls on the disturbance stresslet, and leads to migration away from the walls. In contrast, the $\beta\gamma^{\prime\prime}$ contribution, which characterizes the interaction of the stresslet with the ambient profile curvature, causes migration away from the channel centerline. The Segre-Silberberg equilibria emerge from a balance between these two opposing effects. Substituting $\beta=4(1-2s)$ and $\gamma^{\prime\prime}=-4$ in (49), the final expression for the time-averaged lift velocity for a neutrally buoyant spheroid of an arbitrary aspect ratio $\kappa$, for $Re_{c}\ll 1$, is given by:

[TABLE]

The Fourier integrals in (50) and (51) can be calculated using plane polar coordinates, $k_{1}^{\prime\prime}=k_{\perp}^{\prime\prime}\cos\phi$, $k_{3}^{\prime\prime}=k_{\perp}^{\prime\prime}\sin\phi$, with $d\bm{k}_{\perp}^{\prime\prime}\equiv k_{\perp}^{\prime\prime}dk_{\perp}^{\prime\prime}d\phi$; here $k_{\perp}^{\prime\prime}\in[0,\infty)$ and $\phi\in[0,2\pi]$. The integrals over the transverse coordinate $r_{2}^{\prime\prime}$ and the azimuthal angle $\phi$ may be done analytically, reducing (50) and (51) to the following one-dimensional integrals:

[TABLE]

with $I(k_{\perp}^{\prime\prime},s)$ and $J(k_{\perp}^{\prime\prime},s)$ defined in Appendix B. The $k_{\perp}^{\prime\prime}$-integrals above are evaluated numerically for various $s$, using Gauss-Legendre quadrature, after replacing the infinite interval with a finite one - $(0,K_{\text{max}})$. Numerical convergence, especially near the channel walls, depends sensitively on $K_{\text{max}}$. As the spheroid approaches either wall, the contribution to the $k_{\perp}^{\prime\prime}$-integrals comes from progressively smaller scales (compared to $H$) or, equivalently, from larger and larger $k_{\perp}^{\prime\prime}$. In fact, the relevant scale changes from $H$ ($k_{\perp}^{\prime\prime}\sim O(1)$) to the distance from the wall - either $s$ (lower wall) or $1-s$ (upper wall). To analyze the lift close to the lower wall, for example, one defines a rescaled wavenumber $k_{\perp}^{\prime\prime}=k_{w}/s$; next, expanding the integrands in (53) and (54) for $s\to 0$ with $k_{w}$ fixed, the near-wall lift velocity takes the form:

[TABLE]

The limiting value near the upper wall is equal in magnitude, but of an opposite sign, as must be the case by symmetry. The integral in (55) can be evaluated analytically, and yields

[TABLE]

with the upper and lower signs pertaining to the corresponding channel wall. For a sphere, (56) reduces to $55Re_{p}/36$, a value originally given by Vasseur & Cox (1976), albeit without an accompanying explanation.

The near-wall value above is proportional to $\beta^{2}$, implying that the inertial lift in this limit is generated due to interactions on scales of order the spheroid-wall separation, and the sign points to a repulsion. The curvature-induced $\beta\gamma^{\prime\prime}$ contribution is asymptotically small in this limit, implying that the near-wall lift at leading order may be determined independently from considering a sphere or a spheroid in the vicinity of a single plane wall subject to a linear shearing flow. Another important point is that the near-wall lift velocity approaches a finite value, rather than zero, as one approaches the channel walls. The latter must be the case on account of the eventual diverging lubrication resistance associated with the thin intervening fluid layer between the particle and the wall, and the discrepancy is on account of the regime of validity of the present calculation. As stated in the problem definition in $\S$2, the near-wall limit above pertains to spheroid-wall separations that, although much smaller than $O(H)$ are nevertheless much larger than $L$. We comment further on the nature of the lift force profile for small separations in $\S$4.

Based on the discussion above, it is clear that the choice of any finite $K_{\text{max}}$, however large, only leads to converged lift velocities down to a certain nonzero $s$. In the results shown in the next subsection, we have chosen $K_{\text{max}}=10^{4}$ to ensure accuracy down to $s=0.001$ without the aid of a large-$k_{\perp}^{\prime\prime}$ asymptote; this choice also enables a close approach to the near-wall limiting value obtained above. With this choice, $200$ quadrature points are sufficient to obtain converged results for arbitrary $\kappa$.

3.5 Results and Discussion

In Figure 2, we compare the lift velocity profile for a sphere, obtained by taking $\langle S_{12}\rangle=-\frac{10\pi}{3}$ in (52), with profiles plotted using the data from figure 2 of Ho & Leal (1974) and figure 8 of Vasseur & Cox (1976); the data was obtained by digitizing the latter figure. There is a clear mismatch between our profile and the Ho-Leal one, an aspect that we will comment on in the conclusions. The result of Vasseur & Cox (1976) shows good agreement with the present calculation throughout the channel. The Segre-Silberberg equilibria, corresponding to the zeroes of $\langle V_{p}\rangle$, are located at $s\approx 0.182$ and $0.818$; these are at a (dimensional) distance of $0.636\times H/2$ from the centerline, which agrees well with the intermediate annulus at $\sim 0.6\times$pipe radius in the original experiments(Segre & Silberberg, 1962a). Note that the lift velocity profiles in both our and the Vasseur & Cox (1976) analysis asymptote to the wall values $\pm 55/36$ (the horizontal dashed magenta lines in Figure 2) given in the previous subsection. As already argued therein, $H$ is no longer the relevant scale (at leading order) close to the wall. This implies that $Re_{c}$ should also not be relevant, and the near-wall lift values should therefore be independent of $Re_{c}$. This is validated in the next section where the lift velocity profiles are seen to approach (56) for $s\to 0,1$ even for $Re_{c}\gtrsim O(1)$, although this approach occurs in a shrinking neighborhood of the wall with increasing $Re_{c}$. The inset in Figure 2 plots the $\beta^{2}$ and $\beta\gamma^{\prime\prime}$ contributions to the inertial lift, confirming that the latter curvature-induced contribution becomes vanishingly small at the walls.

As mentioned after (48), while the magnitude of the Jeffery-averaged lift velocity is sensitive to the spheroid aspect ratio, the equilibrium locations remain the same as those for a sphere. This is seen from Figures 3 and 4, which show that the inertial lift profiles for both tumbling prolate spheroids and spinning/tumbling oblate spheroids of different $\kappa$ have the same zero crossings. In Figure 3a, the inertial lift for a tumbling prolate spheroid at a fixed $s$ is seen to decrease with increasing $\kappa$, consistent with the decreasing magnitude of the disturbance velocity field. To better examine the large-$\kappa$ limit, in Figure 3b we plot the lift profiles normalized by the large-$\kappa$ scaling of $\langle S_{12}\rangle$ found earlier. The scaled profiles approach a $\kappa-$independent limiting form for $\kappa\to\infty$, although the approach is non-monotonic - the scaled profile with the maximum amplitude corresponds to $\kappa=5$ - reflecting the non-monotonic approach of the scaled stresslet, $\kappa\ln\kappa\langle S_{12}\rangle$, to its infinite-$\kappa$ limit of $-\frac{2\pi}{3}$. In Figure 4a, the inertial lift velocity of a spinning oblate spheroid is seen to decrease in magnitude as $\kappa$ decreases from $1$ to $0.14$, again consistent with a decrease in the magnitude of the disturbance velocity field. At the latter $\kappa$, a reversal of the inertia-induced orbital drift leads to both tumbling and spinning modes being stabilized, with the respective basins of attraction being demarcated by a pair of unstable limit cycles on the unit sphere (Dabade et al., 2016; Einarsson et al., 2015). In Figure 4b, we therefore plot the inertial lift profiles of both spinning and tumbling oblate spheroids for $0<\kappa<0.14$. The lift profiles for the spinning spheroids approach a finite limiting form for $\kappa\rightarrow 0$, while those for tumbling spheroids are much smaller in magnitude, indicative of the much weaker disturbance field in this case. The latter arises from the much smaller value of the induced stresslet during the prolonged gradient-aligned phase of a thin oblate spheroid in the tumbling mode. The inertial lift for a tumbling oblate spheroid, in fact, goes to zero on account of $\langle S_{12}\rangle$ being $O(\kappa)$ for $\kappa\to 0$, this scaling reflective of the preponderance of orientations with $\bm{p}$ close to the gradient-vorticity plane, as mentioned earlier; the inset in Figure 4b shows that the tumbling-oblate-spheroid lift profiles, when scaled by $\kappa$, approach a finite limiting form for $\kappa\to 0$.

The invariance of the equilibrium locations associated with the lift profiles in Figures 3 and 4, with changing $\kappa$, implies that inertia-induced shape-sorting of spheroidal particles is precluded in sufficiently long channels. All spheroids regardless of $\kappa$ will end up migrating to the same pair of Segre-Silberberg locations in such channels. However, one may achieve separation in channels short enough for the residence time of $O(L_{c}/V_{\text{max}})$ to be comparable or smaller than $t_{\text{lift}}$; this requires $L_{c}\lesssim HRe_{p}^{-1}\lambda^{-1}$, $L_{c}$ being the channel length. Shape-sorting of spheroids would occur in these channels owing to differential rates of migration: for instance, prolate spheroids with $\kappa\sim O(1)$ will migrate to the equilibrium locations relatively rapidly, while those with larger $\kappa$ remain close to their initial positions due to weaker lift forces. A similar scenario would prevail for oblate spheroids provided $\kappa\gtrsim 0.14$. For $\kappa\lesssim 0.14$, the disparity in magnitudes of the lift velocities for spinning and tumbling oblate spheroids would lead instead to sorting of spheroids with the same $\kappa$, but in different (tumbling vis-a-vis spinning) orientation modes.

For even shorter channels with $L_{c}\lesssim O(HRe_{p}^{-1})$, the residence time is no longer enough for a spheroid to be able to migrate to its stable Jeffery orbit. For channel residence times of $O(L_{c}/V_{\text{max}})\lesssim O(H/V_{\text{max}}Re_{p}^{-1})\gg H/V_{\text{max}}$, the cross-stream migration is still determined, at leading order, by a Jeffery-averaged lift velocity, but one pertaining to orbits other than the tumbling or spinning mode. Thus, $\langle V_{p}\rangle$ is now also a function of $C$, being given by (49) with $\langle S_{12}\rangle$ given by its full form involving both $C$ and $\kappa$; see (41) and (42). The $C-$dependent lift velocity profiles for $\kappa=2$ and $\kappa=0.5$ are shown in Figures 5a and b, respectively. Since, for a fixed $\kappa$, the amplitude of the disturbance field is the largest for a tumbling prolate spheroid and a spinning oblate spheroid, the largest amplitude profiles in Figures 5a and b correspond to $C=\infty$ and $C=0$, respectively. Determining single-spheroid trajectories now requires solving the following system of coupled ODEs in $s$ and $C$:

[TABLE]

where $t_{s}=Re_{p}t$ is a slow time variable, and the functions $F_{n}^{f},\,G_{n}^{f},\,I_{n}$ and $J_{n}$ have been defined in Dabade et al. (2016) and Marath & Subramanian (2018).

At higher volume fractions, the trajectory of the spheroid will acquire a stochastic component on account of occasional pair-interactions that may be modeled for via a scattering kernel involving pre- and post-interaction orbit constants (Marath et al., 2017). In this case, the evolution of an initial distribution of spheroids, on time scales long compared to $T_{\text{jeff}}$, is described by a probability density $P(s,C)$ which satisfies a kinetic equation of the form:

[TABLE]

where $P^{\prime}\equiv P(s,C^{\prime})$, $\hat{P}\equiv P(s,\hat{C})$, $\hat{P}\equiv P(s,\hat{C}^{\prime})$; $\dot{s}$ and $\dot{C}$ being given by (57) and (58), respectively. The RHS of (59) denotes a Boltzmann-type kernel involving the pre- ($[\hat{C},\hat{C}^{\prime}]$) and post-($[C,C^{\prime}]$) interaction orbit-constant pairs.

At the end of $\S$3.4, we had highlighted the finite value of the inertial lift velocity attained even as the particle approaches either wall. This runs counter to one’s expectation of the lift velocity vanishing in this limit due to the diverging resistance associated with the thin lubricating layer of fluid between the particle and the wall. As mentioned therein, the discrepancy arises due to the present analysis only being valid for $s,1-s\gg\lambda$, a restriction that comes from treating the spheroid as a point (stresslet) singularity. Our analysis has to be supplemented by one that accounts for the finite size of the spheroid and is therefore valid for $s,1-s\sim O(\lambda)$, which in turn would connect to the lubrication regime corresponding to $s-\lambda,(1-s)+\lambda\ll\lambda$. Such a connection is possible for the case of a sphere based on results available in the literature. The inertial lift on a sphere, in presence of a single plane boundary subject to a linear shearing flow, has been evaluated numerically for $s\sim O(\lambda)$ by Cherukat & Mclaughlin (1994), with the aid of an integral expression obtained using bispherical coordinates. The authors found the lift force to always have a repulsive character (that is, to be directed away from the wall). The numerical value asymptoted to the near-wall limit ($\approx 55\pi Re_{p}/6$) of the two-wall (channel) problem mentioned above for $s\gg\lambda$, while asymptoting to a different finite value ($9.22Re_{p}$) when the sphere touches the wall ($s=\lambda$). The latter value was shown to agree with the inertial lift force acting on a stationary non-rotating sphere, in contact with a plane boundary, evaluated by Leighton & Acrivos (1985). These authors examined a non-rotating sphere in light of results obtained earlier (Goldman et al., 1967) which showed that the angular velocity of a torque-free sphere must decrease to zero logarithmically in the limit of a vanishing sphere-wall separation on account of the lubrication resistance associated with the relative tangential motion in the narrow gap; the finiteness of the force implies that the dominant contributions to the lift arise from the fluid domain outside the thin gap. When combined with the known $O(s-\lambda)^{-1}$ divergence of the translational resistance for normal approach towards the wall, one concludes that the inertial lift velocity for a sphere must start from $55Re_{p}/36$ for $\lambda\ll s\ll 1$, and eventually approach zero linearly in the limit $s-\lambda\ll\lambda$. In the reciprocal theorem formulation used here, the approach to zero would appear via the divergence of the test problem resistance coefficient. To our knowledge, analogous results for a spheroid when $s,1-s\sim O(\lambda)$ are not available, and will involve more effort. This is both because one does not have a spheroidal analog of the bispherical coordinate system used in Cherukat & Mclaughlin (1994), and because the orientation dynamics change from Jeffery rotation to pole-vaulting with approach towards either wall.

4 The inertial lift velocity for $Re_{c}\gtrsim O(1)$

Herein, we calculate the inertial lift profiles for $Re_{c}\gtrsim O(1)$ using a numerical shooting method employed originally by Schonberg & Hinch (1989), and later for higher $Re_{c}$’s by Asmolov (1999), both for the case of a sphere. For $Re_{c}\gtrsim O(1)$, the inertial screening length ($HRe_{c}^{-1/2}$) is of order the channel width or smaller, and one must solve the governing equations (4a,b) with the boundary conditions (5a-c) using a matched asymptotic expansions approach. The inner region is characterized by scales of $O(L)$, and the outer region by scales of $O(HRe_{c}^{-1/2})$, with the channel width $H$ also entering via the wall boundary conditions. Although we use a shooting method below to evaluate the lift velocity, it is worth noting that the reciprocal theorem formulation in section 3.1 remains valid for any $Re_{c}$, provided one uses the finite-$Re_{c}$ disturbance velocity field ($\bm{u}^{\prime}$) in the volume integral. Thus, the scaling arguments given in section 3.2 may be extended to finite $Re_{c}$, and in doing so, one finds that the dominant contribution to the inertial lift continues to come from scales asymptotically larger than $O(L)$. For $Re_{c}\ll 1$ in the previous section, this outer-region dominance meant solving the Stokes equations in a domain confined by plane channel walls with the spheroid approximated as a stresslet singularity. In the present section, this means solving the linearized Navier-Stokes equations at leading order, driven by the same stresslet singularity (since $Re_{p}$ is assumed small). Note that $Re_{c}$, interpreted as the square of the ratio of the two outer-region length scales ($H$ and $HRe_{c}^{-\frac{1}{2}}$) appears as a parameter in the governing equations below, and the inertial lift velocity is therefore now a function of $Re_{c}$.

4.1 The finite-$Re_{c}$ formulation

The Stokesian disturbance field due to a freely suspended particle in an ambient linear flow only decays algebraically, as $1/r^{2}$, at large distances, a fact already used in the scaling arguments in section 3.2. The algebraic decay implies that distinct outer-region expansions for both the velocity and pressure fields become necessary on scales of $O(HRe_{c}^{-1/2})$, with the leading terms satisfying the linearized Navier-Stokes and continuity equations. To write down these equations, one transforms to outer coordinates using $\bm{r}=Re_{p}^{-1/2}\bm{R}$, which corresponds to using $HRe_{c}^{-\frac{1}{2}}$ (rather than $L$ as in section 2) as the relevant length scale. The Stokesian rates of decay in the inner region suggest the scalings $\bm{u}^{\prime}=Re_{p}\bm{U}$ and $p^{\prime}=Re_{p}^{3/2}P$ for the leading order terms in the outer expansions, where $\bm{U}$ and $P$ satisfy:

[TABLE]

The Faxen correction contributes at a higher order in $\lambda$ and therefore, $\bm{u}^{\infty}\approx Re_{p}^{-1/2}(\beta R_{2}+\gamma^{\prime\prime}\,Re_{c}^{-1/2}R_{2}^{2})\delta_{i1}$ has been used in (60a). Equations (60a,b) must be supplemented by the following conditions:

[TABLE]

where (61b) denotes the no-slip conditions on the channel walls, while (61a) is the requirement of matching to the stresslet velocity field that arises as the far-field form of the inner-region Stokesian field. In (60a) and (61a), $\bm{S}$ is the tensorial amplitude defined in (33).

As explained in $\S$3.2, the separation between the migration and drift time scales implies one need only solve for the Jeffery-averaged lift velocity, and towards this end, we average (60a,b) and (61a,b) over a single period of the stable Jeffery orbit. One obtains:

[TABLE]

where $\langle\bm{U}\rangle$ satisfies:

[TABLE]

Here, $\langle S_{12}\rangle(\kappa)$ corresponds to the stresslet averaged over the relevant stable Jeffery orbit, and has been given in (43)-(45). Due to the linearity of (62a,b) and (63a,b), and the fact that the Jeffery-averaged spheroid stresslet tensor differs from that for a sphere only by a scalar multiplicative factor, $\langle\bm{U}\rangle$ and $\langle P\rangle$ differ from their spherical analogs only by $-\frac{3\langle S_{12}\rangle(\kappa)}{10\pi}$, corresponding to the ratio of the aforementioned stresslets. This proportionality relation must hold for any linear functional of the disturbance fields, and in particular, for the lift velocity that is a linear functional of $\langle U_{2}\rangle$; see (67) below. Thus, similar to the case of $Re_{c}\ll 1$, the Jeffery-averaged lift profiles at a given finite $Re_{c}$, for an arbitrary aspect ratio spheroid, have the same shape as those for a sphere at the same $Re_{c}$. It follows that the associated pair of equilibria are identical to those for a sphere regardless of $Re_{c}$, and as for a sphere(Schonberg & Hinch, 1989), must migrate wallward with increasing $Re_{c}$. It is worth reiterating that the requirement for a Jeffery-averaged analysis to remain valid becomes restrictive for extreme-aspect-ratio spheroids. The regime of validity was originally stated after equation (24), and expressed in terms of $Re_{c}$, is given by $Re_{c}\kappa/(\lambda^{2}\ln\kappa)\ll 1$ and $Re_{c}/(\lambda^{2}\kappa^{2})\ll 1$ for $\kappa\gg 1$ and $\kappa\ll 1$, respectively; these point to the restriction becoming more severe with increasing $Re_{c}$. The implications of the finite-$Re_{c}$ Jeffery-averaged analysis above, for shape-sorting, remain the same as those discussed in section 3.5.

For purposes of completeness, we now follow along the lines of Schonberg & Hinch (1989), and briefly present the manner in which inertial lift is determined for $Re_{c}\gtrsim O(1)$. After implementing the partial Fourier transform defined in (29), followed by some algebraic manipulation, one obtains the following coupled ordinary differential equations for $\langle\hat{P}\rangle$ and $\langle\hat{U}_{2}\rangle$:

[TABLE]

with the conditions,

[TABLE]

where $k_{\perp}^{2}=k_{1}^{2}+k_{3}^{2}$ as before. The delta-function forcing in (62a) leads to the following jump conditions across the particle location ($R_{2}=0$):

[TABLE]

where the superscripts ‘+’ and ‘$-$’ denote the limiting values attained on approaching $R_{2}=0$ from the regions $0<R_{2}\leq(1-s)\,Re_{c}^{1/2}$ and $-s\,Re_{c}^{1/2}\leq R_{2}<0$, respectively; these conditions are derived in Appendix C. The limiting form of $\langle\bm{U}\rangle$ in the matching region ($\bm{R}\ll 1$) is the sum of the singular stresslet contribution given in (63a), and a uniform flow along the gradient (cross-stream) direction that is a consequence of fluid inertia. The neutrally buoyant spheroid being force-free is convected by this uniform flow which therefore equals the inertial lift velocity, and may be determined from the limit of the inverse transform for $\bm{R}\to 0$:

[TABLE]

Here, $\Re\{.\}$ denotes taking the real part which eliminates the purely imaginary stresslet contribution. As indicated, one may use either $\hat{\langle U_{2}\rangle}^{-}$ or $\hat{\langle U_{2}\rangle}^{+}$ on account of continuity; see (66c). The governing ODEs (64a,b), the boundary conditions (65a,b) and the jump conditions (66a-d) can be solved using the shooting technique described in Appendix A of Schmid et al. (2002), a brief description of which is given in Appendix D.

The integral in (67) is evaluated numerically using a two-dimensional Gauss-Legendre quadrature over a circle of a sufficiently large radius $K_{m}$ in the $k_{1}\!-\!k_{3}$ plane. An analytical large-$k_{\perp}$ asymptote was calculated using the steps outlined in Hogg (1994), and added to the numerical integral, to improve convergence. The analysis involves expanding $\hat{\langle U_{2}\rangle}$ and $\hat{\langle P\rangle}$ in inverse powers of $k_{\perp}$, an ansatz valid only in an $O(k_{\perp}^{-1})$ neighborhood of $R_{2}=0$. As a result, satisfaction of (65a,b) is replaced by a farfield decay requirement for $R_{2}k_{\perp}\gg 1$. Including only the exponentially decaying solutions on either side of $R_{2}=0$, satisfying the jump conditions (66a-d), and then performing the inverse Fourier transform, one obtains:

[TABLE]

at leading order. As in $\S$3, the choice of $K_{m}$ is dictated by the distance of the particle from the walls. The relevant scale sufficiently close to either wall is still the spheroid-wall separation ($s$ or $1-s$) - the near-wall lift remains the same as that for $Re_{c}\ll 1$, being the farfield lift experienced by a particle in the presence of a single plane boundary. While one needs to keep increasing $K_{m}$ with approach to either wall (that is, for sufficiently small $s$ or $1-s$), to obtain a converged result, this increase is only necessary once $s$ or $(1-s)$ becomes less than $O(Re_{c}^{-1/2})$. Thus, to capture the wall-induced repulsion for finite $Re_{c}$, one needs to ensure $K_{m}Re_{c}^{1/2}s_{\text{min}},K_{m}Re_{c}^{1/2}(1-s)_{\text{min}}\gg 1$. In our calculations, we chose $K_{m}=200$ for $0.1\leq Re_{c}\leq 10$ to ensure accurate lift velocities down to $s,1-s\approx 0.03$. For $Re_{c}>10$, accurate lift profiles were obtained down to $s,1-s\approx 0.05$ for $K_{m}=40$.

4.2 Results and Discussion

We begin with Figure 6 which shows the inertial lift profiles for a sphere for $0.5\leq Re_{c}\leq 75$, along with the limiting small-$Re_{c}$ profile given by (52); only profiles in the half-channel have been plotted owing to their anti-symmetry about the centerline. The profiles for $Re_{c}=1$ and $75$ exhibit good agreement with the data extracted from Schonberg & Hinch (1989). Interestingly, the profiles for $Re_{c}=0.5,1$ and $10$ compare closely with the small-$Re_{c}$ limiting form. The inset shows the approach of the finite-$Re_{c}$ profiles to the near-wall value given by (56). Figure 7 shows that the sphere lift profiles for higher $Re_{c}$ agree well with those extracted from Asmolov’s data(Asmolov, 1999), all the way upto $Re_{c}=3000$; note that our profiles have been continued to smaller $s$ to emphasize the approach to the common near-wall limiting value mentioned above. For $Re_{c}\gtrsim 300$, the lift profiles begin to exhibit an intermediate region of oppositely signed curvature. It has recently been shown that including finite-size effects pushes this intermediate region towards the zero-lift line, eventually leading to the emergence of new equilibria closer to the centerline for sufficiently large $Re_{c}$(Anand & Subramanian, 2022). The inset in the said figure, on a logarithmic ordinate scale, helps highlight the rapid decrease in the lift magnitude for $Re_{c}\gtrsim O(10)$, reflective of weakening particle-wall interactions. Inertia-induced faster decay of the velocity field, on scales larger than $O(HRe_{c}^{-\frac{1}{2}})$, is responsible for the reduced influence of the walls with increasing $Re_{c}$. The rescaled abscissa in the inset highlights the $O(Re_{c}^{-\frac{1}{2}})$ neighborhood of the wall where the lift profiles begin to rise towards to near-wall limit.

Figure 8 shows the magnitude of the sphere lift velocity, as a function of $Re_{c}$, at different locations on either side of the Segre-Silberberg equilibrium ($s_{\text{eq}}=0.182$). In accordance with the above discussion, for all $s$ values considered, the lift velocity starts off on a small-$Re_{c}$ plateau which extends until $Re_{c}\approx 10$. Thus, the small-$Re_{c}$ approximation remains a good approximation well beyond $Re_{c}$’s of order unity. For $s=0.3$ and $0.4$, with increasing $Re_{c}$, the lift velocity directly transitions from the plateau to an eventual algebraic decrease, this being typical for all $s>s_{\text{eq}}$. On the other hand, the lift velocity magnitude for $s=0.1$ (and for all $s<s_{\text{eq}}$) exhibits a non-monotonic variation with an intermediate zero-crossing which, for $s=0.1$, is at $Re_{c}\approx 300$. This is due to the Segre-Silberberg equilibrium crossing the given $s$ in course of its wallward movement (with increasing $Re_{c}$). The lift velocity increases again at larger $Re_{c}$, but to a value smaller than the small-$Re_{c}$ plateau, finally transitioning to a steeper algebraic decrease.

As mentioned earlier, the Jeffery-averaged lift profiles for a spheroid of an arbitrary aspect ratio may simply be obtained by multiplying the sphere lift profile at the same $Re_{c}$ by the ratio of the stresslets. Consequently, features pertaining to the $Re_{c}$-dependence mentioned above, including the range of validity of the small-$Re_{c}$ approximation, remain true for spheroids. Figures 9a and b show the inertial lift profiles for tumbling prolate spheroids over a range of $\kappa$, for $Re_{c}=300$, with and without the large-$\kappa$ stresslet scaling. Figures 10a and b shows the inertial lift profiles for both spinning ($0.14\leq\kappa<1$) and spinning/tumbling ($0<\kappa<0.14)$ oblate spheroids, again for $Re_{c}=300$, with and without the small-$\kappa$ stresslet scaling. Note that the latter scaling is only used for tumbling oblate spheroids in the inset of Figure 10b, since $\langle S_{12}\rangle$ remains of order unity for spinning spheroids even as $\kappa\to 0$. As in $\S$3, the scaled lift profiles in Figure 9b, and in the inset of 10b, approach $\kappa$-independent limiting forms for $\kappa\to\infty$ and [math], respectively. The magnitudes of the lift velocity profiles reflect that of the disturbance field (via the stresslet), and therefore, decrease with increasing (decreasing) $\kappa$ for tumbling prolate (oblate) spheroids. Although not shown, the finite-$Re_{c}$ lift profiles may be evaluated for arbitrary $C$, and the $C$-dependence again correlates to the magnitude of the disturbance velocity field, being the largest for tumbling prolate and spinning oblate spheroids.

5 Conclusions

The primary result of this manuscript is to show that the inertial lift velocity profiles for neutrally buoyant spheroids in plane Poiseuille flow, for sufficiently small $Re_{p}$, can be calculated within a Jeffery-averaged approximation. Within this framework, spheroid lift profiles differ from the ones for a sphere only by a multiplicative factor, regardless of $Re_{c}$, which leads to the original Segre-Silberberg equilibria being dependent on $Re_{c}$ but not on the aspect ratio $\kappa$. The multiplicative factor above is a function of $C$ and $\kappa$ in the general case; for times relevant to cross-stream migration, and that are much longer than those characterizing the inertia-induced orientation drift, it is only a function of $\kappa$. Although the Jeffery-averaged approximation only requires $Re_{p}\ll 1$ for spheroids with $\kappa\sim O(1)$, it becomes increasingly restrictive for spheroids with asymptotically large and small $\kappa$, and with increasing $Re_{c}$. Deviations from this approximation owing to inertia-induced slow down of rotation and eventual arrest, that lead to $\kappa$-dependent equilibria for extreme-aspect-ratio spheroids, will be reported in a future investigation. Herein, we only note that accounting for the rotational slow down leads to the equilibria migrating back towards the centerline beyond a certain $\kappa$-dependent $Re_{c}$ threshold, this being consistent with both experiments(Masaeli et al., 2012) and computations(Chen et al., 2012).

En-route to deriving the Jeffery-averaged lift on a spheroid, we have presented, in some detail, the inertial lift profiles for a sphere in plane Poiseuille flow. While most of these results are known from earlier literature, and are spread across multiple efforts(Ho & Leal, 1974; Vasseur & Cox, 1976; Schonberg & Hinch, 1989; Asmolov, 1999), each of these only pertains to specific ranges of $Re_{c}$. Moreover, there are important features of the lift profile not reported in any of these efforts. We have expressed the small-$Re_{c}$ sphere lift velocity, obtained within the framework of a point-particle approximation, in terms of a one-dimensional Fourier integral, with detailed expressions for the integrands given in Appendix B. These expressions allow one to analytically determine the near-wall lift value. The arguments leading to this limiting value also show that it must be independent of $Re_{c}$, numerical evidence for which is provided by the finite-$Re_{c}$ lift profiles obtained using a shooting technique. Further, the near-wall limiting lift was also identified with the farfield limit of the lift acting on a finite-sized sphere moving parallel to a single plane wall subject to a linear shearing flow. This connection allows one, in principle, to construct a uniformly valid lift profile across the entire channel. We have also established the surprisingly large range of validity (up to $Re_{c}\approx 10$) of the small-$Re_{c}$ approximation derived first by Ho & Leal (1974) and Vasseur & Cox (1976).

The scaling arguments in $\S$3.2 are a crucial element of the overall analysis, and clearly show that the dominant scales contributing to the inertial lift remain much greater than the particle size, regardless of $Re_{c}$. This leads to use of a point-particle framework to obtain the leading order approximation for the inertial lift, with scaling arguments identifying both competing contributions due to profile curvature and wall-shear-induced repulsion. The first lift force calculation for a sphere, by Ho & Leal (1974) for $Re_{c}\ll 1$, used the same reciprocal theorem formulation as that given here. Although in a less transparent form, the said authors, via scaling arguments, did recognize the dominant contribution of the linearized inertial terms, on scales of $O(H)$, to the reciprocal theorem volume integral. The effect of confining plane boundaries was modeled differently. Rather than directly consider the relevant point singularity (Stokeslet or stresslet) between plane parallel walls, the authors used a partial Fourier representation of the full form of the unbounded domain velocity field, including both the finite-size terms and the quadrupolar disturbance induced by the quadratic component of the ambient flow, in order to derive the wall-induced contribution (the first reflection). For evaluation of the final volume integral, the physical space integration appears to have been done analytically, with the partial Fourier integral done numerically. Although more circuitous, carrying out the physical space integration should have led to a residual Fourier integral identical to the one obtained directly here from use of the convolution theorem. While the inaccuracy of the resulting profile, and the errors in the equilibrium locations, have been pointed out in the context of Figure 2, the detailed calculational procedure in Ho & Leal (1974) is correct, and their origin remains uncertain. In contrast to the intuitive scaling arguments here, the later small-$Re_{c}$ lift calculation of Vasseur & Cox (1976) made use of a velocity field that emerged from a formal matched asymptotics expansions approach, developed in earlier articles by Brenner and Cox (Cox & Brenner, 1967, 1968).

The scaling arguments in §3.2 can be generalized in several directions. First, the arguments apply virtually unchanged to pipe Poiseuille flow. Thus, for small pipe Reynolds numbers ($Re$), the dominant contributions to the inertial lift must arise from scales of order the pipe radius, and a calculation of this lift would involve approximating the particle as a stresslet within a cylindrical domain. The required Stokesian velocity fields should either be obtainable from that known for a Stokeslet in this domain(Liron & Shahar, 1978), or be derivable in a manner similar to that in Appendix A (adapted to cylindrical coordinates). Such an analysis would yield the radius of the Segre-Silberberg annulus for $Re\to 0$. Although finite-$Re$ lift force profiles for pipe Poiseuille flow have been computed earlier(Matas et al., 2009), the profiles differ significantly for the two lowest $Re$’s examined ($1$ and $30$), preventing one from inferring the range of validity of a small-$Re$ approximation. The small-$Re$ analysis in $\S$3 can also be generalized to ducts of non-circular cross sections, provided one has the confined Stokeslet field for the relevant cross-sectional geometry. The latter may be derivable for rectangular or elliptical cross sections, using the procedure in Appendix A, owing to availability of the Greens function the (2D) Laplacian. The broken symmetry for the duct case should then allow for the prediction of discrete inertial equilibria in the transverse plane, and their variation with cross-sectional aspect ratio. The requirement that the Stokeslet and stresslet fields be available for a given confined geometry, for a small-$Re$ inertial migration analysis in the same geometry to be possible, calls into question the effort of Hood et al. (2015) who seem to have analyzed migration in rectangular ducts based on the approach of Ho & Leal (1974) which was solely tailored to plane Poiseuille flow.

Importantly, the scaling arguments in $\S$3.2 allow for the systematic incorporation of finite-size effects that modify the leading order point-particle estimate of the inertial lift velocity for both spheres and anisotropic particles in plane Poiseuille flow. As mentioned in §3.2, and discussed in more detail in Anand & Subramanian (2022), the finite-size contributions involve an additional factor of $\lambda$, and pertain to the inner region (scales of order the particle size), implying that they arise independently of the outer-region point-particle contribution. It is shown in Anand & Subramanian (2022) that, for spheres, these contributions invariably become important for sufficiently large $Re_{c}$, leading to the emergence of a pair of new equilibria closer to the channel centerline, consistent with the results of recent experiments. In contrast to the point-particle analysis presented here, calculation of the finite-size contributions involves consideration of the nonlinear terms, and also requires knowing the disturbance velocity field induced by a sphere in an ambient quadratic flow. The analogous calculation for a spheroid will proceed along similar lines. Although more complicated, the disturbance velocity field in an ambient quadratic flow, for a spheroid of an arbitrary aspect ratio and orientation, may be constructed using a superposition of the appropriate vector spheroidal harmonics(Dabade et al., 2015, 2016). Crucially, the nonlinear inertial terms as well as the time dependence of the test velocity field, that arise in the inner region, imply that there is no longer a simple proportionality relationship between the finite-size contribution for a sphere and the Jeffery-averaged version of the same for a spheroid. As a result, one expects the incorporation of finite-size contributions to lead to $\kappa$-dependent equilibria, for neutrally buoyant spheroids, even within a Jeffery-averaged framework. This will be examined in a separate communication.

Appendix A

Herein, we will solve for the disturbance field due to a Stokeslet (point force) confined between plane parallel boundaries. Recall that the partially Fourier transformed Stokeslet, $\hat{\bm{u}}_{\text{St}}$, appears in the final expression for the lift velocity integral viz. (30). In physical space, the disturbance field $\bm{u}_{\text{St}}$ satisfies the following equations:

[TABLE]

with the boundary conditions:

[TABLE]

The problem of a Stokeslet in the vicinity of a single plane wall was solved for by Blake (1971), using the method of images. Liron & Mochon (1976) calculated the disturbance field due to a Stokeslet between two parallel walls, by taking repeated reflections of Blake’s single wall solution and superposing these as an infinite but convergent series. This approach, however, turns out to be quite tedious, and one can instead derive $\bm{u}_{\text{St}}$ using another method described originally by Vasseur & Cox (1976), and later used by Swan & Brady (2010). In this procedure, instead of taking repeated reflections of the single wall solution, one can satisfy the no-slip condition on both channel walls at one go. One writes the velocity field as

[TABLE]

with an analogous decomposition for the pressure field. Here, $\bm{u}^{\infty}_{\text{St}}$ is the Stokeslet velocity field in an unbounded domain, given by:

[TABLE]

where \bm{J}^{\infty}=\frac{1}{8\pi}\big{(}\frac{\bm{I}}{r}+\frac{\bm{rr}}{r^{3}}\big{)} is the Oseen-Burger’s tensor; here $\bm{r}=\bm{x}-\bm{y}$ defines the position of any point $\bm{x}$ in the domain relative to the Stokeslet location $\bm{y}$. The second contribution in (71) is the one that accounts for the no-slip conditions at the channel walls, and satisfies,:

[TABLE]

with no-slip boundary conditions written as:

[TABLE]

where $y_{2}=s\lambda^{-1}$. The solution to (73a,b) is easily obtained by Fourier transforming the flow and vorticity coordinates, with the partial Fourier transform being defined as:

[TABLE]

Fourier-transforming (73a,b) in accordance with (75), one obtains:

[TABLE]

where $i=1,2,3$, and $k_{\perp}^{2}=k_{1}^{2}+k_{3}^{2}$. The partial Fourier transform, $\bm{\hat{u}}^{w}_{\text{St}}$, satisfies the same boundary conditions as the physical space velocity field above except that one uses $\hat{\bm{u}}^{\infty}_{\text{St}}$ instead of $\bm{u}^{\infty}_{\text{St}}$ on the RHS. Here, $\hat{\bm{u}}^{\infty}_{\text{St}}=\hat{\bm{J}}^{\infty}\cdot\bm{1}_{2}$ with the partial Fourier-transform of the Oseen-Burger’s tensor given by:

[TABLE]

It is worth mentioning that, in contrast to the velocity field due to a Stokeslet in an unbounded domain (72) which only depends on $\bm{y}$ via the position vector $\bm{r}=\bm{x}-\bm{y}$, the bounded domain contribution, apart from its dependence on $\bm{r}$, also depends explicitly on the location of the singularity via the wall boundary conditions (74). Thus, $\bm{u}^{w}_{\text{St}}\equiv\bm{u}^{w}_{\text{St}}(\bm{r};\bm{y})$.

To begin with, one derives the equation governing the pressure field by taking the divergence of both sides in (76a) and using the incompressibility condition (76b), leading to,

[TABLE]

This ODE can be solved to give: $\hat{p}^{w}_{\text{St}}=[A_{m}(\bm{k}_{\perp};y_{2})e^{-k_{\perp}r_{2}}+B_{m}(\bm{k}_{\perp};y2)e^{-k_{\perp}r_{2}}]\delta_{m2}$, where $A_{m}$ and $B_{m}$ are unknown vectors. After substituting this solution in (76a), one may use variation of parameters (Arfken & Weber, 1999) to solve the resulting inhomogeneous ODE to obtain,

[TABLE]

where $\hat{\bm{J}}^{w}$, the Fourier transform of the second order tensor $\bm{J}^{w}$ ($\bm{u}^{w}_{\text{St}}=\bm{J}^{w}\cdot\bm{1}_{s}$), is defined as:

[TABLE]

Here, $d_{i}=k_{\perp}\delta_{i2}+\iota(k_{1}\delta_{i1}+k_{3}\delta_{i3})$ and $\bar{d}_{i}=k_{\perp}\delta_{i2}-\iota(k_{1}\delta_{i1}+k_{3}\delta_{i3})$. We will now determine the unknown second order tensors $C_{im}$ and $D_{im}$ and the vectors $A_{m}$ and $B_{m}$, using the no-slip conditions on the walls, which can be written as:

[TABLE]

where the superscripts ‘$L$’ and ‘$U$’ denote the value of the Fourier-transformed Oseen-Burger’s tensor calculated on the lower wall and upper wall, respectively. Using the incompressibility condition (76b),

[TABLE]

The wall boundary conditions (84) and (85) along with the relations (86) and (87) can be solved simultaneously to obtain the following:

[TABLE]

Finally, the partial Fourier transforme of the velocity field due to the Stokeslet confined between plane parallel walls is written as:

[TABLE]

where $\hat{\bm{J}}(k_{1},r_{2},k_{3};y_{2})=\hat{\bm{J}}^{\infty}(k_{1},r_{2},k_{3})+\hat{\bm{J}}^{w}(k_{1},r_{2},k_{3};y_{2})$, with $\hat{\bm{J}}^{\infty}$ and $\hat{\bm{J}}^{w}$ being defined in (80) and (83), respectively.

Appendix B

The functions $I(k_{\perp}^{\prime\prime},s)$ and $J(k_{\perp}^{\prime\prime},s)$ that appear in the integrands, in the expressions for $F(s)$ and $G(s)$ given by (53) and (54) in the main paper, are defined below:

[TABLE]

[TABLE]

Appendix C

Herein, we derive the jump conditions (66a-d). Starting off with the governing equations (64ab),

[TABLE]

where the prime (′) denotes differentiation wrt $R_{2}$. Since the singular forcings on the RHS arise from the highest order derivative, one can postulate the following forms:

[TABLE]

where the superscripts ‘$+$’ and ‘$-$’ denote the function definitions for $R_{2}>0$ and $R_{2}<0$, respectively, and $\mathcal{H}(R_{2})$ is the Heaviside function. Differentiating (100) twice gives,

[TABLE]

Differentiating (101) once gives,

[TABLE]

Integrating (101) gives,

[TABLE]

To obtain the jump in the aforementioned fields across $R_{2}=0$, one can multiply (99a) and (99b) with $R_{2}^{n}$, where $n=0$ and $1$. The resulting expressions can then be integrated with the help of the expressions (100)-(105). Multiplying (99a) with $R_{2}$ and integrating the resulting expression from $R_{2}=0^{-}$ to $0^{+}$ gives,

[TABLE]

Next, (99a) can be integrated across $R_{2}=0$ to yield,

[TABLE]

Equation (99b) can be integrated across $R_{2}=0$ to yield,

[TABLE]

Finally, one can multiply (99b) with $R_{2}$ and integrate from $R_{2}=0^{-}$ to $0^{+}$ to obtain:

[TABLE]

This concludes the derivation of all the jump conditions.

Appendix D

Equations (64a,b), written as a set of four first order ODEs, take the following form:

[TABLE]

where

[TABLE]

The general solution to equation (110) can be written as:

[TABLE]

where each of the $\bm{\Phi}_{i}^{-}$’s and $\bm{\Phi}_{i}^{+}$’s constitute a set of four linearly independent solution vectors, with the $c_{i}^{-}$’s and $c_{i}^{+}$’s being the corresponding integration constants. We choose the following $\bm{\Phi}_{i}^{-}$’s and $\bm{\Phi}_{i}^{+}$’s on the walls:

[TABLE]

[TABLE]

Using the boundary coundition at the lower wall (114) in (112), and the boundary condition on the upper wall (115) in (113), one obtains: $c_{3}^{-}=c_{4}^{-}=c_{3}^{+}=c_{4}^{+}=0$. Therefore,

[TABLE]

The NDSOLVE subroutine in the symbolic computation software Mathematica may be employed to solve for the solutions that appear in (116) and (117). To accomplish this, the entire integration interval in $R_{2}$ can be divided into sub-intervals bounded by ‘orthonormalization’ points $y_{i}$. One may then use NDSOLVE to integrate the relevant solution from $y_{i}$ to $y_{i+1}$, starting from either the upper or lower boundary. At $y_{i+1}$, one may use Gram-Schmidt orthogonalization to orthonormalize the solution vectors $\bm{\Phi}_{1}^{-}$ and $\bm{\Phi}_{2}^{-}$, and $\bm{\Phi}_{1}^{+}$ and $\bm{\Phi}_{2}^{+}$. The orthonormalization at regular intervals is necessary since the initially linearly independent vectors $\bm{\Phi}_{1}^{-}$ and $\bm{\Phi}_{2}^{-}$ (for $R_{2}<0$), and $\bm{\Phi}_{1}^{+}$ and $\bm{\Phi}_{2}^{+}$ (for $R_{2}>0$), become increasingly collinear as the numerical integration progresses, leading to a loss of accuracy.

Following the above procedure, one shoots all the way upto the particle location, $R_{2}=0$, to obtain:

[TABLE]

One may now use the jump conditions (66a-d), which yields:

[TABLE]

where

[TABLE]

The matrix $\bm{C}$ can be inverted to obtain the constants $c_{1}^{+}$, $c_{2}^{+}$, $c_{1}^{-}$ and $c_{2}^{-}$, which can then be used to write:

[TABLE]

Either $\hat{\langle U_{2}\rangle}^{+}(k_{1},0^{+},k_{3})$ or $\hat{\langle U_{2}\rangle}^{-}(k_{1},0^{-},k_{3})$ can be used in (67) to obtain the spheroid migration velocity.