Pairwise interaction of spherical particles aligned in high-frequency oscillatory flow

Abstract

We present a systematic simulation campaign to investigate the pairwise interaction of two mobile, monodisperse particles submerged in a viscous fluid and subjected to monochromatic oscillating flows. To this end, we employ the immersed boundary method to geometrically resolve the flow around the two particles in a non-inertial reference frame. We neglect gravity to focus on fluid–particle interactions associated with particle inertia and consider particles of three different density ratios aligned along the axis of oscillation. We systematically vary the initial particle distance and the frequency based on which the particles show either attractive or repulsive behaviour by approaching or moving away from each other, respectively. This behaviour is consistently confirmed for the three density ratios investigated, although particle inertia dictates the overall magnitude of the particle dynamics. Based on this, threshold conditions for the transition from attraction to repulsion are introduced that obey the same power law for all density ratios investigated. We furthermore analyse the flow patterns by suitable averaging and decomposition of the flow fields and find competing effects of the vorticity induced by the fluid–particle interactions. Based on these flow patterns, we derive a circulation-based criterion that provides a quantitative measure to categorize the different cases. It is shown that such a criterion provides a consistent measure to distinguish the attractive and repulsive arrangements.

1. Introduction

The behaviour of particles in oscillating flows is a widespread phenomenon in many engineering systems and has attracted much attention in recent decades. In industrial applications, oscillating flows are used to improve agglomeration and removal of particles in water treatment facilities (Halfi, Brenner & Katoshevski 2019; Halfi et al. 2020) and in diesel engines (Katoshevski, Dodin & Ziskind 2005; Ruzal-Mendelevich, Katoshevski & Sher 2016; Gupta et al. 2019). Recently, research has also been conducted to investigate how oscillations can be used as a non-invasive biomedical application for drug delivery inside the ear (Sumner, Mestel & Reichenbach 2021; Harte et al. 2023).

Of particular interest for the present study are applications and investigations of fluid–particle interactions in oscillating conditions where gravity is negligible and at the same time the inertia of the particles is of decisive importance. Such conditions exist in space, for instance, and have already been used for a variety of scientific studies. For example, experiments were performed as part of the First International Microgravity Laboratory in $1992$ on board a space shuttle to investigate the growth of crystals in the absence of gravity to gain insights into their growth process (Trolinger et al. 1996a; Trolinger, Rangel & Lal 1996b). Studies of the same kind were recently conducted on board the International Space Station to examine the flocculation behaviour of a variety of cohesive particles (Kleischmann et al. 2021). Such experimental set-ups are subjected to small oscillations, known as g-jitter, caused by the spacecraft's propulsion system and on-board machinery. For obvious reasons, such experimental campaigns involve a great deal of effort as well as expense and are, therefore, rare. Nevertheless, such special conditions can also be found in Earth-bound environments. In the field of microfluidics, the inertial properties of particles can be used to manipulate their dynamics in oscillatory fluid flows (Mutlu, Edd & Toner 2018; Dietsche et al. 2019; Mutlu et al. 2020). As the particles are usually only a few micrometres in size, the settling due to gravity can be largely neglected (Owen et al. 2023).

Regardless of the area of application or the specifications of the set-up, the behaviour of particles in oscillating flows is determined by the fluid–particle interaction of each individual particle. For this purpose, Coimbra & Rangel (2001) proposed an analytical solution for particle motion in an oscillating background flow based on the work of Tchen (1947). The theoretical prediction has subsequently been confirmed by the experiments of Coimbra et al. (2004) and L'Espérance et al. (2005, 2006). In these experimental studies, the authors investigated the response of an individual spherical particle oscillating in a viscous liquid-filled container, neutralizing gravity by tethering the particle to prevent vertical movement based on its density. Similar experiments were performed by Hassan et al. (2006), Hassan & Kawaji (2007), Hassan & Kawaji (2008) and Saadatmand & Kawaji (2010). Simic-Stefani, Hu & Kawaji (2005) extended the idea of studying the particle motion in oscillatory flow without the impact of gravity to experiments in microgravity on parabolic flights. However, for parabolic flights or comparable methods, such as the drop tower, the state of weightlessness only lasts for a very short period of time and long-term studies are not possible. Numerical simulations, on the other hand, easily achieve microgravity by setting gravitational acceleration to zero (e.g. Simic-Stefani et al. 2005; Saadatmand & Kawaji 2010; Satish et al. 2022). Consistent findings across theory, experiments and numerical simulations demonstrate that particle excursion increases with increasing deviation of the density ratio from unity. Hence, particle motion reflects the oscillation frequency with a delay related to particle inertia.

Another important aspect is the flow field generated by the fluid–particle interaction, in which either the fluid or the particle oscillates, characterized by the instantaneous flow field or the flow field averaged over one oscillation period, which is commonly referred to as steady streaming. Due to the fact that the averaged flow field is usually non-zero, it causes a change in particle motion that deviates from the prevailing oscillatory motion (Riley 1966). Analytical models for steady streaming results of flows around a sphere have been developed by Lane (1955) and Riley (1966, 1967), in which inner and outer solutions were obtained by perturbation theory and associated streamlines were described. To this end, vortices near the particle surface (within the oscillating boundary layer) were defined as primary vortices whereas secondary vortices emerge further away from the particle that in turn are driven by the rotation of the primary vortices. In this context, Riley (1966) has elaborated on two distinct cases for small-oscillation amplitudes. For larger inertial forces, the field of the steady streaming shows a pattern consisting of both the primary and secondary vortices. For smaller inertial forces, only one dominant structure of primary vortices exists, which expands over the whole domain (Rott 1964). In the case of a stationary sphere in an oscillating fluid, Li et al. (2023) recently provided a phase diagram that can be utilized to determine the presence of only primary or primary and secondary vortices by means of the non-dimensional frequency and amplitude of the oscillation.

The flow structures around a single spherical particle were extensively investigated by theoretical and numerical studies (Chang & Maxey 1994; Alassar & Badr 1997; Blackburn 2002; Alassar 2008; Fabre et al. 2017; Jalal 2018). Together with experimental measurements (Kotas, Yoda & Rogers 2007; Otto, Riegler & Voth 2008), the size and the rotational directions of the structures of the steady streaming were analysed. According to Voth et al. (2002) and Klotsa et al. (2009), detailed knowledge of the flow structures of an individual particle should provide information about the interaction behaviour of two particles. Hence, experimental and numerical studies explored the interaction of two particles under either vertical (Voth et al. 2002) or horizontal oscillations (Klotsa et al. 2007, 2009; Van Overveld et al. 2022; Van Overveld, Clercx & Duran-Matute 2023). The latter studies, in which the particles were positioned on the bottom plate of a fluid container perpendicular to the oscillation, focused on the equilibrium distance between the particles. A major finding was that the equilibrium is determined by a balance of short-range repulsive and long-range attractive forces induced by the flow structures.

Although these studies provide valuable insights into the particle behaviour subjected to oscillating flows, they are limited to initial particle configurations that are perpendicular with respect to the oscillation direction. Moreover, the particle behaviour is affected by the presence of the bottom plate and the effect of gravity. Despite the fact that Klotsa et al. (2007, 2009) and Van Overveld et al. (2022, 2023) ensured that the friction between the particles and the bottom plate is minimized, the flow field that is formed is nevertheless restricted by the bottom plate. This restriction is reflected in particular by the fact that the flow structures are limited below the particles, which introduces a distortion of the flow field that might alter the system. To avoid this effect, Fabre et al. (2017) and Jalal (2018) numerically investigated two stationary particles in an oscillating domain. The studies provide key results for interpreting the interaction forces that might lead to either attraction or repulsion of the particles. However, for freely moving particles the distance changes over time and, consequently, this may also affect the particle interaction.

Despite the large number of studies on the interaction of two particles in an oscillating flow, no investigations have been conducted in which the particles can move freely, independently of gravity and without spatial constraints. The present study addresses this research gap. We perform particle-resolved direct numerical simulations of two mobile particles in an oscillating box filled with a viscous fluid. The set-up enables us to analyse the behaviour of the freely moving particles as well as the resulting flow structures, with respect to the change of the distance between the particles as a function of time. In order to focus only on the hydrodynamic effects of the fluid–particle interactions, we neglect gravity and ensure that the particles are not affected by short-range wall effects from the outer boundaries of the domain. Based on this, we aim to identify the impact of the particle arrangement, in terms of mutual distance and orientation angle, the oscillation frequency and the particle inertia, characterized by different particle densities and tie these observations to the respective flow characteristics.

The article is structured as follows. First, we introduce the numerical model along with the description of the computational approach and provide a comprehensive validation in § 2. We then analyse the impact of the non-dimensional frequency $S$, the initial inter-particle distance $\zeta _i$ and orientation of the particle alignment $\theta _i$, as well as the density ratio $\rho _s$ on the mutual particle behaviour. This is followed by a detailed analysis of the resulting flow fields in § 4 and, finally, a summary of the conclusions and a brief outlook in § 5.

2. Computational model

2.1. Numerical set-up

