Three-dimensional streaming around an obstacle in a Hele-Shaw cell

Abstract

Driving oscillatory flow around an obstacle generates, due to inertial rectification, a steady ‘streaming’ flow that is useful in a host of microfluidic applications. While theory has focused largely on two-dimensional flows, streaming in many practical microfluidic devices is three-dimensional due to confinement. We develop a three-dimensional streaming theory around an obstacle in a microchannel with a Hele-Shaw-like geometry, where one dimension (depth) is much shorter than the other two dimensions. Utilizing inertial lubrication theory, we demonstrate that the time-averaged streaming flow has a three-dimensional structure. Notably, the flow reverses direction across the depth of the channel, which is a feature not observed in less confined streaming set-ups. This feature is confirmed by our experiments of streaming around a cylinder sandwiched in a microchannel. Our theory also predicts that the streaming velocity decays as the inverse cube of the distance from the cylinder, faster than that expected from previous two-dimensional approaches. We verify this velocity decay quantitatively using particle tracking measurements from experiments of streaming around cylinders with different aspect ratios at different driving frequencies.

1. Introduction

Streaming refers to a secondary steady flow driven by the inertia of primary oscillations of a viscous fluid (Riley 2001; Boluriaan & Morris 2003; Wu 2018). The generation of streaming requires a spatially varying oscillatory flow, which can be established by the propagation of acoustic waves, by oscillations of an interface, or by bulk oscillation of fluid around an obstacle (Eckart 1948; Bolaños-Jiménez et al. 2017; Karlsen et al. 2018). Streaming flows around obstacles in micro-scale confinement are typically characterized by vortical structures and have been used widely in microfluidic applications due to their fast flow speeds, simplicity of generation, and versatility. For instance, eddies of acoustic microstreaming can be used to modify the composition of polydisperse suspensions by trapping smaller microparticles and releasing particles of larger size (Stone & Kim 2001; Beebe, Mensing & Walker 2002; Lutz, Chen & Schwartz 2003; Thameem, Rallabandi & Hilgenfeldt 2017), or to control vesicle deformation and lysis for bioengineering applications (Marmottant & Hilgenfeldt 2003; Marmottant, Biben & Hilgenfeldt 2008; Tandiono et al. 2012; Bhosale et al. 2022). Acoustic streaming induced by surface acoustic waves has been utilized widely to enhance mixing in laminar flow in microchannels (Westerhausen et al. 2016; Ahmed et al. 2019), to precisely control the intercellular distance and spatial arrangement of cells (Guo et al. 2015; Mutafopulos et al. 2019), and in the design of ‘acoustic tweezers’ that can trap suspended microparticles (Shi et al. 2009; Ding et al. 2012; Zhu et al. 2021).

Streaming around cylindrical obstacles (rigid cylinders and bubbles) is well understood under a two-dimensional framework (Holtsmark et al. 1954; Bertelsen, Svardal & Tjøtta 1973; Hamilton, Ilinskii & Zabolotskaya 2003; Chong et al. 2013; Červenka & Bednařík 2016; Lei et al. 2018) and has been verified experimentally (Andrade 1931; Holtsmark et al. 1954; Bertelsen et al. 1973; Lutz, Chen & Schwartz 2005; Rallabandi, Wang & Hilgenfeldt 2014). The flow is organized into four symmetric pairs of vortices around a rigid cylinder; each vortex pair consists of an inner and an outer vortex whose relative size depends on the dimensionless Stokes layer thickness (Holtsmark et al. 1954; Lutz et al. 2005). Vishwanathan & Juarez (2019) recently reported a method to generate and control the size of these vortices in a microchannel using sub-kilohertz oscillations of a loudspeaker. Recent work has also elucidated the role of geometry and topology to control three-dimensional (3-D) streaming flows around isolated curved objects (Chan et al. 2022).

Vertical confinement by walls, for example in microfluidics, can also lead to 3-D streaming. Lutz et al. (2005) reported 3-D streaming near the ends and the surface of an axially confined cylinder, which converged to quasi-two-dimensional flow far away from the cylinder but close to the axially central plane. Marin et al. (2015) showed that semi-cylindrical bubble oscillations under axial confinement in a microfluidic channel also produce 3-D streaming. Measuring the trajectories of tracers using astigmatism particle-tracking velocimetry, they reported toroidal vortex structures separated by two symmetry planes parallel and perpendicular to the bubble axis. The flow was modelled by Rallabandi et al. (2015) using a superposition of Stokes solutions to approximate the effects of lateral confinement, and Volk et al. (2020) showed that such a 3-D streaming flow can trap particles based on their size. In a much more confined (Hele-Shaw-like) experimental set-up, Costalonga, Brunet & Peerhossaini (2015) quantified streaming around vibrating obstacles, concluding that the extent of confinement influenced both the structure and spatial decay of the streaming.

In this paper, we study theoretically and experimentally the effects of vertical confinement on streaming flows, focusing on the limit where vertical length scales are much smaller than lateral scales. We develop a lubrication-like theory describing 3-D streaming around cylindrical obstacles with radius much greater than the depth of cell. In the Hele-Shaw limit, our theory predicts that the flow reverses direction across the depth of the channel, which, to the best of our knowledge, is a feature not observed in less confined microchannels. We confirm this prediction directly using experiments of streaming around cylinders sandwiched in a microchannel. In addition, our theory also predicts that the streaming velocity decays as the inverse cube of the distance from the cylinder, faster than that expected from previous two-dimensional approaches. We verify this decay rate quantitatively using particle tracking measurements from experiments of streaming around cylinders with different aspect ratios at different driving frequencies.

2. Model set-up

We develop a theory for streaming flows in which the characteristic length scale along one dimension (the height of the channel) is much shorter than those along the other two dimensions. Here, we adapt ideas from lubrication theory to the perturbation framework used to analyse streaming flows generated by small-amplitude oscillations with inertia (Riley 2001). The resulting theory follows closely the ideas set out by Shih (1970), although we arrive at different results and conclusions, as we discuss later.

