1. Introduction
The settling dynamics of irregularly shaped particles and particle aggregates is of considerable interest for a variety of environmental transport processes and engineering applications. Examples include the settling of mud, clay and silt (Winterwerp Reference Winterwerp2002; Strom & Keyvani Reference Strom and Keyvani2011; Te Slaa et al. Reference Te Slaa, van Maren, He and Winterwerp2015), as well as marine snow (Alldredge & Gotschalk Reference Alldredge and Gotschalk1988; Diercks & Asper Reference Diercks and Asper1997) and microplastics (Khatmullina & Isachenko Reference Khatmullina and Isachenko2017; Wang et al. Reference Wang, Dou, Ren, Sun, Jia and Zhou2021; Yan et al. Reference Yan, Wang, Dai, Sun and Liu2021). Similar issues arise with regard to the dynamics of ice particles in clouds (Gustavsson et al. Reference Gustavsson, Jucha, Naso, Lévêque, Pumir and Mehlig2017) and in the context of industrial processes involving particle flocculation (Licskó Reference Licskó1997; Kurniawan et al. Reference Kurniawan, Chan, Lo and Babel2006). Environmental systems and engineering applications such as deep-sea mining (Meiburg & Kneller Reference Meiburg and Kneller2010; Peacock, Alford & Stevens Reference Peacock, Alford and Stevens2018; Gillard et al. Reference Gillard, Purkiani, Chatzievangelou, Vink, Iversen and Thomsen2019; Ouillon et al. Reference Ouillon, Kakoutas, Meiburg and Peacock2021; Wells & Dorrell Reference Wells and Dorrell2021) may also be strongly affected by where particles are transported and deposited, so that it is desirable to have accurate predictions of their settling and drift velocities as functions of their geometry.
While the settling of individual spheres across different flow regimes is a well-studied topic that has received much attention (Leal Reference Leal1980; Dietrich Reference Dietrich1982; Brown & Lawler Reference Brown and Lawler2003; Yang et al. Reference Yang, Fan, Liu and Dong2015), more recent studies have addressed the settling dynamics of irregularly shaped objects in inertial flows, such as oblate spheroids (Moriche, Uhlmann & Dušek Reference Moriche, Uhlmann and Dušek2021), cubes and tetrahedra (Rahmani & Wachs Reference Rahmani and Wachs2014), slender bodies (Khayat & Cox Reference Khayat and Cox1989) and disks (Heisinger, Newton & Kanso Reference Heisinger, Newton and Kanso2014). These studies show that the spatial orientation and settling velocity of a body are highly dependent on its geometry. Changes in the size and shape of an aggregate can greatly affect its drag (Loth Reference Loth2008; Li et al. Reference Li, Yu, Gao and Flemming2020), leading to large changes in the settling rate.
An additional key feature of aggregate shape concerns the effective porosity and permeability. The formation of particle aggregates can be triggered by a variety of forces, such as van der Waals forces (Visser Reference Visser1989) or biocohesive forces caused by the bonds between organic molecules (Malarkey et al. Reference Malarkey2015). Recent numerical works have taken first steps towards simulating such cohesive forces, as analysed in the investigation of aggregates with brittle tensile bonds by Langlois, Quiquerez & Allemand (Reference Langlois, Quiquerez and Allemand2015), or in the study of the settling dynamics of cohesive sediment held together by van der Waals forces (Vowinckel et al. Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019). These aggregates, in turn, due to their porosity, allow for fluid to pass through the pore spaces, which influences the settling behaviour (Prairie et al. Reference Prairie, Ziervogel, Camassa, McLaughlin, White, Dewald and Arnosti2015). The size and distribution of these pore spaces are controlled by the geometry of the aggregate, such that the arrangement of particles leads to aggregates with highly variable settling velocities.
In contrast to bodies with fore–aft asymmetry, for non-symmetric bodies the settling motion can become significantly more complex, causing them to settle along complicated trajectories (Johnson, Li & Logan Reference Johnson, Li and Logan1996; Tang, Greenwood & Raper Reference Tang, Greenwood and Raper2002). In the present work, we consider the sedimentation of rigid bodies consisting of two connected spheres of unequal diameters. Except for a simplified model variant (Candelier & Mehlig Reference Candelier and Mehlig2016), the dynamics of this case has not yet received much attention, while the similar situation of particle pairs with equal diameters but unequal densities has been examined in some depth by Nissanka, Ma & Burton (Reference Nissanka, Ma and Burton2023). The existing literature on the topic primarily focuses of the case of highly viscous flows at or just above the Stokes limit. As we will see, the asymmetric particle pair represents a relatively simple case that can exhibit a complicated settling dynamics depending on the ratio of the diameters.
The present study focuses on the dynamics of elementary aggregates consisting of two spherical particles of different sizes. We employ both particle-resolved simulations as well as laboratory experiments, in order to investigate the dynamics of such aggregates settling in a fluid at rest. The problem configuration is described in detail in § 2, including the experimental and numerical approaches, as well as validation results. Section 3 presents comparisons of experimental and simulation results, and it explores the influence of the ratios of the diameters and also of the gravitational and viscous forces in depth. Scaling laws are obtained for the dynamics of the aggregate as a function of these dimensionless parameters, and the accompanying fluid flow field is analysed. Section 4 summarizes the main findings of the study, and discusses potential further extensions.
2. Methods
2.1. Problem definition
We aim to characterize the dynamics of a pair of rigidly connected spherical particles with density 
$\rho _{p}$ that settles in a fluid of density 
$\rho _{f}<\rho _{p}$ and kinematic viscosity 
$\nu$, in a large tank of square cross-section with width 
$W$ and height 
$H$, as illustrated in figure 1. The large and small particles have diameters of 
$D_{L}$ and 
$D_{S}$, respectively, with 
$D_{L} \geqslant D_{S}$, so that the geometry of the aggregate is characterized by the diameter size ratio
\begin{equation} \alpha = D_{L}/D_{S} . \end{equation}
Figure 1. Schematics (a) of the pair of particles considered at 
$t = 0$ (top) and at a later time (bottom), (b) of the experimental set-up and (c) of example aggregates prepared in the laboratory for various values of 
$\alpha$. Straight orange dashed lines in (a) indicate the lines forming the orientation angle 
$\theta$, and the curved orange dashed lines indicate the elliptical shape employed to extract the orientation angle from the experimental data during image processing.
At the initial time, 
$t=0$, the particle centres of gravity are aligned horizontally in the 
$z$-direction, and both the fluid and the particles are at rest, as shown in figure 1(a). Gravity points in the 
$y$-direction. Due to the symmetry of the particle pair and the moderate influence of fluid inertia considered in the present study, the motion of the aggregate in the 
$x$-direction is negligible compared with that in the ( 
$y,z$)-plane. We hence describe the motion of the particle pair in terms of the vertical settling and the lateral drift velocities of its centre of mass, 
$u_y$ and 
$u_z$, respectively, as well as its orientation angle 
$\theta$, defined as the angle between the upward vertical axis and the line connecting the spheres’ centres.
2.2. Experimental methods
The experiments were conducted in a tall tank of height 
$H = 59.4$ cm and square cross-section of dimensions 
$W \times W=14.1 \times 14.1\ {\rm cm}^{2}$. The tank was filled with one of two types of mineral oils: STE Food Grade 200 (dynamic viscosity 
$\mu _{f,1} = 96\ {\rm mPa} \ {\rm s}$ and density 
$\rho _{f,2} = 864 \ {\rm kg} \ {\rm m}^{-3}$ at 
$20\,^\circ {\rm C}$), or STE Food Grade 70 (dynamic viscosity 
$\mu _{f,2} = 23\ {\rm mPa} \ {\rm s}$ and density 
$\rho _{f,1} = 850\ {\rm kg } {\rm m}^{{-3}}$ at 
$20\,^\circ {\rm C}$).
Aggregates made of two spheres were released in a horizontally aligned configuration (see figure 1(a), top). The aggregates were manufactured by gluing together two nylon spheres of density 
$\rho _{p} = 1135\ {\rm kg}\ {\rm m}^{-3}$ (figure 1c). The larger particle diameters were 
$D_{L} = 4.76$, 6.35, 9.53 and 
$11.11\ {\rm mm}$, while the smaller particle size ranged from 
$D_{S} = 1.59 \ {\rm mm}$ to 
$11.11 \ {\rm mm}$. This leads to values of 
$\alpha = D_{L}/D_{S}$ in the range 
$1 \leqslant \alpha \leqslant 7$, with the majority of the aggregates having 
$\alpha \leqslant 4$, which is the range in which the change in dynamics is most pronounced. For the sake of experimental simplicity, some of the experiments at lower viscosity were performed in a smaller tank, without changing quantitatively the results. The parameter ranges for both the experiments and the numerical simulations are indicated in table 1.
Table 1. Key parameter ranges considered in the experiments and numerical simulations performed in the present study, for the diameter aspect ratio 
$\alpha$, density ratio 
$\rho '$, Galileo number 
$Ga$, dimensional domain height 
$H$ and width 
$W$, and initial orientation angle 
$\theta (t=0)$. For the numerical simulations corresponding to the experimental cases, we define 
$\theta _1$ to be the first recorded value of the orientation angle obtained via measurement.
The aggregate was first immersed in a beaker to be coated with the same oil as in the bath to prevent bubble entrainment during its later immersion into the tank. Then, the aggregate was gently submerged in the tank approximately 1 cm below the fluid surface with an adjustable metallic wrench. The wrench rested on a flat horizontal support placed on the top of the tank, and was opened manually to release the aggregate. The experimental set-up was backlit with a LED panel and a Nikon camera recorded the settling of the aggregate in the fluid at 60 frames per second (see video examples in supplementary materials available at https://doi.org/10.1017/jfm.2024.711), so that the instantaneous velocity and angle could be extracted. Three copies of each aggregate geometry (
$D_{L}$, 
$\alpha$) were built from three distinct pairs of spheres, and each aggregate's settling motion was recorded to later average the results over the effective variability in sizes and densities of the nylon spheres.
In most cases, the aggregates reached a quasi-steady state at the end of an experiment, i.e. the aggregate velocity and orientation angle no longer varied with time. Only a few cases did not reach a steady state, due to the limited size of the experimental domain. The value of 
$\theta$ was determined by aligning an elliptical shape onto the images of the settling pair, and calculating the angle formed by its major axis with the vertical direction (see the curved dashed lines in figure 1a).
A few experimental cases, where the aggregate rotated significantly out of the ( 
$y,z$)-plane, were detected via tracking the apparent aggregate shape, and subsequently discarded. The influence of fluid inertia in the present work is not large enough to expect out-of-plane rotation, so that we attribute any rotation in this direction in the experiments to inadvertent angular velocity generated during aggregate release. The numerical simulations support this approach by showing negligible out-of-plane rotation for all cases.
For detailed comparisons with the numerical simulations, a series of settling experiments of single spheres (the same constituting the aggregates) were performed in the same fluids and tanks as the aggregate experiments. The viscosity values at the average room temperature (
$20\,^{\circ }{\rm C}$) for the two fluids, 
$\mu _{f,1}(T = 20\,^{\circ }{\rm C})$ and 
$\mu _{f,2}(T = 20\,^{\circ }{\rm C} )$, were obtained by comparing the experimental terminal settling velocities, 
$u_{term}$, with a classical drag coefficient law (Schiller Reference Schiller1933)
\begin{equation} C_D = \frac{8 F_{drag}}{\rho_{f} \, u_{term}^2 {\rm \pi}D_{L}^2} = \frac{8 \left(m_{p}-m_{f}\right) g}{\rho_{f} \, u_{term}^2 {\rm \pi}D_{L}^2} = \frac{24}{Re}\left(1 + 0.15Re^{0.687}\right),\end{equation}where the Reynolds number is defined as
\begin{equation} Re = \frac{D_{L} \, u_{term}}{\nu}, \end{equation}
with the single (or large) particle diameter 
$D_{L}$. Here, 
$F_{drag}$ denotes the drag force on the spherical particle, g the gravitational acceleration, 
$m_{p}$ its mass, 
$m_{f}$ the mass of fluid in a volume the size of the particle and 
$\nu = \mu _{f}/\rho _{f}$ is the kinematic viscosity.
In order to build a sufficiently large tank without being able to determine the time required to reach a terminal (steady-state) velocity, we took the following approach. To decide on the tank dimensions, oil viscosities and particle sizes, we chose a target range of Reynolds numbers 
$Re \in [1, 100]$, and tank dimensions such that the time to settle over the distance of the tank height 
$H$ is greater than 
$20 \times D_{max}/u_{term}(D_{max})$ and the tank width is such that 
$W> 10 \times D_{max}$, with 
$D_{max} = 11.1\ {\rm mm}$. We determined the corresponding range of 
$u_{term}(D_{max})$ using (2.2) and (2.3). The experimental data presented here are those for which 
$H$ was large enough so that the settling velocity reached its steady state.
To account for slight temperature variations in the room (between 
$17\,^{\circ }{\rm C}$ and 
$23\,^{\circ }{\rm C}$), we measured the temperature dependence of the oil viscosities with a plate–plate geometry on an Anton Paar MCR 302 rheometer. These measurements provided 
$\mu _{{f},i,{rheom}}(T)$ for each of the two oils. We then fitted the measured data (provided in supplementary materials) by the curves 
$\mu _{{f},i}(T) = c_{i} \mu _{{f}, i,{rheom}}(T)$, with 
$c_{\rm 1} = 1.15$ and 
$c_{2} = 1.1$, so that the rheometer measurements agree with our best determination of absolute viscosity values from the settling sphere experiments, obeying (2.2) at 
$20\,^{\circ }{\rm C}$. To summarize, we base the absolute viscosity value on the comparison with (2.2), and the relative variation of the viscosity with temperature on the rheometer measurements.
2.3. Governing equations
The fluid flow is governed by the incompressible continuity equation
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}} = 0, \end{equation}and the unsteady Navier–Stokes equation
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\boldsymbol{u} ={-}\frac{1}{\rho_{f}}\boldsymbol{\nabla} p +\nu\nabla^2\boldsymbol{u} +\frac{1}{\rho_{f}}\boldsymbol{f}_{IBM} , \end{equation}
where 
$\boldsymbol {u}$ represents the fluid velocity, 
$t$ the time, 
$p$ the pressure with the hydrostatic component subtracted out and 
$\boldsymbol {f}_{IBM}$ the distributed force exerted on the fluid by the particles, which is calculated based on an immersed boundary method (IBM) approach, as will be discussed below.
The motion of the finite-size solid particles is governed by
\begin{gather} m_i\frac{{\rm d}\boldsymbol{u}_i}{{\rm d}t} = \oint_{\varGamma_i}\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n}\, {\rm d}A+V_i\left(\rho_{{p}, i}-\rho_{f}\right)\boldsymbol{g}+\boldsymbol{F}_{{b},i} , \end{gather}
\begin{gather}I_i\frac{{\rm d}\boldsymbol{\boldsymbol{\omega}}_i}{{\rm d}t} = \oint_{\varGamma_i}\boldsymbol{r}\times \left(\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n}\right){\rm d}A+\boldsymbol{M}_{{b},i} , \end{gather}
where 
$m_i$ indicates the mass of the 
$i$th particle, 
$\boldsymbol {u}_i$ and 
${\boldsymbol {\omega }}_i$ its translational and angular velocities, respectively, 
$V_i$ its volume, 
$\varGamma _i$ its surface, 
$\rho _{{p}, i}$ its density and 
$I_i$ its moment of inertia. Here, 
$\boldsymbol {F}_{{b},i}$ and 
$\boldsymbol {M}_{{b},i}$ represent the sum of all forces and moments acting on particle 
$i$ via rigid bonds with other particles, 
$\boldsymbol {g}$ denotes the gravity vector and 
$\boldsymbol {\tau }$ indicates the hydrodynamic stress tensor.
We also define 
$\boldsymbol {n}$ to be the outward normal vector on 
$\varGamma _i$, and 
$\boldsymbol {r} = \boldsymbol {x}-\boldsymbol {x}_i$ to be the position vector from the particle centre 
$\boldsymbol {x}_i$ to the surface point 
$\boldsymbol {x}$. Unless otherwise stated, we restrict our study to the case of two particles with varying diameters and identical densities, so that 
$\rho _{{p}, i} = \rho _{p}$.
2.4. Non-dimensionalization
We define characteristic length, velocity, time and pressure scales as
\begin{equation} x_{ref} = D_{L} , \quad u_{ref} = \sqrt{g' \,D_{L}} , \quad t_{ref} = \sqrt{D_{L} / g'} , \quad p_{ref} = \rho_{f} \, u_{ref}^2, \end{equation}
where 
$g'=g (\rho ' - 1)$ is the reduced gravity, with 
$\rho '=\rho _{p}/\rho _{f}$ denoting the ratio of particle-to-fluid density. We remark that we choose the diameter of the larger sphere as our characteristic length scale, in order to focus on how the behaviour of the aggregate deviates from that of a single spherical particle when an additional, smaller sphere is attached to it.
In this way, the governing dimensionless equations take the form
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\boldsymbol{u}} = 0 , \end{gather}
\begin{gather}\frac{\partial \hat{\boldsymbol{u}}}{\partial t} + \left(\hat{\boldsymbol{u}} \boldsymbol{\cdot}\boldsymbol{\nabla}\right)\hat{\boldsymbol{u}} ={-} \boldsymbol{\nabla}\hat{p}+\frac{1}{Ga}\nabla^2\hat{\boldsymbol{u}} +\hat{\boldsymbol{f}}_{IBM} , \end{gather}
\begin{gather}\frac{{\rm d}\hat{\boldsymbol{u}}_i}{{\rm d}\hat{t}} = \frac{1}{\rho'}\oint_{\hat{\varGamma}_i} \hat{\boldsymbol{\tau}}\boldsymbol{\cdot}\boldsymbol{n}\,{\rm d}\hat{A}+\hat{\boldsymbol{F}}_{{b},i} - \frac{1}{\rho'} , \end{gather}
\begin{gather}\frac{{\rm d}\hat{\boldsymbol{\omega}}_i}{{\rm d}\hat{t}} = \frac{1}{\rho'} \oint_{\hat{\varGamma}_i}\hat{\boldsymbol{r}}\times \left(\hat{\boldsymbol{\tau}} \boldsymbol{\cdot}\boldsymbol{n}\right){\rm d}\hat{A}+\hat{\boldsymbol{M}}_{{b},i} , \end{gather}
where 
$\hat {\boldsymbol {u}}_i = \boldsymbol {u}_i/u_{ref}$ and 
$\hat {\boldsymbol {\omega }}_i=\boldsymbol {\omega }_i \, t_{ref}$ are the dimensionless particle translational and angular velocities, respectively. Also, 
$\hat {\boldsymbol {F}}_{{b},i} = \boldsymbol {F}_{{b},i}/(m_ig')$ and 
$\hat {\boldsymbol {M}}_{{b},i} = D_{L}\boldsymbol {M}_{{b},i}/(I_ig')$ denote the dimensionless bond force and moment. The Galileo number,
\begin{equation} Ga =\frac{ D_{L}\, u_{ref}}{\nu} = \frac{D_{L}\sqrt{D_{L}g\left(\rho '-1\right)}}{\nu}, \end{equation}
indicates the ratio of gravitational to viscous forces. In the rest of the manuscript, we will discuss dimensionless quantities only, so that we omit the 
$\, \hat {} \,$ symbol.
2.5. Numerical approach
We solve the governing equations based on the IBM implementation in the code PARTIES, as described in Biegert, Vowinckel & Meiburg (Reference Biegert, Vowinckel and Meiburg2017) and Biegert (Reference Biegert2018). In summary, it employs an Eulerian mesh with uniform spacing 
$h = \Delta x = \Delta y = \Delta z$ for calculating the fluid motion. The viscous and convective terms are discretized by second-order central differences, and a direct solver based on Fast Fourier Transforms is used to obtain the pressure. Time integration is generally performed by a third-order Runge–Kutta scheme, although viscous terms are integrated implicitly with second-order accuracy. To compute 
$\boldsymbol {f}_{IBM}$, we implement the method of Kempe & Fröhlich (Reference Kempe and Fröhlich2012), which employs a mesh of Lagrangian marker points on the surface of each particle. At these marker points, we interpolate between the Eulerian and Lagrangian grids in order to determine the value of 
$\boldsymbol {f}_{IBM}$ so as to ensure the no-slip and no-normal flow conditions at the particle surface. In the bond region we found that, for the range of 
$Ga$ considered, any overlap between marker points did not have a noticeable effect on the settling dynamics, and as such did not deactivate those close to the contact point as done in Biegert (Reference Biegert2018). Additional details regarding the implementation for individual spherical particles, along with validation results, can be found in the references by Vowinckel et al. (Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019) and Zhu et al. (Reference Zhu, He, Zhao, Vowinckel and Meiburg2022).
Those earlier investigations did not consider rigid bonds between the individual particles, which are implemented here for the first time. Hence, we provide in the following a detailed description of the method. First, let us consider the particles 
$P_1$ and 
$P_2$, bonded at a shared contact point 
$\boldsymbol {x}_{c}$. We define the contact plane between these particles as the plane that crosses through 
$\boldsymbol {x}_{c}$ and is orthogonal to the line between the particle centres. The bond is implemented via a corrective force and moment designed to hold the particles together at the contact point. This will ensure that, at all times, 
$\dot {\boldsymbol {x}}_{c,1} = \dot {\boldsymbol {x}}_{c,2}$, i.e. the velocities of the contact point evaluated on the surface of each particle are the same, and 
$\boldsymbol {\omega }_{1}=\boldsymbol {\omega }_{2}$, i.e. the two particles have the same angular velocity so that they rotate as a solid object. The first condition prevents the particles from moving apart or from sliding in opposite directions, while the second prevents the particles from rolling along each other's surface, and from having different angular velocities around the axis connecting them. We note that, in principle, the IBM for complex aggregates could be implemented without the use of bonds, by treating the aggregate as a single rigid object and distributing the Lagrangian marker points over its more complex surface. However, we found the present approach to be more straightforward to parallelize for aggregates of complex shapes, and also to be easily applicable to aggregates that may undergo breakup, which are scenarios that we plan to address in the future. We approximate the rigid bond through the use of spring forces and moments, as illustrated in figure 2, in order to satisfy an approximate form of the above conditions.
Figure 2. Sketch of the forces and torques acting on the two spherical particles 
$P_1$ and 
$P_2$ that form the aggregate, due to the rigid bond connecting them.
To compute the force and moment due to the bond, we define the bond force and moment acting on particle 
$i$ by particle 
$j$ by splitting each into their normal and tangential components
\begin{gather} \boldsymbol{F}_{{b},ij} = \boldsymbol{F}_{{b},ij}^n+\boldsymbol{F}_{{b},ij}^t, \end{gather}
\begin{gather}\boldsymbol{M}_{{b},ij} = \boldsymbol{r}_{{c}, i}\times \boldsymbol{F}_{{b},ij}^t + \boldsymbol{M}_{{b},ij}^n+\boldsymbol{M}_{{b},ij}^t, \end{gather}
with 
$\boldsymbol {F}_{{b},ij}^n$ and 
$\boldsymbol {F}_{{b},ij}^t$ being the components of the bond force normal and tangential to the contact plane between the two particles (with the normal and tangential components marked in superscript as 
$n$ and 
$t$, respectively), 
$\boldsymbol {M}_{{b},ij}^n$ and 
$\boldsymbol {M}_{{b},ij}^t$ being the corresponding components of the moment and 
$\boldsymbol {r}_{{c}, i} = \boldsymbol {x}_{c}-\boldsymbol {x}_{i}$ being the distance between the contact point and the centre of the particle 
$i$.
We compute the normal and tangential forces and moments using a variation of the model described by Potyondy & Cundall (Reference Potyondy and Cundall2004), which considers the particle–particle bond as an area of infinitely small springs along the contact plane. We initialize the force and moment for each bond to be zero at 
$t=0$. Then, at each time step, we compute the dimensionless forces and moments as
\begin{gather} \boldsymbol{F}_{{b},ij}^n\left(t+\Delta t\right) = {\boldsymbol{\mathsf{R}}}\boldsymbol{F}_{{b},ij}^n\left(t\right)-k_{1}\Delta \dot{\boldsymbol{x}}_{{c},ij}^n\Delta t, \end{gather}
\begin{gather}\boldsymbol{F}_{{b},ij}^t\left(t+\Delta t\right) = {\boldsymbol{\mathsf{R}}}\boldsymbol{F}_{{b},ij}^t\left(t\right)-k_{1}\Delta \dot{\boldsymbol{x}}_{{c},ij}^t\Delta t, \end{gather}
\begin{gather}\boldsymbol{M}_{{b},ij}^n\left(t+\Delta t\right) = {\boldsymbol{\mathsf{R}}}\boldsymbol{M}_{{b},ij}^n\left(t\right)-k_{2}\Delta \boldsymbol{\omega}_{ij}^n\Delta t, \end{gather}
\begin{gather}\boldsymbol{M}_{{b},ij}^t\left(t+\Delta t\right) = {\boldsymbol{\mathsf{R}}}\boldsymbol{M}_{{b},ij}^t\left(t\right)-k_{2}\Delta \boldsymbol{\omega}_{ij}^t\Delta t, \end{gather}
with 
${\boldsymbol{\mathsf{R}}}$ a rotation operator that aligns the previous bond force and moment with the contact plane at the current time, 
$\Delta \dot {\boldsymbol {x}}_{{c},ij}$ the translational velocity difference between the particles at the contact point and 
$\Delta \omega _{ij}$ the angular velocity difference between the particles (with the normal component being the rotation aligned with the contact plane, and the tangential component indicating the rotation perpendicular to the plane). Here, 
$k_{1}$ denotes the non-dimensional spring constant for the bond force, and 
$k_{2}$ indicates the corresponding constant for the bond moment. For computational efficiency, we choose a first-order-in-time method in accordance with the original method outlined in Potyondy & Cundall (Reference Potyondy and Cundall2004), and employ a sufficiently small time step to ensure convergence. The advantage of the current bond model lies in the fact that, for sufficiently stiff springs, we approximate a rigid aggregate that is ‘glued together’ and settles as a single body. Further, by modelling a complex body as many elementary spherical shapes, the Lagrangian marker distribution is simplified compared with performing an IBM implementation on a single, irregularly shaped body. In the current implementation we do not consider the breakage of bonds. This allows us to collapse the spring coefficients of the original model, which are based on the material properties and strength of the bonding material, to a single coefficient 
$k = k_{1} = k_{2}$, under the assumption that for rigidity to be ensured, we only need to make the springs sufficiently stiff, as will be discussed in further detail below.
2.6. Initial and boundary conditions
The simulations are initiated with the fluid at rest, and employ triply periodic boundary conditions. The aggregate is released from rest, with a horizontal orientation 
$\theta (t=0) = \theta _{0} = {\rm \pi}/2$. When other initial orientations are used, this will be mentioned explicitly. The computational domain has size 
$W_{c}\times H_{c} \times W_{c} = 10 D_{L} \times 50 D_{L} \times 10 D_{L}$, which ensures that the influence of the boundaries remains minimal throughout the simulations. For some cases with large Galileo numbers, we increased the domain height up to 
$120 D_{L}$, in order to allow the aggregate to reach a steady state without interacting with its own wake via the periodic boundaries in the vertical direction. Additionally, we chose a tall domain to be able to capture the entire wake generated by the aggregate. Due to the periodic boundary conditions in the vertical direction, the downward force of the particle acting on the fluid is not balanced and would lead to a constant downward acceleration of the fluid. To counteract this acceleration, we apply a corrective, distributed, upward force per volume 
$\boldsymbol {f}_{adj}$ to the fluid that is uniformly distributed and equal and opposite to the negative buoyancy of the particles (Höfler & Schwarzer Reference Höfler and Schwarzer2000), which in its dimensionless form is
\begin{equation} \boldsymbol{f}_{adj} ={-}\frac{D_{L}^3}{W_{c}^2H_{c}}\frac{\rm \pi}{6}\left(1+\alpha^{{-}1/3}\right) . \end{equation}2.7. Resolution and validation
Validation and convergence information for the case of individual settling spheres can be found in the supplementary material, as well as in Biegert (Reference Biegert2018). For the current scenario of settling aggregates, we obtained some general guidance for the appropriate time step size 
$\Delta t$ by keeping the Courant–Friedrichs–Lewy (CFL) number
\begin{equation} {\rm CFL} = \frac{u \, \Delta t}{h} \end{equation}
below 0.5, where 
$u$ indicates the maximum vertical velocity value within the domain at any given time. A sufficiently low CFL number ensures the time step is small enough to prevent instability in the numerical solution, and to ensure the discretization converges towards an accurate solution of the case being studied. Beyond that, we further decreased the time step based on convergence studies, in order to ensure that any restrictions imposed by the bond model would not affect the numerical results, such that 
${\rm CFL} = 0.1$. This resulted in typical time steps of order 
$\Delta t \sim 10^{-3}$.
When simulating a pair of bonded particles, the spatial and temporal resolutions require additional consideration. The spatial step size 
$h$ needs to be small enough to resolve the fluid–particle interactions, and the temporal step size 
$\Delta t$ has to be sufficiently small to capture the spring-like motion of the bond, whose dimensional period
\begin{equation} T = 2{\rm \pi} \sqrt{\frac{m_{S}}{k}} \end{equation}
is estimated based on the dimensional spring coefficient 
$k$ and the mass 
$m_{S}$ of the smaller particle.
We follow the approach of Vowinckel et al. (Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019) by taking the grid spacing 
$h$ to be at most one twentieth of the average particle diameter, 
${(D_{L}+D_{S})}/{2} \geqslant 20\,h$. To fully resolve the spring-like behaviour, we limit the time step to be at most 1/32 of the spring period. Tests showed the values to be sufficient to ensure both stability and convergence.
To determine a suitable computational bond strength coefficient 
$k$, we assess the deformation of the aggregate by tracking 
$\delta _{err} = ||\boldsymbol {x}_{L} - \boldsymbol {x}_{S}|-(D_{L}+D_{S})/2|$, for particle centres 
$\boldsymbol {x}_{L}$ and 
$\boldsymbol {x}_{S}$, and 
$\theta _{err} = |\theta _{L}-\theta _{S}|$, for particle orientations 
$\theta _{L}$ and 
$\theta _{S}$, as functions of time. Figures 3(a) and 3(b) show that a value 
$k=1000$ suffices to keep these deformation measures below 
$10^{-4}$, which we determined to be sufficiently small compared with the characteristic length scale that it could be safely neglected. Figure 3(c) shows that the evolution of the aggregate's orientation is fully converged for this value of 
$k$, so that we select it for our simulations.
Figure 3. Convergence of simulation results for increasing bond strength: measure of the error (a) 
$\delta _{err}$ associated with the gap size between the particles, and (b) 
$\theta _{err}$ associated with the bending of the aggregate. These error measures are shown as functions of time for 
$Ga = 13$ and 
$\alpha = 1.25$, and for several values 
$k$. (c) Corresponding orientation angles 
$\theta$.
To validate the implementation of the bond between two spheres, we compare with the recent results of Nissanka et al. (Reference Nissanka, Ma and Burton2023), who considered the settling behaviour of two connected spheres in a fluid of density 
$\rho _{f} = 971\ {\rm kg} \ {\rm m}^{-3}$ and viscosity 
$\nu = 0.01\ {\rm m}^2 \ {\rm s}^{{-1}}$. The spheres have equal diameters 
$D = 0.002\ {\rm m}$, but different densities 
$\rho _1 = 1420\ {\rm kg} \ {\rm m}^{{-3}}$ and 
$\rho _2 = 2790\ {\rm kg} \ {\rm m}^{-3}$. In the following, we use the average particle density 
$(\rho _1+\rho _2)/2 = 2105 \ {\rm kg} \ {\rm m}^{-3}$ for determining a representative Galileo number 
$Ga = 0.03$ based on (2.13), to be consistent with the notation of Nissanka et al. (Reference Nissanka, Ma and Burton2023). In the experiments, the two spheres were glued together and released from rest in a tank of dimensions 
$0.004 \ {\rm m}\times 0.19 \ {\rm m}\times 0.15 \ {\rm m}$. Initially, the denser particle was located approximately above the less dense one, so that the pair will rotate upon release until the denser particle is below the less dense one. In the experiment, the location of the contact point between the spheres, as well as the orientation of the pair, were recorded as functions of time.
For the simulation, we take a computational domain of 
$2D\times 40D \times 10D$, which matches the experiment in the 
$x$-direction and is sufficiently large in the 
$y$- and 
$z$-directions so that the boundaries do not influence the settling behaviour. The boundaries themselves are modelled as no-slip walls. The experiments in Nissanka et al. (Reference Nissanka, Ma and Burton2023) were performed with 
$Ga \sim O(10^{-2})$, which is too small for our code to simulate, as it would require too small a time step. We instead perform several simulations for somewhat larger values of 
$Ga$, to show that we converge to the experimental results of Nissanka et al. (Reference Nissanka, Ma and Burton2023) as 
$Ga$ approaches that of their experiments. We initialize the aggregate at an orientation angle of 
$17{\rm \pi} /18$, so that the lighter particle is slightly displaced from being directly under the denser one.
Figure 4 shows the orientation angle 
$\theta$ of the aggregate, along with the vertical position 
$y/D$ of the geometrical centre of the pair (
$t=0$ is chosen to be where 
$\theta = {\rm \pi}/2$) for the simulations, as well as the corresponding experimental data of Nissanka et al. (Reference Nissanka, Ma and Burton2023). Here, time is scaled by the reference time 
$t_{St} = D/u_{St}$, where 
$u_{St}$ denotes the Stokes settling velocity of the lighter sphere. We find that, even when 
$Ga$ is slightly higher than the value of the experiments, when the system is scaled in reference to the Stokes settling velocity the results collapse onto the experimental results, which further validates the computational bond model used in the following.
Figure 4. Comparison with the experimental results of Nissanka et al. (Reference Nissanka, Ma and Burton2023) (figure 3 in their article), for the settling of an aggregate consisting of two spheres of equal size and different densities: (a) orientation angle 
$\theta$ and (b) vertical position 
$y/D$. All results are non-dimensionalized by using the particle diameter 
$D$ as the length scale and 
$t_{St} = D/u_{St}$ as the time scale, where the Stokes settling velocity 
$u_{St}$ is calculated using the density of the lighter particle.
3. Results and discussion
We begin by exploring the parameter range 
$Ga \in [8,28]$ and 
$\alpha \in [1,2]$, for which the experiments most clearly demonstrate the existence of a steady-state terminal settling velocity and orientation. Following comparisons between experiments and simulation results in this regime, we will expand the discussion to the broader range of 
$Ga = O(1)-O(10^2)$ and 
$\alpha \leqslant 4$.
3.1. Temporal evolution
Figure 5 compares the temporal evolution of the experimental and numerical settling behaviour. In the simulations, the aggregates are released from rest, with an orientation angle 
$\theta _{0}$ equal to the initial value observed in the corresponding experiments (note that, for experiments, there is likely always a small initial rotation and velocity generated during the release of the aggregate). The vertical settling velocity 
$u_{y, {term}}/u_{ref}$ and the horizontal drift velocity 
$u_{ z, {term}}/u_{ref}$ of the centre of mass of the pair of particles are recorded as functions of time, along with the orientation angle 
$\theta$. The figure demonstrates good agreement between simulation results and experimental observations, both during the initial transient stage and for the subsequent steady state. Figure 5(b) shows that the horizontal drift velocity generally points in the direction of the lower, larger sphere. We remark that, for a given pair of 
$\alpha$ and 
$Ga$ values, the initial orientation of the aggregate in the simulation has no influence on the final, steady-state settling velocity or orientation.
Figure 5. Comparison of numerical and experimental data for 
$\alpha = 8/7$, 
$Ga = 15$. (a) Reports the orientation angle 
$\theta$ (blue) and the vertical velocity 
$u_y/u_{ref}$ (red). (b) Shows the horizontal velocity 
$u_z/u_{ref}$ (black). Dashed lines indicate the numerical results, while the shaded region indicates the experimental uncertainties of one standard deviation around the average experimental result.
Figures 6(a) and 6(b) show comparisons between experimental and simulation results for the instantaneous drift velocity as a function of the orientation angle for different sets of parameters 
$Ga$, 
$\alpha$ and 
$\theta _{0}$. The figures demonstrate good agreement between numerical and experimental data, and highlight the fact that the largest horizontal drift velocities are observed for orientation angles near 
$\theta ={\rm \pi} /4$.
Figure 6. Instantaneous rescaled drift velocity 
$u_{z}/u_{ref}$ as a function of the orientation angle 
$\theta$, showing both experimental data (symbols) and simulation results (line). Each plot depicts many time steps of a single experiment and simulation, for (a) 
$Ga = 21$, 
$\alpha = 1.4$, 
$\theta _{0} = 2{\rm \pi} /3$ and (b) 
$Ga = 23$, 
$\alpha = 2$, 
$\theta _{0} = 5{\rm \pi} /9$. The largest drift velocities are generally observed for orientations 
$\theta$ near 
${\rm \pi} /4$, while zero horizontal velocity corresponds to the symmetric orientation 
$\theta = {\rm \pi}/2$.
Figure 7(a,b) compares simulation results for two different values of 
$Ga$ but identical 
$\alpha = 1.5$. For the lower 
$Ga$-value, the orientation angle 
$\theta$ decreases monotonically with time until the aggregate becomes aligned in the vertical direction (figure 7a). For the larger value of 
$Ga$, reported in figure 7(b), the orientation angle no longer evolves monotonically and the terminal orientation is non-vertical. Here, the aggregate behaves like an underdamped oscillator, in that it overshoots its equilibrium orientation and settles into a steady state through damped oscillations due to inertia. Again, we observe that the steady state is reached more quickly for the larger Galileo number. These results suggest that, for larger values of 
$Ga$, the aggregate tends to approach its terminal state in an oscillatory fashion, and the terminal state is less likely to be vertically oriented. We will return to these points in more detail in the following.
Figure 7. Numerical simulation results for the time evolution of the orientation angles 
$\theta$ and the velocity component 
$u_y/u_{ref}$ as functions of time for 
$\alpha = 1.5$ and (a) 
$Ga=4$, and (b) 
$Ga = 42$. For the larger 
$Ga$-value the aggregate behaves as an underdamped oscillator.
Figures 8(a) and 8(b) show the temporal evolution of the orientation 
$\theta$ of the aggregate as 
$\alpha$ is varied, for two different values of the Galileo number. Generally, larger values of 
$\alpha$ result in terminal orientations closer to the vertical direction, although at smaller 
$Ga$ a vertical orientation is achieved for much smaller values of 
$\alpha$ than at larger 
$Ga$. At lower Galileo numbers, the aggregates again take significantly longer to rotate into their terminal configuration than at higher Galileo numbers. Indeed, in the example case at 
$Ga=13$ shown in figure 8(a), the length of the transient phase increases with 
$\alpha$, especially when the final orientation is close to vertical (
$\theta = 0$).
Figure 8. Orientation angle of the aggregate over time, for varying values of 
$\alpha$ and fixed Galileo numbers: (a) 
$Ga = 13$, and (b) 
$Ga = 42$. The time required to converge to a steady state depends on both 
$\alpha$ and 
$Ga$.
The evolution of 
$\theta$ with time as 
$Ga$ is varied is shown in figures 9(a) and 9(b), for two different values of the diameter ratio 
$\alpha$. The time to reach the terminal state usually decreases for larger 
$Ga$, in agreement with the previous paragraph. We note that this observation does not hold for certain cases where the terminal orientation approaches the vertical, such as 
$\alpha =2$ with 
$Ga = 13$. For the same 
$\alpha$, a smaller 
$Ga$ generally leads to a more vertically aligned terminal state.
Figure 9. Orientation angle of the aggregate as a function of time, for various Galileo numbers and (a) 
$\alpha = 1.5$ and (b) 
$\alpha = 2$. For a constant value of 
$Ga$, a larger 
$\alpha$ leads to a more vertical terminal orientation.
Figures 10(a) and 10(b) show the time-dependent horizontal drift velocity 
$u_z/u_{ref}$ as a function of the orientation angle 
$\theta$, for different values of 
$Ga$ and 
$\alpha$. For a constant 
$\alpha$, we find that the maximum drift velocity increases with 
$Ga$, while for constant 
$Ga$, it increases with decreasing 
$\alpha$. The emerging spiral shapes seen for the largest values of 
$Ga$ and the smallest values of 
$\alpha$ reflect the underdamped oscillations described above.
Figure 10. Horizontal drift velocity as a function of 
$\theta$, for (a) different values of 
$Ga$ and 
$\alpha = 1.5$, and (b) 
$Ga = 22$ and various values of 
$\alpha$. Varying 
$Ga$ and 
$\alpha$ lead to significant variation in the magnitude of the maximum drift velocity, as well as the value of 
$\theta$ corresponding to that maximum.
We also vary the density ratio 
$\rho '$, to elucidate its effect on the settling behaviour. In figure 11 we show simulations for aggregates with different 
$\rho '$, while keeping 
$Ga$ fixed. We find that altering 
$\rho '$ only influences the transient behaviour, while the terminal settling properties do not depend on 
$\rho '$.
Figure 11. Orientation angle (a) and horizontal settling velocity (b) for aggregates with different density ratios 
$\rho '$, for 
$Ga = 18$ and 
$\alpha = 1.5$. While 
$\rho '$ modifies the transient dynamics, it does not affect the terminal settling properties.
In summary, all of the above combinations of 
$\alpha$ and 
$Ga$ demonstrate the emergence of a terminal state characterized by steady values of the settling and drift velocities, and of the orientation angle. Lower Galileo numbers 
$Ga$ and larger aspect ratios 
$\alpha$ favour a monotonic evolution of these quantities towards their terminal values, whereas larger 
$Ga$ and smaller values of 
$\alpha$ are seen to promote the emergence of underdamped oscillations. Having considered the transient evolution of the aforementioned values, we consider in the following section the terminal behaviour.
3.2. Terminal behaviour
To compare some aspects of our experimental and simulation results with the existing literature, it is useful to discuss our data in terms of both the Galileo number 
$Ga$, which is formed with the buoyancy velocity (cf. (2.13)), and the Reynolds number 
$Re$, which is based on the terminal settling velocity (cf. (2.3)). Based on a series of experiments for a single spherical particle, Nguyen et al. (Reference Nguyen, Stechemesser, Zobel and Schulze1997) proposed an empirical relationship between 
$Re$ and 
$Ga$
\begin{equation} Re = \left[\frac{18}{Ga^2}+\frac{3}{16}\left(1+0.079Ga^{1.498}\right)^{{-}0.755}\right]^{{-}1}. \end{equation} Figure 12 demonstrates that, across the parameter range that we investigated, this empirical relationship captures qualitatively well our experimental and numerical results, including for different values of 
$\alpha$. We further notice that there seems to be a nearly single-valued relationship between 
$Ga$ and 
$Re$, despite 
$Ga$ not being affected by changes in 
$D_{S}$ (and by extension the overall aggregate mass) while 
$Re$ does vary with 
$D_S$, since the presence of the smaller sphere modifies the settling velocity. This reflects the fact that in the present range 
$Ga \in [4,75]$, the settling velocity is a weak function of 
$\alpha$. Hence, for this range of 
$Ga$ we can use (3.1) to compare the results as functions of 
$Re$ and 
$Ga$, respectively.
Figure 12. Relationship between 
$Ga$ and 
$Re$. The solid line represents equation (3.1) derived by Nguyen et al. (Reference Nguyen, Stechemesser, Zobel and Schulze1997). Filled circles show simulation results. The colour bar indicates the value of 
$\alpha$, and 
$+$ symbols represent experimental data.
We now focus on the terminal values of the orientation angle, 
$\theta _{term}$, as well as the corresponding settling and drift velocity components, 
$u_{y, {term}}$ and 
$u_{ z, {term}}$. In figure 13, we compare the terminal orientation angle and the settling and drift velocities measured experimentally with their numerical counterparts, for varying values of 
$Ga$ and 
$\alpha$.
Figure 13. Comparison between experimental (
$\triangle$ symbols with error bars) and numerical (filled circles) results for (a) the terminal orientation angle, (b) the settling velocity and (c) the drift velocity, for various values of 
$Ga$ and 
$\alpha$. Colours indicate average 
$Ga$ values for experiments with similar 
$Ga$, as the Galileo number could not be fixed to an exact value between experiments for different 
$\alpha$ values. Details regarding the best-fit curves are provided in the text.
Within the parameter range explored here, we observe that all three terminal variables vary strongly with 
$Ga$ and 
$\alpha$. In figure 13(a), both simulations and experiments demonstrate a pronounced decrease of 
$\theta _{term}$ from 
$\theta _{term} = {\rm \pi}/2$ when increasing 
$\alpha$ at approximately (taking into account variations in the experimental viscosity as discussed in § 2.2) constant 
$Ga$, indicating that more unequal sizes of sphere tend to align the aggregates increasingly in the vertical direction (i.e. 
$\theta _{term} = 0$). We also note that the horizontal settling orientation seen for 
$\alpha = 1$ corresponds to the results of Khayat & Cox (Reference Khayat and Cox1989), which show that objects with fore–aft symmetry will align themselves horizontally with the flow. Figure 13(b) suggests that, for 
$Ga \approx 22$, the terminal settling velocity of the aggregates decreases with 
$\alpha$ for small values of 
$\alpha$, reaches a minimum for an intermediate value 
$\alpha \approx 1.4$ and subsequently grows again. For 
$Ga \approx 8$ and 16, the variation of the settling velocity with 
$\alpha$ evolves similarly but is less pronounced than for higher 
$Ga$. Interestingly, figure 13(c) demonstrates that the terminal horizontal drift velocity increases with 
$\alpha$ for small 
$\alpha$, then reaches a maximum and subsequently decreases as 
$\alpha$ grows further. This is consistent with the assumption that as 
$\alpha \rightarrow \infty$ the effect of the smaller particle on the settling velocity of the aggregate becomes negligible so that the aggregate behaves similarly to a single settling sphere, which has no horizontal terminal velocity component.
The solid lines in figure 13 represent least squares fits of the simulation results for the terminal variables as functions of 
$\alpha$ for a fixed value of 
$Ga$, which were obtained as follows. In the experiments, 
$Ga$ cannot be fixed to an exact value as 
$\nu$ varies due to room temperature variations. Hence, no two simulations plotted in figure 13 have the exact same 
$Ga$: the fittings were instead performed for an average 
$Ga$ calculated from nearby values of 
$Ga$. For the terminal orientation angle, we assume that the relationship between 
$\theta _{term}$, 
$\alpha$ and 
$Ga$ can be fitted by the form
\begin{equation} \theta_{term}\left(\alpha,Ga\right) = \left[{\rm \pi}/2-\theta_{\infty} \left(Ga\right)\right] \exp({-A (Ga)(\alpha-1)})+\theta_{\infty} (Ga), \end{equation}
where 
$\theta _{\infty }(Ga)$ indicates the orientation angle as 
$\alpha \rightarrow \infty$ as a function of 
$Ga$, and 
$A (Ga)$ is a fitting parameter different for each value of 
$Ga$. For the terminal settling velocity, 
$u_{y, {term}}/u_{ref}$, we expect that it will converge to the settling velocity of a single sphere of diameter 
$D_{L}$, 
$u_{L, term}$, as 
$\alpha \rightarrow \infty$, so that we fit a polynomial of the form
\begin{equation} \frac{u_{y, {term}}\left(\alpha,Ga\right)}{u_{ref}} = \frac{u_{L, term}}{u_{ref}}+C_1\left(Ga\right)\alpha^{{-}3}+C_2\left(Ga\right)\alpha^{{-}2}+C_3\left(Ga\right)\alpha^{{-}1}\end{equation}
to the data. Here, 
$C_i(Ga)$ denote fitting parameters.
As mentioned above, we expect the terminal drift velocity 
$u_{ z, {term}}$ to vanish for 
$\alpha = 1$ and in the limit 
$\alpha \rightarrow \infty$. Hence, we assume a relationship of the form
\begin{equation} \frac{u_{z, {term}}\left(\alpha,Ga\right)}{u_{ref}} = J_1(Ga) \left(\alpha-1\right)\exp({-J_2(Ga)(\alpha-1)}) , \end{equation}
where 
$J_i(Ga)$ are fitting parameters. While other correlations may be possible, we chose this relation based on our observations of the settling aggregate for the values of 
$Ga$ and 
$\alpha$ considered.
The above comparisons show that the results of the numerical simulation agree well with the experimental data for the parameter regime covered by the experiments. This provides the necessary validation to explore a wider parameter range through the numerical simulations alone, as shown in figures 14(a)–14(d). The results of the numerical simulations shown in figure 14(a) indicate that, for 
$Ga\lesssim 25$, the terminal orientation angle asymptotically approaches the vertical direction as 
$\alpha \rightarrow \infty$. As 
$Ga$ increases within this range, this asymptotic convergence with 
$\alpha$ slows down, and outside of this range, the terminal orientation displays the emergence of a finite-valued plateau for 
$\alpha \approx 4$. The dynamics within this range of relatively small 
$Ga$ is qualitatively similar to that derived theoretically by Candelier & Mehlig (Reference Candelier and Mehlig2016) for two spheres connected by a massless, infinitely thin rod settling at small, but finite 
$Re$. If the rod length is much larger than the particle diameters, the authors find
\begin{equation}
\theta_{term}\left(\alpha,Re\right) = \left\{
\begin{array}{@{}ll} \dfrac{\rm \pi}{2}-{\rm \pi} \, {\rm
arcsin}\left(\dfrac{16\alpha\left(\alpha-1\right)}{3Re}\right)
+
O\left(\left(\dfrac{\alpha^2}{D_{L}}-\dfrac{\alpha}{D_{L}}\right)\delta
x\right) & \alpha < \alpha_{0} \\[10pt] 0 & \alpha
\geqslant \alpha_{0} \end{array}, \right.
\end{equation}where
\begin{equation} \alpha_{0} = \frac{1+\sqrt{1+\dfrac{3}{4}Re}}{2} , \end{equation}
where 
$\delta x$ is the dimensionless distance between particle centres. The equation is valid when 
$\delta x Re/\alpha \ll 1$. While we are able to qualitatively capture the decrease in the orientation angle using (3.5), the prediction overshoots the numerical results for the present Galileo numbers. This overshoot can be attributed to the difference in 
$Re$ values from the region where (3.5) holds, and the absence in our simulations and experiments of the long connecting rod assumed in deriving this equation.
Figure 14. Numerical simulation results for the (a) terminal orientation angle, (b) vertical velocity normalized by the velocity of a sphere of the same volume as the aggregate, (c) vertical velocity normalized by the velocity of the larger sphere alone and (d) horizontal drift velocity normalized by the settling velocity of a sphere with the same volume as the aggregate, for various values of 
$Ga$ and 
$\alpha$. The solid lines represent empirical fits, as explained in the main text.
We now return to the observation that, for 
$Ga \gtrsim 25$, figure 14(a) indicated that the orientation angle of the aggregate no longer approaches the vertical direction 
$\theta _{term} = 0$ for 
$\alpha \gg 1$. Instead, the aggregate takes an inclined orientation that is tilted further away from the vertical as 
$Ga$ increases. We propose the following explanation for these non-zero terminal orientation angles: as 
$Ga$ (or 
$Re$) increases, the flow begins to separate near 
$Ga \approx 28$ (
$Re \approx 20$) for our geometry. Clift, Grace & Weber (Reference Clift, Grace and Weber1978) observe that, for a smooth sphere, the separation angle varies with 
$Re$ as
\begin{equation} \theta_{sep}\left(Re\right) = 0.236{\rm \pi}\ln\left(\frac{Re}{20}\right)^{0.438} , \quad 20 < Re < 400. \end{equation} The dashed horizontal lines in figure 14(a) represent the values of 
$\theta _{sep}$ for 
$Ga =$ 42 and 75, given by (3.1) and (3.7). For sufficiently large 
$Ga$, when the flow separates, and the aggregate no longer takes a vertical orientation for large 
$\alpha$, the aggregate instead aligns itself such that the smaller sphere is located near the separation line.
This is demonstrated in figures 15(a)–15(b), which show the streamline pattern in the reference frame moving with the aggregate, along with the vorticity field, in the symmetry plane for various combinations of 
$Ga$ and 
$\alpha$. In addition, the geometry of the dominant vortical structures is visualized by means of the 
$Q$-criterion (Chakraborty, Balachandar & Adrian Reference Chakraborty, Balachandar and Adrian2005), where the quantity
\begin{equation} Q = \tfrac{1}{2}\left(||{\boldsymbol{\varOmega}}||^2-||{\boldsymbol{S}}||^2\right), \end{equation}
is evaluated from the vorticity vector 
${\boldsymbol {\varOmega }}$ and the strain rate tensor 
${\boldsymbol {S}}$. Regions with 
$Q>0$ are those where the magnitude of the vorticity is greater than that of the rate of strain in the flow.
Figure 15. (a,b) Steady-state vorticity fields and streamline patterns in the reference frame moving with the aggregate, shown in the symmetry plane. (c,d) Contours of 
$Q$, with the value of 
$Q$ chosen to best demonstrate the characteristic vortical structures around the aggregate. The centre of mass of the aggregate is initially positioned at 
$(y_0,z_0)$. Panels (a,c) are intermediate orientations for 
$Ga=9$, 
$\alpha = 1.25$, 
$Q=0.08$, and (b,d) for 
$Ga=75$, 
$\alpha = 4$, 
$Q = 0.1$. The colour bar indicates the vorticity.
The streamlines and vorticity fields reflect the symmetry of the aggregate, as shown in figure 15. As 
$Ga$ increases, the vortical region emanating from the larger sphere extends further into the wake, and for 
$\alpha \gg 1$ it completely engulfs the smaller sphere, as shown in figure 15(b). Interestingly, this figure also shows that, for 
$Ga=75$, the flow separates from the rear of the larger sphere, with the smaller sphere positioning itself at the separation point, as discussed earlier.
Figures 14(b)–14(c) also compare the terminal settling velocity of the aggregate with those of a sphere with the same volume as the aggregate (
$u_{Eq, term}$, figure b) and of the larger sphere alone (
$u_{L, term}$, figure c), with the single sphere values obtained via corresponding simulations. For small-to-moderate 
$\alpha$, the aggregate settles more slowly than a single sphere of equal volume, but faster than the larger sphere alone. When 
$\alpha$ becomes large, both ratios of the velocity asymptotically approach unity, as expected. However, the rate at which the respective ratios converge to unity depends strongly on 
$Ga$. The numerical data for the settling velocity ratios shown in figures 14(b) and 14(c) are captured well by functions of the form
\begin{equation} \frac{u_{y, {term}}\left(\alpha,Ga\right)}{u_{Eq, term}\left(Ga\right)} = 1-C_1\exp({-C_2(\alpha-1)})\end{equation}and
\begin{equation} \frac{u_{y, {term}}\left(\alpha,Ga\right)}{u_{L, term}\left(Ga\right)} = 1+C_3\exp({-C_4(\alpha-1)}) ,\end{equation}
respectively, where the 
$C_{i}(Ga)$ represent fitting parameters.
Figure 14(d) shows the terminal drift velocity, normalized by 
$u_{Eq, term}$. We note that for 
$Ga \gtrsim 25$ the aggregate retains a significant terminal drift velocity even for values of 
$\alpha$ as large as 
$\alpha = 4$. This is consistent with our earlier observation that in this parameter range, the aggregate orients itself such that the smaller particle is located at the separation line, and therefore the terminal orientation remains asymmetric. As 
$\alpha \rightarrow \infty$, however, we expect the terminal drift velocity to decay to zero (given an infinitesimally small sphere will have little effect on the overall aggregate), in spite of the asymmetric orientation, so that it can be captured by a function of the form
\begin{equation} \frac{u_{z, {term}}\left(\alpha,Ga\right)}{u_{Eq, term}\left(Ga\right)} = J_1(\alpha-1)\exp({-J_2(\alpha-1)}),\end{equation}
where 
$J_1(Ga)$ and 
$J_2(Ga)$ are fitting parameters. We also note that (3.11) gives a maximal terminal drift velocity for 
$\alpha _{drift,max} = 1/J_2(Ga)+1$ where it is equal to 
$u_{drift,max}/u_{y,Eq,term} = J_1/(J_2e)$.
We furthermore note that, for aggregates with intermediate orientation angles 
$0 <\theta <{\rm \pi} /2$, a region of larger pressure exists on the side of the larger particle directly below the smaller one. This generates a horizontal pressure gradient across the large sphere, which is the primary reason for the horizontal drift velocity.
Figures 16(a)–16(c) summarize the dependence of 
$\theta _{term}$, 
$u_{y,term}/u_{Eq, term}$ and 
$u_{z, {term}}/u_{Eq, term}$ on 
$\alpha$ and 
$Ga$. To plot these figures, we fit 
$\theta _{\infty }(Ga)$ in (3.2) by means of a hyperbolic tangent, and 
$A(Ga)$ via a power law
\begin{gather} \theta_{\infty}\left(Ga\right) = 0.412\left(1+\tanh\left(0.067\left(Ga-35.6\right)\right)\right)\!, \end{gather}
\begin{gather}A\left(Ga\right) = \frac{196}{\left(Ga+1.46\right)^2} +1. \end{gather}
Figure 16. Summary of the dependence of the (a) terminal orientation angle of the aggregate, (b) terminal settling velocity and (c) terminal drift velocity on a (
$\alpha,Ga$) phase diagram, based on the fits described in the text. The dashed line indicates the location of the maximum drift velocity.
These empirical fits match well with the experiments and the numerical simulations, deviating on average by approximately 
$0.11$ radians (7 % of a full 
${\rm \pi} /2$ rotation). Similarly, we obtain for the fitting parameters in (3.9)
\begin{gather} C_1\left(Ga\right) = 0.001Ga+0.155, \end{gather}
\begin{gather}C_2\left(Ga\right) = 1.71\exp({-0.189Ga+1.79})+1.52. \end{gather}
We find that the average deviation for the fitted model is 
$2\,\%$ of a single settling sphere's velocity 
$u_{Eq, term}$. For the horizontal velocity 
$u_{z, {term}}/u_{Eq, term}$, we obtain fits for the parameters in (3.11) of the form
\begin{gather} J_1\left(Ga\right) = 0.233{\rm e}^{{-}0.189Ga+3.31}+0.371, \end{gather}
\begin{gather}J_2\left(Ga\right) = \frac{J_1\left(Ga\right)}{0.001Ga+0.074}{\rm e}^{{-}1}, \end{gather}
with the average deviation being 
$1\,\%$ of a single settling sphere's velocity 
$u_{Eq, term}$. As such, these fits capture the behaviour of the settling aggregates in the parameter space of 
$Ga$ and 
$\alpha$ considered, allowing us to predict the settling of small, two-particle aggregates settling in a fluid.
4. Conclusion
We have investigated the settling dynamics of a model aggregate made of a pair of rigidly connected spherical particles of unequal size. By means of experiments and particle-resolved simulations, we have obtained a detailed picture of both the transient evolution and the terminal values of the orientation angle, settling and drift velocity of the aggregate, as functions of the aspect ratio 
$\alpha = D_{L}/D_{S}$ and the Galileo number 
$Ga = D_{L} u_{ref}/\nu$. For small values of the Galileo number, the orientation of the aggregate and its velocity are seen to converge to their terminal values in a monotonic fashion, whereas, for larger Galileo numbers, the aggregate tends to behave as an underdamped oscillator. The largest drift velocities are generally observed when the aggregate is tilted at approximately 
${\rm \pi} /4$ radians with respect to the vertical direction for lower 
$Ga$, while for higher 
$Ga$ the largest drift velocities tend to be for orientations closer to horizontal. For large aspect ratios and small Galileo numbers, the terminal orientation of the aggregate tends to be vertical, whereas for smaller aspect ratios and larger Galileo numbers, the terminal orientation is inclined with respect to the vertical, which also results in a non-zero terminal drift velocity. When the Galileo number is sufficiently large for flow separation to occur, aggregates with large aspect ratios orient themselves such that the smaller sphere is located at the separation line. Empirical scaling laws are obtained for the terminal settling velocity and orientation angle as functions of the aspect ratio and Galileo number. An analysis of the accompanying fluid velocity field indicates the formation of a horizontal pressure gradient across the larger sphere, which represents the main reason for the emergence of the drift velocity, and it shows the formation of vortical structures exhibiting complex topologies in the aggregate's wake. We also note that, due to these asymmetries in the aggregate, during sedimentation, the particles will likely disperse greatly from their initial locations over long settling times, indicating the importance of characterizing the relationship that the asymmetry, represented by 
$\alpha$, has with the velocity.
Further work should extend the results of the present study to higher ranges of 
$Ga$ than those considered here. However, based on the present results we can make some predictions as to the expected behaviour under those conditions. First, for higher 
$Ga$ one can expect 
$\theta _{term}$ to eventually cease to converge to a steady value, as a result of vortex shedding and turbulent wake dynamics, so that we would expect periodic or chaotic behaviour. Similarly, one can expect the purely planar behaviour, where the aggregate remains in the 
$(y,z)$-plane, to also end, with the aggregate rotating or tumbling out of plane as a result of the asymmetries caused by turbulent flow past the larger particle. Additional work is required to quantify the dynamics under such conditions.
The present computational framework for establishing rigid bonds between spherical particles will enable us to explore more complex aggregates consisting of many primary particles of potentially different sizes (and potentially in different flow geometries). Simulations of this type, in turn, will allow us to investigate the influence of the porosity and permeability of the aggregate on its transient and terminal settling dynamics, and to assess how its effective settling and drift velocities are related to its geometrical properties, such as the fractal dimension. Efforts in this direction are currently underway.
Supplementary material and movies
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2024.711.
Funding
This work is supported by grant W912HZ22C0037 from the Engineer Research & Development Center at the U.S. Army Corps of Engineers, by grant W911NF-23-2-0046 from the Army Research Office, and by National Science Foundation grant CBET-2138583. The authors gratefully acknowledge several stimulating discussions with J. Smith, K. Fall and D. Tarpley, as well as with B. Vowinckel, A. Leonelli and R. Sharma.
Declaration of interests
The authors report no conflict of interest.
Open access