The scenario under consideration is a cubic container with dimensions of $L_{x,y,z} = L = 10 D_p$, where $D_p$ is the particle diameter. The container is filled with a viscous, incompressible fluid characterized by its density $\rho _f$ and kinematic viscosity $\nu _f$. As illustrated in figure 1, two monodisperse particles are submerged in this container and placed with an initial particle distance $\zeta _i$ as well as an initial orientation $\theta _i$ with respect to the horizontal direction of oscillation (black arrows). Initially, the fluid is at rest and the centre of $\zeta _i$ coincides with the centre of the domain. The origin of the Cartesian coordinate system is located in the back lower left corner and oriented according to the axes in figure 1. The particles are free to move in response to the container that oscillates in the $x$ direction with $u_f =- A_f \varOmega \sin {\varOmega t}$. Here, $A_f$ is the distance amplitude and $A_f \varOmega$ the velocity amplitude, $t$ denotes time, $\varOmega = 2 {\rm \pi}f$ is the angular frequency and $f$ is the frequency of the oscillatory motion of the domain. The walls of the domain are characterized by stress-free boundary conditions, i.e. ${\rm d}u_t / {\rm d}n=0$ and $u_n=0$. Here, $u_t$ and $u_n$ are the tangential and the normal fluid velocity components relative to the wall, respectively, and $n$ is the normal direction on that same wall. A no-slip condition was applied on the particle surface.

Figure 1. Two-dimensional sketch of the cubic numerical domain with size $L_{x,y,z} = L = 10D_p$. Two mobile spherical particles with identical diameters $D_p$ are initially placed so that the centre of the initial particle distance $\zeta _i$ coincides with the centre of the domain. The arrangement is aligned with an initial angle $\theta _i$ with respect to the direction of oscillation. The domain oscillates with velocity $u_f$ in the $x$ direction, as indicated by the arrows.

2.2. Governing equations

We use direct numerical simulations to compute the dynamics inside the oscillating container and apply the immersed boundary method (IBM) to account for the fluid–particle interaction (Uhlmann 2005; Kempe & Fröhlich 2012; Biegert, Vowinckel & Meiburg 2017). The simulations are performed in a non-inertial reference frame, accounting for the acceleration of the applied external oscillation. Since there are no free surfaces or density variations in the fluid, the effect of the frame acceleration is accounted for in the particle motion only. Therefore, we can solve the Navier–Stokes equations and the continuity equation for an incompressible Newtonian fluid given by

(2.1)\begin{gather} \frac{\partial{\boldsymbol{u}}}{\partial{t}}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}\boldsymbol{u}) ={-}\frac{1}{\rho_f}\boldsymbol{\nabla} p + \nu_f \nabla^2 \boldsymbol{u} + \boldsymbol{f}_{\!IBM} , \end{gather}

(2.2)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0 . \end{gather}

Here, $\boldsymbol {u}=(u,v,w)^{{\rm T}}$ represents the fluid velocity vector in Cartesian components, $p$ is the fluid pressure and $\boldsymbol {f}_{\!IBM}$ the volume force based on the IBM. The $\boldsymbol {f}_{\!IBM}$ connects the motion of the particle to the fluid phase and acts in the vicinity of the fluid–particle interface to enforce a no-slip condition on the particle surface. Equations (2.1) and (2.2) are integrated by a third-order low-storage Runge–Kutta scheme in time together with a second-order finite-difference method in space, respectively. The fast Fourier transform is applied to calculate the pressure correction for continuity by means of a direct solver.

In the context of the IBM, the motion of each individual spherical particle is calculated by solving the Newton–Euler equations, where Newton's equation for the translational particle velocity needs to be transformed into the non-inertial frame. In the non-inertial reference frame, which moves with $\boldsymbol {u}_f = (u_f, 0, 0)^{\rm T}$ relative to the inertial frame, we can calculate the non-inertial particle velocity as $\boldsymbol {u}_p' = \boldsymbol {u}_p - \boldsymbol {u}_f$, where ${\boldsymbol {u}_p=(u_p,v_p,w_p)^{{\rm T}}}$ represents the translational velocity vector of the particle in an inertial frame. By rearranging, we can state that

(2.3)\begin{equation} \frac{\text{d}\boldsymbol{u}_p}{\text{d}t} = \frac{\text{d}\boldsymbol{u}_p'}{\text{d}t} + \frac{\text{d}\boldsymbol{u}_f}{\text{d}t}. \end{equation}

If we take Newton's equation in an inertial frame

(2.4)\begin{equation} m_p \frac{\text{d}\boldsymbol{u}_p}{\text{d}t}= \boldsymbol{F}_{h} + \boldsymbol{F}_{c} + \rho_f V_p \frac{\text{d}\boldsymbol{u}_f}{\text{d}t}, \end{equation}

and apply (2.3) to transform (2.4) to the non-inertial frame, we obtain the velocities in the non-inertial frame based on the equations of motion for translation:

(2.5)\begin{equation} m_p \frac{\text{d}\boldsymbol{u}_p'}{\text{d}t}= \boldsymbol{F}_{h} + \boldsymbol{F}_{c} - \left(\rho_p - \rho_f \right) V_p \frac{\text{d}\boldsymbol{u}_f}{\text{d}t} \end{equation}

and rotation:

(2.6)\begin{equation} I_p \frac{\text{d}\boldsymbol{\omega}_p}{\text{d}t}= \boldsymbol{M}_h + \boldsymbol{M}_c. \end{equation}

Here, $m_p$ and $V_p$ are the particle mass and volume, respectively, $\rho _p$ the particle density and $I_p = {\rm \pi}\rho _p D_p ^ 5 / 60$ the moment of inertia with $\boldsymbol {\omega }_p = ( \omega _{p,x}, \omega _{p,y}, \omega _{p,z} )^{{\rm T}}$ the angular velocity vector. The last term of the right-hand side of (2.4) and (2.5) accounts for the pressure gradient of the background acceleration in the inertial and non-inertial frame, respectively. As an important consequence arising from the change of the observational frame, the background acceleration of the non-inertial frame is characterized by the motion of the particle mass minus the displaced fluid mass and, hence, represents the effect of particle inertia. Note that (2.6) remains unchanged for the translation from the inertial to the non-inertial reference frame, because we do not apply any rotational motion. As noted earlier, we neglect gravity in (2.5) to limit the governing mechanism to the fluid–particle and particle–particle interactions. The first terms on the right-hand side of (2.5) and (2.6), $\boldsymbol {F}_{h}$ and $\boldsymbol {M}_{h}$, denote the hydrodynamic forces and torques, respectively. The fluid–particle interactions can be computed via

(2.7a,b) \begin{equation} \boldsymbol{F}_{h} = \oint_{S_p} \boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{n} \,{\rm d}S , \quad \boldsymbol{M}_{h} = \oint_{S_p}\boldsymbol{r} \times (\boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{n}) \,{\rm d}S , \end{equation}

where $\boldsymbol {\tau } = -p \boldsymbol {I} + \mu _f [ \boldsymbol {\nabla } \boldsymbol {u} + (\boldsymbol {\nabla } \boldsymbol {u} )^{\rm T} ]$ represents the hydrodynamic stress tensor that includes the identity tensor $\boldsymbol {I}$ and the dynamic viscosity $\mu _f$ of the fluid. Furthermore, $\boldsymbol {n}$ is the normal vector pointing outward from the particle surface $S_p$, and $\boldsymbol {r}$ is the position vector pointing from the particle centre of mass to a point on $S_p$. The forces and torques due to collision are represented by $\boldsymbol {F}_{c}$ and $\boldsymbol {M}_{c}$, respectively. Here, we account for normal and tangential contact as well as unresolved lubrication forces that emerge when fluid is squeezed out of particle gaps $\zeta _t \leq 2h$, with $\zeta _{t}$ the distance between the particles at any time and $h$ the grid cell size. Full details on the computation of short-range hydrodynamic as well as collision forces and torques can be found in Biegert et al. (2017) and Vowinckel et al. (2019), but we note that those effects are not relevant for this study as no collisions occur in the cases examined.

For our numerical study, we introduce characteristic scales to non-dimensionalize the governing equations where $\ell$ and $\boldsymbol {u}$ are the relevant length and velocity scales, respectively, that appear in (2.1)–(2.6):

(2.8)\begin{equation} \left. \begin{gathered} \ell = D_{p} \tilde{\ell}, \quad \boldsymbol{u} = D_p \varOmega \tilde{\boldsymbol{u}}, \quad t = \frac{\tilde{t}}{\varOmega},\\ \rho = \rho_{f} \tilde{\rho}, \quad p = \rho_f D_p^2 \varOmega^2 \tilde{p}, \quad f_{IBM} = D_p \varOmega^2 \tilde{f}_{IBM},\\ m = m_{f} \tilde{m} = \rho_f V_{p} \tilde{m}, \quad V = D_p^3 \tilde{V}, \quad F = m_f D_p \varOmega^2 \tilde{F}, \\ I = m_f D_p^2 \tilde{I}, \quad \omega = \varOmega \tilde{\omega}, \quad M = m_f D_p^2 \varOmega^2 \tilde{M}. \end{gathered} \right\} \end{equation}

Here, $m_f$ is the mass of fluid of an equivalent particle volume and $F$ and $M$ represent forces or torques on the particles, respectively. The tilde symbol indicates the dimensionless variables. Applying (2.8) to (2.1)–(2.6) yields the dimensionless Navier–Stokes, continuity and Newton–Euler equations with the Reynolds number $Re = u_{f,max} D_p / \nu _f$, where $u_{f,max}=A_f \varOmega$ is the velocity amplitude and $\rho _s = \rho _p / \rho _f$ the density ratio:

(2.9)\begin{gather} \frac{\partial{\tilde{\boldsymbol{u}}}}{\partial{\tilde{t}}} + \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} (\tilde{\boldsymbol{u}}\tilde{\boldsymbol{u}}) ={-}\tilde{\boldsymbol{\nabla}} \tilde{p} + \frac{1}{Re} \tilde{\nabla}^2 \tilde{\boldsymbol{u}} + \tilde{\boldsymbol{f}}_{\!\!IBM}, \end{gather}

