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 Reference Uhlmann2005; Kempe & Fröhlich Reference Kempe and Fröhlich2012; Biegert, Vowinckel & Meiburg Reference Biegert, Vowinckel and Meiburg2017). 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. (Reference Biegert, Vowinckel and Meiburg2017) and Vowinckel et al. (Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019), 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 (Reference Coimbra and Rangel2001), Coimbra et al. (Reference Coimbra, L'esperance, Lambert, Trolinger and Rangel2004) and L'Espérance et al. (Reference L'Espérance, Coimbra, Trolinger and Rangel2005), 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 Reference Riley1966, Reference Riley1967). As stated in § 1, Riley (Reference Riley1966) 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. Reference Biegert, Vowinckel and Meiburg2017; Vowinckel et al. Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019, Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021). As a most relevant case for our scenario, Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017) have obtained excellent agreement for a sphere settling in a large container (Mordant & Pinton Reference Mordant and Pinton2000) and for a settling sphere approaching a wall (Ten Cate et al. Reference Ten Cate, Niewstad, Derksen and Van den Akker2002). 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 Reference Longuet-Higgins1998; Tho, Manasseh & Ooi Reference Tho, Manasseh and Ooi2007). As emphasized by Klotsa (Reference Klotsa2009), 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. (Reference Fabre, Jalal, Leontini and Manasseh2017) and Jalal (Reference Jalal2018) for two stationary spheres and by Tho et al. (Reference Tho, Manasseh and Ooi2007) 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 (Reference Riley1966) (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.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. (Reference Biegert, Vowinckel and Meiburg2017) 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. (Reference Fabre, Jalal, Leontini and Manasseh2017) 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. (Reference Lyubimov, Cherepanov, Lyubimova and Roux2001), Voth et al. (Reference Voth, Bigger, Buckley, Losert, Brenner, Stone and Gollub2002), Wunenburger, Carrier & Garrabos (Reference Wunenburger, Carrier and Garrabos2002) and Klotsa et al. (Reference Klotsa, Swift, Bowley and King2007, Reference Klotsa, Swift, Bowley and King2009) 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. (Reference Klotsa, Swift, Bowley and King2007), Fabre et al. (Reference Fabre, Jalal, Leontini and Manasseh2017), Jalal (Reference Jalal2018) and Van Overveld et al. (Reference Van Overveld, Shajahan, Breugem, Clercx and Duran-Matute2022), 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. (Reference Fabre, Jalal, Leontini and Manasseh2017) 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. (Reference Dietsche, Mutlu, Edd, Koumoutsakos and Toner2019). 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. (Reference L'Espérance, Coimbra, Trolinger and Rangel2005).

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. (Reference Fabre, Jalal, Leontini and Manasseh2017), 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. (Reference Fabre, Jalal, Leontini and Manasseh2017) (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.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. Reference Vowinckel, Nikora, Kempe and Fröhlich2017). Such an averaged flow field is commonly referred to as steady streaming (Riley Reference Riley1966).

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 Reference Riley2001). 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 (Reference Riley1966). 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 (Reference Cerretelli and Williamson2003), 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. (Reference Van Overveld, Shajahan, Breugem, Clercx and Duran-Matute2022).

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$.