3.1. Particle statistics

We start by comparing the single-point flow and particle velocity statistics for the two cases studied, i.e. quiescent and turbulent flow. The results are collected when a statistically steady state is reached. Due to the axial symmetry with respect to the direction of gravity, we consider only two velocity components for both phases, the components parallel and perpendicular to gravity, $V_{{\it\alpha},{\it\tau}}$ and $V_{{\it\alpha},n}$ respectively, where ${\it\alpha}=f,p$ indicates the solid and fluid phases.

In figure 3 we report the p.d.f.s of the particle velocities for both volume fractions investigated here, ${\it\phi}=0.5\,\%$ and 1 %; the moments extracted from these distributions are summarized in table 2. The data in figure 3(a) show the p.d.f. for the component of the velocity aligned with gravity, $V_{p,{\it\tau}}$ , normalized by the settling velocity for ${\it\phi}\rightarrow 0$ in a quiescent environment, $v_{t}$ . This is extracted from the simulation of a very dilute suspension discussed in § 2.3.

In the quiescent cases, the mean settling velocity slightly reduces on increasing the volume fraction ${\it\phi}$ , in agreement with the findings of Richardson & Zaki (Reference Richardson and Zaki1954) and Di Felice (Reference Di Felice1999), among others. The sedimentation velocity decreases to 0.96 at ${\it\phi}=0.5\,\%$ and 0.93 at ${\it\phi}=1\,\%$ . Di Felice (Reference Di Felice1999) investigated experimentally the settling velocity of dilute suspensions of spheres ( ${\it\phi}=0.5\,\%$ ) with density ratio 1.2 in quiescent fluid, for a large range of terminal Reynolds numbers (from 0.01 to 1000). Following the empirical fit proposed in Di Felice (Reference Di Felice1999), we obtain $\langle V_{f,{\it\tau}}\rangle \simeq 0.98$ at ${\it\phi}=0.5\,\%$ , approximately 1 % larger than our result. On the other hand, the formula suggested for the intermediate regime in Di Felice (Reference Di Felice1999) and Yin & Koch (Reference Yin and Koch2007), $Re_{t}<150$ , would give an estimated value of 0.88, 6 % lower than our result.

Figure 3.

The p.d.f.s of the particle velocities along the directions (a) parallel, $V_{p,{\it\tau}}$ , and (b) perpendicular, $V_{p,n}$ , to gravity for ${\it\phi}=0.5\,\%$ and 1 %. The velocities are normalized by $v_{t}$ . Quiescent fluid cases are denoted by a blue curve with squares for ${\it\phi}=0.5\,\%$ and a green curve with triangles for ${\it\phi}=1\,\%$ , and turbulent flow data by a red curve with circles for ${\it\phi}=0.5\,\%$ and a black curve with downward-pointing triangles for ${\it\phi}=1\,\%$ . In the insets we show the p.d.f.s in lin–log scale.

Table 2.

First four central moments of the p.d.f.s of $V_{p,{\it\tau}}$ and $V_{p,n}$ normalized by $v_{t}$ . Here, $S_{V_{p,{\it\tau}}}$ ( $S_{V_{p,n}}$ ) and $K_{V_{p,{\it\tau}}}$ ( $K_{V_{p,n}}$ ) are respectively the skewness and the flatness of the p.d.f.s.

The mean settling velocity for the quiescent case at ${\it\phi}=0.5\,\%$ is instead close to the estimated value from (2.9). This is in agreement with what was reported in Uhlmann & Doychev (Reference Uhlmann and Doychev2014). These authors found that the particle mean settling velocity of a dilute suspension increases above the reference value only when the Archimedes number $Ar$ is larger than approximately 24 000 and clustering occurs. In our case, $Ar\simeq 21\,000$ and no clustering is noticed, according with their findings.

Interestingly, we observe an additional non-negligible reduction when a turbulent background flow is considered, in our opinion the main result of this paper. The reduction of the mean settling velocity $\langle V_{p,{\it\tau}}\rangle$ is 12 % at ${\it\phi}=0.5\,\%$ and 14 % at ${\it\phi}=1\,\%$ , see table 2. This result unequivocally shows that the turbulence reduces the settling velocity of a suspension of finite-size buoyant particles, in agreement with the experimental findings by Variano and co-workers (E. A. Variano, private communication) and Byron (Reference Byron2015). We also note that the reduction of the settling velocity with the volume fraction is less important for the turbulent cases.