(2.10)\begin{gather} \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}=0, \end{gather}

(2.11)\begin{gather} \tilde{m}_p \frac{\text{d}\tilde{\boldsymbol{u}}_p}{\text{d}\tilde{t}}= \tilde{\boldsymbol{F}}_{h} + \tilde{\boldsymbol{F}}_{c} - \left(\rho_s - 1 \right) \tilde{V}_p \frac{\text{d}\tilde{\boldsymbol{u}}_f}{\text{d}\tilde{t}} , \end{gather}

(2.12)\begin{gather} \tilde{I}_p \frac{\text{d}\tilde{\boldsymbol{\omega}}_p}{\text{d}\tilde{t}}= \tilde{\boldsymbol{M}}_h + \tilde{\boldsymbol{M}}_c . \end{gather}

In addition, we introduce a streaming Reynolds number $Re_s$ which is defined by the product of the amplitude ratio $\epsilon = A_f / D_p$ and $Re$ to obtain $Re_s = \epsilon Re = A_f^2\varOmega / \nu$. Following Coimbra & Rangel (2001), Coimbra et al. (2004) and L'Espérance et al. (2005), we also define a Strouhal number $Sl=\varOmega D_p / (9 u_{f,max})$ for spherical objects to obtain the non-dimensional frequency $S=Sl Re /4=D_p^2 \varOmega / (36 \nu _f)$. Note that $1/4$ is a geometric factor that allows for a straightforward comparison with the work of those authors who used the particle radius as a characteristic length scale, while we chose $D_p$ for the present study. As an important consequence, the amplitude $A_f$ that controls $u_{f,max}$ cancels out for the definition of the non-dimensional frequency. Note that $S$ represents the equivalent to the frequency parameter $M^2 = D_p^2 \varOmega / (4 \nu _f) = 9 S$ that is often used in the literature (e.g. Riley 1966, 1967). As stated in § 1, Riley (1966) distinguished two cases for a single particle to characterize the flow pattern of the steady streaming under certain flow parameters: (i) for $ Re \ll 1$ and $S \ll 1/9$ the area over which the primary vortices propagate is much larger than $D_p/2$ and (ii) for $S \gg 1/9$ the primary vortices are confined to the particle surface and adjacent secondary vortices propagate into the far field. Although the focus of this study is on the oscillating behaviour of two particles, this classification is still useful for describing the obtained flow fields.

Moreover, we define a Stokes number $St=\tau _p / \tau _f$ using the characteristic particle time scale $\tau _p = | \rho _s - 1 | D_p^2 / (18 \nu _f)$ and the characteristic flow time scale $\tau _f = 1 / \varOmega$, which yields $St = | \rho _s - 1 | 2 S$. Since $St$ is determined by the linear combination of the density ratio $\rho _s$ and the dimensionless frequency $S$, we focus in the following on the influence of these two non-dimensional parameters on the behaviour of the particles by varying the initial gap size $\zeta _i$.

In what follows, all equations, results and quantities are presented in dimensionless form unless stated otherwise. Therefore, we drop the tilde symbol for the remainder of the article. In addition, to guarantee sufficient spatial resolution for the parameter space presented here, we discretize the domain by a uniform rectangular grid of a grid cell size $h = L / 200$ resulting in $D_p / h = 20$, if not stated otherwise. The time step is chosen as $T / 1000$, where $T$ is the oscillation period, to also guarantee a high temporal resolution of the fluid motion. We verified by means of test simulations that the particle motions and their mutual interaction are not affected by numerical parameters such as the time step, spatial discretization and domain size. Note that preliminary tests have shown that doubling the domain size only has an impact of $\mathcal {O}(10^{-3})$ on the particle motion and is therefore negligible.

2.3. Comparison with benchmark data

The numerical code used in the present study has already been validated and used for detailed studies of fluid–particle and particle–particle interactions, such as the rheological behaviour of sediment beds under shear and particles settling under gravity in quiescent fluids (Biegert et al. 2017; Vowinckel et al. 2019, 2021). As a most relevant case for our scenario, Biegert et al. (2017) have obtained excellent agreement for a sphere settling in a large container (Mordant & Pinton 2000) and for a settling sphere approaching a wall (Ten Cate et al. 2002). Since the set-up under investigation considers two particles in an oscillating domain filled with fluid, we extend our validation to set-ups with oscillating motion patterns. In this regard, we compare our computational results with analytical, experimental and numerical data of single- and two-particle set-ups in which either the fluid or the particles oscillate.

In a first step, we validate the excursions of the particle trajectories as well as the qualitative and quantitative streamline patterns for a set-up with a single particle with benchmark data. The results and thorough descriptions of the validation conditions are provided in Appendix A. In a next step, we examine the characteristics of the flow field generated by two stationary particles of equal diameter arranged in alignment with a unidirectional, monochromatic oscillatory fluid flow. The particles are horizontally aligned, i.e. $\theta _i = 0^{\circ }$, and placed with a surface separation distance of $\zeta _i = \zeta _t = 1.00$ so that the centre of the particle pair coincides with the centre of the numerical domain. The dimensions of the numerical domain are the same as described in § 2.1 with a cubic box of size $L_{x,y,z} = L = 10$. The properties of the applied oscillation are $S = 0.89$, $ Re = 0.32$ and $ Re _s = 0.0032$. The flow characteristics of the steady streaming, which has been computed after $100 T$, are presented by the streamlines and vorticity contours in figure 2(a). The vorticity contours are based on the calculation of the polar component of the vorticity $\omega _z(x,y )$ and the colouring is according to the presented colour bar, which ranges from $-0.01$ (blue) for clockwise rotations to $0.01$ (red) for counterclockwise rotations. The sketch in figure 2(b) represents the fluid–particle interaction in a simplified form, highlighting that each particle is surrounded by four vortex structures, as indicated by the streamlines as well as the vorticity contours. These flow patterns are usually referred to as quadrupoles (Longuet-Higgins 1998; Tho, Manasseh & Ooi 2007). As emphasized by Klotsa (2009), the appearance of these quadrupole structures is due to the fact that the presented illustrations are planes through three-dimensional toroidal vortex structures surrounding each particle. These rotate in such a way that the flow is directed towards the respective particle along the axis of oscillation and away from it perpendicular to the oscillation, as indicated in figure 2. Similar flow properties of such planes have already been reported by Fabre et al. (2017) and Jalal (2018) for two stationary spheres and by Tho et al. (2007) for two bubbles attached to a wall in oscillating fluid flow.

Figure 2. Flow characteristics illustrated by streamlines and vorticity contours for steady streaming of two stationary and axially arranged particles in oscillatory flow. (a) The numerical results of the set-up with ${S = 0.89}$, $ Re = 0.32$ and $ Re _s = 0.0032$. (b) A simplified sketch to highlight the quadrupole structures surrounding each particle as well as the direction of the imposed oscillation by the black arrows. The vorticity colouring ranges from $-0.01$ (blue) for clockwise rotations to $0.01$ (red) for counterclockwise rotations.

With reference to the case distinctions for a single-particle set-up according to Riley (1966) (cf. § 2.2), we note that the present case would represent an intermediate regime, as it could not be assigned to either the first (small $Re$ and small $S$) or the second (large $S$) class. The pattern of the steady streaming, on the other hand, would comply with the conditions for class one, as it apparently consists only of the primary vortices. While the representations of the vorticity contours are only apparent in the vicinity of the respective particles, indicating that the highest vorticity is located close to particle surfaces, the vortex structures presented by the streamlines extend over the entire domain revealing the characteristic quadrupole pattern. The extension as well as the shape of the vortex structures is limited by the walls of the numerical domain. However, since the fluid velocities farther away from the particles rapidly decay to low values, their effects on the particle dynamics become negligible if the particles are relatively close to each other. Nevertheless, depending on the properties of the oscillation and the particles, the characteristics of the vorticity and vortex structure change. A more detailed analysis of the impact of these properties on the flow characteristics as well as an explanation of the division into central and peripheral sections is given in § 4.

In summary, we observe that our simulation approach is capable of reproducing the relevant flow features for two fixed spheres in oscillating fluid motion. Considering the additional validation of the set-ups with a single sphere (Appendix A), we conclude that the computational approach is suitable for detailed analyses of the interaction of two particles in oscillatory fluid flow.

3. Particle dynamics

3.1. Initial motion

When the fluid container starts moving, the first oscillation is to a given direction, in this case to the right as determined by the way the oscillation is performed. For an example of two particles with a density greater than that of the liquid, i.e. $\rho _s > 1$, the inertia of the particles is larger than that of the surrounding fluid. Since we solve for particle motion in a non-inertial frame, both particles shift to the left for this case. At the start of the simulation, the fluid is at rest and the particles do not experience a counterflow at the beginning of their movement. This condition results in a relatively large excursion. After the reverse of the oscillation of the domain in the opposite direction (now to the left), the particles denser than the fluid tend to move again to the opposite direction (to the right). However, after this initial oscillation period they encounter the previously created flow field induced by the fluid–particle interaction, which still points to the preceding direction of movement. As a consequence, the particle motion and, hence, the excursion are slowed down.