The set-up of the theory mirrors the experimental configuration that we discuss in detail in § 4. We consider a cylindrical obstacle that is located in the centre of a microfluidic channel with one inlet and one outlet on the top wall (figure 1). While the inlet stays exposed to the atmosphere, the outlet is connected to a piezo buzzer that oscillates the fluid in the channel. The speed of the oscillatory flow is determined by the voltage applied to the buzzer and the angular frequency of oscillation $\omega =2{\rm \pi} f$. We will assume that the channel's width and length are much greater than both the radius of the cylinder $a$ and the height of the channel $2h$, so that the lateral walls have little influence on the flow. As a non-trivial restriction, we also assume that the radius of the cylinder is much greater than the half-height of the channel ($a\gg h$).

Figure 1. Model set-up. (a) A cylinder sandwiched in the cell with a radius $a$ larger than the height of the cell $2h$. Oscillations of the fluid with angular frequency $\omega$ are driven by a piezo buzzer. (b) Side view, showing Stokes layers with thickness $\sqrt {2 \nu /\omega }$ near the walls and cylinder surface.

We describe the flow around the cylinder by introducing a (dimensional) coordinate system $(X,Y,Z)$ centred at the geometric centre of the cylinder and oriented as shown in figure 1: $X$ is along the length of channel, $Y$ is across its width, and $Z$ is across its height. We also define a cylindrical coordinate system $(R,\theta,Z)$ such that $X = R \cos \theta$ and $Y = R \sin \theta$. The buzzer applies pressure oscillations to drive an oscillatory flow through the channel. The channel-height-averaged fluid velocity $\bar {\boldsymbol {V}} =(1/(2h)) \int _{-h}^h \boldsymbol {V}\,\mbox {d}z$ far from the cylinder is directed along the $X$ axis (unit vector $\boldsymbol {e}_{X}$) and satisfies

(2.1)\begin{equation} \bar{\boldsymbol{V}} (R\rightarrow\infty)\sim \bar{V}_{\infty}\, \boldsymbol{e}_{X} \,\mathrm{e}^{\mathrm{i}\omega T}, \end{equation}

where $\bar {V}_{\infty }$ is assumed known, and $T$ represents dimensional time. Here and below, it will be understood that only the real part of any complex equation is physically relevant. Note that the condition $R \to \infty$ refers to distances far from the cylinder that are not too close to the buzzer.

2.1. Governing equation and boundary conditions

The flow satisfies the incompressible Navier–Stokes equations

(2.2)\begin{equation} \rho\left(\frac{\partial\boldsymbol{V}}{\partial T}+\boldsymbol{V}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{V} \right)={-} \boldsymbol{\nabla} P+\mu\,\nabla^2 \boldsymbol{V}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{V}=0, \end{equation}

where $\boldsymbol {V} (\boldsymbol {X},T)$ is the 3-D fluid velocity, and $P(\boldsymbol {X},T)$ is the pressure. The flow satisfies the no-slip condition ($\boldsymbol {V}=0$) at the cylinder surface and the channel walls.

We non-dimensionalize the equations by choosing $a$ as the characteristic length scale in the $XY$ plane, and $h$ as the characteristic scale in the $Z$ direction. We define dimensionless coordinates (lowercase) according to

(2.3a–d)\begin{equation} X=ax,\quad Y=ay,\quad Z=hz, \quad R=ar. \end{equation}

For later convenience, we write the velocity as $\boldsymbol {V} = \boldsymbol {V}_{\parallel } + V_Z \boldsymbol {e}_Z$, where $\boldsymbol {V}_{\parallel }$ represents velocities in the horizontal ($XY$) plane, and $V_Z$ is the vertical velocity component. The velocity scale $\bar {V}_\infty$ is characteristic of the horizontal velocity $\boldsymbol {V}_{\parallel }$. Continuity then suggests a vertical velocity scale $\bar {V}_\infty h/a$. Substituting these velocity scales into (2.2), and comparing viscous and pressure terms, suggests a pressure scale $\mu \bar {V}_{\infty } a/h^2$. We therefore define the dimensionless velocity components $\boldsymbol {v}_{\parallel }$ and $v_z$, and pressure $p$, according to

(2.4a–c)\begin{equation} \boldsymbol{V}_{{\parallel}}= \bar{V}_\infty \boldsymbol{v}_{{\parallel}},\quad V_Z=\frac{\bar{V}_\infty h}{a}\,v_z, \quad P= \frac{\mu \bar{V}_\infty a}{h^2}\,p. \end{equation}

We also define a rescaled time $t = \omega T$ and introduce the dimensionless parameters

(2.5a–c)\begin{equation} \delta = \sqrt{\frac{2\nu}{h^2 \omega}},\quad \epsilon=\frac{\bar{V}_{\infty}}{a \omega}\quad \mbox{and} \quad \eta=\frac{h}{a}, \end{equation}

where $\nu =\mu /\rho$ is the kinematic viscosity of the fluid. Here, $\epsilon$ is the dimensionless amplitude of the oscillatory flow, and $\delta$ is the ratio of the Stokes layer thickness to the channel height, which characterizes the frequency of the flow, while $\eta$ is the aspect ratio of the cylinder, which characterizes the degree of axial confinement.

The dimensionless 3-D fluid velocity is therefore $\boldsymbol {v} = \boldsymbol {V}/\bar {V}^{\infty }=\boldsymbol {v}_{\parallel } + \eta v_z \boldsymbol {e}_z$. Thus (2.2) rescales (without approximations) as

(2.6a)\begin{gather} \frac{2}{\delta^2} \left(\frac{\partial \boldsymbol{v}_{{\parallel}}}{\partial t}+\epsilon\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_{{\parallel}}\right)={-}\nabla_{{\parallel}} p+\left(\eta^2\,\nabla^2_{{\parallel}} +\frac{\partial ^2}{\partial z^2}\right) \boldsymbol{v}_{{\parallel}}, \end{gather}

(2.6b)\begin{gather}\frac{2}{\delta^2} \left(\frac{\partial v_z}{\partial t}+\epsilon \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} v_z\right)={-}\frac{1}{\eta^2}\,\frac{\partial p}{\partial z}+\left(\eta^2\, \nabla^2_{{\parallel}}+\frac{\partial ^2}{\partial z^2}\right)v_z, \end{gather}

(2.6c)\begin{gather}\nabla_{{\parallel}}\boldsymbol{\cdot}\boldsymbol{v}_{{\parallel}}+\frac{\partial v_z}{\partial z}=0, \end{gather}