Looking more carefully at the p.d.f.s, we note that fluctuations are, as expected, larger in a turbulent environment. In addition, the vertical particle velocity fluctuations are the largest component in a quiescent fluid, whereas in a turbulent flow the fluctuations are largest in the horizontal directions, as summarized in table 2. In the quiescent case, the root mean square (r.m.s.) of the tangential velocity ${\it\sigma}_{V_{p,{\it\tau}}}$ is 0.15 and 0.17 for ${\it\phi}=0.5\,\%$ and ${\it\phi}=1\,\%$ respectively, while it increases up to 0.23 in the corresponding turbulent cases. The difference in the width of the p.d.f. is particularly large in the directions normal to gravity, where the r.m.s. of the variance ${\it\sigma}_{V_{p,n}}$ is 0.3 for both turbulent cases, while it is 5 and 4 times smaller for the quiescent flows at ${\it\phi}=0.5\,\%$ and ${\it\phi}=1\,\%$ . We believe that the interactions among the particle wakes, mainly occurring in the settling direction, promote the higher vertical velocity fluctuations found in the quiescent cases. The shape of the p.d.f. is essentially Gaussian for the turbulent cases, showing vanishing skewness and normal flatness. Interestingly, an intermittent and skewed behaviour is exhibited in the quiescent cases. The flatness $K$ is approximately $9$ for both components at ${\it\phi}=0.5\,\%$ and slightly reduces to 5.5 at ${\it\phi}=1\,\%$ . The settling velocity of the quiescent cases is skewed towards intense fluctuations in the direction of the gravity. The skewness $S$ is higher for the more dilute case, being 1.26 at ${\it\phi}=0.5\,\%$ and 0.7 at ${\it\phi}=1\,\%$ .

We interpret the intermittent behaviour suggested by values of $K>3$ by the collective dynamics of the particle suspension. The significant tails of the p.d.f.s shown in figure 3(a) are indeed associated with a specific behaviour: as particles fall, they tend to be accelerated by the wakes of other particles, before showing drafting–kissing–tumbling behaviour (Fortes, Joseph & Lundgren Reference Fortes, Joseph and Lundgren1987). Snapshots of the drafting–kissing–tumbling behaviour between two spherical particles are shown in figure 4. When this close interaction occurs, particles are found to fall with velocities that can be more than twice the mean settling velocity $\langle V_{p,{\it\tau}}\rangle$ . In the quiescent case, the fluid is still, the wakes are the only perturbation present in the field and are long and intense, so their effect can be felt far away from the reference particle. The more dilute the suspension is the more intermittent the particle velocities are. On the contrary, when the flow is turbulent, the wakes are disrupted quickly and therefore fewer particles feel the presence of a wake. The velocity disturbances are now mainly due to turbulent eddies of different size that interact with the particles to increase the variance of the velocity homogeneously along all directions, leading to the almost perfect normal distributions shown above, with variances similar to those of the turbulent fluctuations. The features of the particle wakes will be further discussed in this paper to support the present explanation.

We also note that the sedimenting speed in the quiescent fluid is determined by two opposite contributions: the excluded volume effect, which contributes to a reduction of the mean settling velocity with respect to an isolated sphere, and the pairwise interactions (the drafting–kissing–tumbling), which increase the mean velocity of the two particles involved in the encounter. To try to identify the importance of the drafting–kissing–tumbling effect, we fit the left part (where no intermittent behaviour is found) of the p.d.f. pertaining to the quiescent case at ${\it\phi}=0.5\,\%$ with a Gaussian function. The mean of $V_{p,{\it\tau}}$ is reduced to approximately 0.93 (value due only to the hindrance effect) instead of 0.96 in the full simulation; thereby the increment in mean settling velocity due to drafting–kissing–tumbling can be estimated to be approximately 3 %.

Figure 4.

Drafting–kissing–tumbling behaviour among two spherical particles in the quiescent case with ${\it\phi}=0.5\,\%$ , at a) $t=1789$ ; (b) $t=1794$ ; (c) $t=1799$ ; (d) $t=1801$ . The particles are coloured with the absolute value of their velocity component in the direction of gravity. Three particles are labelled with ‘A, B, C’ in order to show how accelerated the two interacting particles are compared with the others.

Figure 5.

The p.d.f.s of (a) $|{\it\omega}_{p,z}|$ and (b) $\sqrt{{\it\omega}_{p,x}^{2}+{\it\omega}_{p,y}^{2}}$ for ${\it\phi}=0.5\,\%$ and 1 %. The angular velocities are normalized by $v_{t}/(2a)$ , the settling velocity of a single particle in a quiescent environment and its diameter. Quiescent fluid cases are denoted by a blue curve with squares for ${\it\phi}=0.5\,\%$ and a green curve with triangles for ${\it\phi}=1\,\%$ , and turbulent flow data by a red curve with circles for ${\it\phi}=0.5\,\%$ and a black curve with downward-pointing triangles for ${\it\phi}=1\,\%$ .

The p.d.f.s of the particle angular velocities are also different in quiescent and turbulent flows. These are shown in figure 5(a,b) for rotations about axes parallel to gravity, $|{\it\omega}_{p,{\it\tau}}|=|{\it\omega}_{p,z}|$ , and orthogonal to it, $\sqrt{{\it\omega}_{x}^{2}+{\it\omega}_{y}^{2}}$ . In the settling direction, the peak of the p.d.f. is always at $|{\it\omega}_{p,z}|=0$ . For the translational velocities, the p.d.f.s are broader in the turbulent cases. Due to the interaction with turbulent eddies, particles tend also to spin faster around axes perpendicular to gravity. From figure 5(b) we see that the modal value increases slightly in the quiescent cases on increasing the volume fraction. In the turbulent cases, the modal values are more than three times the values of the quiescent cases and the variance is also increased. Unlike the quiescent cases, the curves almost perfectly overlap for the two different values of ${\it\phi}$ , meaning that turbulent fluctuations dominate the particle dynamics. Turbulence hinders particle hydrodynamic interactions.