As this initial displacement has now established an offset, in which the distance of both particles to the left wall of the domain is slightly smaller than to the right, the system thenceforth recovers by having the two particles drifting back towards the centre of the domain. The particle pair approaches the centre of the domain asymptotically so that for a horizontally aligned arrangement ($\theta _i = 0^{\circ }$) with $\rho _s = 4.68$ and $\zeta _i = 0.75$, the artefact from the initial conditions is compensated after approximately $20$ oscillation periods. This can be seen for the examples of $S = 0.35$ and $1.05$ in figure 3, where figure 3(a) shows the instantaneous location of the centre between the particles in relation to its initial position over time and figure 3(b) illustrates the evolution of the inter-particle centre over time as moving average with an averaging window of one oscillation period.

Figure 3. Location of the midpoint between the particles relative to its initial location over time for (a) the instantaneous location $\Delta x_c$ and (b) the moving average $\langle \Delta x_c\rangle$ with an averaging window equal to the oscillation period $T$.

3.2. Inter-particle distance over time

Since we could see in § 3.1 that the influence of the initial condition of fluid at rest becomes negligible after a few oscillation periods, we can now analyse the near-field and far-field hydrodynamic interactions of two particles aligned with $\theta _i = 0^{\circ }$ for various initial particle distances $\zeta _i$ and oscillation frequencies $S$ over longer simulation times. Moreover, we investigate the particle dynamics as a function of particle inertia expressed by the density ratio $\rho _s$.

Since we neglect gravity and consider two particles with an initial orientation angle $\theta _i = 0^{\circ }$, the particle coordinates remain constant in $y$ and $z$ directions at $L_y/2$ and $L_z/2$, respectively. The initial positions in the $x$ direction are modified to result in $\zeta _i = \{0.10, 0.25, 0.50, 0.75, 1.00, 1.50, 2.00, 3.00\}$. For all these runs, we use a spatial discretization of $D_p / h = 20$, as given in § 2.1, except for $\zeta _i = 0.10$ where we employ an even higher resolution with $D_p / h = 40$ to prevent the lubrication model of Biegert et al. (2017) becoming active for $\zeta _{t}<2h$. These choices also allowed for a proper resolution of the particle gap even for the arrangement with a very small gap width. For the following analyses, oscillations in a range of $S = [0.07, 2.09]$ with a constant amplitude of $\epsilon = 0.1$ were applied for a total number of $100$ oscillation periods.

Figure 4 illustrates the instantaneous evolution of $\zeta$ over time for various choices of $S$ for a representative set-up with initial particle distance $\zeta _i = 0.25$ and $\rho _s = 4.68$. The horizontal dashed line represents the respective state of equilibrium, in which $\zeta$ does not change with time. As can be seen from the results, different $S$ lead to different particle dynamics in terms of the change of the particle distance $\zeta$ despite identical initial conditions. In general, results with increasing values of $\zeta > 0$ over time show that the particles drift apart in a repulsive behaviour. Vice versa, values of $\zeta < 0$ indicate that the particles are approaching each other, representing an attractive motion.

Figure 4. Evolution of instantaneous $\zeta$ over $100$ oscillation periods for $\rho _s = 4.68$, $\zeta _i = 0.25$ and a variety of $S$.

Since we have already seen that $S$ has a direct correlation to the resulting behaviour of the particles, we investigate $\zeta$ more systematically. Here, $\zeta$ represents the change in the particle distance over time, computed as $\zeta = \zeta _{t} - \zeta _i$. Therefore, we define $\zeta _{100}$ as the change of the particle distance $\zeta$ relative to its initial value after $100 T$, i.e. 100 oscillation periods, to provide a quantitative measure for the particle dynamics. In this context, we also compare the results for different $\zeta _i$ to better quantify the effects of this configuration parameter. In figure 5, we present $\zeta _{100}$ for the same values of $S$ as in figure 4. In addition to the results for $\zeta _i = 0.25$ (figure 5a), we also show $\zeta _{100}$ for $\zeta _i = 3.00$ in figure 5(b). Both figures reveal that the repulsive and attractive behaviour becomes more pronounced for larger and smaller values of $S$, respectively. An exception can be seen for $S=0.07$ with $\zeta _i=0.25$ in figure 5(a). This is caused by the work required to move the fluid inside the gap for particle pairs at very small $\zeta _t$, i.e. lubrication forces. Such behaviour indicates that there is a peak for small initial gap sizes where the attracting behaviour is most efficient. For the case of $\zeta _i=0.25$, that peak was determined at about $S = 0.35$. In addition, we found that the dynamics of the particles not only depend on $S$, but also on $\zeta _i$. Comparing the two initial gap sizes $\zeta _i=0.25$ and $\zeta _i=3.00$ in figure 5(a,b), we can conclude that the particle interaction becomes less intense for larger $\zeta _i$. Note the scale difference by one order of magnitude in figure 5(a,b). For $\zeta _i=3.00$ in figure 5(b), the damping of particle motion towards each other by lubrication completely vanishes for low frequencies and the threshold for $S$ for the transition from attractive to repulsive behaviour shifts to smaller values of $S$.

Figure 5. Plots of $\zeta _{100}$ as a function of $S$ with $\rho _s=4.68$ for (a) $\zeta _i = 0.25$ and (b) $\zeta _i = 3.00$. The sketches in the insets illustrate the initial configuration and the direction of oscillation (black arrows). Note the difference in scale by one order of magnitude for the $y$ axes of the two panels. Condition $\zeta _{100} < 0$ represents attraction and $\zeta _{100} > 0$ repulsion.

The analysis for particles of relative density $\rho _s=4.68$ presented so far suggests that there is not only an optimum for attractive particle interaction for $\zeta _i\approx 0.25$, but also a threshold condition for which the interaction transitions from attractive to repulsive behaviour. Both of these phenomena are a function of $S$ and $\zeta _i$. To elucidate the effect of $\zeta _i$ on the particle behaviour further, we evaluate $\zeta _{100}$ for a constant frequency ($S=2.09$) and different values of $\zeta _i$ in figure 6. Again, we use an initial orientation of $\theta _i = 0^{\circ }$ and a relative particle density of $\rho _s=4.68$ for this analysis. For such high frequencies, all presented cases result in repulsion. These results demonstrate that the increase of the particle distance decays exponentially for larger values of $\zeta _i$. Again, the case of $\zeta _i = 0.10$ deviates from this behaviour, because for such small gaps, lubrication forces become an important component of the particle interaction. Nevertheless, for $\zeta _i>0.1$, the intensity of the mutual interaction decreases with increasing $\zeta _i$ and vanishes for large initial particle distances. This is in line with the general assumption by Fabre et al. (2017) that the closer the particles are to each other, the more they are able to influence one another by interacting with each other's flow structures. Conversely, this means that the further apart the particles are, the less intense the interaction becomes and the particle motion resembles that of an individual particle as $\zeta _i$ increases.

Figure 6. Plot of $\zeta _{100}$ as a function of $\zeta _i$ with $\rho _s=4.68$ and for $S = 2.09$.

This behaviour is exemplified by the comparison of the particle amplitudes of the two-particle set-ups $A_{2p}(S,\zeta _i)$, in dependence of $\zeta _i$, and the set-ups with individual particles $A_{1p}(S)$ for the same values of $S$ in the non-inertial frame of reference. The comparison is shown by the ratio $\sigma = A_{2p}(S,\zeta _i) / A_{1p}(S)$ for $T = 11$ in figure 7. The time for the calculation of $\sigma$ was chosen so that the particle arrangements still correspond as closely as possible to their initial state. The decision for this instant is further explained in § 4.1. The comparison reveals that the deviation of the particle amplitudes between the two-particle and the individual-particle set-ups increases with decreasing $\zeta _i$. This is in line with the observation and assumption made above that for the given set-up the interaction of two particles increases with decreasing $\zeta _i$. As an important consequence, this analysis also demonstrates that the closer the particles are, the more amplified their excursion lengths become.

Figure 7. Amplitude ratio $\sigma = A_{2p}(S,\zeta _i) / A_{1p}(S)$ between the amplitudes of the two-particle set-ups $A_{2p}(S,\zeta _i)$ and the respective individual-particle set-ups $A_{1p}(S)$ for $\rho _s = 4.68$ and as a function of $S$ and $\zeta _i$.

3.3. Particle arrangements

In order to illustrate the influence of the arrangement orientation, we first look at the behaviour of two monodisperse particles with $\rho _s = 4.68$, whose initial distance is ${\zeta _i = 0.25}$ and the initial angle of the orientation of the particle alignment with respect to the direction of oscillation $\theta _i$ is varied. We modify $\theta _i$ in steps of $15^{\circ }$ in the range from $0^{\circ }$ to $90^{\circ }$ and analyse the change of $\theta$ and $\zeta$ over a period of $100$ oscillations. Here, $\theta$ represents the alignment angle over time. For all set-ups, the particles are arranged so that the midpoint between the particles coincides with the centre of the domain. For the current analysis, we focus exclusively on the two oscillation scenarios, $S = 0.35$ and $S=2.09$, which result in attraction and repulsion, respectively, for a horizontal alignment (cf. § 3.2). The results for different $\theta _i$ are shown in figures 8 and 9, where the temporal evolution of the alignment angle $\theta$ is shown in figures 8(a) and 9(a) and the change in particle distance $\zeta$ is displayed in figures 8(b) and 9(b).

Figure 8. Evolution of (a) $\theta$ and (b) $\zeta$ over $100$ oscillation periods for $\rho _s = 4.68$, $\zeta _i = 0.25$ and $S = 0.35$.

