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.


Introduction
The behavior 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 et al. 2019(Halfi et al. , 2020) ) and in diesel engines (Katoshevski et al. 2005;Ruzal-Mendelevich et al. 2016;Gupta et al. 2019).Recently, research has also been conducted to investigate how oscillations can be used as a noninvasive biomedical application for drug delivery inside the ear (Sumner et al. 2021;Harte et al. 2023).
Of particular interest for the present study are applications and investigations of fluidparticle 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,b).Studies of the same kind were recently conducted on board the International Space Station to examine the flocculation behavior of a variety of cohesive particles (Kleischmann et al. 2021).Such experimental setups 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 et al. 2018(Mutlu et al. , 2020;;Dietsche et al. 2019).As the particles are usually only a few micrometers 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 setup, the behavior 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. (2005Espérance et al. ( , 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 et al. (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 nonzero, 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 (1966Riley ( , 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, that expands over the whole domain (Rott 1964).In 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 et al. 2007;Otto et al. 2008), the size and the rotational directions of the structures of the steady streaming were analyzed.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 behavior 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(Klotsa et al. , 2009;;Van Overveld et al. 2022, 2023).The latter studies, in which the particles were positioned on the bottom plate of the 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 behavior subjected to oscillating flows, they are limited to initial particle configurations that are perpendicular with respect to the oscillation direction.Moreover, the particle behavior is affected by the presence of the bottom plate and the effect of gravity.Despite the fact that Klotsa et al. (2007Klotsa et al. ( , 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 setup enables us to analyze the behavior 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 analyze the impact of the non dimensional frequency S, the initial interparticle distance ζ i and orientation of the particle alignment θ i , as well as the density Figure 1: Two-dimensional sketch of the cubic numerical domain with size L x,y,z = L = 10 D p .Two mobile spherical particles with identical diameters D p are initially placed so that the center of the initial particle distance ζ i coincides with the center of the domain.The arrangement is aligned with an initial angle θ i with respect to the direction of oscillation.The domain oscillates with velocity u f in x-direction, as indicated by the arrows.ratio ρ s on the mutual particle behavior.This is followed by a detailed analyses of the resulting flow fields in §4 and finally, a summary of the conclusions and a brief outlook in §5.

Numerical setup
The scenario under consideration is a cubic container with dimensions of L x,y,z = L = 10D p , where D p is the particle diameter.The container is filled with a viscous, incompressible fluid characterized by its density ρ f and kinematic viscosity ν f .As illustrated in figure 1, two monodisperse particles are submerged in this container and placed with an initial particle distance ζ i as well as an initial orientation θ i with respect to the horizontal direction of oscillation (black arrows).Initially, the fluid is at rest and the center of ζ i coincides with the center 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 Ω sin Ωt.Here, A f is the distance amplitude and A f Ω the velocity amplitude, t denotes time, Ω = 2πf is the angular frequency and f the frequency of the oscillatory motion of the domain.The walls of the domain are characterized by stress-free boundary conditions, i.e. du t /dn = 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.

Governing equations
We use Direct Numerical Simulations (DNS) 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 et al. 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.2) Here, u = (u, v, w) T represents the fluid velocity vector in Cartesian components, p is the fluid pressure and f IBM the volume force based on the IBM.The 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 u f = (u f , 0, 0) T relative to the inertial frame, we can calculate the non-inertial particle velocity as u ′ p = u p − u f , where u p = (u p , v p , w p ) T presents the translational velocity vector of the particle in an inertial frame.By rearranging, we can state that If we take Newton's equation in an inertial frame 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 and rotation Here, m p and V p are the particle mass and volume, respectively, ρ p the particle density, and I p = πρ p D 5 p /60 the moment of inertia with ω p = (ω p,x , ω p,y , ω p,z ) 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), F h and M h , denote the hydrodynamic forces and torques, respectively.
The fluid-particle interactions can be computed via where τ = −pI + µ f ∇u + (∇u) T represents the hydrodynamic stress tensor that includes the identity tensor I and the dynamic viscosity µ f of the fluid.Furthermore, n is the normal vector pointing outward from the particle surface S p , and r is the position vector pointing from the particle center of mass to a point on S p .The forces and torques due to collision are represented by F c and 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 ζ t 2h, with ζ 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 ℓ and u are the relevant length and velocity scales, respectively, that appear in equations (2.1) -(2.6): 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 /ν f , where u f,max = A f Ω is the velocity amplitude and ρ s = ρ p /ρ f the density ratio (2.12) In addition, we introduce a streaming Reynolds number Re s which is defined by the product of the amplitude ratio ǫ = A f /D p and Re to obtain Re s = ǫRe = A 2 f Ω/ν.Following Coimbra & Rangel (2001), Coimbra et al. (2004) and L Espérance et al. (2005), we also define a Strouhal number Sl = ΩD p /(9u f,max ) for spherical objects to obtain the non-dimensional frequency S = SlRe/4 = D 2 p Ω/(36ν f ).Note that 1/4 is a geometric factor that allows for a straightforward comparison to the work of these 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 2 p Ω/(4ν f ) = 9S that is often used in literature (e.g.Riley 1966Riley , 1967)).As stated in the Introduction, 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 ≪ 1 and S ≪ 1/9 the area over which the primary vortices propagate is much larger than D p /2 and (ii) for S ≫ 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 behavior of two particles, this classification is still useful for describing the obtained flow fields.
Moreover, we define a Stokes number St = τ p /τ f using the characteristic particle time scale τ p = |ρ s − 1|D 2 p /(18ν f ) and the characteristic flow time scale τ f = 1/Ω, which yields St = |ρ s − 1|2S.Since St is determined by the linear combination of the density ratio ρ s and the dimensionless frequency S, we focus in the following on the influence of these two non-dimensional parameters on the behavior of the particles by varying the initial gap size 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 timestep 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 timestep, spatial discretization and domain size.Note that preliminary tests have shown that doubling the domain size only has an impact on the order of O(10 −3 ) on the particle motion and is therefore negligible.

Comparison to 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 behavior of sediment beds under shear and particles settling under gravity in quiescent fluids (Biegert et al. 2017;Vowinckel et al. 2019Vowinckel et al. , 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 setup under investigation considers two particles in an oscillating domain filled with fluid, we extend our validation to setups with oscillating motion patterns.In this regard, we compare our computational results to analytical, experimental and numerical data of single-and two-particle setups 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 setup 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 an unidirectional, monochromatic oscillatory fluid flow.The particles are horizontally aligned, i.e. θ i = 0 • , and placed with a surface separation distance of ζ i = ζ t = 1.00 so that the center of the particle pair coincides with the center 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 100T , 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 ω z (x, y) and the coloring is according the presented color 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 et al. 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 toward 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.
With reference to the case distinctions for a single-particle setup 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 class (large S).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 representation 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 below.
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 setups 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.

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. ρ 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 center of the domain.The particle pair approaches the center of the domain asymptotically so that for a horizontally aligned arrangement (θ i = 0 • ) with ρ s = 4.68 and ζ i = 0.75, the artifact 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 center between the particles in relation to its initial position over time and figure 3(b) illustrates the evolution of the inter-particle center over time as moving average with an averaging window of one oscillation period.

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 analyze the near-field and far-field hydrodynamic interactions of two particles aligned with θ i = 0 • for various initial particle distances ζ 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 ρ s .
Since we neglect gravity and consider two particles with an initial orientation an- gle θ i = 0 • , the particle coordinates remain constant in y-and z-direction at L y /2 and L z /2, respectively.The initial positions in x-direction are modified to result in ζ 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 ζ 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) to become active for ζ 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 ǫ = 0.1 were applied for a total number of 100 oscillation periods.Since we have already seen that S has a direct correlation to the resulting behavior of the particles, we will investigate ζ more systematically.Here, ζ represents the change in the particle distance over time, computed as ζ = ζ t − ζ i .Therefore, we define ζ 100 as the change of the particle distance ζ 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 ζ i to better quantify the effects of this configuration parameter.In Figure 5, we present ζ 100 for the same values of S as in figure 4. In addition to the results for ζ i = 0.25 (figure 5a), we also show ζ 100 for ζ i = 3.00 in figure 5(b).Both figures reveal that the repulsive and attractive behavior becomes more pronounced for larger and smaller values of S, respectively.An exception can be seen for S = 0.07 with ζ 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 ζ t , i.e. lubrication forces.Such behavior indicates that there is a peak for small initial gap sizes where the attracting behavior is most efficient.For the case of ζ 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 ζ i .Comparing the two initial gap sizes ζ i = 0.25 and ζ i = 3.00 in figures 5(a) and (b), we can conclude that the particle-interaction becomes less intense for larger ζ i .Note the scale difference by one order of magnitude in figures 5(a) and (b).For ζ i = 3.00 in Figure 5(b), the damping of particle motion toward each other by lubrication completely vanishes for low frequencies and the threshold for S for the transition from attractive to repulsive behavior shifts to smaller values of S.
The analysis for particles of relative density ρ s = 4.68 presented so far suggests that there is not only an optimum for attractive particle interaction for ζ i ≈ 0.25, but also a threshold condition for which the interaction transitions from attractive to repulsive behavior.Both of these phenomena are a function of S and ζ i .To elucidate the effect of ζ i on the particle behavior further, we evaluate ζ 100 for a constant frequency (S = 2.09) and different values of ζ i in figure 6.Again, we use an initial orientation of θ i = 0 • and a relative particle density of ρ 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 ζ i .Again, the case of ζ i = 0.10 deviates from this behavior, because for such small gaps, lubrication forces become an important component of the particle interaction.Nevertheless, for ζ i > 0.1, the intensity of the mutual interaction decreases with increasing ζ 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 ζ i increases.
This behavior is exemplified by the comparison of the particle amplitudes of the twoparticle setups A 2p (S, ζ i ), in dependence of ζ i , and the setups 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 σ = A 2p (S, ζ i )/A 1p (S) for T = 11 in figure 7. The time for the calculation of σ 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 setups increases with decreasing ζ i .This is in line with the observation and assumption made above that for the given setup the interaction of two particles increases with decreasing ζ i .As an important consequence, this analysis also demonstrates that the closer the particles are, the more amplified their excursion lengths become.

Particle arrangements
In order to illustrate the influence of the arrangement orientation, we first look at the behavior of two monodisperse particles with ρ s = 4.68, whose initial distance is ζ i = 0.25 and the initial angle of the orientation of the particle alignment with respect to the direction of oscillation θ i is varied.We modify θ i in steps of 15 • in the range from 0 • to 90 • and analyze the change of θ and ζ over a period of 100 oscillations.Here, θ represents the alignment angle over time.For all setups, the particles are arranged so that the midpoint between the particles coincides with the center of the domain.For  The results illustrate that θ increases and converges toward a value of 90 • for particle configurations with θ i = 0 • , whereas θ = θ i for θ i = 0 • .The inter-particle distances ζ, on the other hand, develop significantly differently depending on S and θ i .In figure 8 (S = 0.35), the particles for initial arrangements of θ i 30 • approach each other, meaning ζ < 0, while particles move away from each other for θ i 60 • .Interestingly, the particles initially behave attractively and then repulsively for θ i = 45 • .For S = 2.09 (figure 9), similar phenomena occur, whereby in this case smaller angles, i.e. θ i 15 • , show a repulsive and larger angles, i.e. θ i 60 • , an attractive behavior.Here, both θ i = 30 • and 45 • show a reversal, both from repulsive to attractive.An indication of the reversal can already be seen at θ i = 15 • , while the behaviour in the considered time frame is still repulsive.The analysis presented above shows the following.On the one hand, our data confirms the observations by Lyubimov et al. (2001), Voth et al. (2002), Wunenburger et al. (2002), and Klotsa et al. (2007Klotsa et al. ( , 2009) ) that two particles exposed to oscillatory flow align themselves perpendicularly to the direction of oscillation for any perturbation, e.g. in terms of θ i , that breaks the symmetry of the oscillation direction.The well developed stage for the alignment angle θ = 90 • 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 behavior 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 will therefore focus on the specific configuration with an alignment angle θ = 0 • in order to investigate the parameter ranges for which particles can be brought into contact or separated from each other.

Interaction behavior
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 to determine these threshold conditions for a wide range of the parameters S, ζ i and ρ s .To this end, we conducted a set of simulations with ǫ = 0.1 for S = [0.07,2.10] with increments of ∆S = 0.035, ζ i = {0.10,0.25, 0.50, 0.75, 1.00, 1.50, 2.00, 3.00} for the three density ratios considered, ρ s = {0.47,1.78, 4.68}.Recall that these values of ρ 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 (⊲) indicate that the particles are approaching each other, i.e. ζ 100 < 0, while left pointing triangles (⊳) refer to setups where the particles are drifting apart, i.e. ζ 100 > 0. We furthermore categorize the mutual interaction as attractive (◮), transitional (open symbols, ⊳⊲), and repulsive (◭), depending on S and ζ i .This categorization is derived from the analysis of A 2p (S, ζ i ) discussed in §3.2.We define arrangements, in which |ζ 100 | A 2p (S, ζ i ) as transitional, whereas setups with exceeding |ζ 100 | > A 2p (S, ζ i ) are classified as either attractive or repulsive, depending on the sign of ζ 100 .The setups 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 toward 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.
Despite the low values of ζ 100 that are assigned to transitional behavior, a more precise distinction into transitionally attractive and transitionally repulsive remained possible as will be shown in §4.5 below.This is also expressed in the clear definition of the threshold condition shown as a dashed line in figure 10.
For ρ s = 4.68 in figure 10(a), our classification scheme reveals that for all arrangements of ζ i 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 setups classified as transitional.With increasing ζ 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 ζ i = 3.00, no significant attractive or repulsive behavior 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 ρ 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. ρ s = 1.However, comparing them to the case ρ 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 setups lead to transitional results for those two density ratios, where the number of repulsive cases are also less than for ρ s = 4.68.For example, the arrangements with larger initial gaps (ζ i 2.00) do not lead to any significant repulsive behavior, 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 ζ i are increased.
As mentioned above, the dashed line in each regime map marks the transition from attractive to repulsive behavior as indicated by a change of sign in ζ 100 shown in figure 4. For figure 10, this threshold condition was determined by best fit of a power law given by ζ = C 1 S β , for which we obtain an excellent correlation with our data.For the highest particle density ratio ρ s = 4.81 (figure 10a), we determine C 1 = 0.18 and β = −2.73.For the two cases ρ s = 1.78 and ρ s = 0.47, we obtain identical fitting parameters C 1 = 0.15 and β = −2.93.This illustrates two things: (i) an increase of inertia does not substantially change the threshold conditions to transition from attractive to repulsive behavior 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 behavior with respect to ρ 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 (ζ i = [0, 2.0]) and due to the immobilization of the particles, the parameter ρ s did not play a role.In principle, our analyses agree in that small S and small

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 nonzero 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 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 its original arrangement.Second, the effect of the initial displacement of the particles in a fluid at rest (viz.§3.1) must have already decayed and thus be negligible.Consequently, a compromise must be made which best fulfills both requirements.To this end, we have performed a convergence study in which we analyzed 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 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 setup 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 in the upper panels) and repulsive behavior (S = 1.05 in the lower panels) using ζ i = 0.75 and ρ s = 4.68.For this purpose, we show streamlines in the xy-and zy-planes.The planes cut through the center 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 11(b) and 11(d) to indicate the particle positions aligned along the xaxis.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 ω 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 color scheme corresponds to the color scale in figure 2 and is applied for all subsequent analyses.
The comparison of the xy-planes for the repulsive and the attractive case (figures 11a and 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 setup.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 (c.f. 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 setup 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 toward 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 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 toward nor away from the center.The inner part of the circle indicates an convergent inflow, where the streamlines point towards the center between the particles.For the part outside of the stagnation line, the flow is directed outwards away from the particles.For the repulsive behavior (figure 11d), the entire flow of the zy-plane is directed towards the center of the particles.The vorticity in zy-plane is not discernible in neither 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 shown in figures 11 show a strong axisymmetric structure.This is true for both, attraction as well as for repulsion.We therefore conclude that the xz-plane (not shown here) and the xy-plane, both arranged through the center 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 fluidparticle interaction and the resulting particle behavior.

Flow decomposition
The detailed illustrations shown in figure 11 already indicate that there are significant differences in the flow structures of attractive and repulsive setups.However, identifying and quantifying a clear physical mechanism that allows for conclusions regarding attractive or repulsive behavior of the particles is 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 ω z (x, y), hereafter also referred to as total steady streaming, allows to decompose the flow field into a symmetric and an antisymmetric part denoted as ω z,S (x, y) and ω z,A (x, y), respectively, The symmetric part is generated by calculating the symmetry about the center 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 Color scheme as in figure 2.
x c is the x-coordinate at the center between the two particles.Therefore, ω z (x ′ , y) can be written as In figure 12, the streamlines and vorticity contours of the symmetric and antisymmetric components are shown in the panels of the left and the right column, respectively.The setups are the same as given in figure 11.In this regard, figures 11(a) and (c) represent the total steady streaming used for the decomposition (4.1).When considering the streamlines of the symmetric part (left two panels of figure 12), two vortex structures with center points in the far field orthogonal to the particle arrangement and in line with the center of the particle distance can be seen for the attractive as well as repulsive behavior.In both cases, the upper vortex rotates in counterclockwise and the lower one in 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 and 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 higher deviation from the center of the domain at the time of calculating the steady streaming.This results in an stronger motion towards the center 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 behavior regarding attraction or repulsion, it is responsible for the flow field of the symmetric component of the steady streaming.
The antisymmetric fields (figures 12b and 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 behavior 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 behavior 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.

Effect of particle inertia
Since the previous figures have shown results for ρ s = 4.68 only, we now turn our attention to the steady streaming observed for two additionally considered particle density ratios ρ s = 1.78 and 0.47 that were already analyzed in terms of their particle interaction in figure 10.We therefore performed the same analysis of flow decomposition as for ρ s = 4.68 given in the previous section.Figure 13 shows the corresponding results of steady streaming for S = 0.35 and ζ i = 0.75 for these two density ratios.In this way, we can compare the flow structures and vortex contours of these two setups (ρ s = 1.78 in figures 13a, c & e and ρ s = 0.47 in figures 13b, d & f) with the previous one (ρ s = 4.68 in figures 11a and 12a,b) and can thus infer the behavior of two particles in oscillating flow that are both denser and lighter than the surrounding fluid.Since we neglect gravity, the particle dynamics are affected by their mass only in terms of their inertia, but not their weight.
According to the regime maps shown in figure 10, for S = 0.35 and ζ i = 0.75, all ρ s result in an attractive behavior of the particles.Only ρ s = 4.68 is assigned to the category of attraction, while ρ s = 1.78 and 0.47 are associated with transitional behavior according to (3.1).The comparison of the total steady streaming (figures 11a, 13a and 13b) indicates the same formations of vorticity contours with respect to the clockwise and counterclockwise rotation directions.For each ρ s , each respective particle is surrounded by four adjacent and four far-field contours.The intensity of the vorticity decreases with decreasing |ρ 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 ρ s = 4.68 and 1.78 reveal very similar patterns of the respective vortex structures.The vortex structures of ρ s = 0.47 show the same patterns as for ρ 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 ρ 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. ρ 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 ρ s > 1 and ρ 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 noninertial 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 center 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 ρ s so that their contours do not show up in figures 12(a), 13(c) and 13(d).
Looking at the antisymmetric component shown in figures 12(b), 13(e) and 13(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, 13a and 13b), but here the mirrorinverted arrangement was filtered out by the symmetric component.This suggests that even though ρ s = 1.78 and 0.47 are classified as transitional, because they are only marginally attractive, and ρ 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).

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 fluidparticle 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 ρ 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 analyze 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 ζ i = 0.50 and ρ s = 4.68.The resulting streamlines and vorticity contours are depicted in figure 14, with S = 0.35  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 case of attraction, they appear to be in equilibrium, when the transition is established, and become clearly superior in 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 the peripheral decrease.As already indicated by the regime map given by figure 10, the behavior 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 ζ i on the flow characteristics in figure 15 by comparing ζ i = {0.10,0.25, 0.50} for S = 0.94 and ρ s = 4.68.Here, ζ i = 0.10 results in attraction (figure 15a), ζ i = 0.25 leads to a transitional state (figure 15b), and ζ 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 ζ i .This is accompanied by the fact that with increasing ζ i , more fluid can flow into the gap, increasing the strength of the flow into the gap trying to push the particles apart.

A circulation-based criterion to determine the particle interaction
The qualitative examination of the antisymmetric component of the steady streaming shows that the visualization of the flow structures can be helpful in determining the behavior 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 analyzing different density ratios, where we find the same flow patterns for strongly and weakly attractive behavior 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 behavior 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 numbers.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, ω z,α , into account.Here, α is used as an index for either attraction or repulsion.Where, 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 Γ α contributing to attraction or repulsion, respectively, as On this basis, we can compute the ratio of attractive to repulsive circulation φ Γ = Γ attr /Γ rep and, thus, draw conclusions about the behavior of two particles in oscillating fluid flow.A ratio of φ Γ = 1 would represent a balance of the central and peripheral vortices, φ Γ > 1 indicates a predominantly attractive circulation leading to the approach of the two particles, and φ Γ < 1 refers to repulsive behavior due to dominant circulation pushing the particles apart.
Figure 16 presents φ Γ for the considered density ratios ρ s .The same set of S is used as in figures 4. The dashed lines highlight the balance of the circulations in each graph, i.e.Γ α = 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 behavior of the two particles in terms of the qualitative distinction of either attraction or repulsion is small.Consequently, the ratios for the three densities ρ 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 ζ i lead to attractive results is reflected by the fact that these frequencies lead to φ Γ > 1, whereas we find φ Γ < 1 for repulsive setups for large S and large ζ i .
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 behavior of repulsion and attraction.Indeed, the finding in §3.2 that the interaction of the flow fields of the respective particles increases with decreasing ζ i is clearly evident by the increase in φ Γ with decreasing ζ i .For larger ζ i , φ Γ 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 behavior with regard to attraction or repulsion.Thus, this is solely determined by the non-dimensional oscillation frequency S and the initial particle distance ζ i .

Conclusions
To investigate the potential hydrodynamic mechanism that induces flocculation of primary particles in oscillatory flows, a systematic simulation campaign was conducted to analyze the mutual particle behavior 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 setup 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 less 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 behavior.To this end, we analyzed 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 behavior 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 behavior was then further analyzed 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 analyzing the plane through the torroidal 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 behavior.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 behavior 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.
Funding FK and BV gratefully acknowledge support through German Research Foundation (DFG) grant VO2413/2-1.FK also gratefully acknowledges the support of a scholarship from the German-American Fulbright Commission.EM and PLF were supported through NSF grant CBET-1638156.EM furthermore acknowledges support through U.S. Army ERDC grant W912HZ22C0037, U.S. ARO grant W911NF-23-2-0046, and through NSF grant HS EAR 2100691.values of ρ s .Note that (A 1) is only valid for the inertial reference frame.We, therefore, transformed the results by η = 1 − A ′ p /A f , to subtract the amplitude of the domain oscillation.Here, A ′ p denotes the particle motion in a non-inertial reference frame, which consequently represents the relative deviation of the particle motion from the fluid motion.The comparison demonstrates excellent agreement of our results and the experimental as well as the analytical outcomes.Hence, our method is well suited to analyze the fluid-particle interaction for more complex scenarios in oscillatory, low-Reynolds number flows.

A.2. Flow characteristics
In addition to the comparison of particle dynamics presented above, we analyze the flow characteristics of a single particle oscillating in a fluid at rest.This is done both qualitatively (by examining the flow patterns in terms of shapes and flow directions), and quantitatively (by analyzing the distance from the stagnation point of the flow perpendicular to the axis of oscillation to the particle surface).Such detailed analysis is necessary to ensure that the fluid-particle interaction is not affected by numerical effects.For this purpose, we compare the streamlines emerging due to steady streaming, which represents the averaging of the flow field over one oscillation period and which has been thoroughly described in §4.1.
This comparison is based on the work of Kotas et al. (2007) and Chang & Maxey (1994).Kotas et al. (2007) experimentally investigated a single sphere oscillating in an otherwise quiescent fluid.To this end, the sphere was mounted on a steel rod and subjected to oscillatory motion while immersed in a container of viscous fluid.The resulting steady streaming flow field around the sphere was visualized by averaging the records of phase-locked particle pathline images over multiple oscillation periods.The experiment run of Kotas et al. (2007) used for comparison is characterized by the nondimensional numbers S = 9.33, Re = 16.79 and Re s = 0.84.Using the same parameters, Chang & Maxey (1994) conducted DNS to investigate the steady streaming generated by a sphere oscillating with low amplitude.
We chose the same characteristics of the setup as in §2.1 and prescribed a particle velocity to match the non-dimensional parameters given above.Figure 18a shows our simulation results for which we calculated the steady streaming after 100 oscillation periods.The figure represents the part above the symmetry axis, which is simultaneously the oscillation axis of the particle.Two vortex structures evolve close to the particle surface, replicating its spherical structure, followed by two more distant vortex structures.The near-surface structures, often referred to as shear-wave or Stokes layer, rotate counterclockwise in the left and clockwise in the right vortex, respectively.The respective adjacent structures further away rotate in exactly the opposite direction.In consideration of the case distinctions according to Riley (1966), shown in §2.2, the present case can be assigned to case two, which in turn is indicated by the characteristic flow numbers and the illustrated flow patterns consisting of primary and secondary vortices.Kotas et al. (2007) and Chang & Maxey (1994) focused mainly on the shape, size, and rotational direction of the primary vortices, which are very well matched by our simulations.No further information is given about the secondary structures except for the directions of rotation.This is probably due to the fact that the extent and shapes of these structures strongly depend on the domain size and form.A sketch of the close-up of the general streaming structures is given in figure 18b.The shift of the further distant flow structures to the left in figure 18a is related to the transients caused by initial motion of the particle, which is explained in §3.1.
The cross above the particle in figure 18b indicates the stagnation point, which represents the extent of the inner streaming region close to the particle surface.The distance of the stagnation point from the particle surface δ SP was measured experimentally by Kotas et al. (2007) and calculated numerically by Klotsa (2009) in dependence of the oscillation properties.The results of δ SP are shown in figure 19, where we compare our numerical simulations ( * ) with the experimental data of Kotas et al. (2007) ( ) and the numerical outcomes of Klotsa (2009) (•) for a variety of S. The solid and dashed lines represent the best fits of Kotas et al. (2007) and Klotsa (2009), respectively.Both, Kotas et al. (2007) and Klotsa (2009) point out that numerical and experimental results may differ due to transient effects.Such effects are primarily caused by the smaller number of oscillation periods used for numerical simulations to calculate δ SP than for experiments.For example, Kotas et al. (2007) states that the extent of the inner flow region is 10% larger after 10T than after 100T .In this regard, Klotsa (2009) argues that the deviation of their results is caused by the fact that their simulations were limited to 50T .In our case, we analyzed the flow properties after 100T and consequently obtain a better agreement with the experimental data.However, we also want to note that the spatial discretization has an effect on the magnitude of δ SP .Nevertheless, since a higher spatial resolution is associated with a considerably increase of computational cost, we consider our results sufficient to support the accuracy of the computation of the flow characteristics.
We conclude that our results are in close agreement with the experimental and numerical comparison data with respect to the shape and size of the vortices as well as to the size of the inner streaming region for an individual particle oscillating in a quiescent fluid.

Figure 2 :
Figure 2: Flow characteristics illustrated by streamlines and vorticity contours for steady streaming of two stationary and axially arranged particles in oscillatory flow.(a) represents the numerical results of the setup with S = 0.89, Re = 0.32 and Re s = 0.0032 and (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 coloring ranges from −0.01 (blue) for clockwise rotations to 0.01 (red) for counterclockwise rotations.

Figure 3 :
Figure 3: Location of the midpoint between the particles relative to its initial location over time for (a) the instantaneous location ∆x c and (b) the moving average < ∆x c > with an averaging window equal to the oscillation period T .
Figure 4 illustrates the instantaneous evolution of ζ over time for various choices of S for an representative setup with initial particle distance ζ i = 0.25 and ρ s = 4.68.The horizontal dashed line represents the respective state of equilibrium, in which ζ 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 ζ despite identical initial conditions.In general, results with increasing values of ζ > 0 over time show that the

Figure 4 :
Figure 4: Evolution of instantaneous ζ over 100 oscillation periods for ρ s = 4.68, ζ i = 0.25 and a variety of S.

Figure 5 :
Figure 5: ζ 100 as a function of S with ρ s = 4.68 for (a) ζ i = 0.25 and (b) ζ i = 3.00.The sketches in the subfigures 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 figures.ζ 100 < 0 represents attraction and ζ 100 > 0 repulsion.

Figure 6 :
Figure 6: ζ 100 as a function of ζ i with ρ s = 4.68 and for S = 2.09.

Figure 7 :
Figure 7: Amplitude ratio σ = A 2p (S, ζ i )/A 1p (S) between the amplitudes of the two-particle setups A 2p (S, ζ i ) and the respective individual-particle setups A 1p (S) for ρ s = 4.68 and as a function of S and ζ i .

Figure 9 :
Figure 9: Evolution of (a) θ and (b) ζ over 100 oscillation periods for ρ s = 4.68, ζ i = 0.25 and S = 2.09.Note the difference in scale for the y-axis in (b) compared to figure 8.
ζ i generally lead to a decrease of ζ, while large S and large ζ i result in an increase of ζ over time.In fact, one can convert our (S, ζ i )-regime map to the regime map of Fabre et al. (2017) (figure 6(b) in that reference) and we find very close agreement for the transition for S = 0.56 and S = 0.89 (termed Ω = 5 and Ω = 8 in that same reference).

Figure 10 :
Figure 10: Regime maps of the fluid-particle interactions as a function of S and ζ i in log-log plots for different particle density ratios: (a) ρ s = 4.68, (b) 1.78 and (c) 0.47.Right pointing triangles (⊲) indicate ζ 100 < 0 and left pointing triangles (⊳) ζ 100 > 0. The coloring of the symbols is according to the classification stated in (3.1),where ◮ represents attraction, ◭ repulsion and open symbols (⊳⊲) 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.

Figure 11 :
Figure 11: Streamlines and vorticity contours of the steady streaming in xy-(figures a, c) and zy-planes (b, d).The top row shows an attractive case with S = 0.35 and the bottom row a repulsive behavior with S = 1.05.Both are arrangements with ζ i = 0.75 and ρ 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.Color scheme as in figure 2.

Figure 12 :
Figure 12: Streamlines and vorticity contours of the symmetric (a, c) and antiymmetric components (b, d) for ζ i = 0.75 and ρ s = 4.68.The figures in the top row, (a, b), represent S = 0.35 resulting in attraction, and in the bottom row, (c, d), S = 1.05 leading to repulsion.The total steady streaming used for decomposition is given in figures 11(a, c).Color scheme as in figure 2.

Figure 14 :
Figure 14: Streamlines and vorticity contours of the antisimmetric part of the flow field for a constant ζ i = 0.50 and varying S.Figure (a) represents S = 0.35 with attractive behavior, figure (b) shows S = 0.70 (transition), and figure (c) illustrates S = 1.05 (repulsion).Color scheme as in figure 2.
Figure 14: Streamlines and vorticity contours of the antisimmetric part of the flow field for a constant ζ i = 0.50 and varying S.Figure (a) represents S = 0.35 with attractive behavior, figure (b) shows S = 0.70 (transition), and figure (c) illustrates S = 1.05 (repulsion).Color scheme as in figure 2.

Figure 15 :
Figure 15: Streamlines and vorticity contours of the antisimmetric part of the fluid field for constant S = 0.94 and varying ζ i .(a) ζ i = 0.10 resulting in attraction, (b) ζ i = 0.25 leading to transition, and (c) ζ i = 0.50 causing repulsion.Color scheme as in figure 2. The dashed lines in (c) indicate the symmetry lines used for the division into four quadrants to compute ω z,α (4.3).

Figure 17 :
Figure 17: Comparison of the numerical simulation results of the amplitude ratio η for various dimensionless frequencies S with the analytical predictions and experimental results of L Espérance et al. (2005) for a constant amplitude A f = 0.1 and three different density ratios ρ s = 0.47, 1.78, and 4.68.

Figure 18 :
Figure 18: Flow characteristics illustrated by streamlines of the steady streaming of a single particle oscillating horizontally (black arrows) in a fluid at rest.(a) represents the numerical results of the setup with S = 9.33, Re = 16.79 and Re s = 0.84 and (b) a simplified sketch of a close-up to highlight the flow characteristics as well as the stagnation point indicated by the cross.