On the clustering of low-aspect-ratio oblate spheroids settling in ambient fluid

Abstract

We have performed particle-resolved direct numerical simulations of many heavy non-spherical particles settling under gravity in the dilute regime. The particles are oblate spheroids of aspect ratio $1.5$ and density ratio $1.5$. Two Galileo numbers are considered, namely $111$ and $152$, for which a single oblate spheroid follows a steady vertical and a steady oblique path, respectively. In both cases, a strongly inhomogeneous spatial distribution of the disperse phase in the form of columnar clusters is observed, with a significantly enhanced average settling velocity as a consequence. Thus, in contrast to previous results for spheres, the qualitative difference in the single-particle regime does not result in a qualitatively different behaviour of the many-particle cases. In addition, we have carried out an analysis of pairwise interactions of particles in the well-known drafting–kissing–tumbling set-up, for oblate spheroids of aspect ratio $1.5$ and for spheres. We have varied systematically the relative initial position between the particle pair and we have considered free-to-rotate particles and rotationally locked ones. We have found that the region of attraction for both particle shapes, with and without rotation, is very similar. However, significant differences occur during the drafting and tumbling phases. In particular, free-to-rotate spheres present longer drafting phases and separate quickly after the collision. Spheroids remain close to each other for longer times after the collision, and free-to-rotate ones experience two or more collision events. Therefore, we have observed a shape-induced increase in the interaction time which might explain the increased tendency to cluster of the many-particle cases.

1. Introduction

Particle-laden flows play an important role in natural and industrial systems, such as sediment transport in rivers, fluidized beds and pollutants or hydrometeors in the atmospheric boundary layer. It is well known that an isolated particle settling under gravity in an otherwise ambient fluid exhibits a variety of regimes of motion depending on its shape, mass density and size (Ern et al. 2012). Even for spherical particles diverse path regimes have been observed, ranging from steady vertical to chaotic motion, and including various intermediate states of different kinematic complexity (Jenny, Dušek & Bouchet 2004; Zhou & Dušek 2015, cf. also figure 1). The transitions between the distinct regimes of particle motion are the consequence of bifurcations in the flow pattern around the mobile particle, arising as the values of the governing parameters are varied. These parameters in the simplest single-particle case are the solid-to-fluid density ratio, $\tilde {\rho }= {\rho _p}/{\rho _f}$, and the Galileo number, ${Ga}=DU_g/\nu$, where D is the particle diameter (in the case of non-spherical particles D is defined as the diameter of a sphere with the same volume), ${U_g}=(|\tilde {\rho }-1|\|{\boldsymbol {g}}\|{D})^{1/2}$ is a gravitational velocity scale, ${\nu }$ is the kinematic viscosity and ${\boldsymbol {g}}$ is the vector of gravitational acceleration.

Figure 1. Single-particle regimes for heavy ($\tilde {\rho }=1.5$) spheres and oblate spheroids with aspect ratio $\chi =1.5$ as a function of the Galileo number ${Ga}$. Reference data for dilute suspensions of spheres with the same density ratio and Galileo number from Uhlmann & Doychev (2014) are also included. The vortical flow structures in each regime are indicated with the aid of iso-surfaces of ${Q}$, the second invariant of the velocity gradient tensor (Hunt, Wray & Moin 1988).

Whenever the volume fraction ($\varPhi$) of the particulate phase becomes non-negligible, particle–particle interactions can lead to significant collective effects. For dilute systems, these interactions are predominantly of indirect type, acting through long-range hydrodynamic forces. For denser systems, direct contacts between two or more particles come into play more frequently, thereby contributing to the exchange of momentum. In the present work we are interested in dilute systems (with solid volume fractions below one per cent), which is why we will concentrate on this regime in the following literature review. Depending on the values of the parameter triplet $(\varPhi,\tilde {\rho },{Ga})$, dilute suspensions can exhibit spatial particle distributions with significant non-homogeneity. A non-homogeneously distributed particle phase in turn is often accompanied by important macroscopic effects, such as an altered mean settling velocity and a modification of the induced fluid flow features.

At Galileo numbers of ${O}(10)$ and density ratio $\tilde {\rho }=2$ it has been observed that dilute suspensions of spherical particles tend to form horizontally aligned pairs (Yin & Koch 2008), while at larger Galileo numbers (${Ga}\approx 200$) and similar density ratios, where the wake flow features a significant recirculation region, a vertical alignment is more probable (Kajishima & Takiguchi 2002; Uhlmann & Doychev 2014). Kajishima & Takiguchi (2002) were the first to investigate the formation of large, columnar-shaped particle clusters due to wake effects. Their particle-resolved direct numerical simulations (PR-DNS) in triply periodic boxes showed that clusters tend to form for particle Reynolds numbers exceeding $200\ldots 300$ ($Ga\approx 153\ldots 210$), leading to strongly enhanced average settling velocities. In their initial work, the particle rotation was suppressed for computational simplicity, and a density ratio $\tilde {\rho }=8.8$ was considered. In the follow-up work of Kajishima (2004), the effect of particle rotation was taken into account, and it was found that rotational motion leads to a lower concentration of particles in clusters, since particles tend to escape a cluster through a rotation-induced lift effect. Kajishima (2004) also determined a lower limit of the solid volume fraction for the occurrence of clustering. The PR-DNS of Uhlmann & Doychev (2014) were performed at a fixed density ratio $\tilde {\rho }=1.5$, and it was observed that columnar clusters do form for ${Ga}=178$ (which corresponds to an isolated particle in the steady oblique regime), while no clusters are formed at $Ga=121$ (steady vertical regime). The authors suggest that the onset of clustering in dilute suspensions of spherical particles is triggered by the bifurcation of the wake flow from steady axisymmetric to steady oblique, with the consequence that the particles in the latter case drift horizontally (with random azimuthal angle) leading to an enhanced probability of particle–particle encounters. Further data on collective effects upon clustering are available from the PR-DNS studies of Zaidi, Tsuji & Tanaka (2014), Fornari, Picano & Brandt (2016a) and Seyed-Ahmadi & Wachs (2021), which were all performed in similar triply periodic configurations, and from the work of Huisman et al. (2016), who conducted experiments in a settling column with glass beads in water. An overview of the average settling velocities from these various datasets for dilute suspensions of spherical particles in an otherwise ambient fluid is given in figure 2. The PR-DNS results appear to give a consistent trend: for increasing values of the Galileo number the particles tend to settle at an average rate which is enhanced with respect to the velocity of an isolated particle, and there appears indeed to be a cross-over (from a reduced settling velocity to an enhancement) for $Ga\approx 150$. The limited available experimental data in this parameter range are not inconsistent with this picture, however, showing a settling enhancement already at lower Galileo number of $110$ and mild columnar cluster formation. Huisman et al. (2016) have identified the presence of the bounding container walls, which presumably causes a large-scale recirculating flow, as a possible cause for the observed discrepancy.

Figure 2. Mean settling velocity vs Galileo number for spherical and non-spherical suspensions of heavy particles with density ratio ${O}(1)$ in the dilute regime. The velocity data are normalized with the corresponding mean settling velocity of an isolated particle in the asymptotic (long-time) limit. The error bars in the experimental data of Huisman et al. (2016) indicate minimum and maximum values of the repetitions performed by the authors. Present results are included for completeness.

Dilute suspensions of non-spherical particles have received much less attention, mainly due to the complexity associated with the particle shape and the corresponding increase in the size of the parametric space. Spheroids have been the preferred choice of several authors, partly due to their convenient parametrization. Their shape is defined uniquely by the aspect ratio $\chi ={d}/ {a}$, where d and a are the equatorial diameter and the length of the symmetry axis, respectively (see figure 3a,b). The spheroidal shape allows for a smooth transition between flat, disk-like shapes (for very large values of $\chi$), and elongated, fibre-like geometries (for very small values of $\chi$), while a sphere is recovered for $\chi =1$. When oblate spheroids ($\chi >1$) are considered, the regime maps of a single particle are more complex when compared with spheres (Zhou, Chrust & Dušek 2017). Despite this increase in complexity, specific combinations of $\chi$ and $\tilde {\rho }$ result in a single oblate spheroid exhibiting the so-called ‘sphere-like scenario’ in which the first bifurcation for increasing Ga is regular, transitioning from a vertical to an oblique regime (see figure 1, Zhou et al. 2017; Moriche, Uhlmann & Dušek 2021). Fornari, Ardekani & Brandt (2018) studied the settling of dilute suspensions of almost neutrally buoyant ($\tilde {\rho }=1.02$) oblate spheroids with a moderately flat shape $\chi =3$. They considered both dilute and dense regimes at Galileo numbers at which a single particle follows a steady vertical path ($Ga=60$) or an unsteady path which is vertical in the mean ($Ga=140$). For their most dilute cases they observed a large enhancement of the mean settling velocity compared with the single-particle case, with only small differences between the two Galileo numbers (cf. figure 2). Fornari et al. (2018) have also detected non-homogeneous particle distributions in the form of vertical particle trains with lateral dimensions of the order of $10$ times the length of the symmetry axis, based upon visualization and on the analysis of pairwise distribution functions. The authors attributed the enhancement of settling velocity to the occurrence of particle-pair interactions which lead to the formation of piles of particles which practically stick to each other in the case of oblate spheroids with $\chi =3$.