Figure 9. Evolution of (a) $\theta$ and (b) $\zeta$ over $100$ oscillation periods for $\rho _s = 4.68$, $\zeta _i = 0.25$ and $S = 2.09$. Note the difference in scale for the $y$ axis in (b) compared with figure 8.

The results illustrate that $\theta$ increases and converges towards a value of $90^{\circ }$ for particle configurations with $\theta _i \neq 0^{\circ }$, whereas $\theta =\theta _i$ for $\theta _i = 0^{\circ }$. The inter-particle distances $\zeta$, on the other hand, develop significantly differently depending on $S$ and $\theta _i$. In figure 8 ($S = 0.35$), the particles for initial arrangements of $\theta _i \le 30^{\circ }$ approach each other, meaning $\zeta < 0$, while particles move away from each other for $\theta _i \ge 60^{\circ }$. Interestingly, the particles initially behave attractively and then repulsively for $\theta _i = 45^{\circ }$. For $S = 2.09$ (figure 9), similar phenomena occur, whereby in this case smaller angles, i.e. $\theta _i \le 15^{\circ }$, show a repulsive and larger angles, i.e. $\theta _i \ge 60^{\circ }$, an attractive behaviour. Here, both $\theta _i = 30^{\circ }$ and $45^{\circ }$ show a reversal, both from repulsive to attractive. An indication of the reversal can already be seen at $\theta _i = 15^{\circ }$, while the behaviour in the considered time frame is still repulsive.

The analysis presented above shows the following. On the one hand, our data confirm the observations by Lyubimov et al. (2001), Voth et al. (2002), Wunenburger, Carrier & Garrabos (2002) and Klotsa et al. (2007, 2009) that two particles exposed to oscillatory flow align themselves perpendicularly to the direction of oscillation for any perturbation, e.g. in terms of $\theta _i$, that breaks the symmetry of the oscillation direction. The well-developed stage for the alignment angle $\theta = 90^{\circ }$ was further investigated by Klotsa et al. (2007), Fabre et al. (2017), Jalal (2018) and Van Overveld et al. (2022), who all demonstrated that two particles continue to move at a relative distance which remains constant over time. On the other hand, Fabre et al. (2017) were able to prove that this equilibrium only occurs for $S < 2.23$, otherwise the particles approach each other until they come into contact. This is indeed the case for our simulations that do not exceed $S=2.09$. According to our results, the closer the initial alignment angle is to zero, the more the direction of motion tends to be aligned with the oscillation direction. This behaviour can be utilized in microfluidic devices, e.g. to induce targeted particle motion, as was investigated in a recent study by Dietsche et al. (2019). In the following, we therefore focus on the specific configuration with an alignment angle $\theta = 0^{\circ }$ in order to investigate the parameter ranges for which particles can be brought into contact or separated from each other.

3.4. Interaction behaviour

The analysis presented in the previous section shows that there are threshold conditions for the interaction of two particles and the oscillating flow that change the particle dynamics from an attractive to a repulsive behaviour depending on the initial gap size and the dimensionless frequency of the oscillations. It is now interesting to investigate whether a systematic behaviour can be identified that allows one to determine these threshold conditions for a wide range of the parameters $S$, $\zeta _i$ and $\rho _s$. To this end, we conducted a set of simulations with $\epsilon = 0.1$ for $S = [0.07, 2.10]$ with increments of $\Delta S=0.035$, $\zeta _i = \{0.10, 0.25, 0.50, 0.75, 1.00, 1.50, 2.00, 3.00\}$ for the three density ratios considered, $\rho _s = \{0.47, 1.78, 4.68\}$. Recall that these values of $\rho _s$ correspond to the same conditions investigated by L'Espérance et al. (2005).

This campaign allowed us to construct the regime maps shown in figure 10, where right-pointing triangles ($\vartriangleright$) indicate that the particles are approaching each other, i.e. $\zeta _{100} < 0$, while left-pointing triangles ($\vartriangleleft$) refer to set-ups where the particles are drifting apart, i.e. $\zeta _{100} > 0$. We furthermore categorize the mutual interaction as attractive ($\blacktriangleright$, red), transitional (open symbols, $\vartriangleleft \vartriangleright$) and repulsive ($\blacktriangleleft$, blue), depending on $S$ and $\zeta _i$. This categorization is derived from the analysis of $A_{2p}(S,\zeta _i)$ discussed in § 3.2. We define arrangements in which $| \zeta _{100} | \leq A_{2p}(S,\zeta _i)$ as transitional, whereas set-ups with exceeding $| \zeta _{100} | > A_{2p}(S,\zeta _i)$ are classified as either attractive or repulsive, depending on the sign of $\zeta _{100}$. The set-ups classified as transitional represent cases in which the interaction of the flow fields has only minor effects on the respective particles. Although the particles move either towards or away from each other, such cases do not exhibit significant changes of the particle spacing over time and can therefore be described as a state of equilibrium.

(3.1)\begin{equation} | \zeta_{100} | \begin{cases} > A_{2p}(S,\zeta_i) \begin{cases} = \text{attraction}, & \text{for } \zeta_{100} < 0,\\ = \text{repulsion}, & \text{for } \zeta_{100} > 0,\\ \end{cases} \\ \leq A_{2p}(S,\zeta_i) = \text{transition}. \end{cases} \end{equation}

Despite the low values of $\zeta _{100}$ that are assigned to transitional behaviour, a more precise distinction into transitionally attractive and transitionally repulsive remained possible as is shown in § 4.5. This is also expressed in the clear definition of the threshold condition shown as a dashed line in figure 10.

Figure 10. Regime maps of the fluid–particle interactions as a function of $S$ and $\zeta _i$ in log–log plots for different particle density ratios: (a) $\rho _s = 4.68$, (b) $1.78$ and (c) $0.47$. Right-pointing triangles ($\vartriangleright$) indicate $\zeta _{100} < 0$ and left-pointing triangles ($\vartriangleleft$) $\zeta _{100} > 0$. The colouring of the symbols is according to the classification stated in (3.1), where $\blacktriangleright$ (red) represents attraction, $\blacktriangleleft$ (blue) repulsion and open symbols ($\vartriangleleft \vartriangleright$) transition. The dashed line has been determined by best fit of a power law and marks the threshold condition for which the particle interaction transitions from attractive to repulsive.

For $\rho _s = 4.68$ in figure 10(a), our classification scheme reveals that for all arrangements of $\zeta _i \le 1.00$, there is a clear trend that smaller $S$ lead to attraction and larger $S$ to repulsion. The threshold condition from attractive to repulsive results is covered by set-ups classified as transitional. With increasing $\zeta _i$, the threshold conditions shift to lower values of $S$ and the cases classified as transitional cover a wider value range of $S$. For initial particle gaps $\zeta _i = 3.00$, no significant attractive or repulsive behaviour was observed and all cases fall into the transitional regime. Again, this confirms that the effectiveness of the mutual interaction decreases with increasing particle distance.

The regime maps of $\rho _s = 1.78$ and $0.47$ in figures 10(b) and 10(c), respectively, are almost identical with only a few deviations. The similarity is due to the approximately similar deviation of their density from the neutrally buoyant case, i.e. $\rho _s=1$. However, comparing them with the case $\rho _s=4.81$ shown in figure 10(a), we find that the magnitude for which particles approach each other or drift apart decreases substantially for the two cases that are closer to neutrally buoyant conditions. In fact, the small effect of particle inertia diminishes such that the particle interaction yields only marginally attractive cases that we classify as transitional according to the criterion given by (3.1). In general, the majority of set-ups lead to transitional results for those two density ratios, where the number of repulsive cases are also fewer than for $\rho _s = 4.68$. For example, the arrangements with larger initial gaps ($\zeta _i \geq 2.00$) do not lead to any significant repulsive behaviour, so that all arrangements are classified as transitional. Nevertheless, even though the range of intensities of the particle interaction becomes weaker as inertia effects of the particles are decreased, a very similar trend is observed for the region in which the scenario changes from attractive to repulsive as $S$ and $\zeta _i$ are increased.

As mentioned above, the dashed line in each regime map marks the transition from attractive to repulsive behaviour as indicated by a change of sign in $\zeta _{100}$ shown in figure 4. For figure 10, this threshold condition was determined by best fit of a power law given by $\zeta = C_1 S^{\beta }$, for which we obtain an excellent correlation with our data. For the highest particle density ratio $\rho _s=4.81$ (figure 10a), we determine $C_1=0.18$ and $\beta =-2.73$. For the two cases $\rho _s=1.78$ and $\rho _s=0.47$, we obtain identical fitting parameters $C_1=0.15$ and $\beta =-2.93$. This illustrates two things: (i) an increase of inertia does not substantially change the threshold conditions to transition from attractive to repulsive behaviour and (ii) as long as the deviation from the neutrally buoyant condition remains the same in terms of particle inertia, the interaction regime is unaffected, even though we change particle properties to denser or lighter than the ambient fluid.

The deviation in particle behaviour with respect to $\rho _s$ is hence insignificant. Therefore, our results can be directly compared with the force maps presented by Fabre et al. (2017), which correlate forces in the context of two stationary particles in an oscillating flow. In their study, they investigated a wide range of oscillations ($S = [0.01, 11.10]$) and particle distances ($\zeta _i = [0, 2.0]$) and due to the immobilization of the particles, the parameter $\rho _s$ did not play a role. In principle, our analyses agree in that small $S$ and small $\zeta _i$ generally lead to a decrease of $\zeta$, while large $S$ and large $\zeta _i$ result in an increase of $\zeta$ over time. In fact, one can convert our ($S,\zeta _i$) regime map to the regime map of Fabre et al. (2017) (figure 6b in that reference) and we find very close agreement for the transition for $S=0.56$ and $S=0.89$ (termed $\varOmega =5$ and $\varOmega =8$ in that same reference).