where $\nabla _{\parallel }=\boldsymbol {e}_{x}({\partial }/{\partial x})+\boldsymbol {e}_{y}({\partial }/{\partial y})$ is the gradient in the horizontal plane. Using overbars to indicate the average across the height of the channel (for a dimensionless function $f(z)$, $\bar {f} \equiv (1/2) \int _{-1}^1f(z)\,\mbox {d}z$), the height-averaged velocity at infinity satisfies $\bar {\boldsymbol {v}}(r \rightarrow \infty )=\boldsymbol {e}_{x}\,\mathrm {e}^{{\rm i}t}$. The flow also satisfies no-slip conditions on the cylinder surface and channel walls.

3. Streaming theory for narrow channels

We focus on the limit where the channel is narrow ($h \ll a$), so that $\eta \ll 1$. In this limit, the domain of the flow can be split into an outer region where distances from the cylinder surface are much greater than the channel height ($r - 1 \gg \eta$) and an inner region (Shih 1970; Balsa 1998) of thickness comparable to the channel height surrounding the cylinder surface (where $r - 1 = O(\eta )$). As we show below, the flow in the outer region can only approximately satisfy the boundary conditions on the cylinder. The role of the inner region is to compensate for the outer flow and meet the no-slip condition on the cylinder surface. Here, we consider the limit $\eta \ll 1$ (so the inner region is asymptotically thin) and focus mainly on the flow in the outer region, although we discuss some details of the inner region in Appendix A.

In the outer region, viscous stresses are dominated by shear across the height of the channel, and (2.6) are approximately

(3.1a)\begin{gather} \frac{2}{\delta^2}\left(\frac{\partial \boldsymbol{v}_{{\parallel}}}{\partial t}+\epsilon \boldsymbol{v} \boldsymbol{\cdot} \nabla_{{\parallel}}\boldsymbol{v}_{{\parallel}}\right) ={-}\nabla_{{\parallel}} p+\frac{\partial^2\boldsymbol{v}_{{\parallel}}}{\partial z^2}, \end{gather}

(3.1b)\begin{gather}\frac{\partial p}{\partial z}=0, \end{gather}

(3.1c)\begin{gather}\nabla_{{\parallel}}\boldsymbol{\cdot}\boldsymbol{v}_{{\parallel}}+\frac{\partial v_z}{\partial z}=0. \end{gather}

We find that the outer flow satisfies the effective boundary condition of zero channel-height-averaged velocity normal to the cylinder surface, up to corrections of $O(\eta )$; see Appendix A. Denoting the unit vector normal to the cylinder surface by $\boldsymbol {n}$, (3.1) are therefore subject to the boundary conditions

(3.2a)\begin{gather} \bar{\boldsymbol{v}} \rightarrow \boldsymbol{e}_x \, \mathrm{e}^{{\rm i}t} \quad \mbox{as} \ r\rightarrow\infty, \end{gather}

(3.2b)\begin{gather}\boldsymbol{v} =\boldsymbol{0}\quad \mbox{at} \ z={\pm} 1\ \mbox{(top and bottom walls)}, \end{gather}

(3.2c)\begin{gather}\bar{\boldsymbol{v}}\boldsymbol{\cdot} \boldsymbol{n}=0 \quad \mbox{at}\ r=1 \ \mbox{(cylinder surface)}. \end{gather}

We note that the neglected $\eta ^2\,\nabla _{\parallel }^2 \boldsymbol {v}$ viscous terms produce only subdominant corrections of $O(\eta ^2)$ to both oscillatory and steady components of the outer flow.

3.1. Small-amplitude approximation and outer flow

We now solve (3.1) and (3.2) to obtain the outer flow in the practically relevant limit of $\epsilon \ll 1$, but with arbitrary $\delta$. As is typical in many streaming calculations (Riley 2001), we assume that the oscillation amplitude is smaller than other geometric scales, so $\epsilon \ll \delta$ and $\epsilon \ll \eta \ll 1$. For small $\epsilon$, we develop a solution by applying a perturbation expansion

(3.3)\begin{equation} (p, \boldsymbol{v})=(p_1, \boldsymbol{v}_1)+\epsilon\,(p_2, \boldsymbol{v}_2)+\cdots. \end{equation}

The subscript $1$ identifies the primary flow (oscillatory flow), whereas the subscript $2$ corresponds to the secondary steady flow (streaming). The above series expansion, when substituted into (3.1) and separating orders of $\epsilon$, yields the governing equations for the primary and secondary flow.

3.1.1. Primary flow

The primary flow at $O(\epsilon ^0)$ is governed by

(3.4) \begin{equation} \left.\begin{array}{c@{}} \dfrac{2}{\delta^2}\,\dfrac{\partial \boldsymbol{v}_{1 \parallel}}{\partial t}={-}\nabla_\parallel p_1+\dfrac{\partial^2 \boldsymbol{v}_{1 \parallel}}{\partial z^2}, \quad \dfrac{\partial p_1}{\partial z}=0, \quad \nabla_{{\parallel}}\boldsymbol{\cdot}\boldsymbol{v}_{1 \parallel}+\dfrac{\partial v_{1z}}{\partial z}=0,\\ \mbox{with}\quad\bar{\boldsymbol{v}}_1|_{r\rightarrow\infty} =\boldsymbol{e}_{x}\, \mathrm{e}^{{\rm i}t},\quad \boldsymbol{v}_{1 \parallel}|_{z={\pm} 1}=\boldsymbol{0}, \quad v_{1z}|_{z={\pm} 1}=0,\quad \bar{\boldsymbol{v}}_{1 \parallel }\boldsymbol{\cdot} \boldsymbol{n}|_{r=1}=0. \end{array}\right\} \end{equation}

The system (3.4) is linear, and the only non-homogeneous condition is the one on velocity at infinity. We therefore seek separable solutions proportional to $\mathrm {e}^{\textrm {i}t}$. Observing that pressure is independent of $z$, the general expression for the primary velocity in the outer region is

(3.5)\begin{equation} \boldsymbol{v}_1={-}\frac{ \boldsymbol{\nabla} p_1}{\alpha^2} \left(1-\frac{\cosh{\alpha z}}{\cosh{\alpha}}\right),\quad\mbox{where}\ \alpha=\frac{1+\mathrm{i}}{\delta}.\end{equation}