Figure 6 shows the temporal correlations of the particle velocity fluctuations,

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle R_{v_{{\it\tau}}v_{{\it\tau}}}({\rm\Delta}t)=\frac{\langle V_{p,{\it\tau}}^{\prime }(p,t)V_{p,{\it\tau}}^{\prime }(p,t+{\rm\Delta}t)\rangle }{{\it\sigma}_{V_{p,{\it\tau}}}^{2}}, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle R_{v_{n}v_{n}}({\rm\Delta}t)=\frac{\langle V_{p,n}^{\prime }(p,t)V_{p,n}^{\prime }(p,t+{\rm\Delta}t)\rangle }{{\it\sigma}_{V_{p,n}}^{2}}, & \displaystyle\end{eqnarray}$$

${\it\phi}=0.5\,\%$

${\it\phi}=1\,\%$

${\rm\Delta}t\sim 50$

$R_{v_{n}v_{n}}$

Figure 6.

Time correlations of the particle velocity fluctuations. The blue curve with circles represents the correlation of $V_{p,{\it\tau}}^{\prime }$ ( $R_{v_{{\it\tau}}v_{{\it\tau}}}$ ), while the red curve with diamonds is used for the correlation of $V_{p,n}^{\prime }$ ( $R_{v_{n}v_{n}}$ ). (a) Quiescent case with ${\it\phi}=0.5\,\%$ , (b) turbulent case with ${\it\phi}=0.5\,\%$ , (c) quiescent case with ${\it\phi}=1.0\,\%$ and (d) turbulent case with ${\it\phi}=1.0\,\%$ .

The normal velocity correlation $R_{v_{n}v_{n}}$ of the turbulent case oscillates around zero before vanishing at longer times. We attribute this to the presence of the large-scale turbulent eddies. As a settling particle encounters sufficiently strong and large eddies, its trajectory is swept on planes normal to gravity in an oscillatory way. To provide an approximate estimate of this effect, we consider as a first approximation the turbulent flow seen by the particles as frozen, since the particles fall at a higher velocity than the turbulent fluctuations (3.4 times). Since the strongest eddies are of the order of the transverse integral scale we can presume that these structures are responsible for this behaviour. In particular, the transverse integral scale $L_{T}=L_{0}/2\simeq 8$ and $u_{rms}^{\prime }=0.3$ , so we expect a typical period of $t=L_{T}/u_{rms}^{\prime }\simeq 26$ , which is of the order of the oscillations found for both turbulent cases, i.e.  $t\approx 20$ . It should be noted that a similar behaviour has been observed by Wang & Maxey (Reference Wang and Maxey1993) for sufficiently small and heavy particles, termed the preferential sweeping phenomenon.

The same process can be interpreted in terms of crossing trajectories and continuity effects, as described by Csanady (Reference Csanady1963). An inertial particle falling in a turbulent environment changes its fluid–particle neighbourhood continuously. It will fall out from the eddy where it was at an earlier instant and will therefore rapidly decorrelate from the flow. In order to accommodate the back flow necessary to satisfy continuity, the normal correlations must then contain negative loops (like those seen in figure 6 b,d). Following Csanady (Reference Csanady1963), we define the period of oscillation of the fluctuations as the ratio of the typical eddy diameter in the direction of gravity (i.e. the longitudinal integral scale $L_{0}$ ) and the particle terminal velocity $v_{t}$ , obtaining $t=L_{0}/v_{t}\simeq 16$ . This value is similar to the period of oscillations in the correlations in figure 6.

Figure 7.

Mean square particle displacement in the directions parallel and perpendicular to gravity for both quiescent and turbulent cases for (a) ${\it\phi}=0.5\,\%$ and (b) ${\it\phi}=1\,\%$ . The mean particle displacement $\langle V_{p,{\it\tau}}\rangle t$ in the settling direction is subtracted from the instantaneous displacement when computing the statistics. In the inset, we show $\langle {\rm\Delta}x_{i}^{2}\rangle /t$ as a function of time as well as the diffusion coefficients $D_{e}$ for the turbulent cases.

As shown in figure 6(c), the quiescent environment presents a peculiar behaviour of the settling velocity correlation at ${\it\phi}=1\,\%$ . In particular, we observe oscillations around zero of long period, $T=160{-}180$ . From the analysis of particle snapshots at different times (not shown), we observe that these seem to be correlated to the formation of regions of different density of the particle concentration. Hence, a particle crossing regions with different local particle concentrations may experience a varying settling velocity.

To further understand the particle dynamics, we display in figure 7(a,b) the single-particle dispersion, i.e. the mean square displacement, for both quiescent and turbulent cases at ${\it\phi}=0.5\,\%$ and 1.0 %. The mean displacement, $\langle V_{p,{\it\tau}}\rangle t$ , is subtracted from the instantaneous displacement in the settling direction, ${\rm\Delta}z(t)$ , to highlight the fluctuations with respect to the mean motion. For all cases we found initially a quadratic scaling in time ( $\langle {\rm\Delta}x_{i}^{2}\rangle \sim t^{2}$ ) typical of correlated motions, while the linear diffusive behaviour takes over at longer times ( $\langle {\rm\Delta}x_{i}^{2}\rangle \sim 2D_{e}t$ , with $D_{e}$ the diffusion coefficient).