4. Flow field analysis

4.1. Steady streaming

The previous results have shown that two axially arranged particles in an oscillating fluid flow either attract or repel each other, depending on the applied frequency and the initial distance. To better understand the mechanisms that lead to a decrease or increase in the particle distance, we average the flow field over one oscillation period to focus on the time-averaged component of the fluid flow that arises due to the presence of particles. In this regard, we apply a method of intrinsic averaging, where only the volume occupied by fluid is considered (Vowinckel et al. 2017). Such an averaged flow field is commonly referred to as steady streaming (Riley 1966).

As briefly introduced in § 1, it represents the non-zero mean flow of the periodically averaged fluid field generated by the interaction of objects and the surrounding fluid, where either the fluid or the object can induce the oscillations (Riley 2001). Such an analysis is particularly useful for high-frequency oscillations, because it allows one to study the well-developed flow patterns of the steady streaming induced by the inertial particles in the oscillating domains. However, two important aspects have to be taken into account when choosing the time for the calculation of the steady streaming. First, the chosen time should be at the initial phase of the simulation, so that the distance between the particles is as close as possible to the original arrangement. Second, the effect of the initial displacement of the particles in a fluid at rest (see § 3.1) must have already decayed and thus be negligible. Consequently, a compromise must be made which best fulfils both requirements. To this end, we have performed a convergence study in which we analysed the impact of the initial conditions on the flow structures of the steady streaming. This study showed that a well-developed state in the oscillatory motion is reached at time $t/T\geq 11$ for each simulation. This is why all the flow characteristics shown hereafter refer to $t/T = 11$. As mentioned earlier in § 2.2, the scenario under consideration does not clearly fall into the case distinction for a single particle proposed by Riley (1966). It is therefore important to note that all the following flow structures, which consider the set-up parallel to the direction of oscillation, consist exclusively of what has been referred to as primary vortices surrounding the respective particle.

In figure 11, the flow structures of the steady streaming are shown for two representative cases of attractive ($S = 0.35$; figure 11a,b) and repulsive ($S = 1.05$; figure 11c,d) behaviour using $\zeta _i = 0.75$ and $\rho _s = 4.68$. For this purpose, we show streamlines in the $xy$ and $zy$ planes. The planes cut through the centre of the domain in the third respective dimension, i.e. at $z = 5$ for the $xy$ plane and at $x = 5$ for the $zy$ plane. Since the $zy$ plane cuts through a plane in between the particles, we show a hollow black circle at this location in figures 11(b) and 11(d) to indicate the particle positions aligned along the $x$ axis. We note, however, that this circle does not provide an obstacle to the flow for these figures. We also note that in all figures showing streamlines, these lines are uniformly distributed and serve only to qualitatively illustrate the flow patterns and directions. They do not represent a value of a streamfunction to visualize flow intensity. In addition, we calculate the polar component of the vorticity $\omega _z(x,y )$ to provide a quantitative measure. The vorticity contours are shown in blue $(-0.01)$ and red $(0.01)$, representing clockwise and counterclockwise rotation, respectively. This colour scheme corresponds to the colour scale in figure 2 and is applied for all subsequent analyses.

Figure 11. Streamlines and vorticity contours of the steady streaming in (a,c) $xy$ plane and (b,d) $zy$ plane. (a,b) An attractive case with $S = 0.35$ and (c,d) a repulsive behaviour with $S = 1.05$. Both are arrangements with $\zeta _i = 0.75$ and $\rho _s = 4.68$. The black circles in the middle of the $zy$ planes indicate the particle positions, which are arranged in line along the $x$ axis. The outer circle in (b) is a stagnation line of the flow field and results from the averaging process of computing the steady streaming. Colour scheme as in figure 2.

The comparison of the $xy$ planes for the repulsive case and the attractive case (figure 11a,c) shows that in both cases, each individual particle is surrounded by four vortex structures, which are depicted by the tori of the streamlines and referred to as quadrupoles, as described in § 2.3. Consequently, two quadrupole structures are formed in each set-up. However, instead of assigning the vortices to the particles, we examine the structures based on their spatial position and divide them into four central and four peripheral vortices (cf. the sketch in figure 2b). This distinction allows for a more detailed description and analysis of the flow mechanisms in the present and subsequent cases. The central vortices, or more precisely, the vortex cores of the central vortices are located between the two particles. The peripheral vortices, on the other hand, are located on the outward-facing sides. Hence, for the set-up considered here, the central vortices cause fluid to flow into the gap, tending to push the particles away from each other (repulsion), while the peripheral vortices direct the flow along the horizontal axis of symmetry towards the outward-facing sides of the particles, pushing the particles towards each other (attraction). Although the double quadrupole structure of the steady streaming is formed in both the attractive (figure 11a) and repulsive (figure 11c) cases, the streamlines differ significantly in their spatial pattern. Nevertheless, neither of the two cases shows a clear trend with respect to the formation of the central and peripheral vortices. As was shown in § 3.1, this less distinct flow pattern is due to the decaying effects of the initial condition. Qualitatively, however, these figures provide an indication that the size of the peripheral vortices is predominant for attraction and vice versa the central vortices are larger for repulsion. While the far-field flow patterns show the aforementioned differences, the vorticity contours in the near field of the particles show a strong similarity. Next to the particle surfaces, four vorticity contours emerge whose shapes mimic the curvature of the sphere. These are followed by another four contours, each with reversed sign. For each individual particle, the diagonally opposite central as well as peripheral contours have the same direction of rotation.

The $zy$ planes also differ in the structures of their streamlines. In the case of attraction (figure 11b), there is a separation of the flow direction, which is marked by the outer circle. This circle is a result of the calculation of the steady streaming and represents a stagnation line where the flow is directed neither towards nor away from the centre. The inner part of the circle indicates a convergent inflow, where the streamlines point towards the centre between the particles. For the part outside of the stagnation line, the flow is directed outwards away from the particles. For the repulsive behaviour (figure 11d), the entire flow of the $zy$ plane is directed towards the centre of the particles. The vorticity in the $zy$ plane is not discernible in either of the two cases shown, which is related to the fact that the flow components in the $y$ and $z$ directions, i.e. orthogonal to the direction of oscillation, are relatively minor and consequently the vorticity resulting from these components is very small as well.

The flow patterns depicted in figure 11 show a strong axisymmetric structure. This is true both for attraction as well as for repulsion. We therefore conclude that the $xz$ plane (not shown here) and the $xy$ plane, both arranged through the centre of the particles with $y=5$ and $z=5$, respectively, show the same flow properties and structures of the steady streaming. For this reason we will limit our analysis to only one of these two planes, namely the $xy$ plane, as it provides sufficient information to study the fluid–particle interaction and the resulting particle behaviour.

4.2. Flow decomposition

The detailed illustrations shown in figure 11 already indicate that there are significant differences in the flow structures of attractive and repulsive set-ups. However, identifying and quantifying a clear physical mechanism that allows for conclusions regarding attractive or repulsive behaviour of the particles are not possible due to the apparent overlap of different flow structures. For this reason, we perform a more detailed examination in terms of flow decomposition to better evaluate the vorticity generated by the fluid–particle interactions.

Following the argument of Cerretelli & Williamson (2003), the vorticity computed from steady streaming $\omega _z(x,y )$, hereafter also referred to as total steady streaming, allows us to decompose the flow field into a symmetric and an antisymmetric part denoted as $\omega _{z,S}(x,y )$ and $\omega _{z,A}(x,y )$, respectively:

(4.1)\begin{equation} \omega_z \left(x,y \right) = \omega_{z,S} \left(x,y \right) + \omega_{z,A} \left(x,y \right) . \end{equation}

The symmetric part is generated by calculating the symmetry about the centre of the distance between the particles, based on averaging the left and right components of the total steady streaming. Subtracting this part from the total steady streaming yields the antisymmetric component. To this end, we introduce the transform $x' = x - x_c$, where $x_c$ is the $x$ coordinate at the centre between the two particles. Therefore, $\omega _z(x',y )$ can be written as