The continuity equation, along with the normal velocity condition at $r = 1$ (see (3.4)), yields the primary pressure as

(3.6)\begin{equation} p_1= G(\alpha)\left(r+\frac{1}{r}\right)\cos{\theta}\,\mathrm{e}^{{\rm i}t}, \quad\mbox{where}\ G(\alpha) ={-}\frac{\alpha^3}{\alpha-\tanh{(\alpha})},\end{equation}

which agrees with Shih (1970) up to scaling. Note that the result (3.5) is valid in the outer region independent of the details of the primary pressure $p_1$. By contrast, the expression (3.6) explicitly ignores the presence of the inner region. In Appendix A, we will see that a matching condition with the inner region leads to $O(\eta )$ corrections of the outer pressure in (3.4), while leaving (3.5) unchanged. Observe that the primary flow has no vertical velocity component ($v_{1z}=0$), and that the primary pressure varies only in the $r\theta$ plane. The channel-averaged primary (oscillatory) velocity is

(3.7)\begin{equation} \bar{\boldsymbol{v}}_1 = \frac{ \boldsymbol{\nabla} p_1}{G(\alpha)} = \boldsymbol{\nabla} \left[\left(r+\frac{1}{r}\right)\cos{\theta}\right] \mathrm{e}^{{\rm i}t}, \end{equation}

and is used in the calculation of the streaming flow below.

3.1.2. Secondary flow

The convective acceleration term in (3.1a), which is neglected when solving for the primary flow, appears as a body force that drives the secondary flow. It will be understood that all secondary quantities are time-averaged. We assume that $\epsilon \ll \delta$, so that the inertia of the time-averaged secondary flow is negligible Riley (2001), and the governing equations at $O(\epsilon ^1)$ are

(3.8a)\begin{align} \frac{2}{\delta^2}\left<\boldsymbol{v}_{1 \parallel}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}_{1 \parallel}\right>&={-}\nabla_\parallel p_2+\frac{\partial^2 \boldsymbol{v}_{2 \parallel}}{\partial z^2}, \quad \frac{\partial p_2}{\partial z}=0, \quad \nabla_{{\parallel}}\boldsymbol{\cdot}\boldsymbol{v}_{2 \parallel}+\frac{\partial v_{2z}}{\partial z}=0, \end{align}

(3.8b)\begin{align} \mbox{with}\quad p_2|_{r\rightarrow\infty} &= 0,\quad \boldsymbol{v}_{2 \parallel}|_{z={\pm} 1}=\boldsymbol{0}, \quad v_{2z}|_{z={\pm} 1}=0,\quad\bar{\boldsymbol{v}}_{2 \parallel }\boldsymbol{\cdot} \boldsymbol{n}|_{r=1}=0. \end{align}

Here, $\langle {\cdot } \rangle$ denotes the time average over an oscillation cycle. We note the useful averaging rule $\langle \textrm {Re} [ A \,\mathrm {e}^{\mathrm {i} t}]\,\textrm {Re} [ B \,\mathrm {e}^{\mathrm {i} t}]\rangle = (1/2)\,\textrm {Re}[A^* B]$ for any complex quantities $A$ and $B$, where the asterisk denotes the complex conjugate (Longuet-Higgins 1998). We then use (3.5) to evaluate the advective term $\langle \boldsymbol {v}_{1 \parallel }\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_{1 \parallel }\rangle = (1/2) \boldsymbol {v}_{1 \parallel }^*\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_{1 \parallel }$. As mentioned before, only the real part is physically relevant. Integrating the horizontal momentum equation in (3.8a) with respect to $z$, applying no-slip conditions and noting that $\alpha ^*=(1-\mathrm {i})/\delta = -\mathrm {i}\alpha$, we obtain

(3.9a)\begin{equation} \boldsymbol{v}_{2 \parallel}=\tfrac{1}{2}(z^2-1)\, \boldsymbol{\nabla} p_2+ g(z,\alpha)\, \boldsymbol{\nabla} | \boldsymbol{\nabla} p_1|^2, \end{equation}

where

(3.9b)\begin{equation} g(z,\alpha)=\frac{\mathrm{i}}{8\alpha^4}\left((z^2-1) \alpha^2+\frac{2\cos{\alpha z}}{\cos \alpha}-\frac{2\cosh\alpha z}{\cosh\alpha}+\frac{\sin\alpha z \sinh\alpha z}{\cos\alpha \cosh\alpha}-\tan\alpha \tanh\alpha\right). \end{equation}

Next, we integrate the continuity equation in (3.8a), and apply boundary conditions at the top and bottom walls, to find that the horizontal flux is divergence-free: ${\boldsymbol {\nabla } \boldsymbol {\cdot }\bar {\boldsymbol {v}}_{2 \parallel }=0}$. We then substitute (3.9) into the above relation to find that the secondary pressure is governed by

(3.10a)\begin{equation} \nabla_{{\parallel}}^2 {p_2} = 3 \bar{g}\, \nabla_{{\parallel}}^2| \boldsymbol{\nabla} p_1|^2, \end{equation}

subject to

(3.10b)\begin{equation} p_2|_{r \to \infty} = 0 \quad \mbox{and} \quad \boldsymbol{\nabla} (p_2 - 3 \bar{g}\,| \boldsymbol{\nabla} p_1|^2) \boldsymbol{\cdot} \boldsymbol{n}|_{r = 1} = 0. \end{equation}

Note that the condition at $r = 1$ corresponds to no penetration of the depth-averaged flow into the cylinder surface; see (3.8b). Solving (3.10), we obtain the secondary pressure as

(3.11a)\begin{align} p_2&=3\bar{g} (| \boldsymbol{\nabla} p_1|^2-|G|^2) \end{align}

(3.11b)\begin{align} &={-}\frac{\mathrm{i}}{16\alpha^5}\,(4 \alpha^3 + 15 (\tanh \alpha - \tan \alpha) +6\alpha\tan\alpha\,\tanh{\alpha})(| \boldsymbol{\nabla} p_1|^2-|G|^2), \end{align}

where we note that $|G|^2 = | \boldsymbol {\nabla } p_1|^2|_{r \to \infty }$ is a constant.