The turbulent cases show a similar behaviour for both volume fractions. The crossover times when the initial quadratic scaling is lost and linear scaling occurs are approximately $t\simeq 10$ and $t\simeq 50$ for the normal and tangential components respectively. This difference is consistent with the correlation time scales previously discussed. The dispersion rates are similar in all directions in a turbulent environment. The quiescent cases present different features. First of all, dispersion is much more effective in the settling direction than in the normal one. The dispersion rate is smaller than in the turbulent cases in the horizontal directions, while, surprisingly, the mean square displacement in the settling direction is similar to that of the turbulent cases, being even higher at ${\it\phi}=0.5\,\%$ , something we relate to the drafting–kissing–tumbling behaviour discussed above. The crossover time scale is similar to that of turbulent cases, with the exception of the most dilute case which does not reach a fully diffusive behaviour at $t\simeq 500$ . This long correlation time makes the mean square displacement of this case higher than for the corresponding turbulent case at long times.

In the insets of figure 7(a,b) we show $\langle {\rm\Delta}x_{i}^{2}\rangle /t$ as a function of time. In all cases except the quiescent one at ${\it\phi}=0.5\,\%$ , the diffusive regime is reached and it is possible to calculate the diffusion coefficients $D_{e}=\langle {\rm\Delta}x_{i}^{2}\rangle /(2t)$ . For the turbulent case with ${\it\phi}=0.5\,\%$ we obtain $D_{e}=0.52$ and 0.28 in the directions parallel and perpendicular to gravity, while for ${\it\phi}=1\,\%$ we obtain $D_{e}=0.57$ and 0.32. In the quiescent cases, the diffusion coefficients in the horizontal directions are approximately 0.03, whereas the coefficient in the gravity direction at ${\it\phi}=1\,\%$ is approximately 0.40. Csanady (Reference Csanady1963) proposed a theoretical estimate of the diffusion coefficients for pointwise particles. Using these estimates, we obtain approximately $D_{e}=1.4$ and 0.7 in the directions parallel and perpendicular to gravity. These are approximately 2.5 times larger than those found here for finite-size particles.

3.2. Fluid statistics

Table 3 reports the fluctuation intensities of the fluid velocities for all cases considered. These are calculated by excluding the volume occupied by the spheres at each time step and averaging over the number of samples associated with the fluid phase volume. As expected, the fluid velocity fluctuations are smaller in the quiescent cases than in the turbulent regime. In the quiescent case, the r.m.s. of the velocity fluctuations is approximately 50 % larger in the settling direction than in the normal direction because of the long-range disturbance induced by the particle wakes. The increase of the volume fraction enhances the fluctuations in both directions. Fluctuations are always larger in the turbulent case, with the most significant differences compared with the quiescent cases in the normal direction, where the presence of the buoyant solid phase breaks the isotropy of the turbulent velocity fluctuations.

Table 3.

Fluctuation r.m.s. of the fluid velocities parallel, ${\it\sigma}_{V_{f,{\it\tau}}}$ , and perpendicular, ${\it\sigma}_{V_{f,n}}$ , to gravity. The turbulent fluid velocity undisturbed r.m.s. is ${\sim}0.3$ .

Hence, the solid phase clearly affects the turbulent flow field. Although the present study focuses on the settling dynamics, it is interesting to briefly discuss how turbulence is modulated. Modulation of isotropic turbulence by neutrally buoyant particles is examined in Lucci et al. (Reference Lucci, Ferrante and Elghobashi2010); however, the results change due to buoyancy, as investigated here. Typical turbulent quantities are reported in table 4, where they are compared with the unladen case at ${\it\phi}=0$ . The energy dissipation ${\it\epsilon}$ increases with ${\it\phi}$ , becoming almost double at ${\it\phi}=1\,\%$ . This behaviour is expected since the buoyant particles inject energy into the system which is transformed into kinetic energy of the fluid phase that has to be dissipated. The higher energy flux, i.e. dissipation, is reflected in a reduction of the Kolmogorov length ${\it\eta}$ . The particles reduce the velocity fluctuations, decreasing the turbulent kinetic energy level. The combined effect on $k$ and ${\it\epsilon}$ results in a decrease of the Taylor microscale ${\it\lambda}$ and of $Re_{{\it\lambda}}$ ; likewise, the integral length $L_{0}$ and $Re_{L_{0}}$ also decrease. The reduction of the large and small turbulence scales is associated with the additional energy injection from the settling particles. Energy injection occurs at the size of the particles, which is below the unperturbed integral scale $L_{0}$ , explaining the lowering of the effective integral $L_{0}$ and of Taylor ${\it\lambda}$ length scales. This additional energy is transferred to the bulk flow in the particle wake. Associated with this energy input there is a new mechanism for dissipation, which is the interaction of the flow with the no-slip surfaces of the particles. The mean energy dissipation field in the particle reference frame for the turbulent case with ${\it\phi}=0.5\,\%$ is therefore shown in figure 8. After a statistically steady state is reached, the norm of the symmetric part of the velocity gradient tensor $\unicode[STIX]{x1D60C}_{ij}$ and the dissipation ${\it\epsilon}=2{\it\nu}\unicode[STIX]{x1D60C}_{ij}\unicode[STIX]{x1D60C}_{ij}$ are calculated at each time step on a cubic mesh centred around each particle; the dissipation is calculated on the grid points outside the particle volume. The data presented have been averaged over all particles and time to get the mean dissipation field displayed in the figure. The maximum $\langle {\it\varepsilon}\rangle$ is found around the particle surface, with maximum values in the front; the mean dissipation drops down to the values found in the rest of the domain on the particle rear. The overall energy dissipation is therefore made up of two parts, the first associated with the dissipative eddies far from the particle surfaces and the second associated with the mean and fluctuating flow field near the particle surface. To conclude, the settling strongly alters the typical turbulence features via an anisotropic energy injection and dissipation, thus breaking the isotropy of the unladen turbulent flow. The energy is injected by the fluctuations into the particle wake, whereas stronger energy dissipation occurs at the front of each particle. As a consequence, the fluid velocity fluctuations change in the directions parallel and perpendicular to gravity, as shown in table 3.