(4.2)\begin{equation} \omega_z \left(x',y \right) = \underbrace{\frac{1}{2} \left[\omega_z \left(x',y \right) + \omega_z \left({-}x',y \right) \right]}_{=\omega_{z,S}} + \underbrace{\frac{1}{2} \left[\omega_z \left(x',y \right) - \omega_z \left({-}x',y \right) \right]}_{=\omega_{z,A}} . \end{equation}

In figure 12, the streamlines and vorticity contours of the symmetric and antisymmetric components are shown in the panels of the left- and right-hand columns, respectively. The set-ups are the same as given in figure 11. In this regard, figure 11(a,c) represents the total steady streaming used for the decomposition (4.1). When considering the streamlines of the symmetric part (figure 12a,c), two vortex structures with centre points in the far field orthogonal to the particle arrangement and in line with the centre of the particle distance can be seen for the attractive as well as repulsive behaviour. In both cases, the upper vortex rotates in the counterclockwise direction and the lower one in the clockwise direction, but the position of the stagnation points is at smaller distance to the particle surfaces for repulsion (figure 12c) than it is for attraction (figure 12a). The vortex contours also differ in magnitude and spatial extent. While almost no contours are visible for attraction (figure 12a), a distinct layer near the particles is present for repulsion. In general, however, the symmetric component discloses the resulting motion of the particle pairs and, along with it, the basic direction of motion of the flow field, which is from left to right in both arrangements presented (cf. figure 12a,c). The net motion of the particles is due to the transients that are still present, caused by the initial condition of the still fluid, as described in § 3.1. Since the initial displacement is larger for $S=1.05$ than for $S=0.35$, there is a greater deviation from the centre of the domain at the time of calculating the steady streaming. This results in an stronger motion towards the centre of the domain for $S=1.05$, which in turn produces more intense vortex contours. Even though this has no influence on the respective particle behaviour regarding attraction or repulsion, it is responsible for the flow field of the symmetric component of the steady streaming.

Figure 12. Streamlines and vorticity contours of the symmetric (a,c) and antisymmetric (b,d) components for $\zeta _i = 0.75$ and $\rho _s = 4.68$. (a,b) Parameter $S = 0.35$ resulting in attraction and (c,d) $S = 1.05$ leading to repulsion. The total steady streaming used for decomposition is given in figure 11(a,c). Colour scheme as in figure 2.

The antisymmetric fields (figure 12b,d) represent the deviations between the flow fields of the total steady streaming and the symmetric component. It thus illustrates the filtered flow field that is generated solely by the pairwise fluid–particle interaction. This flow component actually dictates the governing motion of the particles with respect to one another, because this is the motion that deviates from the underlying oscillatory dynamics. In this regard, the vorticty contours of the antisymmetric parts are very similar to the contours of the total steady streaming. This is due to the fact that the vorticty of the previously described symmetric component is small for both cases, attraction and repulsion. This analysis already reveals that the particle behaviour of attraction or repulsion is caused by small flow conditions deviating from the oscillating main flow. These small flow deviations are clearly visible in the streamlines of the antisymmetric flow fields. As mentioned above, these flow structures cause either repulsive or attractive behaviour and they compete in size depending on $S$. A qualitative analysis of the shape and size of the vortex structures underlines the findings of § 4.1, in which the peripheral vortex structures dominate over the central ones for attraction and vice versa for repulsion.

4.3. Effect of particle inertia

Since the previous figures have shown results for $\rho _s = 4.68$ only, we now turn our attention to the steady streaming observed for two additionally considered particle density ratios $\rho _s = 1.78$ and $0.47$ that were already analysed in terms of their particle interaction in figure 10. We therefore performed the same analysis of flow decomposition as for $\rho _s = 4.68$ given in the previous section. Figure 13 shows the corresponding results of steady streaming for $S = 0.35$ and $\zeta _i = 0.75$ for these two density ratios. In this way, we can compare the flow structures and vortex contours of these two set-ups ($\rho _s = 1.78$ in figure 13a,c,e and $\rho _s = 0.47$ in figure 13b,d, f) with the previous one ($\rho _s = 4.68$ in figures 11a and 12a,b) and can thus infer the behaviour of two particles in oscillating flow that are both denser and lighter than the surrounding fluid. Since we neglect gravity, the particle dynamics is affected by the particle mass only in terms of particle inertia, but not weight.

Figure 13. Streamlines and vorticity contours of the total steady streaming (a,b), and its symmetric (c,d) and antisymmetric (e, f) components for a set-up with $\zeta _i = 0.75$ and $S = 0.35$. Here, $\rho _s = 1.78$ (a,c,e) and ${\rho _s = 0.47}$ (b,d, f). Colour scheme as in figure 2.

According to the regime maps shown in figure 10, for $S = 0.35$ and $\zeta _i = 0.75$, all $\rho _s$ result in an attractive behaviour of the particles. Only $\rho _s = 4.68$ is assigned to the category of attraction, while $\rho _s = 1.78$ and $0.47$ are associated with transitional behaviour according to (3.1). The comparison of the total steady streaming (figures 11a and 13a,b) indicates the same formations of vorticity contours with respect to the clockwise and counterclockwise rotation directions. For each $\rho _s$, each respective particle is surrounded by four adjacent and four far-field contours. The intensity of the vorticity decreases with decreasing $| \rho _s - 1 |$, which serves as a measure of the difference between the fluid and particle inertia and corresponds to the strength of the fluid–particle interaction. As expected, these results show that the smaller the difference between the particle density and the fluid density, the smaller the inertial effect, and thus the particle and fluid motion are more in agreement. Consequently, less intense flow develops, which is recognizable by lower values of vorticity as indicated by the blue and red contours. Despite the lower vorticity, an examination of the streamlines of the total steady streaming velocity field yields that $\rho _s = 4.68$ and $1.78$ reveal very similar patterns of the respective vortex structures. The vortex structures of $\rho _s = 0.47$ show the same patterns as for $\rho _s = 1.78$, but in a mirror-inverted arrangement.

The same observation holds for the streamlines of the symmetric component for the three density ratios shown in figures 12(a), 13(c) and 13(d), respectively. For $\rho _s = 4.68$ and $1.78$, i.e. particles with higher density than the fluid, the upper vortex rotates counterclockwise and the lower one clockwise. For particles less dense than the fluid, i.e. $\rho _s = 0.47$, it is the other way around. This is due to the fact that the direction of action changes sign for particles with $\rho _s > 1$ and $\rho _s < 1$. Consequently, particles denser than the surrounding fluid and, hence, larger inertia tend to follow the fluid motion in a reduced and delayed manner, resulting in an initial movement to the left in a non-inertial frame for the present scenario. Less dense particles, on the other hand, have a lower inertia than the fluid and accelerate faster than the fluid, so that they overshoot and move to a direction opposite to those of denser particles. After the initial displacement, particles slowly migrate back to the centre of the domain, which causes the background flow in the symmetric component (cf. § 3.1). The vorticity for the symmetric component, however, remains small for all $\rho _s$ so that their contours do not show up in figures 12(a) and 13(c,d).

Looking at the antisymmetric component shown in figures 12(b) and 13(e, f), we find the vorticity contours to be identical for all three density ratios. This was already observed for the total steady streaming (cf. figures 11a and 13a,b), but here the mirror-inverted arrangement was filtered out by the symmetric component. This suggests that even though $\rho _s = 1.78$ and $0.47$ are classified as transitional, because they are only marginally attractive, and $\rho _s=4.68$ is strongly attractive, the mechanism that drives particle motion is dictated by the non-dimensional frequency $S$ and remains unaffected by the particle density ratio. This result is consistent with the findings of Van Overveld et al. (2022).

4.4. Effect of oscillation frequency and initial particle distance

The analysis of the decomposition of the total steady streaming demonstrated that the antisymmetric component governs the flow structures that drive the pairwise fluid–particle interactions. Therefore, all subsequent illustrations and analyses refer to the antisymmetric part only. While the flow structures seem to be less influenced by the density ratio $\rho _s$, the applied frequency $S$ appears to have a major effect. This could already be seen in figures 12 and 13 as well as in the regime maps, presented in figure 10. To analyse whether the change of $S$ has an impact on the resulting flow patterns of the antisymmetric part, we show a comparison of the flow field for three different characteristic frequencies. The chosen frequencies ${S = \{0.35, 0.70, 1.05\}}$ correspond to attraction, transition and repulsion, respectively, for the given $\zeta _i = 0.50$ and $\rho _s = 4.68$. The resulting streamlines and vorticity contours are depicted in figure 14, with $S = 0.35$ in figure 14(a), $S = 0.70$ in figure 14(b) and $S = 1.05$ in figure 14(c). The flow structures of the attractive and repulsive set-ups (figure 14a,c), respectively, are again very similar to the corresponding flow structures observed for larger initial gap size $\zeta _i = 0.75$ (figure 12c, f). Hence, for attraction (repulsion) the peripheral (central) vortices dominate.

Figure 14. Streamlines and vorticity contours of the antisymmetric part of the flow field for a constant $\zeta _i = 0.50$ and varying $S$: (a) $S = 0.35$ with attractive behaviour, (b) $S = 0.70$ (transition) and (c) $S = 1.05$ (repulsion). Colour scheme as in figure 2.

Another interesting aspect arises if one considers the transition of the respective vortices from attraction via the transition to repulsion. While, qualitatively, the central vortex structures are completely suppressed in the case of attraction, they appear to be in equilibrium when the transition is established, and become clearly superior in the case of repulsion. Naturally, we find the corresponding opposite trend for the peripheral vortices. The alteration of the central and peripheral vortex sizes can also be observed on the basis of the vorticity contours. Even though the changes of the contours vary only slightly in the course of the increase of $S$, it can be seen that the size of the central vortices increases while that of the peripheral ones decreases.

As already indicated by the regime map given in figure 10, the behaviour of the two particles in oscillation is dictated by frequency and the initial distance of the two particles. We, therefore, also examine the impact of $\zeta _i$ on the flow characteristics in figure 15 by comparing $\zeta _i = \{0.10, 0.25, 0.50\}$ for $S = 0.94$ and $\rho _s = 4.68$. Here, $\zeta _i = 0.10$ results in attraction (figure 15a), $\zeta _i = 0.25$ leads to a transitional state (figure 15b) and $\zeta _i = 0.50$ yields repulsion (figure 15c). At first glance, the vorticity contours of the three arrangements appear to be almost identical, and the flow patterns also show a close resemblance with predominant central and suppressed peripheral vortex structures, respectively. A closer look at the peripheral flow structures, however, shows a slight decrease of peripheral vortex size with increasing $\zeta _i$. This is accompanied by the fact that with increasing $\zeta _i$, more fluid can flow into the gap, increasing the strength of the flow into the gap trying to push the particles apart.

Figure 15. Streamlines and vorticity contours of the antisymmetric part of the fluid field for constant $S = 0.94$ and varying $\zeta _i$: (a) $\zeta _i = 0.10$ resulting in attraction, (b) $\zeta _i = 0.25$ leading to transition and (c) $\zeta _i = 0.50$ causing repulsion. Colour scheme as in figure 2. The dashed lines in (c) indicate the symmetry lines used for the division into four quadrants to compute $\omega _{z,\alpha }$ (4.3).

4.5. A circulation-based criterion to determine the particle interaction

A qualitative examination of the antisymmetric component of the steady streaming shows that the visualization of the flow structures can be helpful in determining the behaviour of two particles in oscillatory flows. Owing to the low flow intensities, however, it is challenging to make a clear distinction of repulsion and attraction. This was already demonstrated in the context of analysing different density ratios, where we find the same flow patterns for strongly and weakly attractive behaviour for dense and less dense particles, respectively. To come to a more general conclusion, a more precise and quantitative measure for the interactions is, therefore, needed. For this purpose, we compute the circulation of the antisymmetric part of the flow field separately for central and peripheral vortices.

In a first step, we assign each vorticity contour value of the $xy$ plane to contribute to either attractive or repulsive behaviour based on the considerations explained in § 4.1 and sketched in figure 2(b). Since this cannot be done on the basis of the direction of rotation only, we divide the flow field into quadrants as shown in figure 15(c), where the dashed lines separate the quadrants and the quadrant numbering is introduced by the Roman numerals. The intersection of the two lines coincides with the stagnation point between the particles. Starting from this point, the subdivision of the field follows the axes of symmetry. The assignment of the vorticity contours is done for each quadrant separately and takes the contribution of the respective vorticity, $\omega _{z,\alpha }$, into account. Here, $\alpha$ is used as an index for either attraction or repulsion. According to figure 2(b), central vortices push fluid inside the gap and act repulsively whereas peripheral vortices push against the particles from the outside so that particles tend to approach each other. With the help of this assignment, we calculate the circulation $\varGamma _\alpha$ contributing to attraction or repulsion as

(4.3)\begin{equation} \varGamma_\alpha = \iint_{A} \omega_{z, \alpha} \,{\rm d}A . \end{equation}

On this basis, we can compute the ratio of attractive to repulsive circulation ${\phi _\varGamma = \varGamma _{attr} / \varGamma _{rep}}$ and, thus, draw conclusions about the behaviour of two particles in oscillating fluid flow. A ratio of $\phi _\varGamma = 1$ would represent a balance of the central and peripheral vortices, $\phi _\varGamma > 1$ indicates a predominantly attractive circulation leading to the approach of the two particles and $\phi _\varGamma < 1$ refers to repulsive behaviour due to dominant circulation pushing the particles apart.

Figure 16 presents $\phi _\varGamma$ for the considered density ratios $\rho _s$. The same set of $S$ is used as in figure 4. The dashed lines highlight the balance of the circulations in each graph, i.e. $\varGamma _\alpha =1$, and the symbols correspond to the regime map in figure 10. As already shown in § 4.3, the influence of particle inertia on the behaviour of the two particles in terms of the qualitative distinction of either attraction or repulsion is small. Consequently, the ratios for the three densities $\rho _s = 4.68$ (figure 16a), $1.78$ (figure 16b) and $0.47$ (figure 16c) do not differ much, if we choose to normalize the attractive circulation by the respective repulsive circulation of the same case. The trend that small $S$ together with small $\zeta _i$ lead to attractive results is reflected by the fact that these frequencies lead to $\phi _\varGamma > 1$, whereas we find $\phi _\varGamma < 1$ for repulsive set-ups for large $S$ and large $\zeta _i$.

Figure 16. Circulation ratio $\phi _\varGamma = \varGamma _{attr} / \varGamma _{rep}$ in dependence of $\zeta _i$ and $S$. The colour scheme represents the initial particle distances: $\vartriangleleft \vartriangleright$ red, $\zeta _i = 0.10$; $\vartriangleleft \vartriangleright$ green, $\zeta _i = 0.25$; $\vartriangleleft \vartriangleright$ dark blue, $\zeta _i = 0.50$; $\vartriangleleft \vartriangleright$ light blue, $\zeta _i = 0.75$; $\vartriangleleft \vartriangleright$ pink, $\zeta _i = 1.00$; $\vartriangleleft \vartriangleright$ brown, $\zeta _i = 1.50$; $\vartriangleleft \vartriangleright$ violet, $\zeta _i = 2.00$; $\vartriangleleft \vartriangleright$ sky blue, $\zeta _i = 3.00$. Right-pointing triangles ($\vartriangleright$) indicate $\zeta _{100} < 0$ and left-pointing triangles ($\vartriangleleft$) $\zeta _{100} > 0$. According to the classification stated in (3.1), filled symbols represent either attraction or repulsion depending on the orientation of the symbol, and open symbols indicate transition. Density $\rho _s = 4.68$ (a), $1.78$ (b) and $0.47$ (c).

In the context of figure 10, it was already shown that the main difference between the considered density ratios is the assignment of different classes, with density ratios closer to neutrally buoyant conditions leading to more cases that we classify as transitional. The circulation-based criterion depicted in figure 16 now provides a solid measure for distinguishing the different behaviour of repulsion and attraction. Indeed, the finding in § 3.2 that the interaction of the flow fields of the respective particles increases with decreasing $\zeta _i$ is clearly evident by the increase in $\phi _\varGamma$ with decreasing $\zeta _i$. For larger $\zeta _i$, $\phi _\varGamma$ decreases to values below unity. This observation now holds for all density ratios considered, where the different cases shown in figure 16 appear to be almost identical.

By means of this quantitative analysis, the conclusion from § 4.3 can be confirmed that the inertia of the particles, characterized by the density, has no effect on the respective behaviour with regard to attraction or repulsion. Thus, this is solely determined by the non-dimensional oscillation frequency $S$ and the initial particle distance $\zeta _i$.

5. Conclusions

To investigate the potential hydrodynamic mechanism that induces flocculation of primary particles in oscillatory flows, a systematic simulation campaign was conducted to analyse the mutual particle behaviour of two freely moving particles subjected to oscillatory fluid flow. In such flows, particles tend to drift based on their inertia, which can cause particle trajectories to deviate from the background oscillatory motion. To this end, particle-resolved simulations were performed for the set-up of two monodisperse mobile particles submerged in a viscous fluid, where gravity was neglected to focus on particle inertia. Three particle density ratios were considered: a large value with pronounced inertia, and two corresponding smaller inertial density values being denser and lighter than the ambient fluid, but approximately equally far from neutrally buoyancy conditions.

In a first step, different orientations of the particle arrangement in relation to the oscillation direction were investigated. The results showed that any arrangement whose orientation angle deviates from the direction of oscillation will be aligned perpendicular to it over time. Previous studies have shown that in the frequency range we investigated, particles whose arrangement is perpendicular to the oscillation reach a state of equilibrium in which the distance between the particles remains constant (Fabre et al. 2017; Jalal 2018). An arrangement in line with the direction of oscillation, on the other hand, remains in its orientation, in which the particles move horizontally towards or away from each other, respectively. On this basis, we have focused in detail on horizontally aligned arrangements and investigated the impact of the oscillation frequency, the initial distance and the particle density on the mutual particle behaviour. To this end, we analysed the distance between the particles over a period of 100 oscillations and created a regime map that defines threshold conditions to distinguish between attractive and repulsive behaviour in dependence on the above-mentioned parameters. The regime map was found to be self-similar for the three different density ratios, whereby the threshold conditions can be described by a power law of approximately the same coefficients. It was also shown that small distances and low frequencies tend to lead to an attraction of the two particles, while in contrast, large distances and higher frequencies lead to a repulsion.

The mutual particle behaviour was then further analysed by visualizing the flow patterns initiated by the fluid–particle interaction. In this context, we have calculated the steady streaming, i.e. the averaging over an oscillation period, and filtered out the impact of the initial state of a fluid at rest by decomposing the flow field into its symmetric and antisymmetric components. The latter component is decisive for the individual motion of the particles and characterized by a quadrupole flow structure surrounding each particle when analysing the plane through the toroidal vortices. Based on this decomposition, the competing effect of the quadrupole structures acting on the particles was revealed and identified as the decisive factor for the respective particle behaviour. On the one hand, the central vortices act to pump fluid into the gap to drive the particles apart. On the other hand, the peripheral vortices act to push particles together from the outside. We have also demonstrated this effect quantitatively by calculating the circulation of the steady streaming and comparing the respective components, i.e. vortex structures that promote either attraction or repulsion of the particles. Our results show that this effect is the governing mechanism for all density ratios investigated in this study and yielded identical results in terms of the competing circulation patterns in all cases. Given the self-similar results of the circulation produced for the three density ratios considered in the present study, the threshold condition for attractive and repulsive behaviour that can be expressed by a power law appears to be universal and solely dependent on the oscillation frequency and the initial particle distance.

Our analysis could therefore open up opportunities to control flocculation of particles based on hydrodynamic effects in oscillating devices with high precision. On the foundation presented here, this would primarily find application in microfluidics, where the number of particles can be controlled and spatial constraints can induce an orientation of the particles that is approximately aligned with the oscillation direction. However, future analyses of this kind should extend the parameter space to oscillation amplitudes and particles of different size. In addition, a larger container size with more particles would extend the presented framework and could provide insights into particle interactions over both shorter and longer ranges in more complex arrangements.