Substituting the expression above into (3.9) yields $\boldsymbol {v}_{2 \parallel }$, and using the continuity equation yields the vertical component of secondary flow $v_{2z}$. Observing that the secondary pressure and velocity both involve $| \boldsymbol {\nabla } p_1|^2=| G(\alpha ) |^2 |\bar {\boldsymbol {v}}_1|^2$ (see (3.7)), we obtain the secondary (time-averaged) velocity

(3.12a)\begin{equation} \boldsymbol{v}_2 = \boldsymbol{v}_{2 \parallel} + \eta v_{2z}\, \boldsymbol{e}_z = F_{{\parallel}}(z,\alpha)\, \boldsymbol{\nabla} (|\bar{\boldsymbol{v}}_1|^2 ) + \eta\,F_z(z, \alpha)\,\nabla^2 (|\bar{\boldsymbol{v}}_1|^2 )\, \boldsymbol{e}_z, \end{equation}

where

(3.12b)\begin{align} F_{{\parallel}}(z,\alpha)&=\frac{\mathrm{i}\, | G(\alpha) |^2}{32\alpha^5}\left(\frac{8\alpha \cos{\alpha z}}{\cos{\alpha}}-\frac{8\alpha \cosh{\alpha z}}{\cosh{\alpha}}+\frac{4\alpha \sin{\alpha z}\,\sinh{\alpha z}}{\cos{\alpha} \cosh{\alpha}}\right.\nonumber\\ &\quad \left. \vphantom{\frac{8\alpha \cosh{\alpha z}}{\cosh{\alpha}}}+ 15(z^2-1) (\tan \alpha -\tanh \alpha)- 2(3 z^2-1) \alpha\tan\alpha \tanh{\alpha}\right), \end{align}

(3.12c)\begin{align} F_z(z,\alpha)&={-}\frac{\mathrm{i}\,| G(\alpha) |^2}{32\alpha^5} \left(\frac{8\sin{\alpha z}}{\cos{\alpha}}- \frac{8\sinh{\alpha z}}{\cosh{\alpha}} +\frac{2(\sin{\alpha z} \cosh{\alpha z}- \cos\alpha z \sinh{\alpha z})}{\cos{\alpha}\,\cosh{\alpha}}\right.\nonumber\\ &\quad \left.\vphantom{\frac{8\alpha \cosh{\alpha z}}{\cosh{\alpha}}}+5z(z^2-3)(\tan \alpha- \tanh \alpha) - 2 z (z^2-1) \alpha \tan \alpha \tanh \alpha\right). \end{align}

The functions $F_{\parallel }(z,\alpha )$ and $F_z(z,\alpha )$ describe horizontal and vertical velocity profiles across the height of the channel, while horizontal variations are captured by $|\bar {\boldsymbol {v}}_1|^2$. We observe that the secondary flow has zero average across the channel everywhere in the flow, not just at the cylinder surface ($\overline {F}_{\parallel } = 0 \implies \bar {\boldsymbol {v}}_2 = \boldsymbol {0}$).

While $\boldsymbol {v}_2$ describes the Eulerian secondary flow, it is the time-averaged motion of material fluid elements (the Lagrangian secondary flow) that is practically relevant and is normally referred to as streaming. The time-averaged Lagrangian velocity is $\boldsymbol {v}_L=\boldsymbol {v}_2+\boldsymbol {v}_d$, where $\boldsymbol {v}_d$ is the Stokes drift defined by $\boldsymbol {v}_d=\langle \int \boldsymbol {v}_1\,\mbox {d}t \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}_1\rangle$ (Riley 2001). Using (3.5), the Stokes drift in the outer region has the general form (again using the averaging rules for products of complex oscillating quantities)

(3.13)\begin{equation} \boldsymbol{v}_d = \frac{1}{2}\left|\frac{1}{\alpha^2} \left(1-\frac{\cosh{\alpha z}}{\cosh{\alpha}}\right) \right|^2 (\mathrm{i}\,\nabla_{{\parallel}} p_1^* \boldsymbol{\cdot} \nabla_{{\parallel}} \nabla_{{\parallel}} p_1). \end{equation}

For the $p_1$ given in (3.6), the Stokes drift vanishes identically, so the Lagrangian streaming is given simply by $\boldsymbol {v}_2$. Thus the streaming velocity has zero average across the channel height ($\bar {\boldsymbol {v}}_{L} = 0$) up to leading order in $\eta$. We show in Appendix A that the leading contribution to the channel-averaged streaming occurs at $O(\eta )$ as a result of a correction to the primary pressure $p_1$ to match the flow in the inner region.

3.2. Flow structure

The variation of the streaming velocity across the height of the channel, characterized by the functions $F_{\parallel }$ and $F_z$, is shown in figure 2 for different values of $\delta$. The horizontal flow velocity is characterized by multiple flow reversals across the channel height (changes in the sign of $F_{\parallel }$). Horizontal velocity profiles for $\delta \lesssim 0.1$ have qualitatively similar shapes. Near the top and bottom walls, the horizontal flow increases and decreases rapidly on the scale of the Stokes layer thickness $\delta$ before changing direction (figure 2a). At the outer edge of the Stokes layers the flow achieves a maximal velocity (the maximal negative value of $F_{\parallel }$), which we identify with the time-averaged slip over walls at the edge of thin Stokes layers (Longuet-Higgins 1953; Nyborg 1958). The central core outside the Stokes layers is characterized by a more gradual variation of the flow, which reverses direction another time as one moves towards the central plane ($z=0$) of the channel. For $\delta \gtrsim 0.13$, the second flow reversal near $z = 0$ disappears (see profile for $\delta = 0.2$ in figure 2a). For $\delta$ larger than approximately $0.25$, the shape of the flow remains roughly the same and is characterized by a single flow reversal with a maximum at the central plane where $z=0$. For $\delta \ll 1$, the central core (outside Stokes layers) is described by the parabolic profile with slip, $F_{\parallel } \sim (1/16)(1-3z^2)$, whereas for $\delta \gg 1$, the flow profile is described by $F_{\parallel } \sim (3 \delta ^{-2}/560) (z^2 - 1)(7 z^4 - 28 z^2 +5)$.