Figure 3. View of the $i$th spheroid in its body-fixed reference system (a) along the symmetry axis and (b) perpendicular to it. The blue line in (a,b) represents a sphere with the same volume. (c) Sketch of the problem in the global reference system.

Turning now to other shapes, Seyed-Ahmadi & Wachs (2021) more recently studied the behaviour of settling cubes with $\tilde {\rho }=2$, $\varPhi =0.01\ldots 0.2$ and Galileo numbers for which a single cube tends to a steady oblique path with very small tilting angle ($Ga=70$) and to a helical path ($Ga=160$) for large times. For their smallest solid volume fraction these authors report an increase in the relative settling velocity, which is similar to what is observed in suspensions of spheres (cf. figure 2). Interestingly, Seyed-Ahmadi & Wachs (2021) observe that the intensity of columnar cluster formation (measured with the aid of particle-pair distribution functions) is somewhat lower in the case of cubes as compared with the spherical counterpart. This effect is attributed by the authors to higher rotation rates in the former case, leading to an enhanced lift force that increases their probability to escape from an existing cluster.

As noted above in the case of spheres, the settling regime of a single particle appears to be relevant to the macroscopic behaviour of dilute suspensions. As a next step towards the understanding of collective effects, it is useful to investigate the interaction of settling particle pairs, a set-up which may lead to the so-called drafting–kissing–tumbling (DKT) process (Fortes, Joseph & Lundgren 1987). For sufficiently high Galileo number, a trailing sphere initially released inside of the wake region of a leading sphere will approach the latter one during the drafting phase, they will touch (‘kissing’), and then interchange positions (‘tumbling’), since a vertically aligned pair of spheres is unstable, before eventually separating. The reduced size of the problem compared with the many-particle cases makes (DKT) simulations a useful laboratory to understand the underlying physics of many-particle cases. Indeed, Fortes et al. (1987) used auxiliary cases with the DKT set-up to support their results on the spatial structure of many spherical particles in a fluidized bed. Regarding non-spherical particles, Ardekani et al. (2016) observed a suppression of the tumbling phase in pairwise interactions of moderately flat oblate spheroids ($\chi =3$), as well as an increased collision domain. The authors conjectured that the essentially infinite interaction time of piled up particles would lead to enhanced clustering in the corresponding many-particle case. This conjecture was later indeed confirmed by Fornari et al. (2018).

The numerical simulations of DKT cases in the literature use a common configuration: two particles settling in an initially quiescent fluid inside a container, whose vertical size should be large enough to accommodate the different stages of the DKT case, and also minimize the (undesired) influence of the lower boundary of the computational domain. The vertical size of the container varies from $20D$ (Glowinski et al. 2001; Breugem 2012) to approximately $40D$ (Patankar et al. 2000; Ardekani et al. 2016). The main drawback of this configuration is that any modification of the parameters (for example the relative initial position between the particles) which increases the duration of the drafting phase, would require a further increase of the computational domain. As a result, the simplicity of the DKT set-up is somehow diminished. Using periodic boundary conditions in the vertical direction, as some authors have done for single-particle cases (Kajishima & Takiguchi 2002; Doychev 2014), also demands large domains in order to minimize wake effects from periodic repetitions. To overcome these difficulties, in the present work we propose to use inflow/outflow boundary conditions along the vertical direction along with a carefully adjusted vertical velocity at the inflow. Thus, the settling particles are not perturbed by their own wakes, and the domain size becomes less of a critical parameter. To the best of our knowledge, this strategy has not been employed previously. Here, it allows us to employ moderately small computational domains, so that we can perform a campaign of simulations varying the relative initial position of the particle pair and, therefore, we can obtain statistically significant results.

In the present work we analyse the clustering behaviour of dilute suspensions of spheroids with a relatively low aspect ratio $\chi =1.5$, a density ratio $\tilde {\rho }=1.5$ and two Galileo numbers $Ga=\{111, 152\}$. For these spheroids a single particle exhibits the so-called ‘sphere-like scenario’ (Zhou et al. 2017), and for suitable initial conditions a pair of them undergo the three stages in a pairwise interaction (DKT). Our aim is to explore collective effects for non-spherical particle shapes, while remaining relatively close to the well-explored spherical shape. The overarching question is: How do relatively small deviations from the spherical shape affect the settling behaviour of a dilute suspension of rigid particles? In addition to many-particle simulations in large domains, we perform a large set of DKT simulations of the mentioned spheroids as well as for the reference case with spheres in which we analyse the effect of particle shape and of the angular motion by either suppressing or allowing the latter.

The manuscript is structured as follows: in § 2 we present the mathematical model describing the settling of particles in unbounded, otherwise quiescent fluid and, in § 3 we specify the numerical method and the physical and numerical parameters. Results are presented in two steps: in § 4.1 we focus on many-particle cases, which represent the core of this work, and in § 4.2 we analyse a set of auxiliary simulations of a DKT configuration in order to support the observations made in the many-particle cases. In § 5 we discuss the implications of the DKT results for the many-particle case. Finally, a summary and conclusions can be found in § 6.

2. Problem description

We study the settling of particles under the action of gravity in an unbounded, initially quiescent fluid. We consider an incompressible Newtonian fluid of density $\rho _{f}$ and kinematic viscosity $\nu$. Particles are oblate spheroids of equatorial diameter d and aspect ratio $\chi ={d}/a$, where a is the length of their symmetry axis (see figure 3a,b). The particles are assumed to be rigid with homogeneous mass density $\rho _{p}$, which is larger than that of the fluid ($\rho _{p}>\rho _{f}$). The gravitational acceleration $\boldsymbol {g}$ points in the negative $z$-direction, $\boldsymbol {g} = -g\boldsymbol {e}_{z}$, where ${g}$ is the modulus of the gravitational acceleration (see figure 3c).

The fluid velocity $\boldsymbol {u}=({u},v,w)$ is governed by the Navier–Stokes equations for an incompressible, constant density fluid

(2.1a)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\boldsymbol{u} ={-}\frac{\boldsymbol{\nabla} p}{\rho_{f}} +\nu \nabla^2\boldsymbol{u}, \end{gather}$$

(2.1b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0 , \end{gather}$$

where $p$ is the pressure. No-slip and no-penetration boundary conditions are imposed on the surface of each particle. The linear and angular velocities of the particles, $\boldsymbol {u}_{p}=(u_{p},v_{p},w_{p})$ and $\boldsymbol {\omega }_p=(\omega _{px},\omega _{py},\omega _{pz})$, respectively, are governed by the Newton–Euler equations

(2.2a)$$\begin{gather} V_p\rho_p \frac{{\rm d}\boldsymbol{u}_p}{{\rm d}t}=\int_{S} \boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n} \, \text{d} \sigma +\left(\rho_{p}-\rho_{f}\right)V_p\boldsymbol{g} + \boldsymbol{F}, \end{gather}$$

(2.2b)$$\begin{gather}\frac{{\rm d}\left({\boldsymbol{I}_{p}}\boldsymbol{\omega}_{p}\right)}{{\rm d}t}= \int_{S} \boldsymbol{r}\times \left(\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n}\right) \text{d} \sigma + \boldsymbol{T}, \end{gather}$$

