2. Particle-free flow

We begin by computing the laminar particle-free velocity field in a two-dimensional wavy channel that extends infinitely in the out-of-plane direction. Consider fluid flow between two symmetric wavy plates with an average separation of $l$, as shown in figure 3. The wavy plates have a sinusoidal profile with amplitude $\epsilon _{w}h$ and wavelength $\lambda$, where $h$ is the mean half-width of the channel, $h=l/2$. We denote the angular frequency of the waviness as $\omega ^{\prime }$, i.e. $\lambda =2{\rm \pi} /\omega ^{\prime }$; additionally, we assume $\lambda$ is at least of $O(h)$. The laboratory frame is chosen with an origin located at the intersection between the middle plane and a cross-section on which the channel is of minimum width. We use the coordinate $\boldsymbol {\mathcal {R}}^{\prime }=(\mathcal {X}^{\prime },\mathcal {Y}^{\prime },\mathcal {Z}^{\prime })$ for the laboratory frame. The flow is driven by a pressure gradient in the $\mathcal {X}^{\prime }$ direction, and the velocity on the middle plane has an average magnitude of $U_{m}^{\prime }$.

Figure 3. Schematics of the flow in a wavy channel and two sets of coordinate systems. The solid and dashed velocity profiles represent the actual flow at $\mathcal {X}^{\prime }=0$ and the zeroth-order plane Poiseuille flow, respectively.

We non-dimensionalize our variables by normalizing the particle-free velocity $\bar {\boldsymbol {\mathcal {V}}}^{\prime }$ with $U_{m}^{\prime }$ and by introducing $\boldsymbol {\mathcal {R}}=\boldsymbol {\mathcal {R}}^{\prime }/h=(\mathcal {X},\mathcal {Y},\mathcal {Z})$. In the rest of the paper, dimensionless variables are indicated by a lack of primes. The particle-free fluid pressure $\bar {\mathcal {P}}^{\prime }$ is normalized with $\rho U_m^{\prime 2}$, where $\rho$ is the fluid density. The dimensionless governing equations of $\bar {\boldsymbol {\mathcal {V}}}$ can then be written as the dimensionless Navier–Stokes equations and the continuity condition

(2.1a)$$\begin{gather} \bar{\boldsymbol{\mathcal{V}}} \boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\boldsymbol{\mathcal{V}}} ={-}\boldsymbol{\nabla} \bar{\mathcal{P}} + \frac{2}{R_c} \nabla^2\bar{\boldsymbol{\mathcal{V}}}, \end{gather}$$