Figure 2. Analytical description of streaming velocity: (a) horizontal and (b) vertical velocity profiles for different values of Stokes layer thickness $\delta$. The flow reverses direction multiple times across the height of the channel depending on the value of $\delta$. (c) The horizontal flow at the central plane $z=0$ reverses direction across $\delta \approx 0.13$. Dashed curves are asymptotic expansions for small and large $\delta$.

The vertical velocity, shown in figure 2(b), exhibits odd symmetry about $z = 0$ but undergoes similar changes in sign as $\delta$ increases. When $\delta$ is small, $F_z$ is mostly negative in the top half of the channel, with a reversal in direction within the Stokes layer. By contrast, for $\delta$ larger than approximately $0.13$, the vertical velocity reverses and $F_z$ is positive throughout the top half of the channel.

Figure 2(c) shows the behaviour of the horizontal velocity at the central plane $z=0$. As $\delta$ increases, the horizontal velocity decreases and then increases rapidly to reach a maximum (negative) value at approximately $\delta =0.35$, with a direction change occurring at $\delta \approx 0.13$. The orange dashed curve represents asymptotic expansions of horizontal central plane velocity function $F_{\parallel }(z = 0) \sim (2-15\delta )/32$ for small $\delta$, whereas the blue dashed curve shows the behaviour of $F_{\parallel }(z = 0) \sim -3/(112\delta ^2)$ for large $\delta$.

In addition to showing velocity distributions, we visualize the 3-D flow using streamline portraits in two-dimensional sections. Streamline plots in different planes – ‘length–width’ (${xy}$), ‘length–depth’ (${xz}$) and ‘width–depth’ (${yz}$) – are shown in figure 3, with the oscillation direction and coordinates marked.

Figure 3. Theoretical streamline plots for $\delta =1$, with oscillations driven along $x$. (a) Streamlines in the horizontal plane show a vortex structure that reverses direction across $z$ (blue and red streamlines are at $z = 0$ and $z = 0.7$, respectively). The reversal of flow is made clearer by streamlines in the vertical planes (b) and (c).

Figure 3(a) shows streamlines around a cylinder (shaded grey) on the horizontal plane but at different depths for $\delta = 1$. Blue streamlines represent the flow at the axial mid-plane ($z=0$), showing a structure qualitatively similar to two-dimensional streaming. However, streamlines at a depths $z=0.7$, shown in red, indicate flow in the opposite direction. This reversal of flow direction (associated with changes in the sign of $F_{\parallel }$) is a feature that is absent not only from two-dimensional theory (which ignores variations along $z$), but also from previous theory (Shih 1970) and experiments (Costalonga et al. 2015) in the Hele-Shaw limit. It is interesting to note that previous descriptions of 3-D streaming flows under less extreme vertical confinement (of rigid cylinders and semi-cylindrical bubbles) do not report such a flow reversal either (Lutz et al. 2005; Marin et al. 2015). In § 4, we confirm that such a flow reversal with $z$ does indeed occur in experiments in Hele-Shaw geometries.

Figures 3(b) and 3(c) indicate flows in vertical planes and also show that the direction of secondary flow reverses across the depth of the channel. Unlike figure 3(b), we see in figure 3(c) that the flow appears to locally penetrate the cylinder surface, even though the depth-averaged velocity normal to the cylinder surface remains identically zero by construction. This is understood by recalling that the present theory is restricted to the ‘outer region’ where $r - 1 \gg \eta$. To accommodate the local velocity normal to the cylinder at $r=1$, it is necessary to acknowledge the presence of an inner region of width $\eta$ around the cylinder that ensures that the streaming velocity vanishes everywhere on the cylinder surface. Though we do not resolve the details of this inner region here, we provide some estimates of its effects in Appendix A.

4. Experiments

We now show that experiments for streaming around a cylinder sandwiched in a Hele-Shaw channel support our theoretical prediction that the direction of flow reverses across the height of cell. We fabricate devices (see figure 1) out of polydimethylsiloxane (PDMS, Dow Sylgard 184) to obtain transparent microchannels with height $2h$ at 160 $\mathrm {\mu }$m or 480 $\mathrm {\mu }$m, containing cylinders with radius $a=600$ $\mathrm {\mu }$m or 750 $\mathrm {\mu }$m inside. The PDMS mixture is poured onto moulds (made from layers of backing tape) and peeled after curing on a hot plate. Deionized water ($\nu = 10^{-6}$ m$^2$ s$^{-1}$, density 1.00 g cm$^{-3}$) with polystyrene microparticles (Magsphere, diameter $2 a_p = 5.0\ \mathrm {\mu }$m, density 1.05 g cm$^{-3}$) is injected into the microchannel using a syringe, which is disconnected once the microchannel is full, exposing the inlet to the atmosphere (see figure 1). The piezo buzzer connected to the other end (outlet) is then driven by applying sinusoidal signals of various frequencies from a function generator (Rigol DG4062) amplified through a power amplifier (Krohn-hite 7500). The buzzer produces pressure oscillations to drive oscillatory flow in the channel, which in turn results in streaming. A high-speed camera (Photron, Fastcam Nova S6) records images through an inverted microscope (Leica DMi8). The frame rate of the camera is chosen to be a integer divisor of the driving frequency so as to isolate the steady component of the particle motion.

4.1. Visualization

Matlab and ImageJ are used to process the recorded images of the microparticles as they are transported by the flow. The microparticles are expected to behave as passive tracers of the flow due to their small size ($a_p/a \lesssim 0.005$, $a_p/h \lesssim 0.016$) and inertia (Stokes number $a_p^2 \omega /\nu \lesssim 0.01$). We extract particle traces from the experimental images by using an in-house implementation of the Flowtrace algorithm (Gilpin, Prakash & Prakash 2017). Figure 4, with the oscillation direction represented by pink arrows, indicates steady pathlines of tracers around a cylinder of radius $a = 600\ \mathrm {\mu }$m, but with different values of $h$ and driving frequency $f$. The flow structure is robust across our experiments: four steady vortices are created adjacent to the cylinder. Figure 4 shows these vortices for two representative experiments (one pair of vortices per experiment). The flow structure projected in the imaging $xy$ plane resembles two-dimensional streaming around a cylinder in an unbounded or concentrically bounded fluid (Bertelsen et al. 1973), and is qualitatively similar to the theoretical streamline portrait of figure 3(a).