where $V_{p}$ and $S$ are the volume and the surface of the particle, respectively, $\boldsymbol {n}$ is a unit vector normal to $S$ pointing towards the fluid, $\boldsymbol {\tau }$ is the stress tensor ($\boldsymbol {\tau }=-pI+\rho _{f}\nu (\boldsymbol {\nabla }\boldsymbol {u}+\boldsymbol {\nabla }\boldsymbol {u}^T)$), $\boldsymbol {r}$ is a position vector with respect to the centre of gravity of the particle and $\boldsymbol {F}$ and $\boldsymbol {T}$ are the solid–solid contact force and torque, respectively.

The problem is governed by four non-dimensional parameters, namely the density ratio between the particles and the fluid, $\tilde {\rho }=\rho _{p}/ \rho _{f}$, the aspect ratio of the particles, $\chi$, the Galileo number, Ga, and the solid volume fraction, $\varPhi$. The solid volume fraction is defined as $\varPhi =\varSigma _p/( \varSigma _p+\varSigma _f)$ where $\varSigma _p$ and $\varSigma _{f}$ represent the volume occupied by the particles and the fluid, respectively.

In this work we are interested in the dilute regime, and we set the solid volume fraction to $\varPhi =5\times 10^{-3}$. Regarding the shape of the particles, we select oblate spheroids of aspect ratio $\chi =1.5$, which is a moderately small deviation from the spherical reference geometry, that in turn has been extensively investigated in the past. In particular, it is known that the settling regime map for oblate spheroids with $\chi =1.5$ resembles that of a sphere (Moriche et al. 2021). We fix the density ratio at $\tilde {\rho }=1.5$ (which corresponds e.g. to some plastic materials in water), and we select two values of the Galileo number: one for which a single particle follows a steady vertical path ($Ga=110.56$), and another which leads to a single particle following a steady oblique path ($Ga=152.02$).

2.1. Definitions

A velocity scale based on the gravitational acceleration, the density ratio and the size of the particle can be defined as follows:

(2.3)\begin{equation} U_{g}=\sqrt{\left|\tilde{\rho}-1\right|{g}D}, \end{equation}

where $D=d\chi ^{-1/3}$ is the diameter of a sphere with the same volume as the spheroid considered. Based on the velocity scale $U_{g}$ (2.3) the Galileo number is defined as

(2.4)\begin{equation} Ga=\frac{U_{g}D}{\nu} , \end{equation}

and an a priori gravitational time scale can be formed as follows: $\tau _{g}=D/U_{g}$. For future reference let us define the relative velocity $\boldsymbol {u}_{pr}^{(i)}=( u_{pr}^{(i)},v_{pr}^{(i)},w_{pr}^{(i)})$ of the $i$th particle with respect to the mean fluid velocity as

(2.5)\begin{equation} \boldsymbol{u}_{pr}^{(i)}(t)=\boldsymbol{u}_{p}^{(i)}(t)-\langle \boldsymbol{u}\rangle_f(t), \end{equation}

where the average operator $\langle {\cdot } \rangle _f$ indicates spatial averaging over the entire domain occupied by the fluid $\varOmega _f$, viz.

(2.6)\begin{equation} \langle {\cdot} \rangle_f =\frac{1}{\varSigma_f} \int_{\varOmega_f} \left({\cdot} \right)\text{d} {\boldsymbol{x}}. \end{equation}

Let us also introduce the time-dependent settling velocity, averaged over the set of particles as

(2.7)\begin{equation} w_{s}(t) = \langle w_{pr}\rangle_p(t) , \end{equation}

where the average operator $\langle {\cdot } \rangle _p$ indicates the ensemble average over the dispersed phase, which for a set of $N_{p}$ particles is expressed as

(2.8)\begin{equation} \langle {\cdot} \rangle_p = \frac{\displaystyle \sum_{i=1}^{{N_{p}}} \left({\cdot}\right)^{(i)}}{ N_{p}}. \end{equation}

Similarly, the standard deviation of the vertical and horizontal components of the linear and angular particle velocity are defined as

(2.9a)$$\begin{gather} w_s^{std}(t) =\langle w_{pr}^{\prime (i)}w_{pr}^{\prime (i)}\rangle_p^{1/2}, \end{gather}$$

(2.9b)$$\begin{gather}u_h^{std}(t) =((\langle u_{pr}^{\prime (i)}u_{pr}^{\prime (i)}\rangle_p+ \langle v_{pr}^{\prime (i)}v_{pr}^{\prime (i)}\rangle_p)/2)^{(1/2)}, \end{gather}$$

(2.9c)$$\begin{gather}\omega_{v}^{std}(t) =\langle \omega_{pz}^{\prime (i)}\omega_{pz}^{\prime (i)}\rangle_{p}^{1/2}, \end{gather}$$

(2.9d)$$\begin{gather}\omega_{h}^{std}(t) =((\langle\omega_{px}^{\prime (i)} \omega_{px}^{\prime (i)}\rangle_p+ \langle\omega_{py}^{\prime (i)}\omega_{py}^{\prime (i)}\rangle_p)/2)^{(1/2)}, \end{gather}$$

where the prime symbol indicates that the quantity is the fluctuating part of the variable, defined as

(2.10)\begin{equation} \left({\cdot}\right)' = \left({\cdot}\right) - \langle {\cdot} \rangle. \end{equation}

Finally, it should be mentioned that the definition of the Galileo number used in this work is equivalent to that used in Fornari et al. (2016b) and Ardekani et al. (2016) for spheroids, and to that of Seyed-Ahmadi & Wachs (2021) for cubes. We have selected this definition over the one used in some works dealing exclusively with spheroids (Zhou et al. 2017; Moriche et al. 2021) in order to recover the same value of Ga when working with spheres ($\chi =1$).

For a multiparticle case with a given parameter pair ($Ga,\tilde {\rho })$, independently of its solid volume fraction $\varPhi$, we define the reference velocity ${w_{ref}}$ as

(2.11)\begin{equation} {w_{ref}}\left({Ga},\tilde{\rho}\right) = \left| w_{pr}\right| , \end{equation}

where the value of $w_{pr}$ in (2.11) is taken from the corresponding single-particle case. Please note that, for the cases considered here, the single-particle case regime is steady, therefore no (temporal or statistical) averaging is needed to compute $w_{ref}$ in (2.11).

According to the above definitions of $w_{s}$ and $w_{ref}$, we define the average particle Reynolds number

(2.12)\begin{equation} Re_{D}=\frac{\left|\langle \omega_s\rangle_t\right|D}{\nu} , \end{equation}

and the single-particle Reynolds number

(2.13)\begin{equation} Re_D^0=\frac{w_{ref}D}{\nu} , \end{equation}

where the operator $\langle {\cdot } \rangle _t$ indicates temporal averaging.

3. Methodology

3.1. Numerical method

The Navier–Stokes equations (2.1) are integrated in time by means of a three-stage Runge–Kutta scheme in which the viscous term is treated implicitly and the advective term explicitly (Rai & Moin 1991; Verzicco & Orlandi 1996). The fractional step method proposed by Brown, Cortez & Minion (2001) is used to fulfil the continuity constraint. Spatial derivatives are approximated with central finite differences of second order on a staggered, uniform, Cartesian grid.

The presence of the body is modelled by the direct-forcing immersed boundary method proposed by Uhlmann (2005) and later extended to track the motion of non-spherical particles by Moriche et al. (2021), where extensive validation of the method can be found. Collisions are modelled with a repulsive short-range normal force (with a range $\Delta x$) so as to avoid non-physical overlapping of particles (details can be found in Appendix A).

In a triply periodic set-up, gravity continuously accelerates the system. Therefore, in order to allow for a steady state, we add a constant-in-space source term to the vertical momentum equation whose volume integral is equal (and opposite in sign) to the net force exerted by the particles on the fluid (Höfler & Schwarzer 2000).

3.2. Computational set-up

The computational domain is a cuboid with triply periodic boundary conditions, whose size is the result of a compromise between computational resource requirements and physical realism. Based on preliminary tests we choose a domain size of approximately $55D$ in the lateral directions and approximately $220D$ in the vertical direction. Please note that a full decorrelation in the vertical direction is not warranted once clustering sets in, and in the horizontal direction once the clusters grow to a size comparable to the computational domain. This lack of full decorrelation, which was also observed in the case of corresponding spheres (Doychev 2014), is further discussed in Appendix B. We select a spatial resolution of $D/\Delta x\approx 21$ (${d}/ \Delta x=24$), which is supported by the work of Moriche et al. (2021). In their work the authors show that, compared with a spectral/spectral-element solution, the error in the mean settling velocity of a single spheroid of aspect ratio $\chi =1.5$, density ratio $\tilde {\rho }=2.14$ at $Ga=152.02$ is smaller than $2\,\%$. In the present case this results in a grid of $[1152 \times 1152 \times 4608]$ points. Table 1 shows the parameters of the two simulations presented in this work.

