1. Introduction
Particles in oscillatory flows and acoustic fields exhibit complex behaviour due to an interplay between inertial and viscous effects. Fluid oscillations around particles lead to time-averaged flows, known as streaming (Riley Reference Riley2001). When the flow is spatially non-uniform, particles may also experience time-averaged secondary radiation forces (King Reference King1934; Gor’kov Reference Gor’kov1962; Settnes & Bruus Reference Settnes and Bruus2012), which is a phenomenon distinct from streaming but ultimately arising from fluid inertia. The combination of these effects can lead to a steady drift of suspended particles over many oscillation cycles.
These phenomena have been used for different applications, such as particle sorting (Friend & Yeo Reference Friend and Yeo2011) and focusing in microfluidics (Mutlu, Edd & Toner Reference Mutlu, Edd and Toner2018; Rufo et al. Reference Rufo, Cai, Friend, Wiklund and Huang2022; Yang et al. Reference Yang2022). They are also useful in the levitation of particles and droplets (Foresti & Poulikakos Reference Foresti and Poulikakos2014; Andrade, Marzo & Adamowski Reference Andrade, Marzo and Adamowski2020). The levitation of particles also makes it easier to observe phenomena which are hard to study due to gravity; for example, Lee et al. (Reference Lee, James, Waitukaitis and Jaeger2018) use acoustic levitation to study electrostatic charging of fine particles due to collisions. Radiation forces have been used in the non-contact extraction and manipulation of droplets from liquid–liquid interfaces (Lirette, Mobley & Zhang Reference Lirette, Mobley and Zhang2019).
Many examples involve multiple particles, and interactions between these particles are often central to the observed behaviours. Experimental studies show that suspensions of particles exposed to oscillations demonstrate collective motion, organising into chains (Klotsa et al. Reference Klotsa, Swift, Bowley and King2009; van Overveld, Clercx & Duran-Matute Reference van Overveld, Clercx and Duran-Matute2023), lattices with long-range attraction and short-range repulsion (Voth et al. Reference Voth, Bigger, Buckley, Losert, Brenner, Stone and Gollub2002), or clusters (Sazhin et al. Reference Sazhin, Shakked, Sobolev and Katoshevski2008; Lim et al. Reference Lim, Souslov, Vitelli and Jaeger2019).
While reduced-order models have been proposed (Voth et al. Reference Voth, Bigger, Buckley, Losert, Brenner, Stone and Gollub2002; van Overveld et al. Reference van Overveld, Ellenbroek, Meijer, Clercx and Duran-Matute2024), a quantitative understanding from first principles remains challenging as these phenomena are controlled by nonlinear inertial flow features across multiple time and length scales. Recent computational work on the topic has focused on understanding the interactions between a pair of particles. Direct numerical simulations (DNS) resolving all length and time scales have provided insight into particle interactions (Klotsa et al. Reference Klotsa, Swift, Bowley and King2007; van Overveld et al. Reference van Overveld, Shajahan, Breugem, Clercx and Duran-Matute2022a , Reference van Overveld, Shajahan, Breugem, Clercx and Duran-Matuteb ; Kleischmann et al. Reference Kleischmann, Luzzatto-Fegiz, Meiburg and Vowinckel2024), but can be computationally expensive. An alternative approach is to develop a perturbation expansion in small oscillation amplitudes. Although the reduced system of equations still typically requires numerical solutions (Ingber & Vorobieff Reference Ingber and Vorobieff2013; Fabre et al. Reference Fabre, Jalal, Leontini and Manasseh2017), it yields insight into scaling behaviours and symmetries of the time-averaged dynamics.
The focus of the present article is to gain analytic insight into the interactions of two particles in oscillatory flow, filling the gap between past experimental and computational studies. We build on recent progress on the dynamics of single particles in oscillatory flows, accounting for both inertial and viscous effects (Zhang & Marston Reference Zhang and Marston2014; Agarwal et al. Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021, Reference Agarwal, Upadhyay, Bhosale, Gazzola and Hilgenfeldt2024; Zhang, Minten & Rallabandi Reference Zhang, Minten and Rallabandi2024). We focus on the regime of small oscillation amplitudes, which allows us to split the problem into fast (oscillatory) and slow (time-averaged) scales (§ 3). Using this decomposition, we develop a theory that combines multipole expansions for the oscillatory flow and the Lorentz reciprocal theorem (§§ 4 and 5). Together, these techniques lead to a highly efficient semi-analytic framework for the time-averaged interaction forces (§ 6), which we find to be in good agreement with numerical solutions of the small-amplitude scheme (Fabre et al. Reference Fabre, Jalal, Leontini and Manasseh2017). Our framework yields fully analytic results in the limit where the particles are sufficiently far apart so that their Stokes boundary layers do not overlap, providing quantitative insights into the scaling of forces with the inter-particle distance and oscillatory Stokes number. Section 7 extends the framework to the time-averaged drift velocity of freely suspended particles, finding good agreement with the DNS of Kleischmann et al. (Reference Kleischmann, Luzzatto-Fegiz, Meiburg and Vowinckel2024). We discuss possible extensions of the framework before concluding in § 8.
2. Problem set-up
We consider a pair of identical spherical particles of radius
$a$
, suspended in Newtonian fluid of density
$\rho$
and viscosity
$\mu$
. The background flow (defined in the absence of the particles) oscillates with a spatially uniform velocity
$\boldsymbol{V}^{\infty }(t) = V^{\infty } \boldsymbol{e}\cos\, (\omega t )$
, where
$V^{\infty }$
is the maximum speed of the flow,
$\boldsymbol{e}$
is the unit vector along the axis of flow and
$\omega$
is the angular frequency of oscillation. The centres of the particles are separated by a distance
$d \gt 2a$
and lie on an axis defined by the unit vector
$\boldsymbol{e}_{\parallel }$
pointing from particle 1 to particle 2, and oriented at an angle
$\theta$
relative to the flow axis
$\boldsymbol{e}$
(figure 1). We also define a unit vector
$\boldsymbol{e}_{\perp }$
perpendicular to the particle axis in the
$\boldsymbol{e}$
-
$\boldsymbol{e}_{\parallel }$
plane. The presence of the particles leads to oscillatory flow gradients, which in turn excite secondary ‘streaming’ flows with non-zero time average due to inertial effects. These inertial flows exert time-averaged forces on the particles, which are the quantities of interest here.
Figure 1. Sketch showing two identical particles suspended in a uniform ambient oscillatory flow
$\boldsymbol{V}^{\infty }(t)$
. Hydrodynamic interactions between the particles, due to advective nonlinearities, drive secondary time-averaged forces and drift velocities of the particles.
In most practical settings, such a pair of particles would be free to oscillate in response to the flow. In this case, the time-averaged inertial forces are balanced by viscous drag, leading to a drift of the particles over many oscillations. To understand the practically important case of mobile particles, it is convenient to first consider particles that are held stationary. We will show later that the case of mobile particles can be understood through a change in reference frame as long as the particles have identical properties.
We decompose the flow around the particles
$\boldsymbol{v}(\boldsymbol{x},t)$
into the ambient contribution
$\boldsymbol{V}^\infty (t)$
and a disturbance contribution
$\boldsymbol{v}^{d}({\boldsymbol{x}},t)$
, so that
$\boldsymbol{v}(\boldsymbol{x},t) = \boldsymbol{V}^{\infty }(t) + \boldsymbol{v}^{d}(\boldsymbol{x},t)$
. The amplitude of the flow oscillation relative to the particle radius is defined by
$\varepsilon = V^{\infty }/(a \omega )$
. The disturbance velocity is characterised by an oscillatory Stokes number
$\mathcal{S} = \omega a^2/\nu$
, which is the measure of the ratio of inertial to viscous forces over an oscillation period. The vorticity of the oscillatory flow is confined to Stokes layers of dimensionless thickness
$\delta = \sqrt {2/\mathcal{S}}$
(defined in units of particle radius) around each particle, where inertial and viscous effects are comparable.
We use the particle radius
$a$
as a characteristic length, the inverse angular frequency
$\omega ^{-1}$
as characteristic time and the ambient flow speed
$V^{\infty } = \varepsilon a \omega$
as the characteristic velocity. The flow is governed by the dimensionless Navier–Stokes equations
\begin{align} &\quad{{\mathcal{S}}}\left (\frac {\partial \boldsymbol{v}}{\partial {t}}+\varepsilon \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}\right )=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}=0,\quad \textrm {with} \end{align}
\begin{align} &\qquad\qquad\qquad\quad\boldsymbol{v}|_{S^{(1)}}=\boldsymbol{v}|_{S^{(2)}} = \boldsymbol{0}, \end{align}
\begin{align} &\boldsymbol{v}(|\boldsymbol{x}| \rightarrow \infty , t) = \boldsymbol{e} \cos t = (\boldsymbol{e}_{\parallel } \cos \theta + \boldsymbol{e}_{\perp } \sin \theta ) \cos t, \end{align}
where
$S^{(1)}$
and
$S^{(2)}$
represent the surfaces of particles
$1$
and
$2$
, respectively, and
$\boldsymbol{\sigma }=-p\boldsymbol{I}+ (\boldsymbol{\nabla }\boldsymbol{v}+\boldsymbol{\nabla }\boldsymbol{v}^{\unicode{x1D64F}} )$
is the dimensionless Newtonian stress tensor (rescaled with
$\mu V^{\infty }/a$
),
$p(\boldsymbol{x},t)$
being the dimensionless pressure. The goal of the present study is to calculate the time-averaged interaction forces between particles held stationary, and later, time-averaged drift velocities of freely suspended particles.
3. Small-amplitude theory
The full nonlinear flow is complicated and requires numerical methods to resolve. In typical applications (Voth et al. Reference Voth, Bigger, Buckley, Losert, Brenner, Stone and Gollub2002; Klotsa et al. Reference Klotsa, Swift, Bowley and King2009), the dimensionless amplitude of oscillation
$\varepsilon$
is small. To gain analytic insight in this limit, we perform a small-amplitude expansion with the parameter
$\varepsilon$
\begin{align} (\boldsymbol{v}, \boldsymbol{\sigma }) &= (\boldsymbol{v}_1, \boldsymbol{\sigma }_1) + \varepsilon \left(\boldsymbol{v}_2, \boldsymbol{\sigma }_2\right)+ O\left(\varepsilon ^2\right). \end{align}
The primary flow
$\boldsymbol{v}_1({\boldsymbol{x}}, t)$
oscillates with a frequency
$\omega$
and scales with the applied ambient flow. The secondary flow
$\boldsymbol{v}_2({\boldsymbol{x}}, t)$
involves a combination of a time-dependent component (with frequency 2
$\omega$
) and a steady (or time-averaged) component. We are ultimately interested in the time-averaged flow, since it is associated with a time-averaged force on the particles. Here and below, we use subscripts to indicate orders of
$\varepsilon$
, and superscripts to identify particles.
Substituting the above expansion into (2.1) and separating orders of
$\varepsilon$
, we obtain the governing equations for the primary and secondary flows. The primary flow satisfies
\begin{align} {\mathcal{S}}\frac {\partial \boldsymbol{v}_1}{\partial t}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }_1,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_1=0,\quad \textrm {with} \end{align}
\begin{align} \boldsymbol{v}_1|_{S^{(1)}} &= \boldsymbol{v}_1|_{S^{(2)}}=\boldsymbol{0},\quad \boldsymbol{v}_1(|\boldsymbol{x}| \rightarrow \infty , t) = (\boldsymbol{e}_{\parallel } \cos \theta + \boldsymbol{e}_{\perp } \sin \theta ) \cos t. \end{align}
Averaging (2.1) over an oscillation cycle, we see that the secondary flow is governed by
\begin{align} &{\boldsymbol{\nabla }}\boldsymbol{\cdot }\left \langle {\boldsymbol{\sigma }}_2-{\mathcal{S}}{\boldsymbol{v}}_1{\boldsymbol{v}}_1\right \rangle =\boldsymbol{0} ,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\left \langle \boldsymbol{v}_2\right \rangle = 0,\quad \textrm {with} \end{align}
\begin{align} &\left \langle \boldsymbol{v}_2\right \rangle |_{S^{(1)}} = \left \langle \boldsymbol{v}_2\right \rangle |_{S^{(2)}} = \left \langle \boldsymbol{v}_2\right \rangle (|\boldsymbol{x}| \rightarrow \infty ) = \boldsymbol{0}, \end{align}
where angle brackets define a time average over an oscillation according to
$\langle g\rangle (\boldsymbol{x}) = (2 \pi )^{-1} \int _{0}^{2\pi } g(\boldsymbol{x}, t) \textrm{d}t$
.
Since there is no net force on the system and the particles are identical, the time-averaged force on each particle is identical in magnitude but opposite in direction to the other particle; a more careful justification of this fact is given in Appendix A. Using the formulation of Doinikov (Reference Doinikov1994), the time-averaged force on particle 1, non-dimensionalised with
$\mu a \varepsilon V^{\infty }$
, is
\begin{align} \left \langle \boldsymbol{F}\right \rangle = \left \langle \boldsymbol{F}\right \rangle ^{(1)} = -\left \langle \boldsymbol{F}\right \rangle ^{(2)}= \int _{S^{(1)}}\boldsymbol{ n}\boldsymbol{\cdot }\left \langle \boldsymbol{\sigma }_2-{\mathcal{S}}\boldsymbol{v}_1\boldsymbol{v}_1\right \rangle \,{\textrm{d}S}. \end{align}
The ambient flow is uniform and oscillatory, and thus makes no contribution to the time-averaged force. Consequently, we replace the kernel of (3.4) by the disturbance stress
$\langle \boldsymbol{\sigma }_2-{\mathcal{S}}\boldsymbol{v}_1\boldsymbol{v}_1\rangle ^d$
, where
$f^d = f - f^\infty$
denotes the disturbance of any flow quantity
$f$
, with
$f^\infty$
representing that quantity in the absence of the particles.
4. Reciprocal theorem for the time-averaged force
The integrand in (3.4) involves the secondary time-averaged stress
$\langle \boldsymbol{\sigma }_2\rangle$
as well as the Reynolds stress due to the primary flow. Even in the perturbation framework, the linear equations (3.2) and (3.3) must be solved numerically (this was the approach adopted by Fabre et al. (Reference Fabre, Jalal, Leontini and Manasseh2017)), partly due to the complicated two-sphere geometry, and partly due to the potentially large separation of length scales between
$\delta a$
,
$a$
and
$d$
. While the primary flow yields to analytic approximation (§ 5), the secondary flow problem remains analytically intractable as it involves a body force that inherits the structure of
$\boldsymbol{v}_1$
.
In order to gain insight into the time-averaged force without solving the secondary problem, we utilise the Lorentz reciprocal theorem (Stone & Samuel Reference Stone and Samuel1996; Masoud & Stone Reference Masoud and Stone2019). We introduce an auxiliary flow
$(\hat {\boldsymbol{v}},\hat {\boldsymbol{\sigma }})$
, which we define as the steady Stokes flow around two spheres translating with velocity
$\pm \hat {\boldsymbol{V}}$
. These fields satisfy
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\hat {\boldsymbol{\sigma }}=\boldsymbol{0},\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\hat {\boldsymbol{v}}=0,\qquad \text{subject to}\qquad \hat {\boldsymbol{v}}|_{S^{(1)}}= -\hat {\boldsymbol{v}}|_{S^{(2)}}=\hat {\boldsymbol{V}},\quad \hat {\boldsymbol{v}}|_{r\to \infty }=\boldsymbol{0}. \end{align}
Combining the governing equations for the secondary time-averaged flow (3.3) (as noted earlier, we only consider the disturbance contribution) and the auxiliary flow (4.1), we construct a symmetry relation (Zhang et al. Reference Zhang, Minten and Rallabandi2024)
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\left \langle \boldsymbol{\sigma }_2 - {\mathcal{S}}\boldsymbol{v}_1 \boldsymbol{v}_1 \right \rangle ^d\boldsymbol{\cdot }\hat {\boldsymbol{v}} = \boldsymbol{\nabla }\boldsymbol{\cdot }\hat {\boldsymbol{\sigma }} \boldsymbol{\cdot }\left \langle \boldsymbol{v}_2\right \rangle ^d, \end{align}
which we then recast as
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\big [\!\kern-1.5pt\left\langle \boldsymbol{\sigma }_2 - {\mathcal{S}}\boldsymbol{v}_1 \boldsymbol{v}_1 \right \rangle ^d\boldsymbol{\cdot }\hat {\boldsymbol{v}} \big ]-\left \langle \boldsymbol{\sigma }_2 - {\mathcal{S}}\boldsymbol{v}_1 \boldsymbol{v}_1 \right \rangle ^d:\boldsymbol{\nabla }\hat {\boldsymbol{v}}= \boldsymbol{\nabla }\boldsymbol{\cdot }\big(\hat {\boldsymbol{\sigma }} \boldsymbol{\cdot }\left \langle \boldsymbol{v}_2\right \rangle ^d\big)-\hat {\boldsymbol{\sigma }}:\boldsymbol{\nabla }\left \langle {\boldsymbol{v}}_2\right \rangle ^d. \end{align}
Integrating over the fluid volume
$V_f$
and applying Gauss’s theorem yields
\begin{align} \begin{split} &\int _{S^{(1)}+S^{(2)} + S^{\infty }}\boldsymbol{n}\boldsymbol{\cdot }\left \langle \boldsymbol{\sigma }_2 - {\mathcal{S}}\boldsymbol{v}_1 \boldsymbol{v}_1\right \rangle ^d\boldsymbol{\cdot }\boldsymbol{\hat {v}}{\textrm{d}S} - \int _{V_f}\left \langle {\mathcal{S}}\boldsymbol{v}_1 \boldsymbol{v}_1\right \rangle ^d:\boldsymbol{\nabla }{\hat {\boldsymbol{v}}}\,\textrm{d}V\\& =\int _{S^{(1)}+S^{(2)} + S^{\infty }}\boldsymbol{n}\boldsymbol{\cdot }\hat {\boldsymbol{\sigma }}\boldsymbol{\cdot }\left \langle \boldsymbol{v}_2\right \rangle ^d{\textrm{d}S}, \end{split} \end{align}
where
$\boldsymbol{n}$
is the fluid-facing normal vector.
The integrals at the bounding surface at infinity,
$S^{\infty }$
, vanish since the secondary fields only contain disturbance quantities that decay rapidly in the far field. Noting that the auxiliary velocity is constant on the particle surfaces and using (3.4), the first term on the left-hand side of (4.4) simplifies to
$\langle \boldsymbol{F}\rangle ^{(1)}\boldsymbol{\cdot }\hat {\boldsymbol{V}} + \langle \boldsymbol{F}\rangle ^{(2)}\boldsymbol{\cdot }(-\hat {\boldsymbol{V}}) = 2\langle \boldsymbol{F}\rangle \boldsymbol{\cdot }\hat {\boldsymbol{V}}$
. Similarly, the term on the right simplifies to
$2\hat {\boldsymbol{F}}\boldsymbol{\cdot }\langle \boldsymbol{V}_2\rangle$
, where
$\hat {\boldsymbol{F}} = \int _{S^{(1)}} \boldsymbol{n} \boldsymbol{\cdot }\hat {\boldsymbol{\sigma }}\,{\textrm{d}S}$
is the force on particle 1 due to the auxiliary flow (the auxiliary force on particle 2 is equal and opposite), and
$\langle \boldsymbol{V}_2\rangle$
is the time-averaged velocity of particle 1. Note that the time-averaged velocity of particle 2 is equal and opposite due to symmetry; see Appendix A. We note that
$\langle \boldsymbol{V}_2\rangle = \boldsymbol{0}$
for stationary particles, but we retain it in the interest of generality (later we will relax the assumption of stationary particles). Thus, we obtain
\begin{align} \left \langle \boldsymbol{F}\right \rangle \boldsymbol{\cdot }\hat {\boldsymbol{V}} - \hat {\boldsymbol{F}}\boldsymbol{\cdot }\left \langle \boldsymbol{V}_2\right \rangle = \frac {{\mathcal{S}}}{2} \int _{V_f}\left \langle \boldsymbol{v}_1 \boldsymbol{v}_1\right \rangle ^d:\boldsymbol{\nabla }{\hat {\boldsymbol{v}}}\,\textrm{d}V. \end{align}
The above result allows for the computation of the secondary force on either particle along an arbitrary direction
$\hat {\boldsymbol{V}}$
, without requiring a solution to the secondary flow. However, the integrand in (4.5) decays slowly at large distances from the particles, making numerical evaluation (and later, analytic approximations) inconvenient. We therefore recast this integrand in terms of the primary vorticity
$\boldsymbol{\omega }_1 = \boldsymbol{\nabla }\times \boldsymbol{v}_1$
, which decays exponentially away from the particles on length scales of
$O(\delta )$
. To this end, we write the integrand in (4.5) according to the identity
$\boldsymbol{v}_1\boldsymbol{v}_1:\boldsymbol{\nabla }{\hat {\boldsymbol{v}}} =\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{v}_1\boldsymbol{v}_1 \boldsymbol{\cdot }\hat {\boldsymbol{v}})-(\boldsymbol{v}_1\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_1)\boldsymbol{\cdot }\hat {\boldsymbol{v}}-\boldsymbol{v}_1\boldsymbol{\cdot }\hat {\boldsymbol{v}}{}(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_1)$
, where the last term is zero due to incompressibility. We then use the identity
$\boldsymbol{v}_1\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_1=(1/2)\boldsymbol{\nabla }|\boldsymbol{v}_1^2|- (\boldsymbol{v}_1\times \boldsymbol{\omega }_1)$
. Substituting these relations into (4.5) and using Gauss’s theorem, we obtain
\begin{align} \left \langle \boldsymbol{F}\right \rangle \boldsymbol{\cdot }\hat {\boldsymbol{V}} - \hat {\boldsymbol{F}}\boldsymbol{\cdot }\left \langle \boldsymbol{V}_2\right \rangle = \frac {{\mathcal{S}}}{2} \left [\int _{V_f} \!\left \langle \boldsymbol{v}_1 \times \boldsymbol{\omega }_1 \right \rangle ^d\boldsymbol{\cdot }{\hat {\boldsymbol{v}}}\,\textrm{d}V - \int _{S^{(1)} + S^{(2)}} \boldsymbol{n} \boldsymbol{\cdot }\left \langle \!\boldsymbol{v}_1 \boldsymbol{v}_1 - \frac {|\boldsymbol{v}_1|^2}{2} \boldsymbol{I} \right \rangle ^d\!\boldsymbol{\cdot }\hat {\boldsymbol{v}} \, {\textrm{d}S}\right ]\!. \end{align}
The surface integrals go to zero due to the no-slip condition, leading to
\begin{align} \left \langle \boldsymbol{F}\right \rangle \boldsymbol{\cdot }\hat {\boldsymbol{V}} - \hat {\boldsymbol{F}}\boldsymbol{\cdot }\left \langle \boldsymbol{V}_2\right \rangle =\frac {{\mathcal{S}}}{2}\int _{V_f} \left (\boldsymbol{v}_1\times \boldsymbol{\omega }_1\right )\boldsymbol{\cdot }\boldsymbol{\hat {v}}\, \textrm{d}V. \end{align}
Equation (4.7) is an exact writing of (3.4). It directly involves the primary vorticity, which is confined to near-surface regions around the particles (Stokes boundary layers) of thickness
$\delta$
. Outside these regions, the vorticity – and thus the integrand – is exponentially small, making (4.7) much more computationally efficient than (4.5). Later we will see that the `localisation’ of the integrand to the Stokes layer makes the integral amenable to analytic approximations.
5. Primary and auxiliary flows: dual multipole expansions
We have shown that the only pieces of information needed to calculate the time-averaged force are the primary and auxiliary flow. Since exact solutions to these problems are not generally possible in two-sphere geometries, we develop approximations using dual multipole expansions.
5.1. Primary flow
Equation (3.2) governing the primary flow is linear, and its temporal structure is set by the oscillatory far field. Noting that the equation is separable in time, we represent all primary flow quantities as complex phasors
$\propto e^{it}$
, the real parts of which represent the solution to (3.2); see e.g. Agarwal et al. (Reference Agarwal, Chan, Rallabandi, Gazzola and Hilgenfeldt2021) and Zhang et al. (Reference Zhang, Minten and Rallabandi2024). We write the primary velocity as a superposition of the ambient and disturbance flows created by each particle
$\boldsymbol{v}_1= \boldsymbol{V}_1^{\infty }+\boldsymbol{v}_1^{d(1)}+\boldsymbol{v}_1^{d(2)}$
. We approximate the disturbance flow of each particle by a truncated multipole expansion, retaining terms corresponding to translating and straining modes of flow.
To construct a multipole solution, we view each particle as being immersed an ‘effective’ oscillatory background flow that is the combination of the ambient flow
$\boldsymbol{V}^\infty$
and the disturbance created by the other particle (Pozrikidis Reference Pozrikidis1992; Kim & Karrila Reference Kim and Karrila1993). The local background flow can be approximated by a Taylor expansion about the particle’s centre. At linear order, the effective oscillatory background around particle
$j$
is
$(\boldsymbol{\mathcal{V}}^{(j)} + \boldsymbol{r}_j \boldsymbol{\cdot }\boldsymbol{\mathcal{E}}^{(j)}) e^{it}$
, where
$\boldsymbol{r}_j = \boldsymbol{x} - \boldsymbol{x}_j$
is the position vector relative to the centre of particle
$j$
. The quantities
$\boldsymbol{\mathcal{V}}^{(j)}$
and
$\boldsymbol{\mathcal{E}}^{(j)}$
are the ‘effective’ velocity and rate of strain (symmetric and traceless) felt locally by the particle. The particle responds to this background by exciting a disturbance flow, which has the general structure (Agarwal et al. Reference Agarwal, Upadhyay, Bhosale, Gazzola and Hilgenfeldt2024; Zhang et al. Reference Zhang, Minten and Rallabandi2024)
\begin{align} \boldsymbol{v}_1^{d(j)}(\boldsymbol{x},t) = \left (\boldsymbol{D}(\boldsymbol{r}_j,\lambda )\boldsymbol{\cdot }\boldsymbol{\mathcal{V}}^{(j)} + \boldsymbol{Q}(\boldsymbol{r}_j,\lambda ):\boldsymbol{\mathcal{E}}^{(j)} + \ldots \right )e^{it}, \end{align}
where
$\lambda = \sqrt {i {\mathcal{S}}} = (1+i)/\delta$
. Defining
$R=\lambda r$
, the dipole
$\boldsymbol{D}$
and quadrupole
$\boldsymbol{Q}$
tensor fields are solutions of (3.2), given by (Zhang et al. Reference Zhang, Minten and Rallabandi2024)
\begin{align} \begin{split} \boldsymbol{D}(\boldsymbol{r},\lambda )&=d_1\left [\frac {\boldsymbol{I}}{r^3}-3\frac {\boldsymbol{r}\boldsymbol{r}}{r^5}\right ]\\&\quad+d_2\left [e^{-R}\left (\frac {1}{R}+\frac {1}{R^2}+\frac {1}{R^3}\right )\boldsymbol{I}-e^{-R}\left (\frac {1}{R}+\frac {3}{R^2}+\frac {3}{R^3}\right )\frac {\boldsymbol{r}\boldsymbol{r}}{r^2}\right ], \end{split} \end{align}
\begin{align} \begin{split} \boldsymbol{Q}(\boldsymbol{r},\lambda )&=q_1\left [6\frac {\boldsymbol{I}\boldsymbol{r}}{r^5}-15\frac {\boldsymbol{r}\boldsymbol{r}\boldsymbol{r}}{r^7}\right ]\\&\quad+q_2\left [e^{-R}\left (1+\frac {3}{R}+\frac {6}{R^2}+\frac {6}{R^3}\right )\frac {\boldsymbol{I}\boldsymbol{r}}{r^2}-e^{-R}\left (1+\frac {6}{R}+\frac {15}{R^2}+\frac {15}{R^3}\right )\frac {\boldsymbol{r}\boldsymbol{r}\boldsymbol{r}}{r^4}\right ], \end{split} \end{align}
with coefficients
\begin{align} d_1=-\frac {3+3\lambda +\lambda ^2}{2\lambda ^2}, d_2=\frac {3e^{\lambda }\lambda }{2}, q_1=\frac {15+15\lambda +6\lambda ^2+\lambda ^3}{9\lambda ^2(1+\lambda )}, q_2=-\frac {5e^{\lambda }\lambda }{3(1+\lambda )}, \end{align}
chosen to satisfy no slip on
$S^{(j)}$
.
As noted earlier, the ‘effective’ local velocity and strain rate
$\boldsymbol{\mathcal{V}}^{(j)}$
and
$\boldsymbol{\mathcal{E}}^{(j)}$
felt by particle
$j$
are the combined result of the oscillating background flow and the disturbance created by the other particle. We thus obtain the relations
\begin{align} \quad\boldsymbol{\mathcal{V}}^{(1)}=&\boldsymbol{V}_1^{\infty }+\boldsymbol{v}_1^{d(2)}(\boldsymbol{x}_1) ,\end{align}
\begin{align} \quad\boldsymbol{\mathcal{V}}^{(2)}=&\boldsymbol{V}_1^{\infty }+\boldsymbol{v}_1^{d(1)}(\boldsymbol{x}_2) ,\end{align}
where 
represents the symmetric and traceless part of a rank-2 tensor
$\boldsymbol{T}$
. Equation (5.4) is a linear system of equations for the effective quantities
$\boldsymbol{\mathcal{V}}^{(j)}$
and
$\boldsymbol{\mathcal{E}}^{(j)}$
, and can be solved either analytically or numerically for any configuration of the particles. On solving for
$\boldsymbol{\mathcal{V}}^{(j)}$
and
$\boldsymbol{\mathcal{E}}^{(j)}$
, the primary flow is fully determined. We then use (5.1) to calculate the primary vorticity
$\boldsymbol{\omega }_1$
. We note that we have neglected contributions due to the antisymmetric part of the velocity gradient tensor, since these terms decay exponentially away from the particles on length scales of
$\delta$
.
5.2. Auxiliary flow multipole expansion
The auxiliary flow is constructed similarly to the primary flow using a multipole expansion. We decompose the auxiliary velocity as a superposition of disturbance flows created by each particle
$\hat {\boldsymbol{v}}= \hat {\boldsymbol{v}}^{d(1)}+\hat {\boldsymbol{v}}^{d(2)}$
. A particle
$j$
translating with velocity
$\hat {\boldsymbol{V}}^{(j)}$
exposed to a linear flow
$(\hat {\boldsymbol{\mathcal{V}}}^{(j)} + \boldsymbol{r}_j \boldsymbol{\cdot }\hat {\boldsymbol{\mathcal{E}}}^{(j)})$
produces a disturbance flow
\begin{align} \hat {\boldsymbol{v}}^{(j)}= \hat {\boldsymbol{D}}(\boldsymbol{r})\boldsymbol{\cdot }\left (\hat {\boldsymbol{V}}^{(j)}-\hat {\boldsymbol{\mathcal{V}}}^{(j)}\right )+ \hat {\boldsymbol{Q}}(\boldsymbol{r}):\hat {\boldsymbol{\mathcal{E}}}^{(j)} + \cdots , \end{align}
where the tensor fields
$\hat {\boldsymbol{D}}(\boldsymbol{r})$
and
$\hat {\boldsymbol{Q}}(\boldsymbol{r})$
are (Leal Reference Leal2007)
\begin{align} \hat {\boldsymbol{D}}(\boldsymbol{r}) &= \frac {3}{4} \left [ \left ( \frac {\boldsymbol{I}}{r} + \frac {\boldsymbol{r} \boldsymbol{r}}{r^3} \right ) + \frac {1}{3} \left ( \frac {\boldsymbol{I}}{r^3} - \frac {3 \boldsymbol{r} \boldsymbol{r}}{r^5} \right ) \right ], \end{align}
\begin{align} \hat {\boldsymbol{Q}}(\boldsymbol{r}) &= -\frac {5}{2}\left [\frac {\boldsymbol{r} \boldsymbol{r}}{r^5} + \frac {1}{5} \left ( \frac {\boldsymbol{I}}{r^5} - 5 \frac {\boldsymbol{r} \boldsymbol{r}}{r^7} \right ) \right ] \boldsymbol{r}. \\[9pt] \nonumber \end{align}
The effective velocity
$\hat {\boldsymbol{\mathcal{V}}}^{(j)}$
and rate of strain
$\hat {\boldsymbol{\mathcal{E}}}^{(j)}$
are due to the disturbance created from the other particle evaluated at the centre of particle
$j$
. We also include Faxén corrections as is standard in particle hydrodynamics in Stokes flows, making use of standard results for spheres (Batchelor & Green Reference Batchelor and Green1972). Thus, the effective velocities and rates of strain are governed by (Durlofsky, Brady & Bossis Reference Durlofsky, Brady and Bossis1987; Brady & Bossis Reference Brady and Bossis1988)
\begin{align} \hat {\boldsymbol{\mathcal{V}}}^{(1)}&=\left (1+\frac {1}{6}\nabla ^2\right )\hat {{\boldsymbol{v}}}^{d(2)}(\boldsymbol{x}_1) ,\end{align}
\begin{align} \hat {\boldsymbol{\mathcal{V}}}^{(2)}&=\left (1+\frac {1}{6}\nabla ^2\right )\hat {{\boldsymbol{v}}}^{d(1)}(\boldsymbol{x}_2) ,\end{align}
It is useful to note that
$-6 \pi \mu a (\hat {\boldsymbol{V}}^{(j)}-\hat {\boldsymbol{\mathcal{V}}}^{(j)} )$
and
$(20 \pi /3) \mu a^3 \hat {\boldsymbol{\mathcal{E}}}^{(j)}$
are, respectively, the force and stresslet exerted by the auxiliary flow on the particle. Similarly to the primary flow, we solve the linear system (5.7) to calculate the effective quantities
$\hat {\boldsymbol{\mathcal{V}}}^{(j)}$
and
$\hat {\boldsymbol{\mathcal{E}}}^{(j)}$
, thereby determining the auxiliary flow. There are two auxiliary flow configurations of interest: particles moving towards each other, and particles moving normal to their connecting axis in opposite directions (figure 2). These configurations, respectively, let us calculate components of
$\left \langle \boldsymbol{F}\right \rangle$
in the
$\boldsymbol{e}_{\parallel }$
and
$\boldsymbol{e}_{\perp }$
directions, as we will discuss in the following section.
Figure 2. Auxiliary flow: particles translating with opposite velocities
$\pm \hat {\boldsymbol{V}}$
. (a) Particles moving towards each other. (b) Particles moving normal to their line of centres.
6. Time-averaged forces on stationary particles
6.1. Calculation of the force
We first consider fixed particles (
$\langle \boldsymbol{V}_2\rangle = {0}$
) and compute the time-averaged forces. In § 7 we calculate time-averaged drift velocities of freely suspended particles. Both primary and auxiliary flows are now known up to the level of force and mass dipoles and depend on the parameters
$\mathcal{S}$
,
$D = d/a$
, and
$\theta$
. We substitute these fields into the expression (4.7), where all quantities in the integrand are now analytically known. We execute the volume integral numerically since the fluid volume has a somewhat complicated geometry. However, since the primary vorticity decays exponentially outside Stokes layers, the domain of integration only needs to be a few multiples of the Stokes layer thickness
$\delta$
to achieve convergence. This procedure yields the projection of
$\left \langle \boldsymbol{F}\right \rangle$
along an arbitrarily chosen direction
$\hat {\boldsymbol{V}}$
.
To identify the general structure of the force, we appeal to tensor symmetries. The force depends on the ambient velocity
$\boldsymbol{V}^{\infty }$
(which in dimensionless units is just
$\boldsymbol{e}$
), the unit vector
$\boldsymbol{e}_{\parallel }$
connecting the particles, the oscillatory Stokes number
$\mathcal{S}$
and the distance
$D$
. Furthermore, the secondary flow is governed by a linear system (3.3), subject to the body force
$\propto \boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{v}_1 \boldsymbol{v}_1\rangle$
. As we have already seen
$\boldsymbol{v}_1(x,t)$
is linear in the ambient oscillatory flow
$\boldsymbol{V}^\infty$
, so the body force in (3.3) is quadratic in
$\boldsymbol{V}^\infty$
. It follows directly from the linearity of (3.3) that the time-averaged secondary force is quadratic in oscillatory velocity
$\boldsymbol{V}^\infty = V^\infty \boldsymbol{e}$
. The most general expression for a secondary force that satisfies this constraint is (restoring dimensions and using an inertial scale)
\begin{equation} \left \langle \boldsymbol{F}\right \rangle = \rho (V^{\infty })^2 a^2 \big [\alpha \left (\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{e}\right )\boldsymbol{e}_{\parallel }+\beta (\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{e}_{\parallel })^2\boldsymbol{e}_{\parallel }+\gamma (\boldsymbol{e}\boldsymbol{\cdot }\boldsymbol{e}_{\parallel })\boldsymbol{e}\big], \end{equation}
where
$\alpha$
,
$\beta$
and
$\gamma$
are scalar functions of
$D$
and
$\mathcal{S}$
. We substitute
$\boldsymbol{e} = \boldsymbol{e}_{\parallel } \cos \theta + \boldsymbol{e}_{\perp } \sin \theta$
to recast the above expression in terms of
$\theta$
, finding
\begin{align} \left \langle \boldsymbol{F}\right \rangle = \rho (V^{\infty })^2 a^2\big[F_{\textit{AA}} \cos ^2\theta \, \boldsymbol{e}_{\parallel } + F_{\textit{TT}} \sin ^2\theta \, \boldsymbol{e}_{\parallel } + F_{\textit{AT}}\sin {\theta }\cos \theta \, \boldsymbol{e}_{\perp }\big], \end{align}
where
$F_{\textit{AA}}= (\alpha +\beta +\gamma )$
,
$F_{\textit{TT}} = \alpha$
, and
$F_{\textit{AT}} = \gamma$
(all functions of
$\mathcal{S}$
and
$D$
). Equations (6.1) and (6.2), obtained here through straightforward tensorial symmetries, are identical to the result obtained by Fabre et al. (Reference Fabre, Jalal, Leontini and Manasseh2017) through more detailed arguments involving the spatial structure of the secondary flow field. The
$AA$
contribution is quadratic in the axial component of the ambient flow (
$\cos \theta$
), the
$TT$
contribution is quadratic in the transverse velocity component (
$\sin \theta$
), while the
$AT$
contribution involves a product of axial and transverse velocity components. We observe that
$F_{\textit{AA}}$
and
$F_{\textit{TT}}$
contributions are both associated with forces along the axis connecting the particle (positive values indicate attraction), while the
$F_{\textit{AT}}$
term is associated with “reorienting” forces transverse to the connecting axis.
To extract the
$F$
coefficients from (4.7), we fix the particle orientation and make specific choices of
$\theta$
(which sets the ambient flow orientation
$\boldsymbol{e}$
) and the auxiliary velocity
$\hat {\boldsymbol{V}}$
. For example, picking
$\theta = 0$
and
$\hat {\boldsymbol{V}} = \boldsymbol{e}_{\parallel }$
and comparing with (6.2) leads to
$F_{\textit{AA}}$
. Choosing
$\theta = \pi /2$
for the same auxiliary flow yields
$F_{\textit{TT}}$
, whereas picking
$\theta = \pi /4$
and
$\hat {\boldsymbol{V}} = \boldsymbol{e}_{\perp }$
determines
$F_{\textit{AT}}$
.
With these three combinations we identify all three
$F$
coefficients for any
$D$
and
$\mathcal{S}$
. We plot them in figure 3 against
$D$
for different
$\mathcal{S}$
. The results of our semi-analytic theory are shown as solid curves, while the results of Fabre et al. (Reference Fabre, Jalal, Leontini and Manasseh2017), which involve numerical solutions of both primary and secondary flows, are indicated as symbols. The present results are in very good agreement with those of Fabre et al. (Reference Fabre, Jalal, Leontini and Manasseh2017) for all separation distances, despite the fact that we truncated the multipole expansion (formally accurate at large
$D$
) at just two terms for both the primary and auxiliary flows. The two-term theory is essentially indistinguishable from the numerical solutions for
$D \gtrsim 4$
(corresponding to one diameter of surface separation), and the agreement remains reasonably accurate down to contact (
$D = 2$
).
In figure 3(a,b), positive values represent attractive forces while negative values represent repulsion. As we can see from the plots, the direction of the force is determined by a combination of distance, oscillatory Stokes number and configuration. In general, the magnitude of the force drops as distance increases. For both very large and very small
$\mathcal{S}$
values, the forces are either attractive or repulsive regardless of the distance. However, for intermediate
$\mathcal{S}$
values (e.g.
${\mathcal{S}}=5$
; yellow curves), there is a sign change in both the axial and transverse configurations as we move along the distance axis, indicating a reversal in the corresponding force. The locations of the force reversal constitute equilibrium points along either axial or transverse directions. This equilibrium point is unstable when particles are aligned along the oscillation axis, while it is stable for particles aligned transverse to the oscillation axis. For moderately large values of
$\mathcal{S}$
(between around 5 and 15) the equilibrium point moves towards smaller
$D$
, where the two-term multipole theory is less accurate.
Figure 3. Comparison of time-averaged forces in all configurations – (a) axial (
$AA$
), (b) transverse (
$AT$
), (c) reorienting (
$AT$
) – showing the numerical calculations of Fabre et al. (Reference Fabre, Jalal, Leontini and Manasseh2017) (symbols) and the present semi-analytic theory (solid) for different distances
$D$
and oscillatory Stokes numbers
$\mathcal{S}$
. Each curve corresponds to a single value of
$\mathcal{S}$
; curves and circles which share the same colour are of the same value of
$\mathcal{S}$
.
The force also shows significant variation with
$\mathcal{S}$
at fixed distance. For example, the axial configuration in figure 3 shows that the force starts attractive for small
$\mathcal{S}$
, but flips sign and becomes repulsive at large
$\mathcal{S}$
. For large
$\mathcal{S}$
, all the curves appear to approach an ‘envelope’ curve, deviating from it only at small
$D$
.
Figure 3(c) shows the ‘reorienting’ component of the force, which is positive for all
$\mathcal{S}$
for small to moderate
$D$
. At larger
$D$
,
$F_{\textit{AT}}$
decays rapidly while also changing sign, although this reversal is rather subtle. From (6.2) we see that positive
$F_{\textit{AT}}$
corresponds to forces tending to reorient the particles into an orientation that is transverse to the oscillation axis.
6.2. Analytic approximations for
$D \gg \delta$
The advantage of the formulation (4.7) is that it allows us to gain analytic insight into the behaviour of the force. We focus on the limit of large separation between the particles (
$D \gg 1$
) and further consider non-overlapping Stokes layers (
$D \gg \delta$
, or
$D^2 {\mathcal{S}} \gg 1$
). Recalling (4.7), the force integral needs three inputs
\begin{align} \boldsymbol{v}_1&=\boldsymbol{V}^\infty +\boldsymbol{v}^{d(1)}_1+\boldsymbol{v}^{d(2)}_1, \end{align}
\begin{align} \boldsymbol{\omega }_1&=\boldsymbol{\omega }^{d(1)}_1+\boldsymbol{\omega }^{d(2)}_1, \textrm {and} \end{align}
\begin{align} \hat {\boldsymbol{v}}&=\hat {\boldsymbol{v}}^{d(1)}+\hat {\boldsymbol{v}}^{d(2)}. \end{align}
The oscillatory vorticity decays exponentially outside Stokes layers of thickness
$\delta$
around each particle. Consequently, the integral in (4.7) splits into two ‘localised’ integrals around each particle. We take advantage of the symmetry to evaluate the integral in the volume surrounding just one of the particles; the integral in the volume surrounding the other particle is identical. To perform the integral centred at particle 1, we approximate the primary disturbance and auxiliary velocity fields created by particle
$2$
,
$\boldsymbol{v}^{d(2)}$
and
$\hat {\boldsymbol{v}}^{d(2)}$
, by a Taylor expansion around the centre of particle
$1$
, retaining terms up to
$O(D^{-4})$
. Since the vorticity decays exponentially, we do not need to involve
$\boldsymbol{\omega }^{d(2)}$
in the integral around particle 1 in the regime where Stokes layers do not overlap.
With these approximations, we perform the integral in spherical coordinates around particle 1 using Mathematica. This yields analytic expressions for the three
$F$
components in (6.2), valid in the limit
$D \gg 1$
,
$D \gg {\mathcal{S}}^{-1/2}$
. These expressions take the general shape
\begin{align} F_i=\frac {{A}_{2i}}{D^2}+\frac {{A}_{3i}}{D^3}+\frac {{A}_{4i}}{D^4} + O\big(D^{-5}\big)\!, \end{align}
where
$i$
denotes
$AA$
,
$TT$
or
$AT$
, and the
${A}_{ni}$
coefficients (9 in total) depend only on
$\mathcal{S}$
. Detailed expressions for these coefficients are in the supplemental
Mathematica
file. We find that the coefficients are interrelated. In particular,
${A}_{3i} = (3/2){A}_{2i}$
for all configurations. The
$AA$
and
$TT$
coefficients have similar shapes but opposite signs, specifically
${A}_{iTT} = -(1/2) {A}_{iAA}$
for
$i=\{1, 2, 3\}$
. We also find that
${A}_{2AT} = {A}_{3AT} = 0$
, and
${A}_{4\textit{AT}} = -({A}_{4AA} - (9/4)A_{2AA})$
. Introducing the shorthand notation
$A_2 = A_{2AA}$
and
$A_4 = A_{4 AA}$
, (6.4) therefore simplifies as
\begin{align} F_{\textit{AA}}&=\frac {{A}_{2}}{D^2}\left (1+\frac {3}{2D}\right )+\frac {{A}_{4}}{D^4} ,\end{align}
\begin{align} F_{\textit{TT}}&=-\frac {1}{2} F_{\textit{AA}} = -\frac {{A}_{2}}{2 D^2}\left (1+\frac {3}{2D}\right )-\frac {{A}_{4}}{2D^4} ,\end{align}
\begin{align} F_{\textit{AT}}&= \frac {{A}_{4\textit{AT}}}{D^4} = -\frac {{A}_{4} - \dfrac {9}{4} {A}_2}{D^4}. \end{align}
Figure 4. Analytic force coefficients
$A_{2}$
and
$A_{4}$
with respect to
$\delta =\sqrt {2/{\mathcal{S}}}$
. Panel (a) shows
$A_{2AA}$
and its asymptote for large
$\delta$
. At
$\delta \approx 0.3501$
,
$A_{2AA}$
changes sign. Panel (b) shows
$A_{4AA}$
and
$-A_{4\textit{AT}}$
(both are nearly identical, note that
$A_{4\textit{AT}}=-A_{4AA}-(9/4)A_{2AA}$
) and their asymptote for small
$\delta$
.
The simplified formulation (6.5) involves only two distinct coefficients, plotted in figure 4. We see that all
$A$
coefficients reverse sign with
$\mathcal{S}$
. We also see that
${A}_{4AA}$
is nearly identical to
$-A_{4\textit{AT}}$
as their difference scales with
$A_2$
, which is much smaller than
${A}_4$
(figure 4).
The leading contribution scales as
$D^{-2}$
. We interpret it as the drag force on a particle held fixed in the far-field streaming flow created by oscillations around the other particle (Li et al. Reference Li, Collis, Brumley, Schneiders and Sader2023). The
$(1 + 3/2D)$
correction factor is due to hydrodynamic interactions between approaching particles in viscous flows; see (Brenner Reference Brenner1961; Rallabandi et al. Reference Rallabandi, Hilgenfeldt and Stone2017). The
$D^{-4}$
terms are more complicated, and arise from a combination of (i) Faxén corrections and hydrodynamic interactions due to particles suspended in the streaming field, and (ii) a generalisation of secondary radiation forces (accounting for both viscous and inertial contributions) due to a product of the oscillatory flow velocity (dominated by the ambient oscillation) and its gradient (
$\propto D^{-4}$
) felt by each particle (Zhang et al. Reference Zhang, Minten and Rallabandi2024).
Figure 5 shows the analytical large-
$D$
results obtained from (6.5) (curves) against the semi-analytic calculations (symbols). Positive values are indicated as solid curves and filled symbols, while negative values are indicated as dashed curves and open symbols. We see that the analytic theory reproduces the calculations for large
$D$
with quantitative precision. Forces along the axis connecting the particles (attraction/repulsion) decay as
$D^{-2}$
, while the force normal to the axis (
$F_{\textit{AT}}$
) decays more quickly as
$D^{-4}$
.
Figure 5. Comparison between the results of the semi-analytic calculations (symbols) and the analytic expression (6.5) (curves). Solid lines and filled circles represent positive values, while dashed lines and open circles represent negative values.
From figure 4, we can also get a sense of how the force changes sign. At large
$D$
, the
${A}_2$
terms dominate, so the reversal of sign of the force is associated with the change in sign of
${A}_{2}$
, which occurs at
$\delta \approx 0.3501$
(
${\mathcal{S}} \approx 16.317$
). Notably, this is the same value of
$\mathcal{S}$
reported by Li et al. (Reference Li, Collis, Brumley, Schneiders and Sader2023) at which the streaming flow due to oscillations around a single sphere reverses direction in the far field. The agreement of the predicted force reversal with the flow reversal of Li et al. (Reference Li, Collis, Brumley, Schneiders and Sader2023) indicates that for sufficiently large separations, the time-averaged force on one particle is a direct consequence of the streaming due to the other particle. We note that
${A}_{4}$
is approximately two orders of magnitude greater than
${A}_2$
, so even though the streaming term of (6.5)
$\propto A_2 D^{-2}$
is longer ranged for large
$D$
, the
$A_4 D^{-4}$
term due to radiation forces can be numerically larger at moderate distances. The crossover into the streaming-dominated regime occurs when
$D \gtrsim |A_4/A_2|^{1/2}$
. For large
$\delta$
(small
$\mathcal{S}$
), we find the asymptotic behaviour
$A_{2} \sim (27/80)\pi (1-2\delta )$
, whereas in the inviscid limit
$\delta \to 0$
(
${\mathcal{S}} \to \infty$
),
${A}_{2}$
approaches zero. By contrast, the
${A}_{4}$
coefficient asymptotes to
$-(1/350)\pi (1023+3177\delta )$
at small
$\delta$
. Notably,
${A}_{4}$
does not vanish as
$\delta \to 0$
and thus becomes the leading term in this limit.
For smaller distances
$D = O(1)$
, the analytic theory is strained as the Taylor and multipole expansions are less accurate, and because the Stokes layers are more likely to overlap. It nonetheless produces useful insights, particularly at large
$\mathcal{S}$
. Comparisons between the analytic theory (curves) and the computations of Fabre et al. (Reference Fabre, Jalal, Leontini and Manasseh2017) (symbols) are shown for
$D$
between 2 and 6 in figure 6. Since the theory is only valid for
$D \gg \delta$
, it becomes restricted to
$\delta \ll 1$
(or
${\mathcal{S}} \gg 1$
) when separation distances become comparable to the particle radius (
$D = O(1)$
). The theoretical
$F_{\textit{AA}}$
and
$F_{\textit{TT}}$
coefficients corroborate this expectation, becoming more accurate at large
$\mathcal{S}$
(figure 6). In the inviscid limit
${\mathcal{S}} \to \infty$
, the
$A_2$
contribution vanishes while the
$A_4$
contribution survives. The dashed curves in figure 6(a, b) represent this limit of the theory, which appears to approximately define an envelope of the data. The data ‘peel off’ from this envelope at different
$D$
(at smaller
$D$
for larger
$\mathcal{S}$
), which we expect is the result of overlapping Stokes layers. It is interesting to note that the analytic theory captures that transverse force contribution
$F_{\textit{AT}}$
for large
$\mathcal{S}$
all the way down to contact, while it is somewhat less accurate for smaller
$\mathcal{S}$
.
Figure 6. Comparison between the numerical calculations of Fabre et al. (Reference Fabre, Jalal, Leontini and Manasseh2017) (circles) and the analytic result (6.5) (solid curves) for moderate separations
$D$
. As expected, the agreement becomes better for
${\mathcal{S}} \gg 1$
(
$\delta \ll 1$
), where the Stokes layers do not overlap. The dashed curve is the limit
$\delta =0$
, where the
$D^{-4}$
term is dominant.
7. Freely suspended particles
We have so far assumed that the particles are fixed, and calculated the time-averaged inertial forces. We now relax this assumption and consider the case where the particles are freely suspended. They now oscillate in response to the ambient flow, and also acquire time-averaged drift velocities from a balance of time-averaged inertial forces with viscous drag.
We consider particles of identical size and density
$\rho _p$
. Due to symmetry and linearity of the primary flow, they oscillate at identical velocity, in phase with each other, when subject to oscillatory ambient flow; see Appendix A for a careful justification. We denote this particle velocity (in dimensionless units) by
$\boldsymbol{V}_{1} \propto e^{i t}$
, and place the frame of reference at the instantaneous centre of one of the particles (say particle 1) when solving for disturbance flow. In the particle-attached reference frame, the ambient flow becomes effectively replaced by the relative velocity
$\boldsymbol{V}^{\infty } - \boldsymbol{V}_{1}$
. To calculate this relative velocity between the particle and the flow, we invoke conservation of momentum of the oscillating particle. For this purpose we approximate each particle as being isolated, neglecting the
$O(D^{-3})$
corrections to the oscillating velocity due to interactions between the particles. These corrections lead to time-averaged velocities of
$O(D^{-5})$
, which are smaller than the accuracy of the two-term multipole expansion used here. Then, we invoke the known result for the relative oscillation velocity (Settnes & Bruus Reference Settnes and Bruus2012; Zhang et al. Reference Zhang, Minten and Rallabandi2024)
\begin{align} \boldsymbol{V}_{1} - \boldsymbol{V}^{\infty } = \mathcal{R}\boldsymbol{V}^{\infty },\quad \textrm {with} \quad \mathcal{R}=-\frac {2\lambda ^2\left (\tilde {\rho }-1\right )}{\lambda ^2\left (2\tilde {\rho }+1\right )+9\lambda +9}, \end{align}
where
$\tilde {\rho } = \rho _p/\rho$
is the density ratio between the particles and the fluid. The oscillating flow (both ambient and disturbance) in the co-moving reference frame is therefore identical to the one obtained for fixed particles, but for the mapping
$\boldsymbol{V}^{\infty } \mapsto -\mathcal{R}\boldsymbol{V}^{\infty }$
.
In addition to this change in the oscillatory flow, the particles are also able to execute a drift over many cycles. If the particles are allowed to move without the application of an external force (
$\langle \boldsymbol{F}\rangle = 0$
), the particles will drift at equal and opposite velocities
$\langle \boldsymbol{V}_2\rangle = \langle \boldsymbol{V}_2\rangle ^{(1)} = -\langle \boldsymbol{V}_2\rangle ^{(2)}$
by symmetry (Appendix A). Then, the left-hand side of (4.7) becomes
$- \hat {\boldsymbol{F}}\boldsymbol{\cdot }\langle \boldsymbol{V}_2\rangle$
, while the right-hand side is unchanged from before. The force in the auxiliary problem
$\hat {\boldsymbol{F}}$
is related to the auxiliary velocity
$\hat {\boldsymbol{V}}$
according to
$\hat {\boldsymbol{F}}=-\mathbb{R}\boldsymbol{\cdot }\hat {\boldsymbol{V}}$
, where
$\mathbb{R}=\mu a (R_{\parallel }\boldsymbol{e}_{\parallel }\boldsymbol{e}_{\parallel }+R_{\perp }\boldsymbol{e}_{\perp }\boldsymbol{e}_{\perp } )$
is a Stokes resistance matrix. The resistance coefficients
$R_{\parallel }$
and
$R_{\perp }$
correspond to steady motion of particles along and perpendicular to their line of centres, and depend on the separation
$D$
(figure 2). Although they can be computed within the two-term multipole expansion by the Faxén relation
$\hat {\boldsymbol{F}}=-6\pi \mu a (\hat {\boldsymbol{V}}-\hat {\boldsymbol{\mathcal{V}}} )$
, such an expansion does not capture contributions from lubrication theory when the particles are close to contact. Instead, we use known results (Brenner Reference Brenner1961; Jeffrey & Onishi Reference Jeffrey and Onishi1984) for these resistances that are valid for all separations (see Appendix B). The right-hand side of (4.7) remains unchanged from before.
The general reciprocal expression given by (4.7), for freely suspended (force-free) particles, yields the secondary velocity of particle 1 as (in dimensional form)
\begin{align} \left \langle \boldsymbol{V}_2\right \rangle = \varepsilon ^2 a \omega \, {\mathcal{S}} |\mathcal{R}|^2 \left (\frac {F_{\textit{AA}}}{R_{\parallel }}\cos ^2{\theta }\boldsymbol{e}_{\parallel}+\frac {F_{\textit{TT}}}{R_{\parallel }}\sin ^2{\theta }\boldsymbol{e}_{\parallel}+\frac {F_{\textit{AT}}}{R_{\perp }}\sin {\theta }\cos {\theta }\boldsymbol{e}_{\perp}\right )\!, \end{align}
where we recall that
$\varepsilon = V^{\infty }/(a \omega )$
is the dimensionless oscillation amplitude. The velocity of particle
$2$
is equal in magnitude and opposite in sign, and
$F_{\textit{AA}}$
,
$F_{\textit{TT}}$
and
$F_{\textit{AT}}$
are force coefficients defined in (6.2). Thus, we obtain the secondary velocity with very little extra effort. The resistance coefficients
$R_{\perp }$
and
$R_{\parallel }$
depend on
$D$
only, approaching
$6 \pi$
for large
$D$
and diverging at contact. The factor of
$\mathcal{S}$
in (7.2) reflects the balance between the inertial force scale in (6.2) and viscous drag.
Figure 7. Secondary velocities corresponding to different configurations: (a) axial velocity due to axial oscillations (
$AA$
), (b) axial velocities due to transverse oscillations (
$TT$
), (c) transverse velocities (
$AT$
) corresponding to a reorientation of the particles with respect to the flow.
In figure 7 we plot the time-averaged particle velocity (with coefficients calculated by the semi-analytic procedure in § 6.1) against distance under the same values of
$\mathcal{S}$
in figure 3. The magnitudes of the velocities are consistent with the corresponding forces, while the peak values are shifted further from the contact point
$D = 2$
. This is due to the fact that the resistance to motion becomes infinite at contact (where inertial forces remain finite), while the force coefficients decay with distance (with the resistance being finite). The location and stability of the equilibrium points of the motion are inherited from those of the forces. It is possible to study the kinematics of the particle pair over long time scales by integrating (7.2) through time. We leave this task to a future study.
Figure 8. Dimensionless time-averaged velocity of particle 1 comparing value extracted from the DNS trajectory data of Kleischmann et al. (Reference Kleischmann, Luzzatto-Fegiz, Meiburg and Vowinckel2024) (squares) and the present semi-analytic theory (circles).
We close this section by comparing our results with those of Kleischmann et al. (Reference Kleischmann, Luzzatto-Fegiz, Meiburg and Vowinckel2024), where the Navier–Stokes equations were solved numerically to simulate the coupled flow and particle motion over time. From the trajectory data in figure 4 of Kleischmann et al. (Reference Kleischmann, Luzzatto-Fegiz, Meiburg and Vowinckel2024) (particle of density ratio
$\tilde {\rho }=4.68$
in the
$\theta = 0$
configuration), we extract the time-averaged particle velocities at 50 cycles from the initial condition (long enough to avoid transient effects). We analyse 7 sets of data ranging from
${\mathcal{S}}=0.63$
to
${\mathcal{S}}=18.81$
, the entire range considered in Kleischmann et al. (Reference Kleischmann, Luzzatto-Fegiz, Meiburg and Vowinckel2024). Due to the set-up of these DNS, each velocity datum occurs at a slightly different
$D$
. We use these values of
$\mathcal{S}$
,
$D$
and
$\tilde {\rho }$
as direct inputs to the semi-analytic theory (7.2), and plot the resulting velocities in figure 8, without any fitting parameters. The theory is in good agreement with the DNS, capturing both the direction and magnitude of the particle motion across
$\mathcal{S}$
values. The direction of motion reverses with
$\mathcal{S}$
, although the location of the reversal appears to be slightly underpredicted by the theory. Differences between the theory and the DNS may be in part due to presence of walls in the simulations. These walls would add additional viscous resistance to particle motion. Furthermore, Stokes layers are also established at the walls, exciting secondary streaming flows in addition to those generated by the particles. Such flows would cause even an isolated off-centred particle to drift, whereas a single particle in an unbounded domain will not drift due to symmetry. Despite these differences, it is clear that the present theory provides a useful means to predict secondary particle motion with relatively low computational effort. Accounting for external boundaries could form the basis of future work.
8. Conclusions
To conclude, we have explored the interactions between two identical particles in an oscillatory flow accounting for the leading effects of inertial nonlinearities. Through a theoretical approach involving small-amplitude expansions, multipole expansions and the application of the Lorentz reciprocal theorem, we derived a method to calculate the time-averaged secondary forces acting between the particles. The analysis demonstrated that, while the full nonlinear flow is complex, it is possible to reduce the problem to an analytically tractable form by focusing on vorticity effects near the particle surfaces. By using the reciprocal theorem, we circumvented the need for direct computation of secondary stress, leading to an efficient framework for understanding particle interactions in oscillatory flows. By providing new analytic insight, the present framework fills the gap between existing work on direct numerical solutions and numerical solutions of scale-separated equations from perturbation theory. Due to its low computational cost, the present work can cover a broader range of frequency and separation distances.
Our findings not only enhance our understanding of particle dynamics in fluid systems but also provide useful theoretical tools for further exploration in practical applications. Although we focused on a single pair of particles for the theoretical development here, the framework has the potential to be extended to multiple particles of identical size and density. Such an extension would be useful to understand particle pattern formation in oscillatory flows. It would also be interesting to extend the ideas here to systems of particles with different sizes, shapes or densities.
Figure 9. Symmetry argument for zero net time-averaged force on the particles. Rotating by
$180^{\circ }$
causes the forces to flip sign, while preserving the configuration of the system.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2025.10510.
Funding
The authors thank the National Science Foundation for support through award CBET-2143943.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Symmetries of forces and motion
A.1. Symmetry of time-averaged forces
Consider the time-averaged force acting on particle
$i$
,
$\left \langle \boldsymbol{F}\right \rangle ^{(i)}$
, in the configuration shown in the left panel of figure 9. We now rotate the geometry by
$180^{\circ }$
(or alternatively rotate the frame of reference). Then, the force vectors also rotate by
$180^{\circ }$
, transforming to
$-\left \langle \boldsymbol{F}\right \rangle ^{(i)}$
, while the particles simultaneously swap positions, as shown in the right panel of figure 9. We note that the oscillatory velocity also changes sign, but this is inconsequential to the time-averaged forces, which depend quadratically on
$\boldsymbol{\boldsymbol{V}}^\infty$
. If the particles are identical, both panels describe same configuration, so we conclude that
$\left \langle \boldsymbol{F}\right \rangle ^{(1)} = -\left \langle \boldsymbol{F}\right \rangle ^{(2)}$
. If the particles are free to drift, the same argument applies to their time-averaged (drift) velocities.
Figure 10. Symmetry argument to show that particles must oscillate with identical velocity. Flipping the sign of the ambient velocity and also rotating by
$180^{\circ }$
swaps particle velocities but recovers the original configuration.
A.2. Symmetry of oscillatory motion
We show here that both particles must oscillate at identical velocities regardless of their separation distance and orientation relative to the ambient flow. We start with the configuration in the top left panel of figure 10, in which the particles oscillate with generally distinct velocities
$\boldsymbol{V}^{(j)}_1(t)$
, in response to the ambient flow
$\boldsymbol{V}^\infty$
. The primary flow is governed by the linear system (3.2), with boundary conditions on the particle surfaces now accounting for
$\boldsymbol{V}_1^{(j)}$
. Since there are no external (i.e. non-hydrodynamic) oscillatory forces acting on the particles,
$\boldsymbol{V}^{(j)}_1$
are linear in
$\boldsymbol{V}^\infty$
. Furthermore, although (3.2) involves time derivatives, the oscillatory fields depend parametrically on time as
$e^{it}$
, as there is no reference to initial conditions (which are presumed to have occurred in the distant past).
We now apply two different transformations to the situation depicted in the top left panel of figure 10. On the one hand, a
$180^{\circ }$
rotation leads to the bottom left panel, where particle swap places and the velocities change sign. On the other hand, a change in sign of the ambient flow (
$\boldsymbol{V}^\infty \mapsto -\boldsymbol{V}^\infty$
) leads to the top right panel, where the particle velocities change sign (
$\boldsymbol{V}_1^{(j)} \mapsto - \boldsymbol{V}_1^{(j)}$
) due to linearity, but the particles remain in their original positions. Applying both a rotation and a sign flip leads to the bottom right panel.
If the particles are identical (in size and density), the scenarios on either diagonal of figure 10 are physically identical. We thus conclude that
$\boldsymbol{V}^{(1)}_1 = \boldsymbol{V}^{(2)}_1 = \boldsymbol{V}_1$
.
Appendix B. Resistance coefficients
Figure 11. Stokes resistances showing the approximate large-
$D$
result from the two-term multipole expansion with Faxén’s law (solid curves) and uniformly valid results valid for all inter-particle separations (dashed curves). (a) Parallel configuration. The dashed curve denotes the exact result of Brenner (Reference Brenner1961). (b) Perpendicular configuration. The dashed curve denotes a two-term approximation, valid for all
$D$
, due to Jeffrey & Onishi (Reference Jeffrey and Onishi1984). In both configurations we use resistances corresponding to the dashed curves to calculate the velocities in figure 7.
We plot the resistance coefficients for both axis-parallel and axis-perpendicular motion in figure 11. Calculations from Faxén’s law with the two-term multipole expansion are indicated as solid curves. They are accurate at large distances
$(D \gtrsim 3)$
but do not account for lubrication effects when particles are near contact. We remedy this by using a combination of known exact and asymptotic results. For
$R_{\parallel }$
we use the exact result of Brenner (Reference Brenner1961), which diverges as the inverse of surface separation distance,
$R_{\parallel } \sim ({6 \pi }/{4})({1}/{(D/2-1)})$
, as the particles approach contact (figure 11
a). The perpendicular resistance behaves as
$R_{\perp } \sim (6\pi/3)\log [(D/2-1)^{-1} ]$
close to contact, and approaches
$6 \pi$
for large
$D$
. We use the approximation in § 5.3 of Jeffrey & Onishi (Reference Jeffrey and Onishi1984), truncated at two terms, to obtain
$R_{\parallel } = 6 \pi [1 + (1/3)\log (1 + (D/2 - 1)^{-1}) ]$
. Including more terms in the series leads to deviations of at most
$2.5\,\%$
. This expression is accurate for all
$D$
and is plotted as a dashed curve in figure 11(b).
Open access