Figure 4. Pathlines of tracer particles around cylinders of radius $a = 600\ \mathrm {\mu }$m for different experimental conditions; double arrows indicate the oscillation direction. Steady pathlines for two experiments with different values of channel height $h$ and frequency $f$. The qualitative flow pattern is robust across our experiments. (a) $h = 240\ {\rm \mu}{\rm m},\ f = 100\ {\rm Hz}\ (\eta = 0.4,\ \delta = 0.23)$; (b) $h = 80\ {\rm \mu}{\rm m},\ f = 300\ {\rm Hz}\ (\eta = 0.13,\ \delta = 0.41)$.

A closer look at the trajectories for small $\eta$, however, reveals a striking difference from two-dimensional streaming. Figure 5(a) shows trajectories of particles in an experiment with $\eta \approx 0.13$ and $\delta \approx 0.41$ (cf. figure 4b), digitally processed such that the brighter part of each trace represents the current position of a particle (head of the trace), and the darker part represents the past position of the same particle (tail of the trace); see Gilpin et al. (2017). This leads to the appearance of thicker heads and thinner tails of the particle tracers (note that this is not an optical effect but a digital one that accounts for temporal information). We observe pairs of apparently nearby tracers that are clearly travelling in opposite directions. We note that this feature is not localized and can be observed everywhere in the flow; see the supplementary movie available at https://doi.org/10.1017/jfm.2023.276. Zooming into one of the marked regions (insets in figure 5(a)) makes the counter-motion of nearby particles clearer. The theory rationalizes this observation as the consequence of the particles being at different $z$ locations in the channel with a reversed flow direction, as shown in figure 5(b); see also figure 3.

Figure 5. Streaming in a Hele-Shaw cell with $a=600\ \mathrm {\mu }$m, $2h=160 \mathrm {\mu }$m at $f = 300$ Hz ($\eta = 0.13$, $\delta =0.41$). (a) Particle traces, now including temporal information: for each trace, the brighter ‘head’ shows the particle's present position, whereas the darker ‘tail’ represents the particle's past position. The insets detail different parts of the flow where pairs of tracer particles are travelling in opposite directions, but at different depths (yellow arrows indicate the direction of motion to guide the eye). (b) Theoretically predicted streamlines near the central plane ($z=0.19$, blue dashed) and near the boundary wall ($z=0.95$, red solid). (c,d) Experiments with the same parameters as (a) but isolating particles at specific $z$ depths (see main text). Particles at $r\approx 1.6$, $\theta \approx 30^{\circ }$ but at different depths, (c) $z=0.19$ and (d) $z=0.95$, move in opposite directions, consistent with the theoretical streamline portraits in (b).

Figures 5(c) and 5(d) provide direct experimental confirmation of this reversal of flow with $z$. These figures correspond to a single experiment conducted with the same parameters as in figure 5(a) but observed under a higher magnification to obtain a shallower depth of focus, which enables the isolation of particles at different $z$ locations. We also utilize somewhat larger tracer particles ($2 a_p = 10.0\ \mathrm {\mu }$m, density 1.05 g cm$^{-3}$) to aid visualization. The two particles in figures 5(c) and 5(d) are captured at positions $r\approx 1.6$, $\theta \approx 30^{\circ }$ but at different locations along the channel height, $z=0.19$ and $0.95$, respectively. The direction of flow is clearly different at the two $z$ locations, and is consistent with the directions predicted theoretically in figure 5(b) at those $z$ locations. Furthermore, we see experimentally that the flow is faster at $z = 0.19$ (near the channel centre) than at $z = 0.95$ (near the channel wall) by a factor ${\approx }2.8$, which is also consistent with the theory, which predicts a factor ${\approx }2.4$. Thus both the vortex structure and direction of flow are consistent between experiments and theory for narrow channels.

4.2. Velocity decay

Our streaming theory predicts a 3-D velocity field $\boldsymbol {v}_2=\boldsymbol {v}_{2 \parallel }+\eta v_{2z}\, {\boldsymbol {e}}_z$. Substituting $|\bar {\boldsymbol {v}}_1|^2=1+1/r^4-2 \cos {(2 \theta )}/r^2$ from (3.5) into (3.12), we obtain

(4.1)\begin{equation} \boldsymbol{v}_2=(v_{r2}, v_{\theta 2}, \eta v_{2z}) = \left(\frac{4 F_{{\parallel}}}{r^5}\,(r^2\cos{2\theta}-1), \frac{4F_{{\parallel}}}{r^3}\sin{2\theta}, \frac{16 \eta F_z}{r^6}\right). \end{equation}

As is typical for streaming flows, the flow speed decays with the distance from the obstacle. Two-dimensional theory predicts that the velocity decays with the square of the distance to the centre of the cylinder: $|\boldsymbol {v}_2| \propto r^{-2}$ (Bertelsen et al. 1973). However, observe from the above solution that horizontal velocities are predicted to have a dominant far-field decay of $|\boldsymbol {v}_2| \propto r^{-3}$.

We analyse experimental streaming velocities using (an adaptation by Blair & Dufresne (2008) of) the particle-tracking algorithms of Crocker & Grier (1996). From the experimental videos, we calculate the speeds of microparticles, which we approximate by the magnitude of the horizontal velocity, as the vertical velocity is smaller. The results of our analysis are shown in figure 6. We first identify particles in a narrow sector, concentric with the cylinder, of angular width $2\delta \theta$ about a mean angle $\theta$ (we pick $\theta = {\rm \pi}/4$, $\delta \theta =0.1$). We subdivide this sector into multiple radial sections (typically 30) starting from the cylinder surface. We collect hundreds of particle trajectories in the sector of interest, and calculate their positions and speeds as they pass through the angle $\theta$ (light grey markers in figure 6). We then average the particle positions and speeds within each radial section of the sector, which we identify with the mean speed at the projected averaged locations of the particles $v(r, \theta )$ (coloured markers in figure 6). Deviations from the average thus represent the spread of the individual particle positions and speeds. The maximum (bin-averaged) speed for each experiment is denoted $v_{m}$, and its radial location is denoted $r_{m}$.