Table 1. Parameters of the present cases and of those in Uhlmann & Doychev (2014).

3.3. Initialization

In order to initialize the flow around the particles we simulate an initial transient during a time interval of $66.67D/\langle w\rangle _{f}$ in which the particles are fixed. The initial position of the particles follows a random uniform distribution. Their orientation is randomly distributed with a maximum deviation of ${\pm }5^\circ$ tilting angle with respect to the vertical axis. The objective of constraining the angular position is to obtain a slight perturbation of the angular position with respect to the stable position of settling spheroids at moderate Ga. The same initial distribution of particles is used in both cases presented in table 1.

During the fixed-particle transient we impose the Reynolds number of the flow relative to the particles, $Re_{D}^{0}$, as obtained (as an output parameter) in the simulations with an isolated mobile particle at the target value of the respective Galileo number, cf. table 2. This is realized by means of a constant vertical pressure gradient, similar to the body force which counteracts the acceleration of the system when particles are freely mobile (Uhlmann & Doychev 2014). As will be shown later, the flow around the particles during the initial fixed-particle transient mostly resembles that of the analogous single-particle case at the given $Re_{D}^{0}$, due to the low concentration of particles in both cases. It should be mentioned that during this initial transient we obtain good decorrelation of all flow velocity components in all spatial directions (see Appendix B). Once the particles are released, the mean fluid flow relative to the particles is maintained through the above mentioned body force. Therefore, if perfectly adjusted, a single particle would not exhibit any vertical motion in the computational reference system, while the vertical motion of the many-particle ensemble is entirely due to mutual interactions (Uhlmann & Doychev 2014).

Table 2. Single-particle regime and time-averaged results of the present cases and those in the work of Uhlmann & Doychev (2014).

4. Results

In this section we present the results obtained for many-particle cases after particles are released. We also present additional simulations of settling particle pairs in order to analyse their ‘DKT’ dynamics.

4.1. Clustering phenomena in many-particle cases

When particles are released, they start to interact with their neighbours by repeated DKT events. The interaction between particles is intense in both cases G111 and G152, and rather similar (see animations in the supplementary material). In the following we present: (i) the enhanced settling and quantification of clustering, (ii) the angular motion of the particles, (iii) the arrangement of particles’ trajectories and (iv) a visualization of the main features of the flow.

4.1.1. Enhanced settling and quantification of clustering

Figure 4(a) shows the time history of the average particle settling velocity defined in (2.9). It can be seen that, in both present cases, the ensemble of mildly oblate spheroids reaches on average a much enhanced settling velocity magnitude after an initial transient of approximately $400\tau _{g}$, after which the values saturate and continue fluctuating around a value of roughly $-1.25$. This behaviour is quite in contrast to the known results for spheres (also included in the figure for reference) for which a transition from expected settling to enhanced settling occurs in the corresponding range of Galileo numbers, and for which the enhancement of the settling velocity was found to be less pronounced (only roughly half of the present increase). Instead, the present suspension of spheroids with $\chi =1.5$ appears to behave qualitatively similar to the much more flattened ($\chi =3$), almost neutrally buoyant ($\tilde {\rho }=1.02$) spheroids of Fornari et al. (2018), for which the magnitude of the mean collective settling velocity was found to increase by over $30$ %.

Figure 4. Time history of (a) enhancement of the settling velocity ($w_s$), and standard deviation of (c) settling and (d) horizontal velocity, normalized with the reference settling velocity from the single-particle counterpart ($w_{ref}$). (b) Temporal evolution of the standard deviation of Voronoï cell volumes ($\langle \tilde {V}_i^{\prime }\tilde {V}_{i}^{\prime }\rangle ^{1/2}$), normalized with the value obtained for a random Poisson process ($\langle \tilde {V}_i^{\prime }\tilde {V}_{i}^{\prime }\rangle _{rnd}^{1/2}$). Reference data for spheres are from Uhlmann & Doychev (2014).

For completeness, let us report that the amplitude of the particle velocity fluctuations (cf. figure 4c,d) follows a similar trend as the mean settling velocity, with a nearly linear growth during the initial transient and subsequent fluctuations around time-averaged values of approximately $w_s^{std} \approx 0.4w_{ref}$, $u_{h}^{std} \approx 0.2w_{ref}$. It should be noted that a clear signature of the oblique motion of individual particles just after their release is visible in case G152, as can be seen in terms of an initial peak of the intensity of the horizontal particle motion visible in the inset in figure 4(d). After approximately $10\tau _{g}$ this peak disappears, and the two suspensions of spheroids at different Galileo numbers exhibit very similar temporal evolutions.

In order to investigate the spatial structure of the disperse phase, we make use of the normalized Voronoï cell volume

(4.1)\begin{equation} \tilde{V}_i = \frac{V_v^{(i)}}{\langle V_v\rangle_p}, \end{equation}

where $V_{v}^{(i)}$ is the volume of the $i$th cell of the Voronoï diagram obtained from a three-dimensional tessellation of space based on the centroid locations of the particles (Monchaux, Bourgoin & Cartellier 2012; Uhlmann & Doychev 2014). In order to quantify the tendency to cluster we compare the standard deviation of $\tilde {V}_i$ with the same quantity obtained in a random Poisson process (RPP) of particles with the same shape and with the same solid volume fraction. For a mean concentration of particles such that $\varPhi =5\times 10^{-3}$, these RPP reference values are $\langle \tilde {V}_{i}^{\prime }\tilde {V}_{i}^{\prime }\rangle _{rnd}^{1/2}=0.4176$ for oblate spheroids of $\chi =1.5$ (determined in the framework of the present work) and $\langle \tilde {V}_{i}^{\prime }\tilde {V}_{i}^{\prime }\rangle _{rnd}^{1/2}=0.4146$ for spheres (given in Uhlmann 2020). The graph in figure 4(b) shows the temporal evolution of the standard deviation of the Voronoï cell volumes, normalized with the reference value from RPP. A direct qualitative correspondence with the temporal evolution of the (magnitude of the) mean settling velocity in figure 4(a) can be observed. More specifically, after a similar initial transient with an approximately linear growth of $\langle \tilde {V}_{i}^{\prime }\tilde {V}_{i}^{\prime }\rangle ^{1/2}$, the present suspensions of spheroids at both Galileo number values reach very large clustering intensities, which by far exceed those reported for spheres (cf. case M178 from Uhlmann & Doychev 2014). Based on the previous knowledge for settling spheres, and on the fact that the present spheroids with $\chi =1.5$ are not too far from a spherical shape, the strong clustering of the lower Galileo number case G111 is unexpected.

Next, we present time-averaged data of some of the time-dependent quantities discussed above. We discard the initial transient, and we start collecting statistics at $t=500\tau _{g}$, resulting in a sampling time interval of approximately $1000\tau _{g}$. Figure 5(a) shows the time-averaged values of $\langle \tilde {V}_{i}^{\prime }\tilde {V}_{i}^{\prime }\rangle ^{1/2}$ plotted vs the Galileo number. It can be clearly seen that the present suspensions of spheroids feature strong clustering with little effect of varying the value of the Galileo number. Interestingly, it turns out that the magnitude of the mean settling velocity normalized by the single-particle reference value is approximately proportional to the standard deviation of the Voronoï cell volumes normalized by its random value, i.e. to the clustering intensity. This relation is shown in figure 5(b), for the present spheroids and for the sphere suspensions of Uhlmann & Doychev (2014) and of Doychev (2014). This observation should be checked with the help of additional datasets in the future, since – if confirmed – it might open up a possibility to determine the mean settling velocity of a collective from knowledge on the spatial structure of the dispersed phase alone.