Figure 8.

Contours of the mean fluid dissipation, $\langle {\it\varepsilon}\rangle$ ( $v_{t}^{3}/(2a)$ ), averaged in the particle frame of reference.

Table 4.

Turbulent flow parameters in particle units for ${\it\phi}=0$ , ${\it\phi}=0.5\,\%$ and ${\it\phi}=1\,\%$ .

3.3. Relative velocity

An important quantity for understanding and modelling the settling dynamics is the particle to fluid relative motion. Although it is still unclear how to properly calculate the slip velocity between the two phases, we consider spherical shells around each particle, centred on the particle centroids, inspired by the works of Bellani & Variano (Reference Bellani and Variano2012) and Cisse, Homann & Bec (Reference Cisse, Homann and Bec2013). We calculate the mean difference between the particle and fluid velocities in each shell as

(3.3) $$\begin{eqnarray}\langle \boldsymbol{U}_{rel}\rangle _{\boldsymbol{x},t,NP}=\left\langle \boldsymbol{u}_{p}-\frac{1}{{\it\Omega}({\it\Delta})}\int _{{\it\Omega}({\it\Delta})}\boldsymbol{u}_{f}\,\text{d}\mathscr{V}\right\rangle _{t,NP},\end{eqnarray}$$

where ${\it\Omega}({\it\Delta})$ is the volume of a shell of inner radius ${\it\Delta}$ . A parametric study on the slip velocity is performed by changing the inner radii of these spherical shells from ${\it\Delta}=0.75$ particle diameters to ${\it\Delta}=5.0$ particle diameters, while keeping the shell thickness ${\it\delta}$ constant and equal to 0.063 in units of $2a$ . In figure 9 we report the component of $\boldsymbol{U}_{rel}$ parallel to gravity as a function of the inner radii of the shells, ${\it\Delta}/(2a)$ . As the inner radii of the shells increase, $|U_{rel,{\it\tau}}|$ tends exponentially to an asymptotic value which corresponds to the mean particle velocity in the same direction, $\langle V_{p,{\it\tau}}\rangle$ . This is expected since the correlation between the fluid and particle velocities goes to zero at large distances. The quiescent cases still show a 0.5 % difference between $|U_{rel,{\it\tau}}|$ and $\langle V_{p,{\it\tau}}\rangle$ at ${\it\Delta}/(2a)=5$ . This difference is again attributed to the long coherent wakes of the particles.

Figure 9.

The component of $\boldsymbol{U}_{rel}$ aligned with gravity as a function of the inner radii of the shells, ${\it\Delta}/(2a)$ , for the four cases studied. The dashed lines represent instead the mean particle velocities in the direction of gravity, $\langle V_{p,{\it\tau}}\rangle$ .

Figure 10.

Comparison of the p.d.f.s of tangential relative velocity and corresponding particle Reynolds number $Re_{p}$ . The colours and symbols are the same as in the previous figures. In (a) the cases for ${\it\Delta}/(2a)=1.5$ are reported, while in (b) we show the results for ${\it\Delta}/(2a)=2.5$ .

Table 5.

Moments of the $\text{p.d.f.}(U_{rel,{\it\tau}})$ for ${\it\phi}=0.5\,\%$ in the quiescent case. The thickness of the shell used to compute the slip velocity is ${\it\delta}/(2a)=0.063$ .

Table 6.

Moments of the $\text{p.d.f.}(U_{rel,{\it\tau}})$ for ${\it\phi}=0.5\,\%$ in a turbulent environment. The thickness of the shell used to compute the slip velocity is ${\it\delta}/(2a)=0.063$ .