Figure 6. Decay of normalized horizontal speed $v/v_{m}$ with dimensionless distance from the centre of the cylinder $r = R/a$, showing experiments (symbols) and the theoretical prediction (dashed line). The maximal speed for each experiment is $v_{m}$. Light grey symbols are velocity data from individual particles at $\theta = {\rm \pi}/4 \pm 0.1$, whereas coloured symbols are data after a local averaging. (a) Data for $v/v_{m}$ from different experiments decay approximately as $r^{-3}$ but do not overlap due to different aspect ratio $\eta$ and Stokes layer thickness $\delta$. (b) Rescaling velocities as suggested by the present theory causes data from different experiments to collapse onto a universal curve $v/(v_{m} r_{m}^3) \approx 1.1 r^{-3}$. The inset shows the same data on a linear scale.

In figure 6, we plot the ratio of (bin-averaged) speed in each segment and maximal speed among all particles, as a function of the dimensionless distance to the cylinder axis. The experimental data are collected from three characteristic experiments with cylinders confined axially in microchannels: (i) with height $2h=480\ \mathrm {\mu }$m and cylinder radius $a=600\ \mathrm {\mu }$m, operating at frequency $f = 100$ Hz ($\eta =0.4$, $\delta = 0.23$; see panel (a) of figure 4); (ii) $2h=480\ \mathrm {\mu }$m, $a=750\ \mathrm {\mu }$m, driven at $f = 50$ Hz ($\eta =0.32$, $\delta = 0.33$); (iii) $2h=160\ \mathrm {\mu }$m, $a=600\ \mathrm {\mu }$m, operating at $f = 300$ Hz ($\eta =0.13$, $\delta = 0.4$; see panel (b) of figure 4 and figure 5a). As shown in figure 5(a) and discussed earlier, the large depth of field picks up particles at different depths, encompassing a wide range of speeds (including speeds close to zero at depths where flow reverses); our data set includes all of these trajectories.

Figure 6(a) shows the experimental data for $v/v_{m}$ versus $r$. As our theory predicts, the streaming velocity is found to decay as the inverse cube of the distance from the cylinder ($|\boldsymbol {v}_2|\propto r^{-3}$). Parameters $\epsilon$, $\eta$ and $\delta$ vary from experiment to experiment, resulting in different $v_{m}$ and $r_{m}$. We observe that the theory predicts $v_m = r_m = 1$ at $\eta \rightarrow 0$, which is consistent with the experiment for the smallest $\eta$.

Although the theory is formally accurate only for small $\eta$, we find that it nonetheless provides a useful way to organize all of the data onto a universal curve; see figure 6(b). We assume that the maximum binned mean speed from experimental data exists in the outer region, even for moderate $\eta$. The theory predicts a velocity decay behaviour $v\propto r^{-3}$ (with a prefactor that generally depends on $\delta$ and $\eta$), so we anticipate that the maximal velocity satisfies $v_{m}\propto r^{-3}_{m}$, with the same prefactor. Rearranging the equation yields a universal curve $v/(v_{m} r^3_{m}) = 1/r^3$. Applying this relationship, we find that the rescaled experimental data in figure 6(b) closely follow this universal curve across two orders of magnitude, though the agreement is better with an additional prefactor $1.1$.

5. Discussion and conclusions

In this paper, we have developed a 3-D streaming theory around cylindrical obstacles confined axially in Hele-Shaw-like microchannels. We exploited the separation between vertical and horizontal length scales, and utilized lubrication theory under the conditions of small oscillation amplitudes but retaining the inertia of the flow. After resolving the primary (oscillatory) and secondary (steady) flow, we investigated the 3-D streaming structure around the cylinder, which has two unique characteristics compared with previous work in less confined set-ups: (i) the flow direction varies across the depth of the channel; and (ii) the streaming velocity decays with the inverse cube of the distance from the cylinder centre. These findings were validated experimentally by driving oscillations around cylinders of different radii sandwiched inside PDMS microchannels of various aspect ratios. We verified the vertical reversal of the horizontal flow by analysing the trajectories of suspended tracer microparticles and the decay of flow velocity by particle-trackingvelocimetry.

As noted in § 3, the theory developed here focuses on flow in an ‘outer region’ wherein distances from the cylinder surface are much greater than the channel height. This is only part of the picture, as a separate ‘inner region’ surrounding the cylinder surface (and with radial width comparable to the channel height) must also exist for the flow to satisfy the no-slip conditions at the cylinder surface (Balsa 1998). Within the inner region where $r-1=O(\eta )$, radial variations of the flow are as important as axial ones. In Appendix A, we analyse the variation of velocity within this inner region, which in turn produces a small $O(\eta )$ normal velocity that must be met by the outer flow as a matching condition at $r = 1$. We anticipate that the correction to the secondary flow is of similar size and can be computed similarly, though the details are significantly more complex and are left to future work. We find (Appendix A) that one important qualitative feature of accounting for the inner region is the introduction of steady flow components of $O(\eta )$ with non-zero channel-averaged velocity (recall that the steady flow at $O(\eta ^0)$ has zero channel average). Such a channel-averaged flow is likely to be important for transport, in particular when $h/a$ is moderately large, but is beyond the focus of this work.

Although our experiments confirm that the flow does change direction across the channel height, we do not resolve the details of the flow reversal here. Future work may aim to quantify the associated flow structures using 3-D particle tracking, with either defocusing methods (Cierpka et al. 2010; Barnkob, Kähler & Rossi 2015; Marin et al. 2015) or multi-camera set-ups (Klank et al. 2002; Lindken, Westerweel & Wieneke 2006).

Previous work on 3-D streaming focused on relatively deep channels (Marin et al. 2015; Rallabandi et al. 2015; Volk et al. 2020) and reported systematic small axial displacements following quasi-planar orbits on a toroidal surface. By contrast, the streaming that we have studied in relatively shallow channels exhibits a qualitatively different structure involving an additional axial recirculation that is necessary to maintain the near-zero depth-averaged flow predicted by the theory. The ability to generate such recirculating 3-D flows in microconfined geometries is promising for particle manipulation and micromixing applications. For a much wider range of microstreaming set-ups, the tools developed here will be helpful to assess, model, and either minimize or enhance 3-D flow effects in practical applications.