Figure 5. (a) Time-averaged values of the standard deviation of Voronoï cell volumes, $\langle \tilde {V}_{i}^{\prime }\tilde {V}_{i}^{\prime }\rangle _t^{1/2}$, normalized with the RPP reference values vs Ga. (b) Magnitude of the time-averaged mean settling velocity, $w_s$, vs $\langle \tilde {V}_{i}^{\prime }\tilde {V}_{i}^{\prime }\rangle _t^{1/2}$. Reference data for spheres from Uhlmann & Doychev (2014) and Doychev (2014).

The distinct spatial arrangement of the spheroids compared with spheres can be further characterized with the aid of the probability density functions (p.d.f.s) of the normalized Voronoï cell volumes $\tilde {V}$. Figure 6 shows the p.d.f. of $\tilde {V}$ of both cases with spheroids (G111 and G152) together with the p.d.f. of a RPP (Uhlmann 2020) and the p.d.f.s of two cases with spheres taken from Uhlmann & Doychev (2014) (their cases M121 and M178). The same qualitative behaviour is observed in both cases with spheroids: the probability of finding volumes smaller than the average ($\tilde {V}=1$) is high compared with a RPP. This indicates that most of the particles form part of clusters. Similarly, the probability of finding large volumes ($\tilde {V}\gtrsim 2$) is higher than that of a RPP, which implies the presence of particles in large void regions. Contrarily, case M121 from Uhlmann & Doychev (2014) shows a tendency to a more ordered state, where the probability of finding volumes similar to the average value is higher.

Figure 6. Probability density function of the normalized Voronoï cell volumes in (a) linear and (b) logarithmic scale. Reference data for spheres are from Uhlmann & Doychev (2014).

4.1.2. Angular velocity and orientation of particles

Figure 7(a) shows the time history of the standard deviations of the horizontal and vertical components of the angular velocity of the particles normalized with $w_{ref}/D$. Figure 7(b) shows the time history of the averaged and standard deviation of the tilting angle $\varphi _{v}$, which is defined as the angle between the symmetry axis of the spheroid and the vertical direction ($0\leq \varphi _{v}\leq 90^\circ$). The time evolution of these four quantities is remarkably different from the time evolution of the quantities reported in figure 4. Within a few gravitational time units after the particle release all quantities in figure 4 increase from zero to values close to their asymptotic time averages, after which they vary only very mildly. This indicates that the angular motion of the particles shows less sensitivity to collective effects than does the linear motion. Note that the converged values are slightly higher in case G152 when compared with case G111, except for $\langle \varphi _{v}^{\prime }\varphi _{v}^{\prime }\rangle _{p}^{1/2}$, that presents values which are approximately equal in both cases.

Figure 7. Time history of (a) standard deviation of the angular velocity and (b) average and standard deviation of the orientation angle $\varphi _{v}$.

Now let us focus on the p.d.f. of these angular quantities. Figure 8 shows the p.d.f. of the vertical and horizontal components of the angular velocity. We have tried several known p.d.f.s and we found Laplace's distribution

(4.2)\begin{equation} f\left(x;\beta\right) = \left(2\beta\right)^{{-}1}\exp\left(\left|x\right|/\beta\right), \end{equation}

with $\beta$ being a free parameter (cf. fitted numerical values in the figure), the best fit to both components of the angular velocity. This highlights the exponential tails, i.e. the importance of extreme events of the angular particle motion. The vertical component shows, however, an approximately $20\,\%$ smaller value for $\beta$ than the horizontal counterpart, indicating a more localized distribution with higher kurtosis. Similarly, we found the best fit of the estimated p.d.f. of the tilting angle $\varphi _{v}$ with respect to the vertical to be a gamma distribution

(4.3)\begin{equation} f\left(x;k,\theta\right) = \frac{x^{k-1}}{\theta^k\varGamma(k)}\exp\left({-}x/\theta\right) , \end{equation}

where $k$ and $\theta$ are the shape and scale parameters. Figure 9(a) shows the p.d.f. of $\varphi _{v}$ and the fitted Gamma distribution. The obtained fitting shows very good agreement in the range $0^\circ <\varphi _{v}\lesssim 30^\circ$. For higher values $\varphi _v>30^\circ$, the p.d.f. of each case shows a lower decay rate compared with the fitted gamma distribution for both cases. The shape parameter of the fitted gamma distribution in both cases ($k\approx 2.3$) indicates a fast, but non-abrupt approach to zero of the p.d.f. in the limit $\varphi _{v}\rightarrow 0^+$. Please recall that, from the shape parameter of the gamma distribution, we can infer the following: when $k\leq 1$ the maximum probability is located at $\varphi _{v}=0$ and then the p.d.f. decreases monotonically as $\varphi _{v}$ increases, when $k>1$ the limit of the p.d.f. as $\varphi _{v}\rightarrow 0^+$ is zero (with strong gradients in the vicinity of $\varphi _{v}=0$ for values of $k$ closer to unity), and when $k\geq 3$ the distribution shows a slow increase of the probability in the vicinity of zero (see figure 9b). Finally, the non-intuitive zero value of the p.d.f. in the limit $\varphi _{v}\rightarrow 0^+$ can be explained by the circumferential shape of the bins used to generate it. In order to obtain the p.d.f. for a specific value of ${\varphi _{v}}_0$, we define a bin by its edges $[{\varphi _{v}}_a,{\varphi _{v}}_b]$, where ${\varphi _{v}}_a < {\varphi _{v}}_0 < {\varphi _{v}}_b$. In the limit $\varphi _{v}\rightarrow 0^+$ and for increasing resolution of the p.d.f. (${\varphi _{v}}_b - {\varphi _{v}}_a\rightarrow 0$), the area on the surface of a sphere between the edges tends to zero, and as a consequence the probability of finding occurrences of $\varphi _{v}$ such that ${\varphi _{v}}_a < \varphi _{v} < {\varphi _{v}}_b$, also tends to zero.

Figure 8. Probability density functions of the (a) vertical and (b) horizontal components of the angular velocity. The curves are fitted to a Laplace distribution whose parameter $\beta$ is indicated in the legend. The Gaussian curve is shown for comparison purposes.