The p.d.f.s of these relative velocities, $U_{rel,{\it\tau}}$ , and their first four central moments are computed and reported in tables 5 and 6 as a function of ${\it\Delta}/(2a)$ for the quiescent fluid and turbulent cases at ${\it\phi}=0.5\,\%$ . In the turbulent case, the moments approach those of a Gaussian distribution, with vanishing skewness and flatness close to 3, especially at large ${\it\Delta}$ . In the quiescent case, the third and fourth moments display higher values that decrease as ${\it\Delta}$ is increased, tending to the values of the particle velocities. The p.d.f.s pertaining to the four cases considered, calculated in spherical shells with inner radii of 1.5 and 2.5 particle diameters, are compared in figure 10. A second axis, reporting the particle Reynolds number $Re_{p}$ based on $U_{rel,{\it\tau}}$ , is also displayed in each figure. In the former case, ${\it\Delta}/(2a)=1.5$ , the shell radius is of the order of the Taylor scale to highlight the particle dynamics, while the relative velocity is approaching asymptotic values for the larger shell. The p.d.f.s of the relative velocity appear to be narrower than those of the particle absolute velocity, indicating that the particles tend to be transported by the large-scale motions, filtering the smallest scales.

The distributions pertaining to the simulations in a turbulent environment are nearly Gaussian, with modal values well below one. The quiescent cases show skewed distributions with long tails at high velocities, as observed for the particle velocities in figure 3(a). The particles settle on average with a velocity close to that of a single particle, with occasional events of higher velocity due to the drafting–kissing–tumbling dynamics. The lower the volume fraction is the more intermittent the dynamics is.

Figure 11.

Joint probability distributions of $V_{f,{\it\tau}}$ and $V_{p,{\it\tau}}$ evaluated in spherical shells located at 1.75 particle diameters from each particle. We report in (a,c) the quiescent cases for ${\it\phi}=0.5\,\%$ and 1 %, while the respective turbulent cases are reported in (b,d). The continuous lines represent the integrals of the joint probability distributions. In the turbulent cases the dashed lines represent the best fit of these integrals.

Knowing the tangential fluid velocity $\langle V_{f,{\it\tau}}\rangle _{\boldsymbol{r}}(p,t)$ , averaged in each shell and at each time step, and the corresponding tangential particle velocities $V_{p,{\it\tau}}(p,t)$ , it is then possible to find their joint probability distribution $P(V_{f,{\it\tau}},V_{p,{\it\tau}})$ (for the sake of simplicity we write $\langle V_{f,{\it\tau}}\rangle _{\boldsymbol{r}}$ as $V_{f,{\it\tau}}$ ). These are evaluated in shells at ${\it\Delta}/(2a)=1.75$ for each case studied and reported in figure 11. In each case, the integral of $P(V_{f,{\it\tau}},V_{p,{\it\tau}})$ ,

(3.4) $$\begin{eqnarray}\langle V_{p,{\it\tau}}|V_{f,{\it\tau}}\rangle =\int _{-\infty }^{\infty }P(V_{f,{\it\tau}},V_{p,{\it\tau}})V_{f,{\it\tau}}\,\text{d}V_{f,{\it\tau}},\end{eqnarray}$$

is also reported (continuous lines). This represents the most probable particle velocity $V_{p,{\it\tau}}$ given a certain fluid velocity $V_{f,{\it\tau}}$ , or, equivalently, the most probable fluid velocity surrounding a particle settling with velocity $V_{p,{\it\tau}}$ . In the turbulent cases, these integrals are well approximated by straight lines (displayed with dashed lines in the figure),

(3.5) $$\begin{eqnarray}V_{p,{\it\tau}}=C_{1}V_{f,{\it\tau}}+C_{2}.\end{eqnarray}$$

In both cases, $C_{1}$ is approximately 1 while $C_{2}$ is approximately 0.86 for ${\it\phi}=0.5\,\%$ and 0.84 for ${\it\phi}=1\,\%$ . These values are in agreement with the values found for the average relative velocities of shells at ${\it\Delta}/(2a)=1.75$ . In a quiescent flow, conversely, the integral in (3.4) gives a curved line and no best fit is therefore reported. In these cases, the joint probability distribution is broader, particularly in the region of higher particle velocities, $V_{p,{\it\tau}}$ . This is again due to the intense particle interactions and the drafting–kissing–tumbling behaviour described in figure 4, which confirms the high flatness of the p.d.f.s of the relative particle velocities.

Further insight can be obtained by plotting isocontours of the average particle relative velocities and their fluctuations, in both quiescent and turbulent flows. To this end, we follow the approach of Garcia-Villalba, Kidanemariam & Uhlmann (Reference Garcia-Villalba, Kidanemariam and Uhlmann2012). We place a uniform and structured rectangular mesh around each particle, with origin at the particle centre. By means of trilinear interpolations we find the fluid and relative velocities on this local mesh and average over time and the number of particles to obtain the mean relative velocity field and its fluctuations.

Figure 12.

Contour plot of the velocity component in the direction of gravity for a single sphere settling in a quiescent fluid.

The vertical velocity of a single sphere settling in a quiescent fluid is reported in figure 12. The relative normal and tangential velocities and their fluctuations are displayed in figure 13 for suspensions with ${\it\phi}=0.5\,\%$ in both quiescent and turbulent environments. In the quiescent fluid simulations, a long wake forms behind the representative particle and, as seen from the single-point particle velocity correlations, it takes a long time for the velocity fluctuations to decorrelate. In the turbulent case instead, the wakes are disrupted by the background fluctuations.