(2.1b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \bar{\boldsymbol{\mathcal{V}}} = 0, \end{gather}$$

where $\bar {\mathcal {P}}$ is the dimensionless particle-free fluid pressure, $R_c=2U_{m}^{\prime } h/\nu$ is the channel Reynolds number (Schonberg & Hinch Reference Schonberg and Hinch1989; Hogg Reference Hogg1994; Asmolov Reference Asmolov1999) and $\nu$ is the kinematic viscosity of the fluid. The velocity $\bar {\boldsymbol {\mathcal {V}}}$ is subject to a no-slip boundary condition on the walls

(2.2)\begin{equation} \bar{\boldsymbol{\mathcal{V}}} = 0 \quad \mbox{at} \ \mathcal{Z}={\pm}{1} \mp \epsilon_{w} \,{\rm Re}[\mbox{exp}(\mbox{i}\omega\mathcal{X})], \end{equation}

where $\omega = \omega ^{\prime }h$. For a small waviness amplitude, i.e. $\epsilon _{w} \ll 1$, the particle-free velocity and pressure can be expressed as a series of higher-order corrections to Poiseuille flow

(2.3a,b)\begin{equation} \bar{\boldsymbol{\mathcal{V}}} = \bar{\boldsymbol{\mathcal{V}}}^{(0)} + \epsilon_{w}\bar{\boldsymbol{\mathcal{V}}}^{(1)} + O(\epsilon_{w}^{2}) \quad \mbox{and} \quad \bar{\mathcal{P}} = \bar{\mathcal{P}}^{(0)} + \epsilon_{w}\bar{\mathcal{P}}^{(1)} + O(\epsilon_{w}^{2}), \end{equation}

where the $\mathcal {X}$-component of the zeroth-order particle-free velocity is $\bar {\mathcal {V}}_{x}^{(0)}=1-\mathcal {Z}^2$, the $\mathcal {Z}$-component of the zeroth-order particle-free velocity is $\bar {\mathcal {V}}_{z}^{(0)}=0$ and the zeroth-order particle-free pressure $\bar {\mathcal {P}}^{(0)}$ varies linearly in the $\mathcal {X}$ direction. Note that this expansion is distinct from the standard lubrication approach in which $\epsilon _w h/\lambda \ll 1$. Here, we instead allow for the case where the amplitude of the waviness may be of the same order as the wavelength.

To solve the first-order particle-free flow, we employ a strategy introduced by Van Dyke in 1975; instead of enforcing the no-slip condition at the wavy wall, we transfer the boundary condition in (2.2) to the average locations of the walls (Van Dyke Reference Van Dyke1975; Selvarajan, Tulapurkara & Vasanta Ram Reference Selvarajan, Tulapurkara and Vasanta Ram1998, Reference Selvarajan, Tulapurkara and Vasanta Ram1999) such that the no-slip condition is enforced at the wall to first order. To this end, we write (2.2) as a Taylor series around $\mathcal {Z}=1$ and then expand $\bar {\mathcal {V}}$ and its derivatives in orders of $\epsilon _w$ with (2.3a,b), finally arriving at

(2.4a)$$\begin{gather} \bar{\boldsymbol{\mathcal{V}}}_{x}^{(1)} \simeq{\pm} [\bar{\boldsymbol{\mathcal{V}}}_{x}^{(0)}]^{\prime}{\rm Re}[\mbox{exp}(\mbox{i}\omega\mathcal{X})] \quad \mbox{and} \end{gather}$$

(2.4b)$$\begin{gather}\bar{\boldsymbol{\mathcal{V}}}_{z}^{(1)} \simeq{\pm} [\bar{\boldsymbol{\mathcal{V}}}_{z}^{(0)}]^{\prime}{\rm Re}[\mbox{exp}(\mbox{i}\omega\mathcal{X})] \quad \mbox{at}\ \mathcal{Z}=\pm1, \end{gather}$$

where the prime symbols denote derivatives with respect to $\mathcal {Z}$.

We now return to (2.1) to solve for the first-order velocity $\bar {\boldsymbol {\mathcal {V}}}^{(1)}$ and pressure $\bar {\boldsymbol {\mathcal {P}}}^{(1)}$. Inserting (2.3a,b) into (2.1) gives rise to the first-order governing equations, which are then solved with the transferred boundary condition (2.4). The first-order problem is linear in $\bar {\boldsymbol {\mathcal {V}}}^{(1)}$ and $\bar {\boldsymbol {\mathcal {P}}}^{(1)}$; additionally, $\mathcal {X}$ and $\mathcal {Z}$ are separated in the transferred boundary conditions (2.4). Thus, we perform separation of variables and make the ansatz $\bar {\boldsymbol {\mathcal {V}}}^{(1)}={\rm Re}[\hat {\boldsymbol {\mathcal {V}}}(\mathcal {Z}) \exp (\mbox {i}\omega \mathcal {X})]$ and $\bar {\mathcal {P}}^{(1)}={\rm Re}[\hat {\mathcal {P}}(\mathcal {Z}) \exp (\mbox {i}\omega \mathcal {X})]$ following the structure of (2.4). Herein, $\hat {\boldsymbol {\mathcal {V}}}(\mathcal {Z})$ and $\hat {\mathcal {P}}(\mathcal {Z})$ are assumed to be complex functions, and the Re symbols ensure that the computed $\bar {\boldsymbol {\mathcal {V}}}^{(1)}$ and $\bar {\boldsymbol {\mathcal {P}}}^{(1)}$ are real-valued functions. The first-order variables $\hat {\boldsymbol {\mathcal {V}}}(\mathcal {Z})$ and $\hat {\mathcal {P}}(\mathcal {Z})$ are then subject to the following differential equations and boundary conditions:

(2.5a)$$\begin{gather} \mbox{i}\omega \hat{\mathcal{V}}_{x} + \frac{ \mbox{d}\hat{\mathcal{V}}_{z}}{\mbox{d}\mathcal{Z}} = 0, \end{gather}$$

(2.5b)$$\begin{gather}\mbox{i}\omega (1-\mathcal{Z}^2)\hat{\mathcal{V}}_{x} - 2\mathcal{Z}\hat{\mathcal{V}}_{z} ={-}\mbox{i}\omega \hat{\mathcal{P}} + \frac{2}{R_c} \left(-\omega^2 + \frac{\mbox{d}}{\mbox{d}\mathcal{Z}^2} \right)\hat{\mathcal{V}}_{x}, \end{gather}$$

(2.5c)$$\begin{gather}\mbox{i}\omega (1-\mathcal{Z}^2)\hat{\mathcal{V}}_{z} ={-}\frac{\mbox{d}\hat{\mathcal{P}}}{\mbox{d}\mathcal{Z}} + \frac{2}{R_c} \left(-\omega^2 + \frac{\mbox{d}}{\mbox{d}\mathcal{Z}^2} \right)\hat{\mathcal{V}}_{z}, \end{gather}$$

(2.5d)$$\begin{gather}\hat{\mathcal{V}}_{x} ={-}2 \quad \mbox{and} \quad \hat{\mathcal{V}}_{z} = 0 \quad \mbox{at} \ \mathcal{Z}={\pm} 1. \end{gather}$$

Note that we drop the Re symbols and the sinusoidal term $\exp (\mathrm {i}\omega \mathcal {X})$ in obtaining (2.5). Equation (2.5) can be solved numerically with the orthonormalization method (Conte Reference Conte1966; Scott & Watts Reference Scott and Watts1977). To avoid numerical instability, we perform integration from both boundaries and match the results at $\mathcal {Z}=0$. In short, we solve (2.5) on two intervals, $\mathcal {Z}\in [-1,0]$ and $\mathcal {Z}\in [0,1]$. For each interval, we write the solution as the summation of a particular solution and a linear combination of orthonormal homogeneous solutions, where the pre-factors for the homogeneous solutions are the unknowns. We integrate (2.5) from both boundaries $\mathcal {Z}=\pm 1$, and finally determine the pre-factors by utilizing the continuity of the functional values at $\mathcal {Z}=0$.

Figure 4 shows the particle-free velocity field in a wavy channel. Unlike the zeroth-order flow which is independent of the $\mathcal {X}$-coordinate (figure 4a), the first-order velocity varies periodically along the $\mathcal {X}$-direction (figure 4b). Due to the waviness, a pair of counter-rotating vortices appear in both the $\mathcal {X}$- and $\mathcal {Z}$-directions in each period of the channel. To better illustrate the distributions of $\bar {\mathcal {V}}_{x}^{(1)}$ and $\bar {\mathcal {V}}_{z}^{(1)}$ across the channel, we select the cross-section $\mathcal {X}=0$ and show $\bar {\mathcal {V}}_{x}^{(1)}(\mathcal {X}=0)$ and $\bar {\mathcal {V}}_{z}^{(1)}(\mathcal {X}=0)$ as functions of $\mathcal {Z}$ at varying channel Reynolds number $R_c$ (figure 4c) and frequency $\omega$ (figure 4d). The first-order $x$-velocity correction, $\bar {\mathcal {V}}_{x}^{(1)}(\mathcal {X}=0,\mathcal {Z})$ is symmetric about $\mathcal {Z}=0$, whereas $\bar {\mathcal {V}}_{z}^{(1)}(\mathcal {X}=0,\mathcal {Z})$ is anti-symmetric. As $R_c$ or $\omega$ increases, $\bar {\mathcal {V}}_{x}^{(1)}$ becomes more uniform, and the magnitude in the middle of the channel decreases; by contrast, $\bar {\mathcal {V}}_{x}^{(1)}$ varies more sharply close to the boundaries. The magnitude of $\bar {\mathcal {V}}_{z}^{(1)}$ does not vary monotonically with $R_c$ (figure 4c). By contrast, the magnitude of $\bar {\mathcal {V}}_{z}^{(1)}$ increases with $\omega$ (figure 4d). The locations of the maximum and minimum $\bar {\mathcal {V}}_{z}^{(1)}$ move toward the boundary of the channel as $\omega$ increases. Figures 4(c) and 4(d) indicate that the shear rate close to the boundaries is positively related to $R_c$ and $\omega$.

Figure 4. Particle-free flow in the wavy channel. The flow field can be approximated as the sum of (a) a zeroth-order Poiseuille flow $\bar {\boldsymbol {\mathcal {V}}}^{(0)}$ and (b) a first-order flow $\bar {\boldsymbol {\mathcal {V}}}^{(1)}$ which introduces eddies. In (b), $R_c=10$ and $\omega =2$. (c) Variations of the first-order velocity components $\bar {\mathcal {V}}_{x}^{(1)}$ and $\bar {\mathcal {V}}_{z}^{(1)}$ across the channel at $\mathcal {X}=0$ for varying channel Reynolds numbers $R_c$, where $\omega =2$. (d) Variations of the first-order velocity components $\bar {\mathcal {V}}_{x}^{(1)}$ and $\bar {\mathcal {V}}_{z}^{(1)}$ across the channel at $\mathcal {X}=0$ for varying frequencies $\omega$ of the waviness, where $R_c=250$. Dashed lines in (c,d) represent analytic solutions in the limits $R_c\rightarrow 0$ and $\omega \rightarrow 0$, respectively.

Figure 4 also includes the analytic solutions at the asymptotic limits $R_c\rightarrow 0$ and $\omega \rightarrow 0$ in 4(c) and 4(d), respectively. In the Stokes-flow limit, the first-order velocity component $\bar {\mathcal {V}}_{x}^{( 1 )}(\mathcal {X} =0,\mathcal {Z} )$ is given by

(2.6)\begin{align} \bar{\mathcal{V}}_{x}^{( 1 )}( \mathcal{X} =0,\mathcal{Z} ) =\frac{4\cosh ( \omega \mathcal{Z} ) \left[ -\omega \cosh ( \omega ) +\sinh ( \omega ) \right] +4\omega \mathcal{Z} \sinh ( \omega ) \sinh ( \omega \mathcal{Z} )}{2\omega -\sinh ( 2\omega )}, \end{align}

and $\bar {\mathcal {V}}_{z}^{( 1 )}(\mathcal {X} =0,\mathcal {Z} )$ is calculated to be zero. As $R_c$ decreases, our numerical computations asymptote to these analytic solutions, as shown in figure 4(c). In the limit of an infinitesimal waviness frequency, i.e. $\omega \rightarrow 0$, our numerical computations also asymptote to the analytic solutions (figure 4d), where

(2.7a,b)\begin{equation} \bar{\mathcal{V}}_{x}^{( 1 )}( \mathcal{X} =0,\mathcal{Z} ) = 1-3\mathcal{Z} ^2 \quad\mathrm{and}\quad\bar{\mathcal{V}}_{z}^{( 1 )}( \mathcal{X} =0,\mathcal{Z} ) =0. \end{equation}

The detailed derivations for the Stokes-flow limit $R_c \rightarrow 0$ and the lubrication limit $\omega \rightarrow 0$ are reported in Appendices A and B.

3. Particle perturbation analysis formulation

We now introduce a spherical particle that gives rise to a perturbation to our particle-free channel flow. Hence, we investigate the perturbation in three dimensions. As shown in figure 3, the particle is located at (dimensional) coordinates $(\mathcal {X}_{p}^{\prime },\mathcal {Y}_{p}^{\prime },\mathcal {Z}_{p}^{\prime })$ and is moving with a (dimensional) translational velocity $\boldsymbol U_{p}^{\prime }=U_{p,x}^{\prime } \boldsymbol {e}_{x}+U_{p,z}^{\prime } \boldsymbol {e}_{z}$, where $\boldsymbol {e}_x$ and $\boldsymbol {e}_z$ are unit vectors in the $x$- and $z$-directions, respectively. We define $\varphi _p=\omega ^{\prime }\mathcal {X}_p^{\prime }$ as the phase angle corresponding to the $\mathcal {X}$-location of the particle (figure 3). To understand the perturbation induced by the particle, we choose a non-inertial frame of reference that is instantaneously located at the particle centre and translates at the particle velocity $\boldsymbol {U}_{p}^{\prime }$ (Hogg Reference Hogg1994; Asmolov Reference Asmolov1999; Hood et al. Reference Hood, Lee and Roper2015). Since the particle is free, the net force (including the virtual force due to the acceleration of the frame) acting on the particle is zero in this frame of reference. We use the coordinate $\boldsymbol {r}^{\prime }=(x^{\prime },y^{\prime },z^{\prime })$ for the particle frame (figure 3). Let $\bar {\boldsymbol {v}}^{\prime }$ denote the particle-free velocity and $\bar {p}^{\prime }$ the particle-free pressure in the particle frame; namely, the velocity and pressure fields computed in the previous section are shifted by the (currently unknown) velocity $\boldsymbol {U}_p^{\prime }$. The particle-free equations of motion can be written in terms of $\bar {\boldsymbol {v}}^{\prime }$ as

(3.1a)$$\begin{gather} \rho \left( \frac{\partial \boldsymbol{\bar{v}}^{\prime}}{\partial t^{\prime}}+\boldsymbol{\bar{v}}^{\prime}\boldsymbol{\cdot} \boldsymbol{\nabla} \mathrm{'}\boldsymbol{\bar{v}}^{\prime}+\boldsymbol{a}_{p}^{\prime} \right) ={-}\boldsymbol{\nabla} '\bar{p}^{\prime}+\mu \nabla'^2\boldsymbol{\bar{v}}^{\prime}, \end{gather}$$

(3.1b)$$\begin{gather}\boldsymbol{\nabla} '\boldsymbol{\cdot} \boldsymbol{\bar{v}}^{\prime} = 0, \end{gather}$$

(3.1c)$$\begin{gather}\boldsymbol{\bar{v}}^{\prime} ={-}\boldsymbol{U}_{p}^{\prime} \quad \mbox{on the walls}, \end{gather}$$

where $\boldsymbol {a}_{p}^{\prime }$ is the particle acceleration measured in the laboratory frame (Kundu, Cohen & Dowling Reference Kundu, Cohen and Dowling2016). Although (3.1) has been solved in the previous section, we rewrite it here in the particle frame as they will be used below to derive the governing equations for the difference between the disturbed and particle-free velocity fields.

We now consider the equations of motion for the flow field which is perturbed by the presence of the particle. Let $\boldsymbol {v}'$ (note no overbar) denote the disturbed fluid velocity and $p'$ the disturbed pressure in the particle frame. The equations for the disturbed flow take the same form of the Navier–Stokes equations and the continuity condition as (3.1)

(3.2a)$$\begin{gather} \rho \left( \frac{\partial \boldsymbol{v}^{\prime}}{\partial t^{\prime}}+\boldsymbol{v}^{\prime}\boldsymbol{\cdot} \boldsymbol{\nabla} {'}\boldsymbol{v}^{\prime}+\boldsymbol{a}_{p}^{\prime} \right) ={-}\boldsymbol{\nabla} {'}p^{\prime}+\mu \nabla {'}^2\boldsymbol{v}^{\prime}, \end{gather}$$

(3.2b)$$\begin{gather}\boldsymbol{\nabla} '\boldsymbol{\cdot} \boldsymbol{v}^{\prime}=0. \end{gather}$$

The boundary conditions include a no-slip boundary condition on the walls, a no-slip boundary condition on the surface of the particle and the absence of perturbation far away from the particle

(3.3a)$$\begin{gather} \boldsymbol{v}^{\prime} ={-}\boldsymbol{U}_{p}^{\prime}\quad \mbox{on the walls}, \end{gather}$$

(3.3b)$$\begin{gather}\boldsymbol{v}^{\prime} = \boldsymbol{\varOmega }_{p}^{\prime}\times \boldsymbol{r}^{\prime}\boldsymbol \quad \mbox{at}\boldsymbol\ \left| \boldsymbol{r}^{\prime} \right|=a, \end{gather}$$

(3.3c)$$\begin{gather}\boldsymbol{v}^{\prime} = \boldsymbol{\bar{v}}^{\prime}\quad\mbox{at} \left| x^{\prime} \right|\rightarrow \infty, \end{gather}$$

where $a$ is the particle radius, and $\boldsymbol {\varOmega }_p^{\prime }$ is the particle's rotational velocity. Next, we define the perturbation velocity $\boldsymbol {u}^{\prime }$ as the difference between the disturbed and the particle-free velocities, i.e. $\boldsymbol {u}^{\prime }=\boldsymbol {v}{'}-\boldsymbol {\bar {v}}{'}$. The governing equations of $\boldsymbol {u}^{\prime }$ can be obtained by subtracting (3.1) from (3.2)

(3.4a)$$\begin{gather} \frac{\partial \boldsymbol{u}^{\prime}}{\partial t^{\prime}}+\boldsymbol{u}^{\prime}\boldsymbol{\cdot} \boldsymbol{\nabla} {'}\boldsymbol{u}^{\prime}+\boldsymbol{\bar{v}}^{\prime}\boldsymbol{\cdot} \boldsymbol{\nabla} {'}\boldsymbol{u}^{\prime}+\boldsymbol{u}^{\prime}\boldsymbol{\cdot} \boldsymbol{\nabla} {'}\boldsymbol{\bar{v}}^{\prime} ={-}\frac{1}{\rho}\boldsymbol{\nabla} {'}(p^{\prime}-\bar{p}^{\prime}) +\nu \nabla{'}^2\boldsymbol{u}^{\prime}, \end{gather}$$

(3.4b)$$\begin{gather}\boldsymbol{\nabla} {'}\boldsymbol{\cdot} \boldsymbol{u}^{\prime} = 0. \end{gather}$$

Note that the particle acceleration $\boldsymbol {a}_p^{\prime }$ in (3.1) and (3.2) cancels out during subtraction. The boundary conditions of (3.4) are

(3.5a) $$\begin{gather} \boldsymbol{u}^{\prime}=0 \quad \mbox{at} \ z^{\prime}=\left\{\begin{array}{@{}l}-\mathcal{Z}_p^{\prime}-h+\epsilon_{w}h\,{\rm Re}\{ \exp [ \mathrm{i}(\omega^{\prime} x^{\prime}+\varphi_p ) ] \} \\ h-\mathcal{Z}_p^{\prime}-\epsilon_{w}h\,{\rm Re} \{ \exp [ \mathrm{i}(\omega^{\prime} x^{\prime}+\varphi_p ) ]\} \\ \end{array},\right. \end{gather}$$

(3.5b)$$\begin{gather}\boldsymbol{u}^{\prime}=\boldsymbol{\varOmega }_{p}^{\prime}\times \boldsymbol{r}^{\prime}-\boldsymbol{\bar{v}}^{\prime}\quad\mbox{at}\ \left| \boldsymbol{r}^{\prime} \right|=a, \end{gather}$$

(3.5c)$$\begin{gather}\boldsymbol{u}^{\prime}=\boldsymbol{0}\quad \mbox{at}\ \left| x^{\prime} \right|\rightarrow \infty. \end{gather}$$

We non-dimensionalize (3.5) by introducing

(3.6)\begin{equation} \left.\begin{gathered} \boldsymbol{r}=\frac{\boldsymbol{r}^{\prime}}{a},\quad \boldsymbol{u}=\frac{\boldsymbol{u}^{\prime}}{U_{m}^{\prime}},\quad \bar{\boldsymbol{v}}=\frac{\bar{\boldsymbol{v}}^{\prime}}{U_{m}^{\prime}},\quad \omega =\omega ^{\prime}h,\\ \boldsymbol{\varOmega }_p=\frac{2h\boldsymbol{\varOmega }_{p}^{\prime}}{U_{m}^{\prime}},\quad \mbox{and}\quad q=\left( \frac{\mu U_{m}^{\prime}}{a} \right) ^{{-}1}\left( p^{\prime}-\bar{p}^{\prime} \right), \end{gathered}\right\}\end{equation}

where $\bar {\boldsymbol {v}}$ is the dimensionless particle-free velocity in the particle frame, $q$ is the dimensionless perturbation pressure and a lack of prime denotes dimensionless variables. Additionally, we introduce a particle Reynolds number $R_p=U_{m}^{\prime }a{{^2}/{( l\nu )}}$ (Schonberg & Hinch Reference Schonberg and Hinch1989; Hogg Reference Hogg1994; Asmolov Reference Asmolov1999), where $\nu$ is the kinematic viscosity of the fluid. We define $\epsilon =R_{p}^{{{1}/{2}}}=R_{c}^{{{1}/{2}}}{{a}/{l}}$ (Asmolov Reference Asmolov1999; Matas et al. Reference Matas, Morris and Guazzelli2009). This small dimensionless parameter $\epsilon$ is taken to be the asymptotic parameter that we will use to analyse the perturbation introduced by the particle, and we assume $\epsilon \ll 1$. Hence, we have the following dimensionless governing equations of $\boldsymbol {u}$:

(3.7a)$$\begin{gather} {\epsilon R}_{c}^{{{1}/{2}}}\left( \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}+\boldsymbol{\bar{v}}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\bar{v}} \right) ={-}\boldsymbol{\nabla} q+\nabla ^2\boldsymbol{u}, \end{gather}$$

(3.7b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$

with boundary conditions

(3.8a)$$\begin{gather} \boldsymbol{u}= \boldsymbol{0}\quad \mbox{at}\ z = \tfrac{1}{2}\epsilon ^{{-}1}R_{c}^{{{1}/{2}}}({\mp} 1 - \mathcal{Z}_p \pm \epsilon_w \,{\rm Re}\{ \exp [ \mathrm{i}(\varphi_p +x\epsilon \kappa )] \}), \end{gather}$$

(3.8b)$$\begin{gather}\boldsymbol{u}= {\epsilon R}_{c}^{-{{1}/{2}}}\boldsymbol{\varOmega }_p\times \boldsymbol{r}-\boldsymbol{\bar{v}}\quad\mbox{at}\ \left| \boldsymbol{r} \right|=1, \end{gather}$$

(3.8c)$$\begin{gather}\boldsymbol{u}=0\quad\mbox{at}\ \left| x \right|\rightarrow \infty, \end{gather}$$

where we replace $2\omega {R_c^{-{1}/{2}}}$ with $\kappa$, and $\bar {\boldsymbol {v}}$ is the non-dimensional particle-free flow field in the frame of the particle. Additionally, we factor out $\epsilon$ in the exponential term of (3.8a) for mathematical convenience in § 4.2.

Finally, we express the particle-frame-based particle-free velocity $\boldsymbol {\bar {v}}$ in (3.7) using quantities that have been calculated in § 2. In terms of the dimensionless particle-free velocity $\boldsymbol {\bar {\mathcal {V}}}$ in the laboratory frame, and the dimensionless particle velocity $\boldsymbol {U}_p$ we have

(3.9)\begin{equation} \boldsymbol{\bar{v}}=\boldsymbol{\bar{\mathcal{V}}}-\boldsymbol{U}_p. \end{equation}

We further rewrite the arguments of $\boldsymbol {\bar {\mathcal {V}}}$, $\boldsymbol {\bar {\mathcal {V}}}^{(0)}$, $\boldsymbol {\bar {\mathcal {V}}}^{(1)}$ and $\hat {\boldsymbol {\mathcal {V}}}$ in terms of the dimensionless coordinates in the particle frame $x$, $y$ and $z$ with the relationships $\mathcal {X}=ax/h+\mathcal {X}_p$ and $\mathcal {Z}=az/h+\mathcal {Z}_p$, e.g. $\hat {\mathcal {V}}_x( \mathcal {Z} ) =\hat {\mathcal {V}}_x( 2\epsilon zR_{c}^{-{{1}/{2}}}+\mathcal {Z}_p )$. Note that a conversion between $\mathcal {Y}$ and $y$ is not necessary because the $y$-components of both $\bar {\boldsymbol {\mathcal {V}}}$ and $\boldsymbol {U}_{p}$ are equal to zero, leading to $\bar {v}_y = 0$ (3.9); similarly, all $y$-components in the rest of this section vanish. Moreover, since the perturbation velocity $\boldsymbol {u}$ is influenced by the lag between the particle velocity $\boldsymbol {U}_p$ and the particle-free velocity at the instantaneous location of the particle centre $\boldsymbol {\bar {\mathcal {V}}}( \boldsymbol {\mathcal {R}} _p( t ) )$, we define a dimensionless particle slip velocity, $\boldsymbol {V}_s=\boldsymbol {\bar {\mathcal {V}}}( \boldsymbol {\mathcal {R}} _p( t ) ) -\boldsymbol {U}_p$. Inserting (2.3a,b), the asymptotic expansions for $\bar {\boldsymbol {\mathcal {V}}}$ and $\bar {\mathcal {P}}$ into (3.9), we find

(3.10a,b)\begin{equation} \bar{v}_x=V_{s,x}+\Delta \bar{\mathcal{V}}_{x}^{(0)}+\epsilon _w\Delta \bar{\mathcal{V}}_{x}^{(1)}+O ( \epsilon _{w}^{2} ) ,\quad \bar{v}_z=V_{s,z}+\epsilon _w\Delta \bar{\mathcal{V}}_{z}^{(1)}+O ( \epsilon _{w}^{2} ), \end{equation}

where $\Delta \bar {\mathcal {V}}_{x}^{(i)}$ represents the difference between the $i$th order particle-free velocity at location $\boldsymbol {r}=( x,y,z )$ and the $i$th-order particle-free velocity at the particle centre $(0,0,0)$. Here, $\Delta \bar {\mathcal {V}}_{x}^{(0)}=\epsilon \gamma zR_{c}^{-{{1}/{2}}}-4\epsilon ^2z^2R_{c}^{-1}$ and $\Delta \bar {\mathcal {V}}_{z}^{(0)}=0$, where $\gamma =-4\mathcal {Z}_p$ is the zeroth-order shear rate at the particle centre normalized with a characteristic shear rate, $U_{m}^{\prime }/{l}$. Additionally,

(3.11a)$$\begin{gather} \Delta \bar{\mathcal{V}}_{x}^{(1)}={\rm Re}\{ \hat{\mathcal{V}}_x( 2\epsilon z{R}_{c}^{-{1}/{2}}+ \mathcal{Z}_p ) \exp [\mathrm{i}( \varphi_p+\epsilon x\kappa )] -\hat{\mathcal{V}}_x( \mathcal{Z}_p ) \exp ( \mathrm{i}\varphi_p ) \}, \end{gather}$$

(3.11b)$$\begin{gather}\Delta \bar{\mathcal{V}}_{z}^{(1)}={\rm Re}\{ \hat{\mathcal{V}}_z( 2\epsilon z{R}_{c}^{-{1}/{2}}+\mathcal{Z}_p ) \exp [\mathrm{i}( \varphi_p+\epsilon x\kappa )] -\hat{\mathcal{V}}_z( \mathcal{Z}_p ) \exp ( \mathrm{i}\varphi_p ) \}. \end{gather}$$

To determine the orders of $\Delta \bar {\mathcal {V}}_{x}^{(1)}$ and $\Delta \bar {\mathcal {V}}_{z}^{(1)}$ in $\epsilon$ for later analyses in § 4, we first write the Taylor series expansions of $\hat {\mathcal {V}}_x( 2{\epsilon z}R_{c}^{-{{1}/{2}}}+\mathcal {Z}_p )$ and $\hat {\mathcal {V}}_z ( 2{\epsilon z}R_{c}^{-{{1}/{2}}}+\mathcal {Z}_p )$ at $z=0$ under the condition $2{\epsilon z}R_{c}^{-{{1}/{2}}} \ll 1$

(3.12a)$$\begin{gather} \hat{\mathcal{V}}_x( 2{\epsilon z}R_{c}^{-{1}/{2}}+\mathcal{Z}_p ) =\hat{\mathcal{V}}_{x0}+\alpha _x{\epsilon z}R_{c}^{-{1}/{2}}+\tfrac{1}{2}\beta _x\epsilon ^2z^2{R}_c^{{-}1}+O ( \epsilon ^3 ), \end{gather}$$

(3.12b)$$\begin{gather}\hat{\mathcal{V}}_z( 2{\epsilon z}R_{c}^{-{1}/{2}}+\mathcal{Z}_p ) =\hat{\mathcal{V}}_{z0}+\alpha _z{\epsilon z}R_{c}^{-{1}/{2}}+\tfrac{1}{2}\beta _z\epsilon ^2z^2{R}_c^{{-}1}+O ( \epsilon ^3 ), \end{gather}$$

where $\hat {\mathcal {V}}_{x0}=\hat {\mathcal {V}}_x( \mathcal {Z}_p )$, $\hat {\mathcal {V}}_{z0}=\hat {\mathcal {V}}_z( \mathcal {Z}_p )$ and $\alpha _x$, $\alpha _z$, $\beta _x$ and $\beta _z$ are prefactors of the Taylor series, which can be determined numerically from the solutions to (2.5). Specifically, the continuity condition gives rise to the relation ${\mathrm {i}\kappa }R_c^{{{1}/{2}}}\hat {\mathcal {V}}_{x0}+\alpha _z=0$. Given (3.12), we further expand $\Delta \bar {\mathcal {V}}_{x}^{(1)}$ and $\Delta \bar {\mathcal {V}}_{z}^{(1)}$ in Taylor series with respect to both $x$ and $z$, and find that $\Delta \bar {\mathcal {V}}_{x}^{(1)}$ and $\Delta \bar {\mathcal {V}}_{z}^{(1)}$ are of $O ( \epsilon )$, e.g. $\Delta \bar {\mathcal {V}}_{x}^{( 1 )}=\mathrm {Re}[ ( \mathrm {i}x\kappa \hat {\mathcal {V}} _{x0}+\alpha _{x} z R_{c}^{-1/2} ) \exp \mathrm {(i}\varphi _p) ] \epsilon +O(\epsilon ^2)$. These expressions will be used to express the boundary conditions order by order at the surface of a particle in the following section. Furthermore, we can analytically obtain the variables in (3.11) and (3.12) in the lubrication limit where the angular frequency $\omega \rightarrow 0$ (Appendix B)

(3.13a,b)\begin{equation} \Delta \bar{\mathcal{V}}_{x}^{( 1 )}\rightarrow -12z\epsilon (\mathcal{Z} _pR_{c}^{-{{1}/{2}}}+z\epsilon R_{c}^{{-}1})\cos(\varphi _p) \quad \textrm{and} \quad \Delta \bar{\mathcal{V}}_{z}^{(1)}\rightarrow 0, \end{equation}

which are of $O(\epsilon )$. We further find in this limit:$\hat {\mathcal {V}}_{x0}\rightarrow 1-3\mathcal {Z} _{p}^{2}$, $\hat {\mathcal {V}}_{z0}\rightarrow 0$, $\alpha _x \rightarrow -12\mathcal {Z} _p$, $\alpha _z \rightarrow 0$, $\beta _x\rightarrow -24$ and $\beta _z\rightarrow 0$.

4. Inertial lift force acting on a neutrally buoyant particle

The particle induces a small perturbation to the steady channel flow and hence experiences an inertial lift force in the lateral direction. In this section, we calculate the inertial lift force on a neutrally buoyant particle (recall that we seek to determine whether particle focusing in wavy channels arises from this lift force). The problem formulated in (3.7) is solved with the method of matched asymptotic expansions (Schonberg & Hinch Reference Schonberg and Hinch1989; Hogg Reference Hogg1994; Asmolov Reference Asmolov1999). We first solve an inner problem to obtain the flow field surrounding the particle. We then regard the particle as a singular point that exerts a force on the channel flow in an outer problem. Finally, matching the inner and outer solutions gives rise to a lift force on the particle.

4.1. Inner solution

We introduce an inner expansion of the perturbation velocity and pressure

(4.1a,b)\begin{equation} \boldsymbol{u}=\boldsymbol{u}_0+\epsilon \boldsymbol{u}_1+O ( \epsilon ^2 ) ,\quad q=q_0+\epsilon q_1+O ( \epsilon ^2 ). \end{equation}

Inserting (4.1a,b) into (3.7), the zeroth-order governing equations represent a Stokes-flow problem

(4.2a)$$\begin{gather} \nabla ^2\boldsymbol{u}_0-\boldsymbol{\nabla} q_0=\boldsymbol{0}, \end{gather}$$

(4.2b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_0=0, \end{gather}$$

(4.2c)$$\begin{gather}\boldsymbol{u}_0={-}V_{s,x}\boldsymbol{e}_x-V_{s,z}\boldsymbol{e}_z,\quad \left| \boldsymbol{r} \right|=1, \end{gather}$$

(4.2d)$$\begin{gather}\boldsymbol{u}_0 \rightarrow \boldsymbol{0},\quad \boldsymbol{r}\rightarrow \infty. \end{gather}$$

In the absence of fluid inertia, we apply Faxén's law and find that $V_{s,x}$ and $V_{s,z}$ are of $O(\epsilon ^2{R}_c^{-1})$. Thus, (4.2c) reduces to $\boldsymbol {u}_0=\boldsymbol {0}$ on the particle surface, and $\boldsymbol {u}_0$ vanishes everywhere in the channel. The first-order governing equations are then

(4.3a)$$\begin{gather} \nabla ^2\boldsymbol{u}_{1}-\boldsymbol{\nabla} q_{1}=\boldsymbol{0}, \end{gather}$$

(4.3b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_1=0, \end{gather}$$

(4.3c)$$\begin{gather}\boldsymbol{u}_1={R}_c^{-{{1}/{2}}}\left( \boldsymbol{\varOmega }_{p}\times \boldsymbol{r}-\gamma z\boldsymbol{e}_{x} \right) - \epsilon _w\epsilon ^{{-}1}\Delta \bar{\mathcal{V}}_{x}^{\left( 1 \right)}\boldsymbol{e}_x - \epsilon _w\epsilon ^{{-}1}\Delta \bar{\mathcal{V}}_{z}^{\left( 1 \right)}\boldsymbol{e}_z,\quad\left| \boldsymbol{r} \right|=1, \end{gather}$$

(4.3d)$$\begin{gather}\boldsymbol{u}_1\rightarrow \boldsymbol{0},\quad\boldsymbol{r}\rightarrow \infty, \end{gather}$$

where both $\epsilon ^{-1}\Delta \bar {\mathcal {V}}_{x}^{(1)}$ and $\epsilon ^{-1}\Delta \bar {\mathcal {V}}_{z}^{(1)}$ are of $O ( 1 )$ (see the end of § 3). We also assume that $\epsilon _w$ and $\epsilon$ are of the same order of magnitude such that the last two terms in (4.3c) are bounded. The first-order problem is still a Stokes-flow problem. Due to its linearity, (4.3) can be decomposed into several sub-problems, which we solve separately based on the fundamental solutions for a sphere immersed in Stokes flow derived by Guazzelli & Morris (Reference Guazzelli and Morris2011). For a rotational boundary condition with dimensionless rotation rate $\boldsymbol {\omega }$ at the particle surface, i.e. $\boldsymbol {u}_1=\boldsymbol {\omega } \times \boldsymbol {r}$ at $| \boldsymbol {r} |=1$, the perturbation velocity $\boldsymbol {u}_1$ is given by

(4.4)\begin{equation} \boldsymbol{u}_1=\boldsymbol{\omega }\times \frac{\boldsymbol{r}}{r^3}, \end{equation}

which is adapted from (2.7) in Guazzelli & Morris (Reference Guazzelli and Morris2011). For a straining boundary condition with dimensionless strain-rate tensor $\boldsymbol{\mathsf{E}}$ at the particle surface, i.e. $\boldsymbol {u}_1=\boldsymbol{\mathsf{E}}\boldsymbol {r}$ at $| \boldsymbol {r} |=1$, the components of the perturbation velocity $\boldsymbol {u}_1$ are given by

(4.5)\begin{equation} u_{1i}=\frac{5}{2}\frac{x_i\left( x_j {\mathsf{E}}_{jk}x_k \right)}{r^5}+\frac{1}{2}{\mathsf{E}}_{jk}\left[ \frac{\delta _{ij}x_k+\delta _{ik}x_j}{r^5}-\frac{5x_ix_jx_k}{r^7} \right], \end{equation}

which is adapted from (2.17) in Guazzelli & Morris (Reference Guazzelli and Morris2011). Based on these fundamental solutions, we decompose $\boldsymbol {u}_1( \boldsymbol {r} )$ into three parts: $\boldsymbol {u}_1( \boldsymbol {r} ) =\boldsymbol {u}_{1A}( \boldsymbol {r} ) +\boldsymbol {u}_{1B}( \boldsymbol {r} ) +\boldsymbol {u}_{1C}( \boldsymbol {r} )$ to satisfy the different components in boundary condition (4.3c). In (4.3c), the term ${R}_c^{-{{1}/{2}}}( \boldsymbol {\varOmega }_{p}\times \boldsymbol {r} )$ at the boundary corresponds to a pure rotational boundary condition with a dimensionless rotation rate $\boldsymbol {\omega }={R}_c^{-{1}/{2}} \boldsymbol {\varOmega }_p$, which gives rise to the solution

(4.6)\begin{equation} \boldsymbol{u}_{1A}( \boldsymbol{r} ) ={R}_c^{-{1}/{2}} \boldsymbol{\varOmega }_p\times \frac{\boldsymbol{r}}{r^3}. \end{equation}

The term $-{R}_c^{-{{1}/{2}}}\gamma z\boldsymbol {e}_{x}$ in (4.3c) represents a simple-shear boundary condition, which can be treated as the combination of a rotational boundary condition with a dimensionless rotation rate $\boldsymbol {\omega }=-{R}_c^{-{1}/{2}} (\gamma /2) \boldsymbol {e}_y$ and a straining boundary condition with a dimensionless strain-rate tensor

(4.7)\begin{equation} \boldsymbol{\mathsf{E}}={-}\tfrac{1}{2} {R}_c^{-{1}/{2}} \left[ \begin{matrix} 0 & 0 & \gamma\\ 0 & 0 & 0\\ \gamma & 0 & 0\\ \end{matrix} \right]. \end{equation}

We then construct the following solution:

(4.8)\begin{equation} \boldsymbol{u}_{1B}( \boldsymbol{r} ) ={R}_c^{-{1}/{2}}\left[ -\frac{\gamma}{2}\boldsymbol{e}_{y}\times \frac{\boldsymbol{r}}{r^3}+\frac{5}{2}\gamma \boldsymbol{r}xz\left( \frac{1}{r^7}-\frac{1}{r^5} \right) -\frac{\gamma}{2r^5}\left( \boldsymbol{e}_{x}z+\boldsymbol{e}_{z}x \right) \right]. \end{equation}

The terms $-\epsilon _w\epsilon ^{-1}\Delta \bar {\mathcal {V}}_{x}^{( 1 )}\boldsymbol {e}_{x}-\epsilon _w\epsilon ^{-1}\Delta \bar {\mathcal {V}}_{z}^{( 1 )}\boldsymbol {e}_{z}$ in (4.3c) can be treated as the combination of a rotational and a straining boundary condition. Given (3.11), we find the Taylor series for these two terms at the origin and only collect the leading-order terms in $\epsilon$. The corresponding strain-rate tensor $\boldsymbol{\mathsf{E}}$ and rotation vector $\boldsymbol {\omega }$ at the origin can be written as

(4.9)\begin{equation} \left.\begin{gathered} \boldsymbol{\mathsf{E}}={-}\epsilon _w{R}_c^{-{1}/{2}}\,{\rm Re}\!\left\{ \exp ( \mathrm{i} \varphi_p ) \left[ \begin{matrix} -\alpha _z & 0 & \dfrac{\zeta}{2}\\[4pt] 0 & 0 & 0\\ \dfrac{\zeta}{2} & 0 & \alpha _z \end{matrix} \right] \right\} \\ \mbox{and} \\ \boldsymbol{\omega }={-}\epsilon _w{R}_c^{-{1}/{2}}\,{\rm Re}\!\left\{ \exp ( \mathrm{i}\varphi_p ) \left[ \begin{array}{@{}c@{}} 0\\ \dfrac{\sigma}{2}\\[4pt] 0 \end{array} \right] \right\}, \end{gathered}\right\} \end{equation}

where $\zeta := \alpha _x+\mathrm {i}\kappa R_{c}^{1/2}\hat {\mathcal {V}}_{z0}$ and $\sigma := \alpha _x-{\mathrm {i}\kappa R}_c^{{{1}/{2}}}\hat {\mathcal {V}}_{z0}$, which reduces to $\zeta = \sigma = -12\mathcal {Z}_p$ for an infinitesimal angular frequency $\omega$. The third component of the first-order velocity field $\boldsymbol {u}_{1C}( \boldsymbol {r} )$ is thus given by

(4.10)\begin{align} \boldsymbol{u}_{1C}( \boldsymbol{r} ) &=\epsilon _w{R}_c^{-{1}/{2}}\,{\rm Re}\!\left\{ \exp( \mathrm{i}\varphi_p ) \left[ \left( -\frac{\sigma}{2}\boldsymbol{e}_{y} \right) \times \frac{\boldsymbol{r}}{r^3} +\frac{5}{2}\boldsymbol{r}\left( \zeta xz-\alpha _zx^2+\alpha _zz^2 \right) \left( \frac{1}{r^7}-\frac{1}{r^5} \right)\right.\right. \nonumber\\ &\quad \left.\left. -\frac{\zeta}{2}\left( \frac{z\boldsymbol{e}_{x}+x\boldsymbol{e}_{z}}{r^5} \right) +\alpha _z\left( \frac{x\boldsymbol{e}_{x}-z\boldsymbol{e}_{z}}{r^5} \right) \right] \right\}. \end{align}

We now determine the rotation velocity $\boldsymbol {\varOmega }_p$ by establishing the balance of moment on the particle. According to Faxén's second law, the leading-order dimensional hydrodynamic torque on the particle $\boldsymbol {T}^{\prime }$ is proportional to the relative rotation velocity between the undisturbed flow and the particle (Guazzelli & Morris Reference Guazzelli and Morris2011)

(4.11)\begin{equation} \boldsymbol{T}^{\prime}=8{\rm \pi} \mu a^3( \boldsymbol{\bar{\varOmega}}^{\prime}(\boldsymbol{r}^{\prime}=\boldsymbol{0})-\boldsymbol{\varOmega}_{p}^{\prime} ), \end{equation}

where $\boldsymbol {\bar {\varOmega }}^{\prime }( \boldsymbol {r}^{\prime }=\boldsymbol {0} ) =\boldsymbol {\nabla } ^{\prime }\times {{\bar {\boldsymbol {v}}^{\prime }}/{2}}$ is the rotation vector of the undisturbed flow evaluated at the particle centre. The torque leads to the angular acceleration of the particle, i.e. $\boldsymbol {T}^{\prime }=I_p\,\mathrm {d}\boldsymbol {\varOmega }_{p}^{\prime }/{\mathrm {d}}t^{\prime }$, where $I_p=(2/5) a^2 m_p$ is the moment of inertia of the particle. Since we assume that the wavelength of the channel $\lambda$ is at least of $O(l)$, the characteristic time scale $t^{\prime }$ for the variation of $\boldsymbol {\varOmega }_{p}^{\prime }$ is taken to be $l/U_m^{\prime }$. Further normalizing the rotation vectors with $U_m^{\prime }/l$ (3.6), we can establish the following moment balance on the neutrally buoyant particle:

(4.12)\begin{equation} \frac{1}{15}\epsilon ^2\frac{\mathrm{d}\boldsymbol{\varOmega }_p}{\mathrm{d}t}=\boldsymbol{\bar{\varOmega}}(\boldsymbol{r}=\boldsymbol{0})-\boldsymbol{\varOmega }_p. \end{equation}

Expanding $\boldsymbol {\varOmega }_p$ into orders of $\epsilon$, we find $\boldsymbol {\varOmega }_p^{(0)}=\boldsymbol {\bar {\varOmega }}(\boldsymbol {r}=\boldsymbol {0})$. Namely, in the framework of the first-order perturbation analysis in $\epsilon$, the particle follows the rotation of the undisturbed flow at the instantaneous particle centre, exerting no net hydrodynamic torque on the fluid. Therefore, the summation of all rotational components in $\boldsymbol {u}_1( \boldsymbol {r} )$ must vanish, and $\boldsymbol {\varOmega }_p^{(0)}$ equals $\{ \gamma +\epsilon _w\,\textrm {Re}[ \exp ( \mathrm {i}\varphi _p ) \sigma ] \} \boldsymbol {e}{{_y}/{2}}$. Instead of a lift force, the perturbation velocity $\boldsymbol {u}_1( \boldsymbol {r} )$ leads to a strainlet (or force dipole) on the particle. The truncated term $\boldsymbol {\varOmega }_p - \boldsymbol {\varOmega }_p^{(0)}$ generates a rotlet, but its strength is negligible compared with the strainlet. To give a sense of the perturbation velocity field associated with the strainlet, figure 5 shows the distributions of the $x$ and $z$ components of $\boldsymbol {u}_{1}(\boldsymbol {r})$ on the plane $y=0$ in the particle frame. Furthermore, at the limit $\epsilon _w \rightarrow 0$, our analytic expression of $\boldsymbol {u}_{1}(\boldsymbol {r})$ recovers the corresponding expression in a straight channel (see (5.3) in Asmolov Reference Asmolov1999).

Figure 5. Distributions of the (a) $x$-component and (b) $z$-component of the first-order perturbation velocity $\boldsymbol {u}_{1}(\boldsymbol {r})$ on the plane $y=0$, where $\boldsymbol {r}=(x,y,z)$ denotes the coordinate system in the particle frame. Arrows represent the vector field $\boldsymbol {u}_{1}(\boldsymbol {r})$, and the circle represents the particle surface. In this example, $R_c=100$, $\omega =2$, $\mathcal {Z}_p=0.5$, $\varphi _p=0.3$ and $\epsilon _w=0.05$.

4.2. Outer solution

To find the lift force, we now turn to the outer region in which the balance between the viscous and inertial effects are important. In the outer problem, we introduce outer coordinates $\boldsymbol {R}=( X,Y,Z )$ that are related to the inner variables via $\boldsymbol {R}=\epsilon \boldsymbol {r}$. The perturbation velocity and pressure fields are then expressed as $\boldsymbol {U}(\boldsymbol {R})$ and $P(\boldsymbol {R})$, respectively.

To match the inner and outer problems, we regard the particle as a singular point that exerts a strainlet on the flow field. We will introduce the strainlet as a forcing term at the origin in the momentum equations of the outer problem. To this end, we first convert $\boldsymbol {u}_1( \boldsymbol {r} )$ (see (4.6)–(4.10)) to a function of outer coordinates

(4.13) \begin{align} \boldsymbol{u}_1(\boldsymbol{R} ) =\frac{5}{2}\epsilon ^2{R}_c^{-{1}/{2}}\frac{\boldsymbol{R}}{R^5}\{ -\gamma XZ+\epsilon _w\,{\rm Re}[ \exp ( \mathrm{i}\varphi_p ) ( -\zeta XZ+\alpha _zX^2-\alpha _zZ^2)]\} +O( \epsilon ^4 ). \end{align}

Given (4.1a,b) and (4.13), we now propose an outer expansion of the perturbation field

(4.14a,b)\begin{equation} \boldsymbol{u}=\epsilon ^3\boldsymbol{U}+O ( \epsilon ^4 ) ,\quad q=\epsilon ^4Q+O ( \epsilon ^5 ). \end{equation}

For a general strain tensor $\boldsymbol{\mathsf{E}}$ and the corresponding strainlet velocity field $\boldsymbol {u}=5x_i {\mathsf{E}}_{ij}x_j{{\boldsymbol {r}}/{( 2r^5 )}}$, the singular forcing term on the right-hand side of the outer momentum equations takes the form of $20{\rm \pi} {\mathsf{E}}_{ij}\partial \delta {{( \boldsymbol {r} )}/{\partial x_j}}$ (Hogg Reference Hogg1994). Hence, the governing equations and boundary conditions for our outer problem are

(4.15a)\begin{align} R_{c}^{\frac{1}{2}}\left(\boldsymbol{\bar{v}}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U} + \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\bar{v}} \right) &=\boldsymbol{\nabla} Q+\nabla ^2\boldsymbol{U}-\frac{10}{3}{\rm \pi} R_{c}^{-{1}/{2}}\gamma \left[ \frac{\partial \delta( \boldsymbol{R} )}{\partial Z}\boldsymbol{e}_x+\frac{\partial \delta( \boldsymbol{R} )}{\partial X}\boldsymbol{e}_z \right]\nonumber\\ &\quad +\epsilon _w\frac{10}{3}{\rm \pi} {R}_c^{-{1}/{2}}\,{\rm Re}\!\left[ \exp( \mathrm{i}\varphi_p ) \left\{ -\zeta \left[ \frac{\partial \delta( \boldsymbol{R} )}{\partial Z}\boldsymbol{e}_x+\frac{\partial \delta( \boldsymbol{R} )}{\partial X}\boldsymbol{e}_z \right]\right.\right.\nonumber\\ &\quad \left.\left. +2\alpha _z\left[ \frac{\partial \delta ( \boldsymbol{R} )}{\partial X}\boldsymbol{e}_x-\frac{\partial \delta ( \boldsymbol{R} )}{\partial Z}\boldsymbol{e}_z \right] \right\} \right] \end{align}

(4.15b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0 \end{equation}

(4.15c)\begin{equation} \boldsymbol{U} = \boldsymbol{0} \quad \textrm{at} \ Z =\tfrac{1}{2}R_{c}^{{{1}/{2}}}\left({\mp} 1 - \mathcal{Z}_p \pm \epsilon_w \,{\rm Re}\{\exp [ {\mathrm{i}}\left(\varphi_p +X \kappa \right)] \} \right) .\end{equation}

In obtaining (4.15a), we neglect the time-derivative term. This is because the time scale for the establishment of the disturbance flow, $l^2/\nu$, is much smaller than the time scale over which the particle migrates laterally relative to the particle-free flow, ${{l}/{U_{lateral}}}$ (Hogg Reference Hogg1994). Here, $U_{lateral}$ denotes the lateral migration velocity of the particle. We also neglect a term $\epsilon ^3\boldsymbol {U}\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {U}$ in (4.15a), as the interaction of the perturbation flow with itself is negligible compared with other inertial terms. Next, we expand ${R}_c^{{{1}/{2}}}\boldsymbol {\bar {v}}$, $\boldsymbol {U}$, and $Q$ in powers of $\epsilon _w$

(4.16a)$$\begin{gather} {R}_c^{{1}/{2}}\bar{v}_x=\bar{V}_{x}^{\left( 0 \right)}+\epsilon _w\bar{V}_{x}^{\left( 1 \right)}+O ( \epsilon _{w}^{2} ) ,\quad {R}_c^{{1}/{2}}\bar{v}_z=\epsilon _w\bar{V}_{z}^{\left( 1 \right)}+O ( \epsilon _{w}^{2} ), \end{gather}$$

(4.16b)$$\begin{gather}\boldsymbol{U}=\boldsymbol{U}^{(0)}+\epsilon _w\boldsymbol{U}^{(1)}+O ( \epsilon _{w}^{2} ) , \quad Q=Q^{(0)}+\epsilon _wQ^{(1)}+O ( \epsilon _{w}^{2} ), \end{gather}$$

where $\bar {V}_{x}^{(0)}=\gamma Z-4Z^2{R}_c^{-{{1}/{2}}}$, $\bar {V}_{x}^{(1)}={R}_c^{{{1}/{2}}}\Delta \bar {\mathcal {V}}_{x}^{(1)}$ and $\bar {V}_{z}^{( 1 )}={R}_c^{{{1}/{2}}}\Delta \bar {\mathcal {V}}_{z}^{(1)}$ (3.10a,b), which reduce to $\bar {V}_{x}^{( 1 )}\rightarrow -12z\epsilon (\mathcal {Z} _p+z\epsilon R_{c}^{-{{1}/{2}}})\cos (\varphi _p)$ and $\bar {V}_{z}^{( 1 )}\rightarrow 0$ in the limit $\omega \rightarrow 0$. Expanding (4.15) and collecting like terms in $\epsilon _w$, we obtain the zeroth-order governing equations and boundary conditions

(4.17a)\begin{align} &\nabla ^2\boldsymbol{U}^{\left( 0 \right)}-\boldsymbol{\nabla} Q^{\left( 0 \right)}-\bar{V}_{x}^{\left( 0 \right)}\frac{\partial \boldsymbol{U}^{\left( 0 \right)}}{\partial X}-U_{z}^{\left( 0 \right)}\frac{\mbox{d}\bar{V}_{x}^{\left( 0 \right)}}{\mbox{d}Z}\boldsymbol{e}_x \nonumber\\ &\quad={-}\frac{10}{3}{{\rm \pi} R}_c^{-{1}/{2}}\gamma \left[ \frac{\partial \delta( \boldsymbol{R} )}{\partial Z}\boldsymbol{e}_x+\frac{\partial \delta( \boldsymbol{R} )}{\partial X}\boldsymbol{e}_z \right], \end{align}

(4.17b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U}^{\left( 0 \right)}=0, \end{equation}

(4.17c)\begin{equation} \boldsymbol{U}^{\left( 0 \right)}=\boldsymbol{0},\quad \mbox{at}\ Z_{-}={-}\tfrac{1}{2}{R}_c^{{1}/{2}}( 1+\mathcal{Z}_p ) \quad \mbox{and} \quad Z_+{=}\tfrac{1}{2}{R}_c^{{1}/{2}}( 1-\mathcal{Z}_p ), \end{equation}

(4.17d)\begin{equation} \boldsymbol{U}^{\left( 0 \right)}=\boldsymbol{0},\quad X\rightarrow \infty, \end{equation}

where $Z_{-}$ and $Z_+$ denote the $Z$-coordinates of the average locations of the walls. The first-order governing equations in $\epsilon _w$ take the form

(4.18a)\begin{align} &\nabla ^2\boldsymbol{U}^{(1)}-\boldsymbol{\nabla} Q^{\left( 1 \right)}-\bar{V}_{x}^{\left( 0 \right)}\frac{\partial \boldsymbol{U}^{\left( 1 \right)}}{\partial X}-U_{z}^{\left( 1 \right)}\frac{\mbox{d}\bar{V}_{x}^{\left( 0 \right)}}{\mbox{d}Z}\boldsymbol{e}_{x}=\bar{V}_{x}^{\left( 1 \right)}\frac{\partial \boldsymbol{U}^{\left( 0 \right)}}{\partial X}+\bar{V}_{z}^{\left( 1 \right)}\frac{\partial \boldsymbol{U}^{\left( 0 \right)}}{\partial Z} \nonumber\\ &\quad +\left( U_{x}^{\left( 0 \right)}\frac{\partial \bar{V}_{x}^{\left( 1 \right)}}{\partial X}+U_{z}^{\left( 0 \right)}\frac{\partial \bar{V}_{x}^{\left( 1 \right)}}{\partial Z} \right) \boldsymbol{e}_x+\left( U_{x}^{\left( 0 \right)}\frac{\partial \bar{V}_{z}^{\left( 1 \right)}}{\partial X}+U_{z}^{\left( 0 \right)}\frac{\partial \bar{V}_{z}^{\left( 1 \right)}}{\partial Z} \right) \boldsymbol{e}_z \nonumber\\ &\quad+\frac{10}{3}{{\rm \pi} R}_c^{-{1}/{2}}\,{\rm Re}\!\left[ \exp ( \mathrm{i}\varphi_p ) \left\{ \begin{array}{c} -\zeta \left[ \dfrac{\partial \delta ( \boldsymbol{R} )}{\partial Z}\boldsymbol{e}_x+\dfrac{\partial \delta ( \boldsymbol{R} )}{\partial X}\boldsymbol{e}_z \right] \\[10pt] +2\alpha _z\left[ \dfrac{\partial \delta ( \boldsymbol{R} )}{\partial X}\boldsymbol{e}_x-\dfrac{\partial \delta ( \boldsymbol{R} )}{\partial Z}\boldsymbol{e}_z \right] \end{array} \right\} \right] , \end{align}

(4.18b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U}^{\left( 1 \right)}=0. \end{equation}

Similar to our strategy in (2.4), we transfer the boundary condition (4.15c) to the average locations of the walls instead of enforcing $\boldsymbol {U}^{(1)}=\boldsymbol {0}$ at the wavy walls, such that $\boldsymbol {U}^{(1)}=\boldsymbol {0}$ is enforced at the wavy wall to the first order (Van Dyke Reference Van Dyke1975)

(4.19a)$$\begin{gather} \boldsymbol{U}^{\left( 1 \right)}={\pm} \frac{1}{2}\frac{\partial \boldsymbol{U}^{\left( 0 \right)}}{\partial Z}{R}_c^{{1}/{2}}\,{\rm Re}\{ \exp [\mathrm{i}( \varphi_p+\kappa X )] \} ,\quad\mbox{at}\ \boldsymbol{R}=\left( X,Y,Z_{{\pm}} \right) \end{gather}$$

(4.19b)$$\begin{gather}\boldsymbol{U}^{\left( 1 \right)}=\boldsymbol {0},\quad X\rightarrow \infty , \end{gather}$$

where $Z_{\pm }=\{ Z_{-}, Z_+ \}$ denotes the $Z$-coordinates of the average locations of the walls.

4.3. Matching inner and outer solutions

Matching the inner and outer solutions enables us to calculate the lift force on the particle in the $z$-direction. From (4.1a,b) and (4.14a,b), we know that the perturbation velocity component at the origin $u_{z}(\boldsymbol {R}=\boldsymbol {0})$, i.e. the lateral migration velocity of the particle relative to the undisturbed flow, can be determined from the difference between $\epsilon ^3 U_z( \boldsymbol {R}=\boldsymbol {0} )$ and $\epsilon u_{1,z}( \boldsymbol {R}=\boldsymbol {0} )$. The migration velocity $u_{z}(\boldsymbol {R}=\boldsymbol {0})$ gives rise to an inertial lift force acting on the particle in the $z$ direction (Schonberg & Hinch Reference Schonberg and Hinch1989; Hogg Reference Hogg1994; Asmolov Reference Asmolov1999) (see derivations of the force balance in Appendix C). Given the characteristic force scale $\mu aU_{m}^{\prime }$, the dimensionless lift force on the particle $F_{L,z}$ is then

(4.20)\begin{equation} F_{L,z}=6{\rm \pi} u_z =6{\rm \pi} \epsilon ^3\lim_{\boldsymbol {R}\rightarrow \boldsymbol {0}} \!( U_z- \epsilon^{{-}2}u_{1,z} ). \end{equation}

Since $u_{1,z}(\boldsymbol {R})$ has already been obtained in (4.13), we now provide a strategy to numerically solve $U_z$. We first convert the partial differential equations (PDEs) (4.17) and (4.18) with Fourier transforms, obtaining governing equations for $\boldsymbol {\tilde {U}}^{(0)}(k_x,k_y,Z)$ and $\boldsymbol {\tilde {U}}^{(1)}(k_x,k_y,Z)$, where the tildes denote quantities in the frequency domain, and $k_x$ and $k_y$ are frequencies. We further convert these PDEs into ordinary differential equations (ODEs) for $\tilde {U}_z^{(0)}$ and $\tilde {U}_z^{(1)}$, respectively. Numerically solving the ODEs for $\tilde {U}_z^{(0)}$ and $\tilde {U}_z^{(1)}$ and performing inverse Fourier transforms will enable us to obtain $F_{L,z}$ at different orders of $\epsilon _w$ based on (4.20). The detailed derivations for the ODEs are documented in Appendix D.1.

After solving $\tilde {U}_{z}^{( 1a )}$ and $\tilde {U}_{z}^{( 1b )}$ for arbitrary values of $k_x$ and $k_y$, we further evaluate the normalized lift force by performing the following inverse Fourier transform:

(4.21a)$$\begin{gather} F_{L,z} = F_{L,z}^{(0)} + \epsilon _w F_{L,z}^{\left( 1 \right)} + O(\epsilon _w^{2}), \end{gather}$$

(4.21b)$$\begin{gather}F_{L,z}^{(0)}=6{\rm \pi} \epsilon ^3\,{\rm Re}\!\left[ \int_{-\infty}^{+\infty}{\int_{-\infty}^{+\infty}\tilde{U}_{z}^{\left( 0 \right)}( k_x,k_y,0 )\,\mathrm{d}k_x\,\mathrm{d}k_y} \right], \end{gather}$$

(4.21c)$$\begin{gather}F_{L,z}^{(1)}=6{\rm \pi} \epsilon ^3\,{\rm Re}\!\left[ \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \tilde{U}_{z}^{\left( 1 \right)}( k_x,k_y,0 ) \,\mathrm{d}k_x\,\mathrm{d}k_y\right]. \end{gather}$$

According to (D6)–(D7), $\tilde {U}_{z}^{( 1 )}$ can be written in a form with the underlying parameters $\mathcal {X}_p$ and $\mathcal {Z}_p$ separated

(4.22)\begin{align} {\tilde{U}_{z}^{\left( 1 \right)}( k_x,k_y,0 )}=\tfrac{1}{2} [ \exp( \mathrm{i}\varphi_p ) \tilde{U}_{z}^{\left( 1a \right)}( \check{k}_x,k_y,0 ) +\exp ( -\mathrm{i}\varphi_p ) \tilde{U}_{z}^{\left( 1b \right)}( \hat{k}_x,k_y,0 )], \end{align}

where $\check {k}_x = k_x - \kappa$, $\hat {k}_x = k_x + \kappa$ and the particle's coordinate components $\mathcal {X}_p$ and $\mathcal {Z}_p$ are embedded in $\varphi _p$ and $\tilde {U}_{z}^{( 1a/b )}$, respectively. Compared with (4.20), (4.21) does not contain any terms from $u_{1,z}( \boldsymbol {R}=\boldsymbol {0} )$ (4.13), because the two-dimensional Fourier transform of a strainlet velocity field is imaginary at the origin (Schonberg & Hinch Reference Schonberg and Hinch1989; Hogg Reference Hogg1994; Asmolov Reference Asmolov1999). Thus, $F_{L,z}$ is different from the outer solution $U_z(\boldsymbol {R}=\boldsymbol {0})$ only by a factor of $6{\rm \pi} \epsilon ^3$. Notably, $F_{L,z}$ is a function of $\mathcal {X}_p$, $\mathcal {Z}_p$, $\omega$, $\epsilon$ and $\epsilon _w$.

It is important to note that the ODEs for $\tilde {U}_{z}^{( 1a )}$ and $\tilde {U}_{z}^{( 1b )}$ (D8)–(D12) reduce to the same form in the limit of $\omega \rightarrow 0$, giving rise to $\tilde {U}_{z}^{( 1a )}( \check {k}_x,k_y,0 ) = \tilde {U}_{z}^{( 1b )}( \hat {k}_x,k_y,0 )$. In this case, one can show that $F_{L,z}^{( 1 )}$ is in phase with the sinusoidal profile of the wavy channel (4.22).

We compute $F_{L,z}^{(0)}$, $F_{L,z}^{(1)}$ and $F_{L,z}$ across the channel based on (4.21) and present the results in figures 6–8, respectively. Figure 6 shows $\epsilon ^{-3}R_{c}^{{{1}/{2}}} F_{L,z}^{(0)}$ as a function of $\mathcal {Z}_p$ at varying $R_c$. The zeroth-order lift-force component $F_{L,z}^{(0)}$ points towards the wall, or outwards, when the particle is in the inner part of the channel; $F_{L,z}^{(0)}$ points towards the centreline, or inwards, when the particle gets close to the wall. As $R_c$ increases, the equilibrium position of $F_{L,z}^{(0)}$ shifts outwards, and the magnitude of $\epsilon^{-3} R_c^{1/2}F_{L,z}^{(0)}$ decreases in the inner part of the channel. Our computed zeroth-order lift force is consistent with the lift force reported in Asmolov (Reference Asmolov1999).

Figure 6. Normalized zeroth-order lift force $\epsilon ^{-3}R_{c}^{{{1}/{2}}}F_{L,z}^{(0)}$ as a function of the $\mathcal {Z}$-coordinate of the particle $\mathcal {Z}_p$ at varying channel Reynolds numbers $R_c$. For varying particle $z$-location $\mathcal {Z}_p$, the outer solution $U_z^{(0)}(\boldsymbol {R}=\boldsymbol {0})$ is different from $F_{L,z}^{(0)}$ only by a factor of $6{\rm \pi} \epsilon ^3$ (4.21b).

Figure 7. Normalized first-order (in $\epsilon _w$) lift force. (a,b) Value of $\epsilon ^{-3}R_c^{1/2}F_{L,z}^{(1)}$ as a function of the $\mathcal {Z}$-coordinate of the particle $\mathcal {Z}_p$ and the particle's phase angle in the x-direction $\varphi _p$ for (a) $R_c=10$, $\omega =2$ and (b) $R_c=100$, $\omega =0.02$. Results for $|\mathcal {Z}_p|>1$ are obtained by extrapolation. (c) Value of $\epsilon ^{-3}R_{c}^{{{1}/{2}}}F_{L,z}^{(1)}( \varphi _p=0 )$ and $\epsilon ^{-3}R_{c}^{{{1}/{2}}}F_{L,z}^{(1)}( \varphi _p={{{\rm \pi} }/{2}} )$ as a function of $\mathcal {Z}_p$ at varying $R_c$ and $\omega =2$. Here, a scaling factor $R_c^{1/2}$ makes the order of magnitude consistent for varying $R_c$. (d) Value of $\epsilon ^{-3}R_{c}^{{{1}/{2}}}F_{L,z}^{(1)}( \varphi _p=0 )$ and $\epsilon ^{-3}R_{c}^{{{1}/{2}}}F_{L,z}^{(1)}( \varphi _p={{{\rm \pi} }/{2}} )$ as a function of $\mathcal {Z}_p$ at varying $\omega$ and $R_c=100$. Dashed lines in (d) represent numerical solutions for $\omega \rightarrow 0$. For varying particle $z$-location $\mathcal {Z}_p$, the outer solution $U_z^{(1)}(\boldsymbol {R}=\boldsymbol {0})$ is different from $F_{L,z}^{(1)}$ only by a factor of $6{\rm \pi} \epsilon ^3$ (D6a) and (4.21c).

Figure 8. Normalized lift force $\epsilon ^{-3}F_{L,z}$ as a function of the $\mathcal {Z}$-coordinate of the particle $\mathcal {Z}_p$ and the phase angle of the particle in the $x$-direction $\varphi _p$ for (a) $R_c=10$, $\omega =2$ and (b) $R_c=100$, $\omega =0.02$. (c) Zero-lift-force location $\mathcal {Z}_p|_{F_{L,z}=0}$ as a function of $\varphi _p$ for varying $R_c$ and $\omega =2$. (d) Zero-lift-force location $\mathcal {Z}_p|_{F_{L,z}=0}$ as a function of $\varphi _p$ for varying $\omega$ and $R_c=100$. $\epsilon _w=0.05$ for all panels.

The first-order lift force $F_{L,z}^{(1)}$ varies in both the $x$- and $z$-directions and is a function of the channel Reynolds number $R_c$ and the frequency $\omega$. Figures 7(a) and 7(b) show the variation of the normalized first-order lift force $\epsilon ^{-3}R_c^{1/2}F_{L,z}^{(1)}$ in a period of the channel for two pairs of $R_c$ and $\omega$; due to the sharp variations close to the walls, we adopt a moderate range of $[-200,200]$ for $\epsilon ^{-3}R_c^{1/2}F_{L,z}^{(1)}$ to better illustrate the variations in the bulk of the channel. The lift force $F_{L,z}^{(1)}$ has a relatively low magnitude close to the centre of the channel, and its magnitude increases sharply close to the walls. The sign of $F_{L,z}^{(1)}$ also varies with the phase angle of the particle in the $x$-direction $\varphi _p$, where $\varphi _p=\omega \mathcal {X}_p$. According to (4.21c), for any particle location,

(4.23)\begin{equation} F_{L,z}^{(1)}(\mathcal{X}_p,\mathcal{Z}_p) = \sin ( \varphi_p ) F_{L,z}^{(1)}( \mathcal{X}_p=0,\mathcal{Z}_p ) + \cos ( \varphi_p ) F_{L,z}^{(1)}\left( \mathcal{X}_p=\frac{\rm \pi}{2 \omega} ,\mathcal{Z}_p \right). \end{equation}

Thus, we present $F_{L,z}^{(1)}$ as a function of $\mathcal {Z}_p$ at two phase angles in the $x$-direction, i.e. $\varphi _p=0$ and $\varphi _p={\rm \pi} /2$ in figures 7(c) and 7(d), respectively. Figure 7(c) shows the influence of $R_c$ at $\omega =2$, and figure 7(d) shows the influence of $\omega$ at $R_c=100$. Notably, increasing $R_c$ and decreasing $\omega$ have a similar effect on the $\mathcal {Z}_p$–$F_{L,z}^{(1)}$ relationship, as illustrated in the resemblance in the variation of the curves between figures 7(c) and 7(d). Figure 7(d) also shows a vanishing $F_{L,z}^{(1)}$ for $\varphi _p={\rm \pi} /2$ in the limit of $\omega \rightarrow 0$, which accords with our observation that $F_{L,z}^{(1)}$ is in phase with the wavy channel.

By adding $F_{L,z}^{(0)}$ and $F_{L,z}^{(1)}$, we show the variation of $\epsilon ^{-3} F_{L,z}$ across the channel for two pairs of $R_c$ and $\omega$ in figures 8(a) and 8(b). We are interested in the locations where $F_{L,z}=0$, i.e. the zero contour lines. Particles are unstable on the zero contour line at the centre because the nearby lift force points away from the centreline. In figure 8(a), we also observe a pair of unstable zero contour lines along part of the channel close to the walls (mostly within $\varphi _p \in [{\rm \pi},2{\rm \pi} ]$). These near-wall zero contour lines arise because $F_{L,z}^{(0)}$ and $F_{L,z}^{(1)}$ have different signs in these regions and $F_{L,z}^{(1)}$ has a large magnitude close to the walls for $R_c=10,\omega =2$ (figure 7c). Since previous studies report that the lift force directs particles away from the walls (Asmolov Reference Asmolov1999; Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009), it is unlikely that particles may be pushed toward the walls as they approach these unstable contour lines. Note that we compute the undisturbed flow and the perturbation velocity up to the first order in $\epsilon _w$. Due to the sharp variations in the first-order flow and perturbation fields near the walls, higher-order computations may be necessary for a better understanding of the near-wall lift force, which is, however, beyond the scope of our current analyses. The critical stable zero contour lines are close to $\mathcal {Z}_p=\pm 0.6$. The shapes of these stable zero contour lines do not always follow the configuration of the channel. For a large $R_c$ and a small $\omega$ (figure 8b), the shape of the stable zero contour line is in phase with the configuration of the wavy channel. As $R_c$ decreases or $\omega$ increases, the shape of the stable zero contour line gradually becomes out of phase, as shown in figures 8(c) and 8(d).

Finally, it is important to note that the zero-lift-force locations are generally different from the focusing locations of particles in wavy channels. Unlike the situation in straight channels where both types of locations coincide, the secondary flow in wavy channels makes particles move in the $z$-direction constantly. Thus, even if a particle finally stabilizes in a wavy channel, it never truly reaches equilibrium. To resolve the focusing location, it is necessary to analyse the particle trajectory, which we address in § 5.

5. Particle focusing locations

After obtaining the lift force on the particle, we here provide a strategy to predict the focusing locations of the particle in the wavy channel. For this, we come back to the laboratory frame (figure 3). We use $\boldsymbol {\mathcal {R}}_p^{\prime}( t ) =( \mathcal {X} _p^{\prime}( t ) ,\mathcal {Z} _p^{\prime}( t ) )$ to denote the instantaneous location of the particle centre. Instead of solving the fluid stress over the surface of the particle and using Newton's second law, we employ the modified Maxey–Riley equation to analyse the trajectory of the particle (Maxey & Riley Reference Maxey and Riley1983; Ferry & Balachandar Reference Ferry and Balachandar2001)

(5.1) \begin{align} m_p \frac{\mbox{d}\boldsymbol{U}_{p}^{\prime}}{\mbox{d}t^{\prime}} &= 6{\rm \pi} a\mu ( \bar{\boldsymbol{\mathcal{V}}}'-\boldsymbol{U}_{p}^{\prime} ) +m_f\frac{\mbox{D}\bar{\boldsymbol{\mathcal{V}}}'}{\mbox{D}t^{\prime}}+\frac{1}{2}m_f\left( \frac{\mbox{D}\bar{\boldsymbol{\mathcal{V}}}'}{\mbox{D}t^{\prime}}-\frac{\mbox{d}\boldsymbol{U}_{p}^{\prime}}{\mbox{d}t^{\prime}} \right) +( m_p-m_f ) \boldsymbol{g}^{\prime} \nonumber\\ &\quad +\frac{6{\rm \pi} a^2\mu}{\sqrt{\nu}}\frac{\mbox{d}^{{{1}/{2}}}}{\mbox{d}t^{\prime {{1}/{2}}}}( \bar{\boldsymbol{\mathcal{V}}}'-\boldsymbol{U}_{p}^{\prime} ) +\boldsymbol{F}_{L}^{\prime}, \end{align}

where $m_p$ is the particle mass, $m_f$ is the fluid mass displaced by the particle, $\bar {\boldsymbol {\mathcal {V}}}^{\prime }=\bar {\boldsymbol {\mathcal {V}}}^{\prime }( \bar {\boldsymbol {\mathcal {R}}}_{p}^{\prime }( t ) )$ is the particle-free fluid velocity at the instantaneous particle centre, ${{\mathrm {D}}/{\mathrm {D}}}t^{\prime }$ is the material derivative following a fluid parcel, $\boldsymbol{g}'$ is the gravitational acceleration and $\mu$ is the dynamic viscosity of the fluid. Additionally, ${\mbox {d}\boldsymbol {U}_{p}^{\prime }}/{\mbox {d}t^{\prime }}=\boldsymbol {a}_p^{\prime }$ is the particle acceleration defined in § 3. The terms on the right-hand side of (5.1) represent the Stokes drag, fluid acceleration, added mass, gravity, Basset force and lift force on the particle. For the scaling analysis, we choose $l$ as the characteristic length, $U_{m}^{\prime }$ as the characteristic velocity, ${{l}/{U_{m}^{\prime }}}$ as the characteristic time and $\mu aU_{m}^{\prime }$ as the characteristic lift force. By comparing the forcing terms, we neglect the Basset force and write (5.1) in the following dimensionless form:

(5.2)\begin{align} \left( 1+\frac{1}{2}\frac{\rho _f}{\rho _p} \right) Stk \frac{\mathrm{D}\!\left( \boldsymbol{U}_p-\bar{\boldsymbol{\mathcal{V}}} \right)}{\mathrm{D}t}={-}\left( \boldsymbol{U}_{p}-\bar{\boldsymbol{\mathcal{V}}} \right) +\left( \frac{\rho _f}{\rho _p}-1 \right) Stk \frac{\mathrm{D}\bar{\boldsymbol{\mathcal{V}}}}{\mathrm{D}t}+B\boldsymbol{e}_g+\frac{1}{6{\rm \pi}}\boldsymbol{F}_L, \end{align}

where $\rho _p$ is the particle density, $B=2a^2( \rho _p-\rho _f ) \boldsymbol {g}{{^{\prime }}/{( 9\mu U_{m}^{\prime } )}}$ is the buoyancy number, $\boldsymbol {e}_g$ is a unit vector in the direction of the gravitational force and $Stk$ is the Stokes number. Specifically, $Stk$ is the ratio between the particle relaxation time $t_p$ and the characteristic time of the flow $t_f$, where $t_p=2\rho _p a^2/{( 9\mu )}$ and $t_f={{l}/{U_{m}^{\prime }}}$. Since our asymptotic parameter $\epsilon =R_{c}^{1/2}a/l={(\rho _f U_m^{\prime }l/\mu )}^{1/2}a/l$, we have

(5.3)\begin{equation} Stk=\frac{2}{9}\frac{\rho _p}{\rho _f}\epsilon ^2\ll 1. \end{equation}

For a neutrally buoyant particle, (5.2) can be converted into a differential equation of the particle location

(5.4)\begin{equation} \frac{1}{3}\epsilon ^2\ddot{\boldsymbol{\mathcal{R}}}_p+\dot{\boldsymbol{\mathcal{R}}}_p=\bar{\boldsymbol{\mathcal{V}}}+\frac{1}{3}\epsilon ^2 \frac{\mathrm{D}\bar{\boldsymbol{\mathcal{V}}}}{\mathrm{D}t}+\frac{1}{6{\rm \pi}}\boldsymbol{F}_L, \end{equation}

where $\dot {\boldsymbol {\mathcal {R}}}_p$ and $\ddot {\boldsymbol {\mathcal {R}}}_p$ denote the velocity and acceleration of the particle, respectively. Here, $\dot {\boldsymbol {\mathcal {R}}}_p=\boldsymbol {U}_p$. With $\bar {\boldsymbol {\mathcal {V}}}$ and $\boldsymbol {F}_L$ solved in previous sections, we can simulate the trajectory of the particle based on (5.4). By holding $R_c$, $\omega$, and $\epsilon _w$ constant, we find that the simulated trajectory gradually stabilizes over time for any initial particle location. We regard the stabilized trajectory in a period of the channel as particle focusing locations. Additionally, the stabilized trajectory asymptotes to a constant shape as $\epsilon$ decreases; as indicated by the stabilized range of the particle's $z$-location (inset of figure 9b), the stabilized trajectory for $\epsilon = 0.2$ is not significantly different from that for $\epsilon = 0.1$.

Figure 9. (a) Simulated trajectories of a particle initially located at ${\mathcal {Z}}_p=0.1$ in channels with varying amplitude $\epsilon _w$, where the lift forces are the zeroth- and first-order lift forces (a-1) and zeroth-order lift force only (a-2). Here, $R_c=10$, $\omega =2$ and $\epsilon =0.2$. (b) Influence of asymptotic parameter $\epsilon$ on the simulated trajectory of a particle initially located at ${\mathcal {Z}}_p=0.1$, where $R_c=10$, $\omega =2$, $\epsilon _w=0.05$ and the lift force contains both the zeroth- and first-order lift forces. (c) Average focusing location $\bar {\mathcal {Z}}_p$ as a function of channel Reynolds number $R_c$ for varying $\epsilon _w$, where $\omega =2$. (d) Average focusing location $\bar {\mathcal {Z}}_p$ as a function of channel angular frequency $\omega$ for varying $\epsilon _w$, where $R_c=100$.

In figure 9(a), we show two sets of simulated particle trajectories in channels with amplitudes $\epsilon _w=0,0.05$ and 0.1 for $R_c=10,\omega =2$ and $\epsilon =0.2$. The results in figure 9(a-1) are based on both the zeroth-order and first-order lift forces, i.e. $F_{L,z}$. For all three values of $\epsilon _w$, the particle is initially at $\mathcal {Z}_p=0.1$ and gradually migrates towards the stable focusing locations. When $\mathcal {X}_p=2000$, the movement of the particle is stabilized: for the straight channel ($\epsilon _w=0$), $\mathcal {Z}_p$ does not change over time and agrees with the asymptotic solution based on Asmolov (Reference Asmolov1999); for the wavy channels ($\epsilon _w=0.05$ and 0.1), $\mathcal {Z}_p$ oscillates sinusoidally about an average level. As $\epsilon _w$ increases, the amplitude of $\mathcal {Z}_p$ also increases but is always lower than $\epsilon _w$; additionally, the average focusing location shifts slightly inwards. Interestingly, the channel length for particle focusing is independent of $\epsilon _w$, but is significantly impacted by $\epsilon$ (figure 9b). We notice that when $\epsilon$ becomes very small, e.g. $\epsilon =0.01$, a significantly longer channel is required to achieve particle focusing. Because $\boldsymbol {F}_L\sim O(\epsilon ^3 )$, a smaller $\epsilon$ gives rise to a smaller lift force, which causes the particle to migrate more slowly.

To understand which part of the lift force contributes to the oscillation and shifting of the focusing locations, we only keep the zeroth-order lift force $F_{L,z}^{(0)}$ and repeat the simulations for the particle trajectories. By comparing the insets of figures 9(a-1) and 9(a-2), we find only slight differences in the average focusing locations for $\epsilon _w=0.05$ and 0.1. All other features of figures 9(a-1) and 9(a-2) are indistinguishable. The zeroth-order lift force is thus the governing force for particle focusing in wavy channels. Despite its sinusoidal characteristics, the first-order lift force $F_{L,z}^{(1)}$ plays an almost negligible role in the oscillation of the particle in the steady state. Rather, the oscillation of the particle stems from the interaction between $F_{L,z}^{(0)}$ and the periodic flow $\bar {\boldsymbol {\mathcal {V}}}$. Furthermore, due to its sinusoidal characteristics, $F_{L,z}^{(1)}$ does not induce any variation in the length of the channel for focusing in figure 9(a). $F_{L,z}^{(1)}$ mainly contributes to the slight variation in the average focusing location in the $z$-direction $\bar {\mathcal {Z}}_p$.

We now mathematically understand the role of $F_{L,z}^{(1)}$ in influencing $\bar {\mathcal {Z}}_p$ in the asymptotic limit $\epsilon _w \ll 1$ and $\epsilon \ll 1$. Further neglecting the $O(\epsilon ^2)$ terms in (5.4), we obtain the following equation for a neutrally buoyant particle:

(5.5)\begin{equation} -\boldsymbol{V}_s=\frac{1}{6{\rm \pi}}\boldsymbol{F}_L, \end{equation}

where $\boldsymbol {V}_s=\bar {\boldsymbol {\mathcal {V}}}-\dot {\boldsymbol {\mathcal {R}}}_p=\bar {\boldsymbol {\mathcal {V}}}-\boldsymbol {U}_p$ is the slip velocity defined in § 3, and $-V_{s,z}$ corresponds to the lateral migration velocity of the particle relative to the particle-free flow $U_{lateral}$. Denoting the average slip velocity in the $z$-direction over time as $\overline {V_{s,z}}$, then $\overline {V_{s,z}( \mathcal {Z} _p( t ) ) }=0$ represents the stable state. We have the stability condition

(5.6)\begin{equation} \overline{F_{L,z}( \mathcal{Z} _p( t ) ) }=0. \end{equation}

Because $\epsilon _w$ is small, the particle does not deviate significantly from $\bar {\mathcal {Z}} _p$. Thus, we can expand $F_{L,z}( \mathcal {Z} _p( t ) )$ into a Taylor series about $\bar {\mathcal {Z}} _p$

(5.7)\begin{equation} F_{L,z}( \mathcal{Z} _p ) =F_{L,z}( \bar{\mathcal{Z}} _p ) +\left( \mathcal{Z} _p-\bar{\mathcal{Z}} _p \right) \left. \frac{\mathrm{d}F_{L,z}}{\mathrm{d}\mathcal{Z} _p} \right|_{\mathcal{Z} _p=\bar{\mathcal{Z}} _p}+O ( ( \mathcal{Z} _p-\bar{\mathcal{Z}} _p ) ^2 ). \end{equation}

Neglecting high-order terms, the stability condition (5.6) becomes

(5.8)\begin{equation} \overline{F_{L,z}( \bar{\mathcal{Z}} _p ) }+\overline{\left( \mathcal{Z} _p-\bar{\mathcal{Z}} _p \right) \left. \frac{\mathrm{d}F_{L,z}}{\mathrm{d}\mathcal{Z} _p} \right|_{\mathcal{Z} _p=\bar{\mathcal{Z}} _p}}=0, \end{equation}

where $\mathcal {Z} _p-\bar {\mathcal {Z}} _p=\int ( \bar {\mathcal {V}}{{_z}/{\bar {\mathcal {V}}_x}} )\,\mathrm {d}\mathcal {X} _p$. Evaluating (5.8) with results from (2.5) and (4.21) and only keeping the leading-order terms in $\epsilon _w$, we finally arrive at the following implicit expression for $\bar {\mathcal {Z}} _p$:

(5.9) \begin{align} & F_{L,z}^{\left( 0 \right)}( \bar{\mathcal{Z}} _p ) -\frac{\epsilon _w}{\omega ( 1-\bar{\mathcal{Z}} _p^{2})}\left. \left\{ {\rm Re}[ \hat{\mathcal{V}}_z( \bar{\mathcal{Z}} _p ) ] \frac{\mathrm{d}F_{L,z}^{\left( 0 \right)}}{\mathrm{d}\mathcal{Z} _p} \right\} \right|_{\mathcal{Z} _p=\bar{\mathcal{Z}} _p} \nonumber\\ &\quad+\frac{\epsilon _{w}^{2}}{2\omega (1-\bar{\mathcal{Z}} _p^{2})}\left\{ {\rm Re}[ \hat{\mathcal{V}}_z( \bar{\mathcal{Z}} _p ) ] \left. \frac{\partial F_{L,z}^{\left( 1 \right)}}{\partial \mathcal{Z} _p} \right|_{\mathcal{Z} _p=\bar{\mathcal{Z}} _p,\varphi_p={\rm \pi}/{2}}\right. \nonumber\\ &\quad +\left.{\rm Im}[ \hat{\mathcal{V}}_z( \bar{\mathcal{Z}} _p ) ] \left. \frac{\partial F_{L,z}^{\left( 1 \right)}}{\partial \mathcal{Z} _p} \right|_{\mathcal{Z} _p=\bar{\mathcal{Z}} _p,\varphi_p=0} \right\} =0. \end{align}

Clearly, the terms related to $F_{L,z}^{(1)}$ are of $O(\epsilon _w^2)$, which explains the limited contribution of $F_{L,z}^{(1)}$ to $\bar {\mathcal {Z}}_p$. To interpret this result physically, the interaction between the first-order lift force and the first-order particle-free flow is not comparable to the interaction between the zeroth-order lift force and the first-order particle-free flow.

By repeating the simulations of particle trajectories with (5.4), we investigate the influence of channel Reynolds number $R_c$ and channel angular frequency $\omega$ on the average focusing location $\bar {\mathcal {Z}}_p$. Figure 9(c) shows $\bar {\mathcal {Z}}_p$ vs $R_c$ for $\epsilon _w=0$, 0.05 and 0.1, where $\omega =2$. As $R_c$ increases, $\bar {\mathcal {Z}}_p$ shifts outwards for all three values of $\epsilon _w$. For a constant $R_c$, $\bar {\mathcal {Z}}_p$ gets closer to the centre of the channel as $\epsilon _w$ increases. In figure 9(d), we show the relationship between $\bar {\mathcal {Z}}_p$ and $\omega$ for $\omega \in [ 0.02,5 ]$ and $R_c=100$. $\bar {\mathcal {Z}}_p$ does not vary monotonically with $\omega$. When $\omega$ approaches zero, particles behave as if in a straight channel; thus, $\bar {\mathcal {Z}}_p$ of varying $\epsilon _w$ also approaches $\bar {\mathcal {Z}}_p$ in the straight channel. For $\omega >2$, $\bar {\mathcal {Z}}_p$ slightly increases with $\omega$. However, the influence of $\epsilon _w$ on $\bar {\mathcal {Z}}_p$ is much less obvious than that of $R_c$. These results suggest that adjusting $\omega$ is not an effective strategy for altering the particle focusing locations.

6. Experiments

In this section, we experimentally investigate particle focusing in wavy channels and compare the results with those obtained from a straight channel. All channels have the same width $l$, where $l=2h=1.18$ mm. The experimental set-up and procedures are described in Appendix E. Briefly, we use two sizes of particles with average radii $a=200\ \mathrm {\mu }\textrm {m}$ and $60\ \mathrm {\mu }\textrm {m}$, which correspond to $a/h=0.17$ and $0.05$, respectively. The particles are neutrally buoyant in our water–glycerol solution. We inject the particle suspensions in the channels and observe particle behaviours at a fixed location. The distance from the channel inlet and the observation location is denoted as $L$, where $L=550$ mm (figure 10a). It is important to note that our theoretical analyses are based on ideal assumptions $\epsilon =R_{c}^{{{1}/{2}}}{{a}/{l}}\ll 1$, ${{a}/{h}}\ll 1$, ${\epsilon _w}\ll 1$, and infinite channel depth $w$. Since these parameters are finite in practice, in addition to testing our asymptotic theory, we aim to experimentally characterize particle focusing when the real set-up deviates from our ideal assumptions.

Figure 10. (a) Schematics of the experimental set-up, where $L$ denotes the distance between the channel inlet and the fixed observation location. (b) Images of a polystyrene particle in channels with varying amplitudes. All four channels have average width $l=1.18$ mm. The scale bars denote $500\ \mathrm {\mu } \textrm {m}$.

In figure 11, we present particle distributions at $L/l=466$ in channels with varying amplitude $\epsilon _w$ at varying $R_c$, where $a/h=0.17$. For a constant $\epsilon _w$, particles gradually get focused as $R_c$ increases. As predicted by our asymptotic theory, we observe that particles are focused on two symmetric lines, which is a hallmark of a strong inertial lift effect. At a small $\epsilon _w$, e.g. $\epsilon _w=0.093$, the focusing locations are along two curves that are almost parallel to the channel profiles. For a relatively large $\epsilon _w$, e.g. $\epsilon _w=0.249$ or 0.490, the focusing curves deviate from the channel profiles and flatten out as $R_c$ increases. When $\epsilon _w=0.490$ and $R_c=250$, some particles even get trapped in the expanded sections of the channels. To understand these variations, we compute the streamlines in the wavy channels with Ansys Fluent. For all values of $\epsilon _w$ and $R_c$, the focusing curves overlap two symmetric streamlines (figure 11), indicating that the particle trajectories follow the particle-free streamlines. For a small $\epsilon _w$ or $R_c$, all streamlines mostly follow the shape of the channels; however, when both $\epsilon _w$ and $R_c$ are large enough, vortices appear in the expanded section, and the continuous streamlines only appear in the middle part of the channel. Consequently, most particles are focused along the flattened streamlines close to the centre, but some particles are occasionally trapped in the vortices. Furthermore, to probe the capability of our asymptotic theory in predicting particle focusing locations, we present in figure 11 the simulated results based on the asymptotic lift force $F_{L,z}^{(0)}$ and (5.4). Remarkably, although the assumptions $\epsilon \ll 1$ and $a/h \ll 1$ are not strictly satisfied in the experiments, the asymptotic theory provides fairly close estimates of the experimental focusing locations.

Figure 11. Experimental and asymptotic particle focusing locations in channels with different amplitudes $\epsilon _w$ at varying channel Reynolds numbers $R_c$. Asymptotic focusing locations are based on simulations in § 5, and experimental focusing locations are measured at $L/l=466$. (*) Asterisks denote little or no focusing at a low $R_c$. Here, $a/h=0.17$.

Another takeaway from figure 11 is that wavy channels can help particles focus at a slightly lower channel Reynolds number $R_c$ than the straight channel while keeping the average channel width $l$ and observation position $L/l$ constant. For $R_c=10$, for instance, particles are almost fully focused in the wavy channels, whereas little focusing is observed in the straight channel. To rationalize these observations, we present in figure 12(a) the distributions of two sizes of particles at increasing $R_c$ in the straight channel and the wavy channel with $\epsilon _w=0.093$. The critical $R_c$ for the focusing transition is of $O(10)$ for $a/h=0.17$, and the critical $R_c$ is of $O(10^2\unicode{x2013} 10^3)$ for $a/h=0.05$. The variation in the critical $R_c$ with $a/h$ can be explained by the asymptotic results: in figures 6 and 7, we show that $F_{L,z}\sim O(\epsilon ^3 R_c^{-1/2} )$. We can bring $F_{L,z}$ back into the dimensional format $F_{L,z}^{\prime }=C \mu U_m^{\prime } R_c a^4/l^3$ with $\epsilon =R_c^{1/2} a/l$ and $F_{L,z}^{\prime }=\mu aU_m^{\prime } F_{L,z}$, where $C$ is a constant. Moreover, the lift force must be in balance with a drag force $F_D^{\prime }$ during lateral migration, where $F_D^{\prime }=6 {\rm \pi}\mu a U_{lateral}^{\prime }$. Since $F_{L,z}^{\prime }=F_D^{\prime }$, we obtain $L_{focus}/l=U_m^{\prime }/U_{lateral}^{\prime }=(3{\rm \pi} /4C) R_c^{-1} h^{3}/a^3$, where $L_{focus}$ is the length of the channel required for particle focusing. For our fixed length $L/l$, the critical $R_c$ for particle focusing thus follows a power law with $a/h$

(6.1)\begin{equation} R_{c, critical}=\frac{3{\rm \pi}}{4C}\frac{l}{L}\left( \frac{a}{h} \right) ^{{-}3}, \end{equation}

which coincides with the theories for a straight channel (Matas et al. Reference Matas, Morris and Guazzelli2004). We select $C=0.1\unicode{x2013} 1$ and show the critical $R_c$ as a function of $a/h$ in figure 12(a).

Figure 12. (a) Comparison of particle distributions in the wavy and straight channels at increasing channel Reynolds number $R_c$, where $\epsilon _w=0.093$ for the wavy channel, and the normalized particle sizes $a/h$ are $0.17$ (left) and $0.05$ (right). (b) Experimentally obtained average focusing location $\bar {\mathcal {Z}}_p$ as a function of $R_c$ compared with results based on the asymptotic theory. (c) Transition of the focusing states occurring as the channel Reynolds number $R_c$ increases from 5 to 25 at time $t=0$. The upper panel shows the time required for the particles to be focused in channels with varying amplitudes $\epsilon _w$. This time is converted into the equivalent normalized channel lengths $L/l$ in the lower panel. Here, $a/h=0.17$. The online supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.558 shows particle focusing in a wavy channel occurring as the channel Reynolds number increases from 5 to 25 at time $t = 0$.

The insets of figure 12(a) show that particles are focused at a slightly lower $R_c$ in the wavy channel, especially for $a/h=0.17$; for the same $R_c$ and $l$, it requires a shorter distance for particles to be focused in a wavy channel. This finding, however, does not contradict the results in figure 9(a) based on the asymptotic theory. First, the experimentally obtained $R_{c,critical}$ is of the same order of magnitude in both types of channels. Second, the difference in $L_{focus}$ is due to the deviation from the asymptotic assumption $a/h \ll 1$ in our experiments. For a finite-sized particle, we can still decompose $F_{L,z}$ into $F_{L,z}^{( 0 )}+ \epsilon _w F_{L,z}^{( 1 )}$. However, unlike the asymptotic case in which $F_{L,z}^{(1)}$ varies sinusoidally along the $x$-direction at a constant $\mathcal {Z}_p$ (4.21c), the symmetry is broken for a finite $a/h$; $F_{L,z}^{(1)}$ of a finite-sized particle has a higher magnitude in the constricted sections of the wavy channel than in the expanded sections. This is because the distance to the wall is now finite relative to the particle size, and a stronger hydrodynamic interaction with the wall gives rise to greater repulsion effects in the constricted section, leading to slightly faster particle focusing in wavy channels. This effect becomes weaker as the particle size approaches the asymptotic limit, e.g. $a/h=0.05$ in figure 12(a). Third, when $\epsilon _w$ is even larger, e.g. $\epsilon _w=0.249$ and 0.490 in figure 11, the flattening of the continuous streamlines also causes the particles to experience stronger repulsion effects in the constricted sections.

Figure 12(b) shows the relationship between the average focusing location $\bar {\mathcal {Z}}_p$ and $R_c$ for the channels with $\epsilon _w=0$ and 0.093. As predicted by the asymptotic theory, the measured $\bar {\mathcal {Z}}_p$ is positively related to $R_c$ for $a/h=0.05$. However, as the particle size increases, the measured $\bar {\mathcal {Z}}_p$ gradually deviates from the asymptotic predictions and moves towards the centre of the channel. This behaviour is consistent with previous experimental observations in circular pipes (Matas et al. Reference Matas, Morris and Guazzelli2004) and square cross-section channels (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009). Because the motion of a finite-size particle causes a significant disturbance to the particle-free flow, it is well understood that the lift force is subject to a different scaling law than $O( \epsilon ^3R_{c}^{-{{1}/{2}}} )$ and that the equilibrium positions move closer to the centreline (Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Toner2009; Hood et al. Reference Hood, Lee and Roper2015). Additionally, the measured $\bar {\mathcal {Z}}_p$ in the wavy channel is consistently lower than that measured in the straight channel, which agrees with the asymptotic predictions (figures 9(c) and 11). Our experiments show that $\bar {\mathcal {Z}}_p$ starts to decrease for $R_c = 250$ and $\epsilon _w=0.093$; by contrast, our asymptotic theory demonstrates this phenomenon only at a higher $\epsilon _w$ (figure 11). This is because a finite-sized particle moves closer to the channel centreline, and $\bar {\mathcal {Z}}_p$ decreases at $R_c = 250$ even for a relatively small $\epsilon _w$.

Finally, because the difference in the critical $R_c$ for the focusing transition can be subtle in figure 11, we here provide a second strategy to demonstrate that particles get focused faster in wavy channels than in straight channels under the same conditions. Figure 12(c) shows the response of the particle distribution at $L/l= 466$ in channels with varying $\epsilon _w$ as $R_c$ increases from 5 to 25 at time $t = 0$. The detailed experimental method is described in Appendix E. When $t<0$, we set $R_c=5$ and particles are unfocused in all four channels. As $R_c$ increases to 25, the particles become focused after a certain time that decreases with increasing $\epsilon _w$. Since $R_c=U_m^{\prime } l/\nu$, where $U_m^{\prime }$ is twice the average velocity in the channel $U_a^{\prime }$, we can estimate the equivalent length required for the focusing transition $L_{focus}$ as $L_{focus}=U_a^{\prime } t$. $L_{focus}$ decreases with increasing channel amplitude.