Figure 9. (a) The p.d.f. of the tilting angle $\varphi _{v}$ of cases G111 and G152 (the same information with the $y$ axis in logarithmic scale is shown in (b). A fitted gamma distribution is included (parameters from fitting included in the legend).

4.1.3. Trajectories

Now we turn our attention to the trajectories of the particles. Figures 10, 11 and 12 show the trajectories followed by the centre of gravity of the particles during a time span of $0.54\tau _{g}$ leading up to three successive time instants in two perpendicular projections. With this representation, considering the initialization procedure (§ 3.3) and if, hypothetically, collective effects were absent, particles would be represented by single points in the case G111 (single particle follows steady vertical trajectory) and by straight horizontal lines in the case G152 (single particle follows steady oblique trajectory). The time instants selected are $t/\tau _{g}\approx 10\,400$ and $1500$, which correspond to: (i) shortly after the particles are released, (ii) the end of the linear growth of $w_s$ and $\langle \tilde {V}_{i}^{\prime }\tilde {V}_{i}^{\prime }\rangle ^{1/2}$ (see figure 4) and (iii) the asymptotic state close to the end of the simulated time, respectively. Shortly after particles are released ($t/\tau _{g}\approx 10$), there is almost negligible motion of the particles in case G111 (figures 10(a) and 11(a) show point-like trajectories) and a small lateral motion in case G152 (figures 10(b) and 12(a) show trajectories with a noticeable horizontal component). As both cases evolve, the trajectories show elongated shapes along the vertical direction, indicating an enhancement of the settling velocity compared with the single-particle configuration. On average, the length of the trajectories in the vertical direction keeps growing during the constant acceleration phase ($t\lesssim 400\tau _{g}$) forming columnar clusters whose size is comparable to that of the computational domain (figures 10c and 10d, 11b and 12b). Again, the similarity of both cases is clearly noticeable. There is, however, a tendency of the clusters in the case of lower Galileo number to be more stable compared with the case of higher Galileo number when $t\approx 400\tau _{g}$ (see figures 11b and 12b). We believe that the smaller flow disturbances of the low Galileo number case allow the presence of clusters with a cross-section of the order of tens of $D$ (see figure 10c) to fill the entire domain in the vertical direction, whereas clusters of the same size in the horizontal direction of the higher Galileo number case (see figure 10d) are not stable and thus, appear to be more localized. This could explain the higher enhancement of the settling velocity with respect to the single particle case observed in the case at lower Galileo number (figure 4a) but it should be noted that this is only possible due to the periodic configuration. Therefore, the larger enhancement in $w_s$ observed in figure 4(a) for case G111 should be interpreted with care. In the horizontal direction, clusters show a continuous growth until they reach a size comparable to the computational domain (figure 10e, f).

Figure 10. Top view of trajectories of the particles’ centre of gravity during a time span of $[t_0-T_t,t_0]$, where $T_t=0.54\tau _{g}$ for cases G111 (a,c,e) G152 (b,d, f). Trajectories are coloured according to the particle's velocity relative to the mean velocity of the mixture (red downwards, blue upwards).

Figure 11. As in 10 but for case G111 only and viewed from the side.

Figure 12. As in 10 but for case G152 only and viewed from the side.

4.1.4. Flow visualization

Figure 13 shows visualizations of the cases G111 and G152 just before releasing the particles (panels a,b,e, f) and two converged states (panels c,d,g,h). In the figure the particles are represented together with isocontours of the second invariant (Q) of the velocity gradient tensor (Hunt et al. 1988) and isocontours of the filtered vertical velocity $\tilde {w}$. The filter applied to ${w}$ is a Gaussian filter of width $2.3D$ used to visualize large-scale velocity fluctuations. Just before particles are released ($t/\tau _{g}=0$) we observe that in case G111 ($Re_{D}^{0}=105$) the most common flow structure is a toroidal vortex around each particle (figure 13e), whereas in case G152 ($Re_{D}^{0}=158$) we also frequently observe a double-threaded wake (figure 13f). The similarity of these vortical structures with the analogous single-particle case at the given $Re_{D}^{0}$ is due to the low concentration of particles in both cases. After convergence has been reached ($t/\tau _{g}>1000$) both cases show large regions of high-speed downward flow, which are absent when particles are released ($t=0$). The strong clustering of both cases discussed above is clearly seen when comparing the initial and converged snapshots of both cases.

Figure 13. Visualization of isocontours of $QD^2/U_{g}^2=0.7$ (case G111) and $0.83$ (case G152) and $\tilde {w}=\langle w\rangle _f-0.5U_{g}$. Particles are represented in pink, isocontours of $Q$ with grey-coloured surfaces and isocontours of $\tilde {w}$ with yellow surfaces. (a–d) Show the whole domain ($Q$ and $\tilde {w}$), (e–h) show a part of the domain (only Q). First and second columns correspond to the instant before the release of the particles of cases G111 and G152, respectively. Third and fourth columns correspond to a converged state of each case.

4.2. Drafting–kissing–tumbling

Since animations show that interactions between particle pairs occur frequently in the many-particle cases, we now proceed to an analysis of the interaction of such pairs in isolation. As we will see, these interactions are quite sensitive to the particle shape, and we believe that pairwise interactions are the key to understanding the tendency of oblate spheroids to cluster at Galileo numbers for which spheres do not exhibit clustering. We restrict the analysis to one combination of Galileo number and density ratio ($Ga=111$, $\tilde {\rho }=1.5$) for which a single spheroid of aspect ratio $\chi =1.5$ and a single sphere result in a steady vertical regime (Jenny et al. 2004; Moriche et al. 2021). As mentioned above, many-particle cases in the dilute regime ($\varPhi = 5\times 10^{-3}$) with this combination of Ga and $\tilde {\rho }$ present strong clustering in the case of oblate spheroids with $\chi =1.5$ and no clustering in the case of spheres (Uhlmann & Doychev 2014). We include an additional series in this set of DKT cases, in which the angular motion is suppressed. Therefore, we have these four configurations:

1. (i) Free-to-rotate spheres (angular motion enabled).

2. (ii) Rotationally locked spheres (angular motion suppressed).

3. (iii) Free-to-rotate spheroids (angular motion enabled).

4. (iv) Rotationally locked spheroids (angular motion suppressed).

This is motivated by the study of Kajishima (2004) who demonstrated that suppressing the rotation of spheres leads to an enhanced clustering tendency. We perform a parametric sweep of $72$ initial relative particle positions for each configuration, resulting in a total number of $288$ simulations (see figure 14i). The problem set-up is analogous to the set-up of the many-particle cases presented in § 3, except that only two particles are present, and that we use an inflow/outflow configuration in the vertical direction as described in detail in Appendix C.

Figure 14. Trajectories of the trailing and leading particles for selected initial positions of spheres (a–d) and spheroids with $\chi =1.5$ (e–h), all with $\tilde {\rho }=1.5$ at $Ga=111$. The reference frame is translating downwards at a constant speed slightly smaller than the settling velocity of a single particle ($0.975w_{ref}$). Each panel contains the data of the cases with angular motion enabled and suppressed for a single initial condition and particle shape (see legend). (i) A sketch of the problem and the coordinates used is presented, in which the leading particle is represented with its actual shape and the trailing particle with a marker. The point markers correspond to all the initial conditions which we have computed, and the symbol markers to those initial conditions which are shown in panels a–h. (j) Sketch of the $x',y'$ coordinates.

Figure 14 shows the trajectories of $16$ selected cases projected onto a vertical plane intersecting the particles at their initial position. These cases correspond to different initial conditions of the four selected configurations. Please note that in these four cases particles do collide at least once, while this is not the case for all initial separations (as will be shown in figure 16 below). If we focus on the trajectories of the trailing particles, there is a clear similarity in free-to-rotate and rotationally locked spheroids: for these two configurations the trailing particle of the four initial conditions shown in the figure moves laterally until its centre is vertically aligned with the leading particle. Then, the trailing particle drifts towards the leading particle following an almost vertical path. On the contrary, spheres present a different path along which the trailing particle approaches the leading particle. The trailing particle of both free-to-rotate and rotationally locked spheres presented in the figure drifts towards the leading one following an oblique path for a significant time during the final approach. Interestingly, free-to-rotate spheres present a lateral shift (by roughly $D/4$) in their trajectory compared with their rotationally locked counterparts. This is clearly evident in figure 14(a), where the trailing particle of the free-to-rotate pair first moves away laterally from the leading particle, then makes its way in a straight, oblique path towards it.

In figure 15 we show the trajectories of the trailing particle relative to the centre of the leading particle for the same configurations as presented in figure 14. It can be seen how the rapid alignment along the vertical axis exhibited by the spheroidal particles results in the possibility of finding the trailing particle on either side of the vertical axis after the initial contact. Spheres, on the other hand, which do not fully align vertically, do not cross the vertical axis through the leading particle during the entire interaction. A very interesting result is the robustness of the lateral shift of the trailing particle for free-to-rotate spheres. Independently of the four initial conditions shown in the figure, the path followed by the trailing particle in the free-to-rotate sphere configuration converges to a single master curve. This feature has been observed for a number of other initial conditions. We attribute this lateral shift to a rotation-induced lift force resultant from the finite angular velocity (around the horizontal axis perpendicular to the plane shown in the figure) reached by these particles (see inset of figure 15a). The origin of the rotation of the particles ($\omega _{p,y'}<0$) is the shear seen by the trailing particle because of the deficit in vertical velocity in the wake of the leading particle. A detailed analysis of the forces acting on a particle in motion in the wake of another particle, however, is outside the scope of the present contribution, and it is left as a worthwhile topic for future studies.

Figure 15. Trajectory of the trailing particle relative to the leading one for (a) spheres and (b) spheroids. The close-up trajectories shown in the insets in (a,b) are coloured with the angular velocity perpendicular to the plane shown (see legend). Line colour and marker type follow the same convention as in figure 14.

Figure 16. Maps of (a,b) time to first collision $t_{cI}$ and (c,d) interaction time of the DKT cases as a function of the initial condition of the trailing particle. The $x'$ axis for the rotationally locked cases is flipped to facilitate the comparison. Cases in which no interaction occurred in the evaluated time are represented with black dots and interacting cases are represented with coloured markers. The red line in panels (a,b) is an isocontour of $t_{cI}=100\tau _{g}$. Panels (e, f) contain the ratio of $t_{cI}$ of free-to-rotate cases with respect to their rotationally locked counterparts ($t_{cI}^{FTR}/t_{cI}^{RL}$).

Next we consider on the complete set of initial conditions of all four configurations. Figure 16 shows maps of the time to first collision, $t_{cI}$, and the interaction time, $t_{cR}$, for the cases showing at least one DKT event. The interaction time, $t_{cR}$, measures the time between the first collision and the smallest time after which the distance between the particles grows monotonically (see precise definitions in Appendix C), evaluated as a function of the initial position of the trailing particle. The curve which delimits those initial conditions that lead to particle contact from those which do not is very similar in all four cases (figure 16a,b). This means that the region of attraction in the explored range of relative positions is essentially insensitive to the precise particle shape (sphere vs mildly oblate spheroid) and to their ability to rotate. Both spheres and the spheroids considered, independently of angular motion being suppressed or not, present a roughly cylindrical region of radius approximately $3D$ located downstream of the leading particle in which particles will interact. Now let us consider the time to first collision, $t_{cI}$, which can give us insight into the particles’ response to wake attraction mechanisms in an integral sense, as a function of the initial relative position of the particles (figure 16a,b). We find similar values of $t_{cI}$ for free-to-rotate spheroids and rotationally locked spheroids and spheres, whereas free-to-rotate spheres present somewhat larger values for the same initial conditions. This result suggests that the angular motion of spheres causes these particle pairs to approach more slowly than corresponding non-rotating spheres and oblate spheroids. The latter do not rotate continuously, instead the angular motion during the approach phase is characterized by small oscillations around the equilibrium position (with the symmetry axis vertically aligned). For the sake of clarity we include the ratio of $t_{cI}$ of free-to-rotate particles with respect to their rotationally locked counterparts in the auxiliary panels (e, f) for each initial condition. It can be seen that free-to-rotate spheres with a small horizontal shift in their relative position show ratios as large as $2.5$. Furthermore, this ratio increases with the vertical distance within the range of parameters evaluated. These differences disappear for the cases whose initial relative position has almost no horizontal shift, indicating that a larger horizontal shift is required to trigger the rotation of the trailing sphere. Next let us focus on the interaction of particle pairs after the first collision. Only a maximum of one collision is observed in the entire datasets for rotationally locked spheres and spheroids and free-to-rotate spheres. Free-to-rotate spheroids, however, present two or three collision events per parametric point for most of the chosen initial conditions, except for a few of them in which a single, or even four collisions occur. Figure 16(c,d) shows the interaction time $t_{cR}$, where a clear ascending order of configurations with respect to the duration of particle-pair interactions is identified: (i) free-to-rotate spheres present an almost negligible interaction time, (ii) rotationally locked spheres interact during approximately $5\tau _{g}$, the value of $t_{cR}$ increases to approximately $20\unicode{x2013}40\tau _{g}$ for the (iii) rotationally locked spheroids and it reaches values above $100\tau _{g}$ for the (iv) free-to-rotate spheroids.

In the following we will analyse the temporal evolution of the distance separating the two interacting particles after the first contact. For this purpose we set the time of the first contact arbitrarily to zero, and then average the data over the ensemble of realizations with different initial separations. We define the average post-collisional distance as follows:

(4.4)\begin{equation} L\left(\tilde{t}\right)=\frac{\displaystyle\sum_i^{N} \lVert \boldsymbol{x}^i_{r,trail}\left( \tilde{t}\right)\rVert}{N} , \end{equation}

where $\boldsymbol {x}^i_{r,trail}$ represents, for each case $i$, the relative position of the trailing particle with respect to the leading particle and $\tilde {t}$ measures the time after the first collision (recall from figure 16(c,d) that the value of $t_{cI}$ is case dependent). Figure 17 shows the average post-collisional distance, $L(\tilde {t})$, for the four investigated configurations. First, let us compare rotationally locked spheres vs rotationally locked spheroids. On average, rotationally locked spheres drift away from each other a distance of approximately $2.5D$ in the first $10\tau _{g}$ after the initial contact (figure 17b). The scenario for rotationally locked spheroids is, on average, slightly more complex. There is a local maximum of the inter-particle distance at $\tilde {t}\approx 20\tau _{g}$, after which particles start to approximate each other again, and a local minimum of $L$ at $\tilde {t}\approx 30\tau _{g}$, after which particles drift away from each other. Second, let us focus on the effect that angular motion has on both particle shapes considered. For the spheres, the angular motion results in particles staying much further away from each other after the collision, with an approximately constant offset of around $1D$ for $\tilde {t} > 10\tau _{g}$. This considerable effect due to the spheres’ rotation can be understood from figure 18, which depicts the relative motion of particle pairs for exemplary cases, as will be discussed in more detail below. Conversely, when angular motion is allowed for spheroids, particles remain on average closer to each other at times larger than approximately $10\tau _{g}$ after the first collision, and the local maximum of $L$ observed for the rotationally locked counterparts is more pronounced.

Figure 17. (a) Time history of average distance (4.4) after the first contact for the four configurations considered in the DKT configuration (see legend). Non-colliding cases are excluded from the plot. (b) Zoom of panel (a) (see dashed rectangle in (a)).

Figure 18. Overlay of consecutive snapshots of the different DKT configurations for the cases whose trailing particle starts at $\boldsymbol {x}_r=(0.625,7.5)D$. The initial and final snapshots are indicated by highlighting the particles’ contour. The time interval selected is such that it starts at the first contact ($t=t_{cI}$) and ends after $16\tau _{g}$ for spheres and $32\tau _{g}$ for spheroids, sampling $8$ equispaced time instants. The time between consecutive snapshots $\Delta t$ is indicated in the figure. The reference frame is translating downwards at a speed slightly smaller than the settling velocity of a single particle ($0.975w_{ref}$). Particles are identified by colour (trailing: green, leading: purple). Time and angular position are indicated with a small mark whose colour changes with time.

Let us now identify possible phenomena that lead to the strong differences in the pairwise interaction times after the first collision. First, the larger separation observed in spheres as compared with spheroids just after the first collision may be attributed to the mechanism by which pairs of particles start to tumble after they contact each other. According to Fortes et al. (1987), vertically aligned spheres which are in contact form an unstable system which is the actual cause of the tumbling phase. This instability is due to the fact that two spheres vertically aligned can be seen as an elongated body. Such body would tilt so that it settles maximizing drag. Since the pairs of particles considered are not connected, they separate from each other shortly after they start to tilt. The ratio between height and width of the body formed by two vertically aligned and touching spheres is equal to 2, while the same quantity for spheroids (with their axis of symmetry aligned with the vertical) equals $2/\chi$, which amounts to $4/3$ in the present case. Therefore, we can expect a more abrupt initial tumbling motion in the case of spheres. Please note that if the oblate spheroids considered had an aspect ratio of $\chi =3$, as considered by Fornari et al. (2018), the ratio between height and width becomes $2/3$, which is smaller than unity. In such a situation the tandem of vertically aligned particles would be even more stable. Indeed, Fornari et al. (2018) reported a stick mechanism in which the oblate spheroids stay together, suppressing the tumbling phase. Regarding the angular motion, two different effects are considered as candidates to explain the lower values of the particle-pair separation distance $L$ for spheroids and higher values of $L$ for spheres, when compared with their rotationally locked counterparts. In the case of spheroids the angular motion allows the particles to have a stronger rocking motion around a horizontal axis after their first collision. This rocking motion is a consequence of the leading particle tilting when pushed downwards by the trailing particle. Figure 18(d) shows the tilting of the leading particle and the consequent rocking motion after the first collision. This results in a zig-zag trajectory and, as a consequence, particles stay, on average, closer to each other, thereby increasing the probability for repeated collisions. The subsequent collisions show similar features, but the amplitude of their rocking motion is decreased, and the distance between the particles after the collisions is increased. In the case of spheres, the distance between centres of free-to-rotate spheres is approximately $1.5$ times larger than that of the rotationally locked counterparts. We attribute this to a rotation-induced lift force on the trailing particle that increases the distance between the two particles. This rotation can be appreciated in figure 18(b), and its effect on the relative trajectory between the particles is visible in figure 15(a).

5. Discussion

On the one hand, our present many-particle simulations show that mildly oblate spheroids (with an aspect ratio of $1.5$) form strong clusters at a comparably low value of the Galileo number, for which an isolated particle in ambient fluid settles in the steady vertical regime. This observation is in contrast to the known behaviour of ensembles of spheres, for which the Galileo number has been identified as a critical parameter with respect to the onset of wake-induced clustering (Uhlmann & Doychev 2014). This clustering transition in the case of spheres occurs in the range $Ga=121\ldots 178$ which encloses the value $Ga=155.8$ at which an isolated particle's path regime bifurcates from steady vertical to steady oblique (Zhou & Dušek 2015). The explanation proposed by Uhlmann & Doychev (2014) links the onset of clustering in the case of spheres to the enhanced horizontal mobility of individual particles above the oblique-path threshold, which causes an increase of the frequency of particle–particle encounters. In order to make this point clearer, let us consider the fact that, in the wake of a given particle, a ‘region of attraction’ can be defined, i.e. a spatial region of limited extent within which a trailing particle near equilibrium will experience a modified hydrodynamic force which has the effect that – without further perturbations – it will continuously approach the leading particle until contact. Now let us consider the hypothetical case of a set of mono-disperse particles which are initially all located outside of any other particles’ region of attraction. If these particles are settling in the steady vertical regime in an unbounded ambient fluid, they simply move as a group without changing their relative positions, and they will, therefore, never come into contact. If, however, each individual particle of this ensemble settles in a steady oblique regime (with a random azimuthal angle), there exists a finite probability for individual particles to encounter another particle's region of attraction, and to come into direct contact. If the conditions are such that an interacting particle pair (on average) encounters one (or more) additional particle(s) before separating again, initial seeds will eventually lead to large-scale cluster formation through accretion. Hence, we can conclude that horizontal particle mobility can have the effect of enhancing the tendency to cluster in this scenario. However, as we have observed in the present many-particle simulations, oblate spheroids with aspect ratio $\chi =1.5$ apparently do not require this mechanism in order to form columnar clusters. This is in accordance with the results of Fornari et al. (2018) for oblate spheroids at larger aspect ratio ($\chi =3$), lower Galileo number ($Ga=60$) and smaller density ratio ($\tilde {\rho }=1.02$) which likewise appear to form clusters.

On the other hand, we have seen in the present series of DKT simulations that pairs of mildly oblate spheroids, which are initially positioned such that the trailing particle is located inside of the leading particle's region of attraction, interact for a significantly longer time than spherical counterparts at corresponding parameter points. We have observed that the difference in interaction time is caused by two factors: first, spheroids approach each other at a somewhat faster rate than spheres during the drafting phase; second, spheroids remain in close proximity over a significantly longer duration after the first contact. This is qualitatively in line with the results of Ardekani et al. (2016) in their DKT simulation of spheroids with larger aspect ratio ($\chi =3$), smaller Galileo number ($Ga=80$) and smaller density ratio ($\tilde {\rho }=1.14$). These authors observe that those flatter objects do not tumble, but instead remain locked in position as a stack after the initial contact. Translating our result on the isolated pairwise interaction to the (dilute) many-particle configuration implies that the observed increase in the interaction time between pairs of particles can also enhance the tendency to form large-scale clusters. This is due to the fact that once a pair of particles gets into contact, it will have a higher probability of attracting an additional particle before separating, if the original pair has an extended contact duration. Hence, this argument suggests that an increase in the interaction time between pairs of particles can lead to large-scale clustering, which is a mechanism that is different from (and complementary to) the above process by enhancement of lateral mobility.

Figure 19 shows a sketch which is intended to contrast the above two routes of cluster formation in the dilute regime for sphere-like particles. Starting from the (non-clustering) baseline case with spherical particles at a Galileo number of ${O}(100)$ in the steady vertical regime (lower left corner of the $(Ga,\chi )$-plane in figure 19), clustering can be enabled by either one of the following two options: (i) by increasing the probability of particles entering their peers’ attractive wake region (which is achieved by increasing the horizontal particle mobility, i.e. through increasing the spheres’ Galileo number Ga); (ii) by increasing the temporal interval over which particles remain close to each other during wake-induced encounters (which can be achieved by replacing spheres with oblate spheroids, i.e. by way of increasing the particles’ aspect ratio $\chi$). Please note that additional triggers of cluster formation (besides lateral particle mobility and close interaction duration) might be at play, such as the growth of the spatial extent of the wake with increasing Galileo number. In addition, it remains to be understood how the density ratio influences the mechanisms mentioned above. Further research would be necessary to properly answer these questions.

Figure 19. Summary sketch of clustering mechanisms analysed in this work (intense DKT interactions) and from the reference work of Uhlmann & Doychev (2014) (promoted particle encounters by horizontal motion).

6. Conclusions

We have performed PR-DNS of many heavy non-spherical particles settling under gravity. The particles are oblate spheroids of aspect ratio $1.5$ (which represent a modest deviation from a spherical shape) and density ratio $\tilde {\rho }=1.5$; the global solid volume fraction is $0.005$ such that the suspension can be considered as dilute. Two Galileo numbers are considered, namely $111$ and $152$ for which a single oblate spheroid follows a steady vertical and a steady oblique path, respectively. In contrast to previous results for spheres (Uhlmann & Doychev 2014) we have found that the qualitative difference in the single particle regime does not result in a qualitatively different behaviour of the multiparticle cases: at both Galileo number values a strongly inhomogeneous spatial distribution of the disperse phase in the form of columnar clusters is observed, with a significantly enhanced average settling velocity as a consequence. A similar result has previously been reported for PR-DNS of significantly flatter spheroids (with $\chi =3$) and for lower density ratio ($\tilde {\rho }=1.02$) by Fornari et al. (2018). Here, we have used Voronoï tessellation as a basis for the analysis of the structure of the particulate phase. The intensity of clustering has been measured with the aid of the standard deviation of the Voronoï cells’ volume, normalized with the value obtained from a RPP (Monchaux, Bourgoin & Cartellier 2010). It turns out that the enhancement of the average settling speed is approximately proportional to the standard deviation of the Voronoï cell volumes when considering both the present spheroids as well as the spheres of Uhlmann & Doychev (2014) as a joint dataset. This result, however, requires further confirmation through additional data points before its potential implications can be evaluated. Note that the amount of enhancement of the settling velocity may still depend on the domain size, since the size of particle clusters approaches the former.

Motivated by the lack of influence of the single-particle regime upon the statistical features of the multi-particle settling we have carried out a thorough analysis of pairwise interactions of particles in the well-known DKT set-up, conducted in a computational domain with inflow/outflow boundary conditions in the vertical direction. We have considered four configurations, namely oblate spheroids of aspect ratio $1.5$ and spheres, with and without suppression of the angular motion, with density ratio $\tilde {\rho }=1.5$ and a Galileo number such that a single particle would follow a steady vertical path ($Ga=111$). Through systematic variation of the particle pair's relative initial position we have found that the region of attraction for both particle shapes, with and without rotation, is very similar. However, in the case of free-to-rotate spheres the trailing particle's trajectory is horizontally shifted towards larger radial distances, resulting in a prolonged drafting phase. Regarding the particles’ tumbling phase, we have shown that spheres and spheroids behave in a qualitatively different manner. Spheres undergo at most a single collision, and they quickly separate afterwards. Rotationally locked spheroids also experience a maximum of one collision, but they remain close to each other for relatively long times. Finally, free-to-rotate spheroids exhibit two collision events in most of the cases in which a DKT event is observed, and, consequently, their average interaction time is the maximum out of the four investigated configurations.

To summarize, we observe a shape-induced increase in the interaction time when two particles happen to ‘meet’ (i.e. when one of them enters the other particle's wake region), such that the probability of additional particles joining the initial pair is increased with respect to the baseline case of spheres. Hence, the tendency to form a large-scale cluster increases. This is in contrast to the mechanism for spheres (as proposed by Uhlmann & Doychev 2014), where the clustering transition is believed to be triggered by the primary bifurcation of the isolated particle's wake flow (from axisymmetric vertical to planar oblique) that leads to lateral mobility, hence increasing the probability of mutual particle encounters. As a consequence of these two observations, we conclude that the mechanism for the initiation of columnar clusters in the case of a dilute suspension of modestly oblate spheroids is in a sense orthogonal to the mechanism that is believed to be at work in the counterpart with spherical particles.

Since the qualitatively different behaviour of spheroids (as compared with spheres) has now been established both for moderately flat geometries ($\chi =3$, Fornari et al. 2018) and for modestly flat ones ($\chi =1.5$, present work), the question of a possibly finite critical aspect ratio poses itself naturally. Future work should be aimed in this direction. Furthermore, a thorough analysis of collective effects in prolate spheroids in a comparable parameter range is still lacking.