Interesting observations can be drawn from the relative velocity fluctuation fields. Comparing figures 13(c) and 13(d) we note that intense vortex shedding occurs around the particles in the turbulent case, with important fluctuations of $U_{rel,{\it\tau}}$ . From figure 13(e,f) we also see that the relative velocity fluctuations are drastically increased in the horizontal directions in a turbulent environment. It is noteworthy that vortex shedding occurs at particle Reynolds numbers below the critical value above which this is usually observed (Bouchet, Mebarek & Dušek Reference Bouchet, Mebarek and Dušek2006). Vortex shedding is unsteady in nature, and unsteady effects may therefore play an important role in the increase of the overall drag, as further discussed below. Lower fluctuation intensities are found on the front parts of the particles, where the energy dissipation is highest, and in the immediate wake in the recirculating region where the instability is found to develop.

Figure 13.

Fields of $\langle U_{rel,{\it\tau}}\rangle$ , $U_{rel,{\it\tau}}^{\prime }$ and $U_{rel,n}^{\prime }$ for the quiescent (a,c,e) and turbulent (b,d,f) cases.

3.4. Drag analysis

As in Maxey & Riley (Reference Maxey and Riley1983), we write the balance of the forces acting on a single sphere settling through a turbulent flow. The equation of motion for a spherical particle reads

(3.6) $$\begin{eqnarray}\frac{4}{3}{\rm\pi}a^{3}{\it\rho}_{p}\frac{\text{d}\boldsymbol{V}_{p}}{\text{d}t}=\frac{4}{3}{\rm\pi}a^{3}({\it\rho}_{p}-{\it\rho}_{f})\boldsymbol{g}+\oint _{\partial \mathscr{V}_{p}}{\bf\tau}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S,\end{eqnarray}$$

where the integral is over the surface of the sphere $\partial \mathscr{V}_{p}$ , $\boldsymbol{n}$ is the outward normal and ${\bf\tau}=-p\unicode[STIX]{x1D644}+2{\it\mu}\unicode[STIX]{x1D640}$ is the fluid stress. As commonly done in aerodynamics, we replace the last term of (3.6) with a term depending on the relative velocity $\boldsymbol{U}_{rel}$ and a drag coefficient $C_{D}$ ,

(3.7) $$\begin{eqnarray}\frac{4}{3}{\rm\pi}a^{3}{\it\rho}_{p}\frac{\text{d}\boldsymbol{V}_{p}}{\text{d}t}=\frac{4}{3}{\rm\pi}a^{3}({\it\rho}_{p}-{\it\rho}_{f})\boldsymbol{g}-\frac{1}{2}{\it\rho}_{f}{\rm\pi}a^{2}|\boldsymbol{U}_{rel}|\boldsymbol{U}_{rel}C_{D},\end{eqnarray}$$

with ${\rm\pi}a^{2}$ the reference area. Generally, the drag coefficient $C_{D}$ is a function of a Reynolds number and a Strouhal number which accounts for unsteady effects. In the present case, we consider it to be a function of the Reynolds number based on the relative velocity $Re_{p}=2a|\boldsymbol{U}_{rel}|/{\it\nu}$ (in a turbulent field it is proper to define $Re_{p}$ in terms of the relative velocity between the particle and the fluid) and of the Strouhal number defined as

(3.8) $$\begin{eqnarray}St=\frac{\displaystyle \frac{\text{d}\boldsymbol{V}_{p}}{\text{d}t}(2a)}{|\boldsymbol{U}_{rel}|^{2}}.\end{eqnarray}$$

The drag on the particle depends on both nonlinear and unsteady effects (such as the Basset history force and the added-mass contribution) through these two non-dimensional numbers.

Commonly, the unsteady contribution is neglected and $C_{D}$ is assumed to depend only on the Reynolds number. Since we want to investigate both nonlinear and unsteady effects, we decide to express the drag coefficient as $C_{D}(Re_{p},St)=C_{D_{0}}(Re_{p})[1+{\it\psi}(Re_{p},St)]$ , yielding

(3.9) $$\begin{eqnarray}\frac{4}{3}{\rm\pi}a^{3}{\it\rho}_{p}\frac{\text{d}\boldsymbol{V}_{p}}{\text{d}t}=\frac{4}{3}{\rm\pi}a^{3}({\it\rho}_{p}-{\it\rho}_{f})\boldsymbol{g}-\frac{1}{2}{\it\rho}_{f}{\rm\pi}a^{2}|\boldsymbol{U}_{rel}|\boldsymbol{U}_{rel}C_{D_{0}}(Re_{p})[1+{\it\psi}(Re_{p},St)],\end{eqnarray}$$

where ${\it\psi}=0\,\forall \,Re_{p}$ if $St=0$ (steady motion). We can therefore identify a quasi-steady term and the extra term that accounts for unsteady phenomena.

By ensemble averaging equation (3.9) over time and the number of particles, and assuming the settling process to be at statistically steady state, we can find the most important contributions to the overall drag. The steady-state average equation reads

(3.10) $$\begin{eqnarray}0={\textstyle \frac{4}{3}}{\rm\pi}a^{3}({\it\rho}_{p}-{\it\rho}_{f})\boldsymbol{g}-{\textstyle \frac{1}{2}}{\it\rho}_{f}{\rm\pi}a^{2}\langle |\boldsymbol{U}_{rel}|\boldsymbol{U}_{rel}C_{D_{0}}(Re_{p})[1+{\it\psi}(Re_{p},St)]\rangle .\end{eqnarray}$$

Denoting the entire time-dependent term simply as ${\it\Psi}(t)$ and rearranging, we obtain the following balance:

(3.11) $$\begin{eqnarray}\frac{4}{3}{\rm\pi}a^{3}\left(\frac{{\it\rho}_{p}}{{\it\rho}_{f}}-1\right)\boldsymbol{g}=\frac{1}{2}{\rm\pi}a^{2}\langle |\boldsymbol{U}_{rel}|\boldsymbol{U}_{rel}C_{D_{0}}(Re_{p})\rangle +{\it\Psi}(t).\end{eqnarray}$$

The term on the left-hand side is known, whereas the time-dependent term ${\it\Psi}(t)$ is of difficult evaluation. The nonlinear steady term can be calculated using the approach described in § 2.2. At each time step we calculate the relative velocity in the spherical shells surrounding each particle. From these we compute the particle Reynolds number $Re_{p}=2a|U_{rel}|/{\it\nu}$ , and using (2.8) (where we replace the terminal Reynolds number with the new particle Reynolds number) we obtain the drag coefficient. The first term on the right-hand side can therefore be evaluated by averaging over the number of particles and time steps considered. Finally, we divide everything by the buoyancy acceleration to estimate the relative importance of each term,

(3.12) $$\begin{eqnarray}100\,\%=\mathscr{S}(Re_{p})+{\it\Psi}(t),\end{eqnarray}$$

where $\mathscr{S}(Re_{p})$ represents the nonlinear steady term, while the unsteady term has been denoted again as ${\it\Psi}(t)$ for simplicity.

This approach is applied to the results of the simulations of a single sphere and to the quiescent and turbulent cases with ${\it\phi}=0.5\,\%$ and 1 %. The inner radius of the sampling shells is chosen to be 5 particle radii, and the results obtained are reported in table 7. The single-sphere simulation provides an estimate of the error of the method. Since our terminal Reynolds number is smaller than the critical Reynolds number above which unsteady effects become important, and our velocity field is indeed steady, the term ${\it\Psi}(t)$ should be negligible. The critical Reynolds number $Re_{c}$ has been found to be approximately 274 by Bouchet et al. (Reference Bouchet, Mebarek and Dušek2006), while we recall that our terminal Reynolds number for the isolated particle is 200. The nonlinear steady term provides in this case, however, an overestimated value of the drag, with a percentage error of approximately $+$ 3 % with respect to the buoyancy term. The possible causes of this error are the long wake and the fact that (2.8) is empirical. The data from the single-particle simulations are used to correct the results from the other runs, i.e. the data are normalized with the total drag obtained in this case. In the quiescent case at ${\it\phi}=0.5\,\%$ , unsteady effects are negligible (approximately 0.5 % of the total drag), while they increase to approximately 6 % on increasing the particle evolve fraction to 1 %. In a turbulent flow, importantly, we notice that the contribution of ${\it\Psi}(t)$ adds up to approximately 10 % of the total at ${\it\phi}=0.5\,\%$ and to approximately 12 % of the total at ${\it\phi}=1\,\%$ .

Table 7.

Percentage contributions of the nonlinear and unsteady terms for the quiescent and turbulent cases with ${\it\phi}=0.5\,\%$ . The data are normalized by the mean drag from the simulation of a single sphere in quiescent fluid.

It should be noted that one can write the steady drag as mean and fluctuating components,

(3.13) $$\begin{eqnarray}\mathscr{S}(Re_{p})=\langle |\boldsymbol{U}_{rel}|\boldsymbol{U}_{rel}C_{D_{0}}(Re_{p})\rangle =\langle U_{rel}\rangle ^{2}C_{D_{0}}(\langle Re_{p}\rangle )+\mathscr{S}^{\prime }(\langle Re_{p}\rangle ).\end{eqnarray}$$

The fluctuations $\mathscr{S}^{\prime }(\langle Re_{p}\rangle )$ would be responsible for the reduction of the settling velocity if this were to be attributed to nonlinear drag effects only (see also Wang & Maxey Reference Wang and Maxey1993). We verified that for our results the total and mean components differ by approximately 2 %, $\langle |\boldsymbol{U}_{rel}|\boldsymbol{U}_{rel}C_{D_{0}}(Re_{p})\rangle \approx \langle U_{rel}\rangle ^{2}C_{D_{0}}(\langle Re_{p}\rangle )$ . This leads us to the conclusion that the main contribution to the overall drag is due to the steady term, but the reduction of the mean settling velocity in a turbulent environment is almost entirely due to the various unsteady effects. These can be related to unsteady vortex shedding, see figure 13, as in the experiments on a single sphere by Mordant & Pinton (Reference Mordant and Pinton2000). These observations are also in agreement with the results of Homann et al. (Reference Homann, Bec and Grauer2013). These authors observed that the enhancement of the drag of a sphere towed in a turbulent environment can be explained by the modification of the mean velocity and pressure profile, and thus of the boundary layer around the sphere, by the turbulent fluctuations.