3.1 Particle fields

In this section we explore the spatial structure of the inertial particle fields and the length scales over which clustering occurs. In the literature this has been characterized by two main tools: the radial distribution function (RDF) and, more recently, Voronoï tessellation. Both methods provide different and complementary information – we apply them both for a comprehensive description.

The RDF describes the scale-by-scale concentration in the space surrounding a generic particle, compared to a uniform distribution. For 2-D distributions such as those obtained by planar imaging, this is defined as:

(3.1) $$\begin{eqnarray}g(r)=\frac{N_{r}/A_{r}}{N/A_{tot}},\end{eqnarray}$$

where $N_{r}$ represents the number of particles within an annulus of area $A_{r}$ , while $N$ is the total number of particles within the planar domain of area $A_{tot}$ . In the presence of clustering, the RDF is expected to increase for decreasing $r$ , and the range over which it remains greater than unity indicates the length scale over which clustering occurs (e.g. Sundaram & Collins Reference Sundaram and Collins1997; Reade & Collins Reference Reade and Collins2000; Wood, Hwang & Eaton Reference Wood, Hwang and Eaton2005; Saw et al. Reference Saw, Shaw, Ayyalasomayajula, Chuang and Gylfason2008; de Jong et al. Reference de Jong, Salazar, Woodward, Collins and Meng2010; Ireland et al. Reference Ireland, Bragg and Collins2016a ,Reference Ireland, Bragg and Collins b ). We compute RDFs by binning particle pairs based on their separation distance. To avoid projection biases at separations below the illuminated volume thickness (Holtzer & Collins Reference Holtzer and Collins2002), we only calculate $g(r)$ for $r>1.5~\text{mm}$ . As noted by de Jong et al. (Reference de Jong, Salazar, Woodward, Collins and Meng2010), imaging-based RDF measurements are sensitive to the size and shape of the observation region, and some sort of edge-correction strategy is needed for particles near the image boundaries. One can omit statistics for radial annuli that cross the image boundary, but this approach has two shortcomings: the maximum separation becomes limited to the radius of the domain-inscribed circle; and the number of particle pairs per unit area used to calculate the RDF decreases as the separation increases. Both effects combine to thwart the reliable assessment of large-scale clustering. Indeed, past RDF measurements at distances $O(L)$ in flows with wide scale separation were obtained using single-point probes and invoking Taylor’s hypothesis (Saw et al. Reference Saw, Shaw, Ayyalasomayajula, Chuang and Gylfason2008, Reference Saw, Shaw, Salazar and Collins2012; Bateson & Aliseda Reference Bateson and Aliseda2012). Here, following de Jong et al. (Reference de Jong, Salazar, Woodward, Collins and Meng2010), we leverage the spatial homogeneity of our fields and apply a periodic-domain correction: the particle field is mirrored across the image boundaries, so that the same number of radial annuli can be used for each particle location, yielding a maximum separation equal to the full size of the FOV. Although this assumption introduces unphysical correlations between particles near the reflected boundaries, numerical experiments using DNS showed the associated error to be small (Salazar et al. Reference Salazar, de Jong, Cao and Woodward2008).

Figure 4.

(a) Example RDFs where symbols represent the calculated RDFs and dashed lines their fits according to (3.2). Panels (b,c) show the fitted values versus $St_{\unicode[STIX]{x1D702}}$ , and (d) shows the value of the correlation dimension.

Due to the large range of scales separation ( $L/\unicode[STIX]{x1D702}\sim O(10^{3})$ , $L/d_{p}\sim O(10^{4}{-}10^{5})$ ) it is not feasible to simultaneously resolve all scales at play, and thus the small-FOV and large-FOV imaging will suffer from large-scale and small-scale cutoffs, respectively. Comparing the different FOVs, however, provides quantitative information over a wide range of scales. Figure 5(a) shows examples from the small-FOV measurements. At small separations, several authors have found satisfactory fit to the data using a power law, which indicates a self-similar spatial distribution (Chun et al. Reference Chun, Koch, Rani, Ahluwalia and Collins2005; Salazar et al. Reference Salazar, de Jong, Cao and Woodward2008; Zaichik & Alipchenkov Reference Zaichik and Alipchenkov2009; Ireland et al. Reference Ireland, Bragg and Collins2016a ,Reference Ireland, Bragg and Collins b ):

(3.2) $$\begin{eqnarray}g(r/\unicode[STIX]{x1D702})=c_{0}(r/\unicode[STIX]{x1D702})^{-c_{1}},\end{eqnarray}$$

where $c_{0}$ and $c_{1}$ are coefficients dependent on $St_{\unicode[STIX]{x1D702}}$ (and, in presence of gravity, $Sv_{\unicode[STIX]{x1D702}}$ ). While theoretical arguments consistent with this formulation strictly apply for dissipative separations ( $r/\unicode[STIX]{x1D702}<1$ , Chun et al. (Reference Chun, Koch, Rani, Ahluwalia and Collins2005) and Saw et al. (Reference Saw, Shaw, Ayyalasomayajula, Chuang and Gylfason2008) argued that the power-law form should continue into the correlation scale of the velocity gradients ( $r/\unicode[STIX]{x1D702}=O(10)$ ). In figure 4(a) we see indeed that RDFs closely follow a power-law decay up to $r/\unicode[STIX]{x1D702}\approx 40$ for $St_{\unicode[STIX]{x1D702}}$ close to unity. The departure from the power law at larger separations indicates the particle set is not self-similar at those scales (Bragg, Ireland & Collins Reference Bragg, Ireland and Collins2015). We evaluate the coefficients $c_{0}$ and $c_{1}$ by least-square fits over the range $10<r/\unicode[STIX]{x1D702}<30$ , and plot them in figure 4(b,c) as a function of $St_{\unicode[STIX]{x1D702}}$ . The error bars for the coefficients comes from the covariance matrix of the fit. The results, which are only weakly sensitive to varying the fit upper bound between $r/\unicode[STIX]{x1D702}=20$ and 40, display the higher values in the approximate range $1.5<St_{\unicode[STIX]{x1D702}}<4$ . This confirms that particles with Stokes number $O(1)$ display the stronger degree of clustering over the near-dissipative range. However, the most intense clustering occurs for $St_{\unicode[STIX]{x1D702}}>1$ , possibly because of the significant effect of gravitational settling as discussed below. The trend and values are in fair quantitative agreement with the DNS of Ireland et al. (Reference Ireland, Bragg and Collins2016b ) in similar conditions. The $c_{1}$ coefficient is related to the correlation dimension used in dynamical system theory (Bec et al. Reference Bec, Cencini, Hillerbrand and Turitsyn2008), $D_{2}=n-c_{1}$ , where the number of spatial dimensions is $n=2$ for our planar realizations. Figure 5(d) plots $D_{2}$ for the different cases, showing trends and values consistent with the channel flow experiments by Fessler et al. (Reference Fessler, Kulick and Eaton1994) and the grid turbulence experiments by Monchaux et al. (Reference Monchaux, Bourgoin and Cartellier2010, Reference Monchaux, Bourgoin and Cartellier2012). For increasing particle inertia, one expects a loss of spatial correlation as the particle response time grows beyond the fine turbulent scales. Although the present range does not extend to very large $St_{\unicode[STIX]{x1D702}}$ , we note that the return to a homogeneous distribution appears slow. Recent numerical studies compared settling and non-settling conditions, and concluded that gravity hinders clustering for $St_{\unicode[STIX]{x1D702}}<1$ but enhances it for $St_{\unicode[STIX]{x1D702}}>1$ , resulting in significant clustering over a wide range of Stokes numbers (Bec, Homann & Ray Reference Bec, Homann and Ray2014b ; Gustavsson, Vajedi & Mehlig Reference Gustavsson, Vajedi and Mehlig2014; Ireland et al. Reference Ireland, Bragg and Collins2016b ; Baker et al. Reference Baker, Frankel, Mani and Coletti2017; Matsude, Onishi & Takahashi Reference Matsude, Onishi and Takahashi2017). Ireland et al. (Reference Ireland, Bragg and Collins2016b ) attributed this behaviour to the competing effects of the particle path history and preferential flow sampling. Sahu et al. (Reference Sahu, Hardalupas and Taylor2016) measured RDFs for spray droplets and also noticed an increasing tendency to cluster for increasing $St_{\unicode[STIX]{x1D702}}$ (although their range was very close to unity).

The large-FOV measurements allow us to probe the spatial distribution over much greater scales. Figure 5(a) clearly indicates that considerable clustering occurs over lengths $O(L)$ . For a quantitative assessment, we consider the original power-law model proposed by Reade & Collins (Reference Reade and Collins2000):

(3.3) $$\begin{eqnarray}g(r/\unicode[STIX]{x1D702})-1=c_{0}^{\ast }(r/\unicode[STIX]{x1D702})^{-c_{1}^{\ast }}\text{e}^{-c_{2}^{\ast }r/\unicode[STIX]{x1D702}},\end{eqnarray}$$

which, unlike (3.2), does recover the return to unity at large separations. The excellent fit to the data over the entire window confirms the observation made by Reade & Collins (Reference Reade and Collins2000), that the RDFs of preferentially concentrated particles have a power-law decay at small scales and an exponential tail at large scales. Figure 5(b–d) displays the least-square-fit coefficients as a function of $St_{\unicode[STIX]{x1D702}}$ . The $c_{2}^{\ast }$ coefficient decreases as $St_{\unicode[STIX]{x1D702}}$ increases, implying a greater spatial extent of clustering for the more inertial particles. This may be due to the more inertial particles responding to larger time scales of the turbulence, and to the influence of gravitational settling as mentioned above. The length scale of the large-scale clustering can be estimated from the exponential decay as $\unicode[STIX]{x1D702}/c_{2}$ , which for $St_{\unicode[STIX]{x1D702}}>1$ is approximately $300\unicode[STIX]{x1D702}{-}400\unicode[STIX]{x1D702}$ , close to the integral scale of the turbulence. Taken together, these results confirm that clustering can extend over larger scales for heavier particles. This is in agreement with the conceptual picture of Goto & Vassilicos (Reference Goto and Vassilicos2006) and Yoshimoto & Goto (Reference Yoshimoto and Goto2007) and the simulations of Bec et al. (Reference Bec, Biferale, Cencini and Lanotte2010) and Ireland et al. (Reference Ireland, Bragg and Collins2016a ,Reference Ireland, Bragg and Collins b ), which showed that particles of $St_{\unicode[STIX]{x1D702}}>1$ respond to eddies in the inertial range. However, as we will reiterate in the next sub-section, the present results indicate that some level of clustering may extend even beyond, approaching the integral scales.

Figure 5.

(a) Example RDFs where symbols represent the calculated RDFs and dashed lines their fits according to (3.3). Panels (b–d) show the fitted coefficients versus $St_{\unicode[STIX]{x1D702}}$ .

The degree of clustering can also be evaluated from the Voronoï diagrams. We remark that the type of information provided by this method is somewhat different than the RDF. The latter is strictly a two-particle quantity, while the shape and size of the Voronoï cells result from the mutual position of multiple particles. Therefore, we look for an insight complementary to our RDF results. In figure 6(a) we plot the standard deviation of the Voronoï cell areas $\unicode[STIX]{x1D70E}_{A}$ as a function of the Stokes number, normalizing it by the expected value for particles distributed according to a random Poisson process, $\unicode[STIX]{x1D70E}_{RPP}\approx 0.53$ (Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2010). As a general trend, clustering is most pronounced for particles of $St_{\unicode[STIX]{x1D702}}\approx 1$ , in agreement with previous studies (Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2010; Tagawa et al. Reference Tagawa, Mercado, Prakash, Calzavarini, Sun and Lohse2012; Dejoan & Monchaux Reference Dejoan and Monchaux2013; Monchaux & Dejoan Reference Monchaux and Dejoan2017). However, the significant scatter suggests that other parameters may also play a role. Indeed, in their grid turbulence study, Sumbekova et al. (Reference Sumbekova, Cartellier, Aliseda and Bourgoin2017) found that $\unicode[STIX]{x1D70E}_{A}$ was strongly dependent on $Re_{\unicode[STIX]{x1D706}}$ , moderately on $\unicode[STIX]{x1D719}_{v}$ and negligibly on $St_{\unicode[STIX]{x1D702}}$ . While the considerable degree of polydispersity in their experiments may have influenced such conclusion, their results convincingly indicated that clustering is affected by a range of turbulent scales, whose breadth is controlled by $Re_{\unicode[STIX]{x1D706}}$ . Moreover, as pointed out by Baker et al. (Reference Baker, Frankel, Mani and Coletti2017), $\unicode[STIX]{x1D70E}_{A}$ is not only a function of the concentration of the clustered particles, but also of the size and distribution of the voids, and the latter are strongly influenced by the inertial and integral scales of the turbulence (Yoshimoto & Goto Reference Yoshimoto and Goto2007).

The decrease in $\unicode[STIX]{x1D70E}_{A}$ for $St_{\unicode[STIX]{x1D702}}>1$ is mild. Since increasing $St_{\unicode[STIX]{x1D702}}$ also implies increasing $Sv_{\unicode[STIX]{x1D702}}$ , this again suggests that, in this range, gravitational settling may enhance clustering. This idea is supported by figure 6(b), showing the fraction of particles belonging to coherent clusters (according to the definition in § 2.3) plotted versus $\unicode[STIX]{x1D70E}_{A}$ . A clear correlation is visible, indicating that the number of clustered particles, $N_{c}$ , is similarly affected by the physical parameters. The values are possibly underestimated, because particles in highly concentrated regions are more likely to be overshadowed by neighbouring particles and go undetected. Still, the results are in fair agreement with the DNS of Baker et al. (Reference Baker, Frankel, Mani and Coletti2017), where less than 3 % of the particles with $St_{\unicode[STIX]{x1D702}}<1$ belonged to coherent clusters, with the percentage increasing up to 14 % for $St_{\unicode[STIX]{x1D702}}=O(10)$ . Therefore, the more inertial (and faster falling) particles are more likely to belong to clusters; or, equivalently, they tend to form clusters that are more numerous, larger, or denser. The question of the cluster size and the concentration within them is addressed in the following section.

Figure 6.

Two qualifications of clustering intensity using the Voronoï tessellation method: (a) $\unicode[STIX]{x1D70E}_{A}/\unicode[STIX]{x1D70E}_{RPP}$ representing departure from a random particle distribution versus Stokes number. (b) Fraction of inertial particles belonging to coherent clusters as a function of normalized $\unicode[STIX]{x1D70E}_{A}/\unicode[STIX]{x1D70E}_{RPP}$ .

3.2 Individual clusters

The multi-scale nature of the clustering process is reflected in the features of the individual clusters. Figure 7 shows several sample clusters from various instantaneous realizations, as captured by the large-FOV measurements and identified by the Voronoï diagram method, illustrating the variety of sizes and complex shapes of these objects. Some of them are even larger than the integral scales of the flow, often exceeding the limits of the imaging window. Their borders are jagged and convoluted, and their bodies non-simply connected. In the following, we provide quantitative support to these observations. We stress that the objects captured by 2-D imaging are cross-sections of 3-D clusters; this naturally conditions our ability to assess their topology. Such limitation, however, is not expected to overshadow the main conclusions of the analysis.

Figure 7.

Example clusters for the cases $St_{\unicode[STIX]{x1D702}}=1.6$ and 4.6 imaged in the large FOV, highlighting the wide variety of sizes and shapes.

Figure 8 shows the PDFs of the areas $A_{C}$ of the clusters, coherent and not, distinguishing between measurements obtained over the small, intermediate and large FOV. In agreement with Baker et al. (Reference Baker, Frankel, Mani and Coletti2017) and Sumbekova et al. (Reference Sumbekova, Cartellier, Aliseda and Bourgoin2017), most cases display a power-law behaviour over several decades, suggesting a self-similar hierarchy of structures, possibly associated with the scale-invariant properties of the underlying turbulent field (Moisy & Jiménez Reference Moisy and Jiménez2004; Goto & Vassilicos Reference Goto and Vassilicos2006). The data are consistent with the previously suggested values of $-2$ and $-5/3$ for the power-law exponent for planar measurements (Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2010; Obligado et al. Reference Obligado, Teitelbaum, Cartellier, Mininni and Bourgoin2014; Sumbekova et al. Reference Sumbekova, Cartellier, Aliseda and Bourgoin2017).

As expected, the spatial resolution influences the size distributions. The small FOV is affected by a cutoff at large scales. At small scales, the limited resolution in the large FOV makes particles more likely to go undetected due to the glare of their neighbours, reducing the probability of finding small clusters. The latter effect can partly explain why the area threshold for self-similarity (see vertical dashed line in figure 3 d) varies significantly between cases, while this was found to be very consistent in the simulations of Baker et al. (Reference Baker, Frankel, Mani and Coletti2017). Another factor influencing this threshold is the particle volume fraction. Although the value of $A^{\ast }$ was shown to be robust to particle sub-sampling (Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2012; Baker et al. Reference Baker, Frankel, Mani and Coletti2017), varying the number of particles in the domain results in a shift of the cluster area distribution (figure 9), which in turn affects the number of detected coherent clusters above the self-similar threshold. Finally, as $\unicode[STIX]{x1D719}_{v}$ increases, the possibility of significant two-way coupling effects also increases, which may alter the turbulence structure and consequently the clustering process. This aspect will be discussed in § 5.5.

Figure 8.

PDFs of cluster area normalized by $\unicode[STIX]{x1D702}$ and separated by FOV size. Largest FOV ( $30\times 30~\text{cm}^{2}$ ) is shown in (a), medium field ( $14.5\times 14.5~\text{cm}^{2}$ ) in (b) and smallest field in (c) ( $4.5\times 4.5~\text{cm}^{2}$ ).

Figure 9.

PDFs of cluster area for experiments with $St_{\unicode[STIX]{x1D702}}$ between 1 and 2.5 ( $30~\unicode[STIX]{x03BC}\text{m}$ glass spheres) in the medium FOV showing how increasing volume fraction shifts PDFs to smaller scales.

A remarkable aspect of the distributions in figure 8 is the non-negligible probability of finding clusters of size comparable to the inertial scales of the turbulence. Since our definition of coherent clusters entails a power-law size distribution, and this is found to have an exponent close to $-2$ , the mean area of the coherent clusters is ill defined. In order to compare with past studies, we calculate the mean area of all clusters $\langle A_{C}\rangle$ , below and above the self-similarity threshold, and plot its square root in figure 10. The majority of cases display mean sizes between $10\unicode[STIX]{x1D702}$ and $40\unicode[STIX]{x1D702}$ , with a generally increasing trend with $St_{\unicode[STIX]{x1D702}}$ . Several previous studies have reported mean cluster sizes around $10\unicode[STIX]{x1D702}$ (Aliseda et al. Reference Aliseda, Cartellier, Hainux and Lasheras2002; Wood et al. Reference Wood, Hwang and Eaton2005; Dejoan & Monchaux Reference Dejoan and Monchaux2013). Most of these studies, however, considered turbulent flows with relatively low $Re_{\unicode[STIX]{x1D706}}$ and thus limited scale separation. Recently, Sumbekova et al. (Reference Sumbekova, Cartellier, Aliseda and Bourgoin2017) investigated droplets in grid turbulence at $Re_{\unicode[STIX]{x1D706}}$ approaching 500, and found cluster size distributions and averages comparable with ours. At large Reynolds numbers the spectrum of temporal scales widens, and particles with a broad range of response times become susceptible to clustering mechanisms (Yoshimoto & Goto Reference Yoshimoto and Goto2007). As noted in § 3.2, the more inertial particles respond to larger eddies and therefore can agglomerate in larger sets. The observed dependence of the cluster size with $St_{\unicode[STIX]{x1D702}}$ is consistent with this view.

Figure 10.

Mean cluster area versus particle inertia. Marker size is proportional to $Re_{\unicode[STIX]{x1D706}}$ of the underlying turbulent flow.

To investigate the degree of self-similarity exhibited by individual clusters, and to provide a descriptor of their complex shape, we calculate their box-counting dimension. This has been widely used to characterize the topology of both turbulent structures (Moisy & Jiménez Reference Moisy and Jiménez2004; Lozano-Durán, Flores & Jiménez Reference Lozano-Durán, Flores and Jiménez2012; Carter & Coletti Reference Carter and Coletti2018) and particle clusters in turbulence (Baker et al. Reference Baker, Frankel, Mani and Coletti2017). The domain is partitioned into non-overlapping square boxes of side length $r$ , and for each cluster we count the number of boxes $N_{B}$ containing at least one particle. If $N_{B}(r)$ follows a power law, i.e. $N_{B}\sim r^{-D}$ , over a sizable range of scales, the exponent $D$ is taken as the box-counting dimension of the object, which is in turn a measure of its fractal dimension. (Several other definitions of fractal dimension exist, and typically they only coincide for mathematical constructs, Falconer (Reference Falconer2003).) Relatively large objects are needed for a robust estimate of $D$ over a wide range of scales, and we thus consider only clusters of area larger than $10^{4}\unicode[STIX]{x1D702}^{2}$ . Additionally, we neglect clusters touching the image boundary, as their silhouette would include spurious straight segments. Figure 11 shows, for three sample cases, $N_{B}(r)$ normalized by the maximum number of boxes for each cluster (corresponding to the smallest box size, $r=\unicode[STIX]{x1D702}$ ). For each case, curves for only 20 example clusters are shown for clarity. These reveal a remarkably consistent box-counting dimension $D\approx 1.6$ over at least a decade of scales; the same trend is followed by all other cases. Baker et al. (Reference Baker, Frankel, Mani and Coletti2017) found $D\approx 1.9$ for 3-D clusters. Relating the box-counting dimension of 3-D objects and their 2-D cross-sections is not straightforward (Tang & Marangoni Reference Tang and Marangoni2006; Carter & Coletti Reference Carter and Coletti2018). Rather, the present result may be compared with that of Carter & Coletti (Reference Carter and Coletti2018) who evaluated the box-counting dimension of turbulent coherent structures using 2-D PIV in the same facility. They found $D\approx 1.5$ , which suggests a strong link between the particle cluster topology and the underlying turbulent flow. Beside the precise value of the box-counting dimension, the main observation is that large clusters of inertial particles do exhibit a scale-invariant shape in the present range of $St_{\unicode[STIX]{x1D702}}$ and $Sv_{\unicode[STIX]{x1D702}}$ .

Figure 11.

$N_{B}(r)$ for three cases with $St_{\unicode[STIX]{x1D702}}=4.6$ in (a), 3.6 in (b) and 2.0 in (c), normalized by the respective maximum. The box-counting dimension of 1.6 is consistent for all cases.

In order to characterize the spatial distribution of particles within each cluster, we use the singular value decomposition (SVD) method introduced by Baker et al. (Reference Baker, Frankel, Mani and Coletti2017). The SVD provides the principal axes and corresponding singular values for a particle set. In two dimensions, the primary axis lies along the direction of greatest particle spread from the cluster centroid, the secondary axis being orthogonal to it (figure 12 a). The corresponding singular values $s_{1}$ and $s_{2}$ measure the particle spread along the respective principal axes, and can be used as simple shape descriptors through the aspect ratio $s_{2}/s_{1}$ : the limit $s_{2}/s_{1}=0$ corresponds to particles arranged in a straight line, whereas $s_{2}/s_{1}=1$ corresponds to a perfect circle. The PDF of the aspect ratio for all considered cases (figure 12 b) shows that clusters are likely to exhibit aspect ratios between 0.4 and 0.5, reflecting a tendency to form somewhat elongated objects. Furthermore, the distribution has a positive skew, indicating that globular shapes are more common than extremely long streaks. While these observations are influenced by the 2-D nature of the technique, they are consistent with the results of Baker et al. (Reference Baker, Frankel, Mani and Coletti2017) who found that 3-D clusters had $s_{2}/s_{1}$ distributions peaking around 0.5, and were positively skewed.

Figure 12.

(a) A single coherent cluster, with primary and secondary axes as computed through SVD. The length of each arrow is proportional to the corresponding singular value. (b) PDFs of the aspect ratios of clusters for each case.

The orientation of the primary axis from the SVD analysis provides information on the cluster orientation in space. In figure 13(a) we plot the PDF of the cosine of $\unicode[STIX]{x1D703}_{g}$ , i.e. the angle between the cluster primary axis and the vertical, evidencing a strong preference for the clusters to align with gravity. That particles tend to agglomerate along their falling direction was found in past one-way-coupled simulations (Woittiez, Jonker & Portela Reference Woittiez, Jonker and Portela2009; Dejoan & Monchaux Reference Dejoan and Monchaux2013; Bec et al. Reference Bec, Homann and Ray2014a ; Ireland et al. Reference Ireland, Bragg and Collins2016b ; Baker et al. Reference Baker, Frankel, Mani and Coletti2017), indicating the mechanism is not necessarily related to the particle backreaction on the flow. Baker et al. (Reference Baker, Frankel, Mani and Coletti2017) reasoned that, especially for cases with high $St_{\unicode[STIX]{x1D702}}$ and high $Sv_{\unicode[STIX]{x1D702}}$ , particles are influenced by intermittent downward gusts that add to their fall speed, channelling them and creating elongated quasi-vertical structures. The joint probability distribution of cos( $\unicode[STIX]{x1D703}_{g}$ ) and $A_{C}$ (figure 13 b) supports this view, showing that a vertical alignment corresponds to generally larger clusters. Further studies, possibly including time-resolved information, are needed to gain a mechanistic understanding of the cluster formation process.

Figure 13.

(a) PDFs of the cosine of the angle between clusters’ primary axis and the vertical. (b) Joint PDF of angle from vertical and size of cluster.

We finally consider the concentration of particles within each coherent cluster, $C_{C}=N_{PC}/A_{C}$ , where $N_{PC}$ is the number of particles in each cluster. Figure 14(a,b) shows scatter plots of cluster areas and number of particles for two representative cases. The excellent fit using a power law of exponent close to unity indicates that the relationship is approximately linear, i.e. the concentration within each cluster is approximately the same for a given case. This trend is recovered for all cases. Considering the wide range of sizes, this result (reported by Baker et al. (Reference Baker, Frankel, Mani and Coletti2017) at a much lower $Re_{\unicode[STIX]{x1D706}}$ ) indicates again that the clusters display scale-invariant features. Figure 14(c) illustrates the average in-cluster concentration as a function of $St_{\unicode[STIX]{x1D702}}$ . Despite the scatter (which points to the concurrent effects of the multiple parameters at play), one notices an increase up to $St_{\unicode[STIX]{x1D702}}\approx 2$ , followed by a plateau. The concentration within clusters can be up to an order of magnitude higher than the average over the whole particle field ( $C_{0}=1/\langle A\rangle$ ); these values are likely underestimated as particles may shadow each other at high local concentration. The present results are comparable to those from the experiments by Monchaux et al. (Reference Monchaux, Bourgoin and Cartellier2010) and the DNS by Baker et al. (Reference Baker, Frankel, Mani and Coletti2017).

Figure 14.

(a,b) Scatter plots of cluster area (normalized by the Kolmogorov scale) versus the number of particles in each cluster (normalized by the total number of particles in the domain) for two sample cases. (c) Average particle concentration in the clusters, normalized by the global concentration, plotted versus $St_{\unicode[STIX]{x1D702}}$ and coloured by overall solid-volume fraction.

4.1 Mean settling velocity

In this section we present and discuss the settling velocity measurements $W_{s}$ , obtained ensemble averaging over all particles and realizations for each case. In figure 15(a) this is normalized by $W_{0}$ (so that values greater and smaller than one indicate turbulence-enhanced and turbulence-inhibited settling, respectively) and plotted against $St_{\unicode[STIX]{x1D702}}$ . The main contribution to the error is the uncertainty in $\unicode[STIX]{x1D70F}_{p}$ due to the particle size variance, and the non-zero vertical air velocity measured at the same time as the settling. We only plot cases in which the mean vertical fluid velocity is smaller than $0.25W_{s}$ , and in fact in most cases it is $0.01{-}0.05W_{s}$ . The first observation from the plot is that the vast majority of cases display strong settling enhancement, especially for $St_{\unicode[STIX]{x1D702}}\approx 1$ , which is consistent with most previous numerical (Wang & Maxey Reference Wang and Maxey1993; Bosse & Kleiser Reference Bosse and Kleiser2006; Dejoan & Monchaux Reference Dejoan and Monchaux2013; Bec et al. Reference Bec, Homann and Ray2014a ; Ireland et al. Reference Ireland, Bragg and Collins2016b ; Rosa et al. Reference Rosa, Parishani, Ayala and Wang2016) and experimental studies (Aliseda et al. Reference Aliseda, Cartellier, Hainux and Lasheras2002; Yang & Shy Reference Yang and Shy2003, Reference Yang and Shy2005; Good et al. Reference Good, Ireland, Bewley and Bodenschatz2014). The amount of such increase is more remarkable, with the settling velocity being enhanced by a factor 2.6 for $St_{\unicode[STIX]{x1D702}}\approx 1{-}2$ . As mentioned in the Introduction, most numerical studies reported maximum increase in fall speed between approximately 10 and 90 %; some experiments (Yang & Shy Reference Yang and Shy2003, Reference Yang and Shy2005) found even smaller values. The present results instead indicate that turbulence can lead to a multi-fold increase in settling rate, which agrees with the conclusions from the field study of Nemes et al. (Reference Nemes, Dasari, Hong, Guala and Coletti2017).

Similar levels of settling enhancement could also be deduced from the data of Aliseda et al. (Reference Aliseda, Cartellier, Hainux and Lasheras2002) and Good et al. (Reference Good, Ireland, Bewley and Bodenschatz2014); and while the former used concentrations where collective effects are expected ( $\unicode[STIX]{x1D719}_{v}\geqslant 10^{-5}$ ), the latter used particle loadings small enough to neglect two-way coupling ( $\unicode[STIX]{x1D719}_{v}\approx 10^{-6}$ ). In fact, those authors did not explicitly mention a multi-fold increase in vertical velocity, as they mostly plotted their data as $(W_{s}-W_{0})/u^{\prime }$ . We present this scaling in figure 15(b), where a maximum settling enhancement of $0.28u^{\prime }$ is found for $St_{\unicode[STIX]{x1D702}}\approx 1$ , again in good agreement with those authors. (Note that the vertical velocity is positive when downward, hence negative values imply settling enhancement and vice versa.) The scatter and the superposition of multiple factors prevent distinguishing a clear trend at the larger $St_{\unicode[STIX]{x1D702}}$ . Those data points are also at relatively high $\unicode[STIX]{x1D719}_{v}$ , which may have a non-trivial influence on settling, as we discuss later. Thus, the reduced settling exhibited by some of the most inertial cases, while it might appear consistent with recent results (Good et al. Reference Good, Ireland, Bewley and Bodenschatz2014; Rosa et al. Reference Rosa, Parishani, Ayala and Wang2016), should be considered with caution.

Figure 15.

(a) Measured settling velocity normalized by the still-fluid Stokes value versus $St_{\unicode[STIX]{x1D702}}$ and coloured by solid-volume fraction. (b) Measured settling velocity variation from the still-fluid Stokes value normalized by the fluid r.m.s. fluctuations.

Figure 16 shows the results in the $St_{\unicode[STIX]{x1D702}}{-}Sv_{\unicode[STIX]{x1D702}}$ plane. This provides a clearer view of the data, as both parameters are expected to have significant influence on the dynamics. The maximum enhancement of settling rate occurs when both $St_{\unicode[STIX]{x1D702}}$ and $Sv_{\unicode[STIX]{x1D702}}$ are close to unity, in broad agreement with Good et al. (Reference Good, Ireland, Bewley and Bodenschatz2014) and Rosa et al. (Reference Rosa, Parishani, Ayala and Wang2016). The distribution of values suggests that a dependence with $StSv$ may capture the observed trend. Figure 16(b) shows the settling enhancement ratio against the group $St_{\unicode[STIX]{x1D702}}Sv_{L}$ , displaying a significantly improved collapse of the data. This scaling follows the argument of Nemes et al. (Reference Nemes, Dasari, Hong, Guala and Coletti2017) that $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ and $u^{\prime }$ are the main time and velocity scales, respectively, determining increase of fall speed by turbulence. That both the small and large eddies impact the settling process has been acknowledged (Good et al. Reference Good, Ireland, Bewley and Bodenschatz2014), and already Wang & Maxey (Reference Wang and Maxey1993) favoured $u^{\prime }$ over $u_{\unicode[STIX]{x1D702}}$ as the driving parameter. Yang & Lei (Reference Yang and Lei1998) explicitly indicated $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ and $u^{\prime }$ as the correct flow scales, reasoning that the former controlled clustering and the latter controlled the drag experienced by the particles. The group $St_{\unicode[STIX]{x1D702}}Sv_{L}$ can be interpreted as the ratio of the particle stopping distance ( $\unicode[STIX]{x1D70F}_{p}^{2}g$ ) and a mixed length scale ( $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}u^{\prime }$ ); settling enhancement appears most effective when this ratio is $O(0.1)$ . This is approximately the condition at which Nemes et al. (Reference Nemes, Dasari, Hong, Guala and Coletti2017) reported turbulence-augmented fall speeds of snowflakes in the atmospheric surface layer ( $Re_{\unicode[STIX]{x1D706}}\approx 10^{3}$ ). Mixed-scaling arguments have been successfully used in various turbulent flows (e.g. in boundary layers, De Graff & Eaton Reference De Graff and Eaton2000) but their theoretical underpinning poses issues which are beyond the scope of the present study.

Figure 16.

(a) Settling increase plotted as a function of both $St_{\unicode[STIX]{x1D702}}$ and $Sv_{\unicode[STIX]{x1D702}}$ . (b) Settling increase as a function of the mixed scale, $St_{\unicode[STIX]{x1D702}}Sv_{L}$ .

Overall, the results presented in this section indicate that turbulence greatly enhances the settling velocity of sub-Kolmogorov particles with Stokes number around unity, which is consistent with the preferential sweeping mechanism proposed by Wang & Maxey (Reference Wang and Maxey1993). However, a full demonstration of this view requires the simultaneous measurements of particle and fluid velocities. These will be presented in § 5.2.

4.2 Settling velocity conditioned on particle concentration

To explore the interplay between the particle accumulation and settling mechanisms, we consider the fall speed associated with individual coherent clusters. In figure 17 we plot the cluster settling velocity $W_{C}$ , obtained by averaging the vertical velocity of all particles belonging to a given clustered set. This is normalized by the mean settling velocity $W_{s}$ , and plotted against the cluster area. Overall, clusters settle significantly faster than the mean, and there is an apparent trend of increasing fall speed with cluster size, especially for the larger objects. There are two possible interpretations for this result. On the one hand, clustered particles may affect the flow by virtue of their elevated concentration, exerting a collective drag on the surrounding fluid that results in increased settling velocity. This view reflects the argument proposed by Bosse & Kleiser (Reference Bosse and Kleiser2006) in interpreting their two-way-coupled DNS study. On the other hand, particles may be merely oversampling downward regions of flow according to the preferential sweeping mechanism, and therefore cluster in such regions, leading to the observed trend. This latter interpretation, which does not require any significant two-way coupling between the dispersed and continuous phase, is consistent with the results of Baker et al. (Reference Baker, Frankel, Mani and Coletti2017), who reported cluster fall speeds up to twice the mean particle settling velocity in their one-way-coupled DNS.

Figure 17.

Average settling velocity of all clusters as a function of cluster size. Error bars indicate $\pm$ the standard deviation.

Contrasting the effect of local and global concentration may provide further hints. In figure 18, the particle settling velocity of all particles $W_{s}$ (normalized by the still-air fall speed $W_{0}$ ) is plotted against the local relative concentration $C/C_{0}$ (which is readily available for each particle from the Voronoï diagrams). As expected, $W_{s}/W_{0}$ increases monotonically with $C/C_{0}$ , in agreement with the trends reported by Aliseda et al. (Reference Aliseda, Cartellier, Hainux and Lasheras2002). Indeed, particles residing in regions of low concentration are often associated with an upward velocity. We also observe, although with some scatter, the beginning of a plateau in the settling enhancement around $C/C_{0}\approx 5$ . More importantly, the plot clearly indicates that the fall speed dependence with concentration is strongly mitigated at larger global volume fractions, $\unicode[STIX]{x1D719}_{v}$ . If the high fall speed of the clusters was mainly due to a collective effect of the particles on the fluid, we would expect such a speed to be further enhanced with increasing $\unicode[STIX]{x1D719}_{v}$ . The fact that the opposite is true rather suggests that the augmented cluster settling is mainly caused by preferential sweeping (or other mechanisms not depending on the mass loading). In fact, figure 18 suggests that two-way coupling may be significant over the considered range of $\unicode[STIX]{x1D719}_{v}$ , but its effect may be subtle: if the particles are altering the turbulence structure, this backreaction can have a non-trivial effect on the settling rate. In general, it should be remarked that the simultaneous variation of multiple physical parameters between the considered cases (in this as in other studies) is a confounding factor in determining the role of two-way coupling, and one cannot rule out the influence of collective drag on the settling velocity (as argued by Huck et al. (Reference Huck, Bateson, Volk, Cartellier, Bourgoin and Aliseda2018)). Future dedicated studies, in which the global volume fraction is systematically varied while keeping all other parameters constant, may help shed light on this point.

Figure 18.

Normalized vertical particle velocity conditioned on the local particle concentration.

5.1 Preferential concentration

The fact that inertial particles oversample high-strain/low-vorticity regions, as theorized by Maxey (Reference Maxey1987) and demonstrated numerically by Squires & Eaton (Reference Squires and Eaton1991b ), was confirmed by several later DNS studies of homogeneous turbulence, at least for $St_{\unicode[STIX]{x1D702}}\leqslant 1$ (Chun et al. Reference Chun, Koch, Rani, Ahluwalia and Collins2005; Bec et al. Reference Bec, Biferale, Boffetta, Celani, Lanotte, Musaccio and Toschi2006; Cencini et al. Reference Cencini, Bec, Biferale, Celani, Lanotte, Musacchio and Toschi2006; Coleman & Vassilicos Reference Coleman and Vassilicos2009; Salazar & Collins Reference Salazar and Collins2012; Esmaily-Moghadam & Mani Reference Esmaily-Moghadam and Mani2016; Ireland et al. Reference Ireland, Bragg and Collins2016a ; Baker et al. Reference Baker, Frankel, Mani and Coletti2017). To our knowledge, this prediction has not been directly verified by experiments in fully turbulent flows. Indeed, most previous laboratory studies on this topic only captured the dispersed phase (Fessler et al. Reference Fessler, Kulick and Eaton1994; Aliseda et al. Reference Aliseda, Cartellier, Hainux and Lasheras2002; Wood et al. Reference Wood, Hwang and Eaton2005; Salazar et al. Reference Salazar, de Jong, Cao and Woodward2008; Saw et al. Reference Saw, Shaw, Ayyalasomayajula, Chuang and Gylfason2008; Gibert, Xu & Bodenschatz Reference Gibert, Xu and Bodenschatz2012), and as such could only provide results consistent with a certain picture of preferential concentration, rather than demonstrating it. We characterize the local balance of strain-rate versus rotation in the particle-laden air flow measured by PIV, using the second invariant of the velocity gradient tensor $Q=1/2(\unicode[STIX]{x1D734}^{2}-\boldsymbol{S}^{2})$ , where $\boldsymbol{S}$ and $\unicode[STIX]{x1D734}$ are the symmetric and anti-symmetric parts of the velocity gradient tensor (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988). To this end, we calculate spatial derivatives using a second-order central difference scheme on the small and medium-FOV fields, where our resolution is sufficient to capture the Kolmogorov scales (Worth, Nickels & Swaminathan Reference Worth, Nickels and Swaminathan2010; Hearst et al. Reference Hearst, Buxton, Ganapathisubramani and Lavoie2012). From the planar data we can only determine the four components in the upper-left $2\times 2$ block of the full $3\times 3$ velocity gradient tensor. This limitation needs to be kept in mind, because 2-D sections of 3-D flows can sometimes be misleading (Perry & Chong Reference Perry and Chong1994). However, several studies showed how high-resolution 2-D imaging of homogeneous turbulence yields features of the coherent structures and high-order statistics in quantitative agreement with 3-D imaging and DNS (Fiscaletti, Westerweel & Elsinga Reference Fiscaletti, Westerweel and Elsinga2014; Carter & Coletti Reference Carter and Coletti2018; Saw et al. Reference Saw, Debue, Kuzzay, Daviaud and Dubrulle2018). Therefore, we do not expect the qualitative results of the present analysis to be biased by the nature of the measurements.

Figure 19 shows the fraction of inertial particles found in regions where $Q<0$ . As expected, this fraction is larger than 50 % for all cases, confirming that the particles are more likely to be found in strain-dominated regions than rotation-dominated ones. The figure also presents the percentage of clustered particles found in $Q<0$ regions. Interestingly, the fraction is systematically lower compared to the entire particle set. This suggests that the preferential sampling of high-strain regions might not be the main factor (or at least not the only one) for the formation of clusters over the considered parameter space. Indeed, most of our cases feature particles with $St_{\unicode[STIX]{x1D702}}>1$ , and in this regime several numerical studies indicate that the nature of the clustering mechanism is different compared to weakly inertial particles (Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Masucchio and Toschi2007; Coleman & Vassilicos Reference Coleman and Vassilicos2009; Bragg & Collins Reference Bragg and Collins2014; Bragg et al. Reference Bragg, Ireland and Collins2015). In particular, Bragg & Collins (Reference Bragg and Collins2014) argued for the importance of path-history effects (i.e. particles retaining memory of the velocity fluctuations they experienced), while Vassilicos and coworkers (Chen, Goto & Vassilicos Reference Chen, Goto and Vassilicos2006; Goto & Vassilicos Reference Goto and Vassilicos2008) proposed that clustering in this range is due to a sweep-stick mechanism (i.e. particles sticking to zero-acceleration points which are swept and clustered by large-scale motions). A critical discussion of these and other possible explanations is beyond the scope of this work. In fact, while instantaneous realizations and velocity statistics may provide support to a given theory (see Obligado et al. Reference Obligado, Teitelbaum, Cartellier, Mininni and Bourgoin2014; Sumbekova et al. Reference Sumbekova, Cartellier, Aliseda and Bourgoin2017), time-resolved measurements would be better suited to inform a mechanistic understanding of the process.

Figure 19.

Fraction of particles in rotation-dominated regions. Black circles represent all particles in the field, while the purple diamonds represent only particles in coherent clusters.

5.2 Preferential sweeping

As discussed in §§ 1 and 4.1, preferential sweeping is considered the most impactful mechanism by which turbulence affects the fall speed of sub-Kolmogorov particles. Its main manifestation is the tendency of particles with Stokes number of order one to oversample regions of downward velocity fluctuations. This was first theorized by Maxey & Corrsin (Reference Maxey and Corrsin1986) and Maxey (Reference Maxey1987), demonstrated numerically by Wang & Maxey (Reference Wang and Maxey1993) and confirmed by several other analytical and computational studies (Yang & Lei Reference Yang and Lei1998; Dávila & Hunt Reference Dávila and Hunt2001; Dejoan & Monchaux Reference Dejoan and Monchaux2013; Frankel et al. Reference Frankel, Pouransari, Coletti and Mani2016; Baker et al. Reference Baker, Frankel, Mani and Coletti2017). While laboratory studies (Aliseda et al. Reference Aliseda, Cartellier, Hainux and Lasheras2002; Yang & Shy Reference Yang and Shy2003, Reference Yang and Shy2005; Good et al. Reference Good, Ireland, Bewley and Bodenschatz2014) and field observations (Nemes et al. Reference Nemes, Dasari, Hong, Guala and Coletti2017) showed results consistent with this picture, no direct experimental verification has been reported. Similar to preferential concentration, the challenges associated with two-phase measurements may be responsible.

Figure 20.

Example case of the relative particle concentration conditioned on the vertical fluid velocity (normalized by the still-air settling velocity.) Black circles indicate the full particle set, while purple diamonds represent only particles in clusters.

We provide such verification first by considering the particle concentration conditionally averaged on the local fluid velocity. This is obtained by counting the number of inertial particles in each PIV interrogation window, and binning the results by the value of $W_{f}/W_{0}$ (because the mean vertical fluid velocity is negligibly small, total and fluctuating components coincide). The relative concentration is calculated as the number of particles in each bin, divided by the sum of window areas associated with that bin and finally normalized by the global concentration. The procedure is equivalent to that originally adopted by Wang & Maxey (Reference Wang and Maxey1993) and later by Baker et al. (Reference Baker, Frankel, Mani and Coletti2017) to analyse DNS results, with the PIV interrogation windows in lieu of the computational cells. In figure 20 we show the result for a representative case, clearly indicating downward fluid velocity corresponding to higher local concentration. When the process is repeated only considering particles belonging to clusters, the trend is significantly more pronounced. This is consistent with the result that clusters fall at faster speeds than the rest of the particles (figure 17). At the same time, it also supports the idea that preferential sweeping plays an important role in the clustering of settling particles.

To quantify the impact on the settling rate, we consider the vertical component of the fluid velocity at the particle location, $W_{f}(\boldsymbol{x}_{\boldsymbol{p}})$ (figure 21). The latter is approximated via a piecewise linear interpolant between the particle position and the four closest fluid velocity vectors; tests with other schemes indicate only a weak dependence with the interpolation method. Error analysis based on the fluid velocity gradient statistics (see Carter et al. Reference Carter, Petersen, Amili and Coletti2016; Carter & Coletti Reference Carter and Coletti2017) yields uncertainty on $W_{f}(\boldsymbol{x}_{\boldsymbol{p}})$ around 2–7 % $u^{\prime }$ . Despite some scatter partly attributable to the several factors at play, the results indicate that preferential sweeping is important for most considered regimes, being the strongest for $St_{\unicode[STIX]{x1D702}}=O(1)$ .

Figure 21.

Vertical fluid velocity at particle location, normalized by the fluid r.m.s. velocity, plotted versus particle Stokes number.

Comparing figures 15(b) and 21, the oversampling of downward fluid velocity regions seems to account for a large part of the settling enhancement. A more quantitative account can be given in the framework of the point-particle approximation. Retaining only drag and gravity in the particle equation of motion, the fall speed can be approximated as (Wang & Maxey Reference Wang and Maxey1993):

(5.1) $$\begin{eqnarray}\langle W_{s}\rangle \approx \langle \boldsymbol{u}_{\boldsymbol{f}}(\boldsymbol{x}_{\boldsymbol{p}},t)\rangle +\frac{\unicode[STIX]{x1D70F}_{p}g}{\langle f\rangle },\end{eqnarray}$$

where $\boldsymbol{u}_{\boldsymbol{f}}(\boldsymbol{x}_{\boldsymbol{p}})$ is the fluid velocity vector at the particle location obtained via the piecewise-linear interpolation, and $\langle f\rangle$ is the ensemble average of the Schiller & Naumann correction factor in (2.1), $f=1+0.15Re_{p}^{0.687}$ . We can directly verify this approximation using the instantaneous $Re_{p}$ measured from the simultaneous PIV/PTV measurements (which will be discussed further in the next section). Figure 22 shows the ratio between the fall speed calculated from (5.1) and the measured values. This formulation consistently underpredicts the measured fall speed. Such a discrepancy between experiments and theory suggests that the one-way-coupled point-particle approach, while providing the correct qualitative trend, is missing significant aspects of the particle–fluid interaction. We investigate possible sources of the mismatch in § 5.4, where we consider the instantaneous slip velocity.

Figure 22.

Prediction of settling velocity based on (5.1) normalized by the measured settling velocities, versus particle Stokes number.

5.3 Crossing trajectory and continuity effects

Since Yudine (Reference Yudine1959), it has been recognized that the drift induced by body forces such as gravity causes heavy particles to decorrelate from their past velocity faster than a fluid element. This so-called crossing trajectory effect can be quantified by the Lagrangian autocorrelation of the particle velocity (Elghobashi & Truesdell Reference Elghobashi and Truesdell1992). For large drift velocities, this reduces to the fluid space–time correlation in an Eulerian frame (Csanady Reference Csanady1963; Squires & Eaton Reference Squires and Eaton1991a ). In this limit, as the longitudinal integral scale is twice the transverse one, the particle dispersion parallel to the drift direction is double the dispersion perpendicular to it (the so-called continuity effect; Csanady Reference Csanady1963; Wang & Stock Reference Wang and Stock1993). The footprints of these effects are visible in the Eulerian particle velocities. In figure 23 we present a scatter plot of the vertical and horizontal r.m.s. particle velocity fluctuations ( $W_{p,rms}$ and $U_{p,rms}$ , respectively), normalized by the r.m.s. fluid fluctuations in the respective directions ( $u_{z}^{\prime }$ and $u_{x}^{\prime }$ ) to account for the anisotropy in our facility. The normalized vertical fluctuations of the particles exceed those in the horizontal direction, the disparity being more substantial for larger $Sv_{\unicode[STIX]{x1D702}}$ . This trend is consistent with the continuity effect, and in line with previous analysis by Wang & Stock (Reference Wang and Stock1993) and measurements by Good et al. (Reference Good, Ireland, Bewley and Bodenschatz2014): the falling particles have more time to respond to the vertical fluid fluctuations, due to the larger longitudinal integral scale compared to the transverse one.

Figure 23.

Ratio of vertical particle to fluid r.m.s. velocity versus ratio of the same horizontal velocities.

5.4 Particle–fluid relative velocity

The relative (slip) velocity between suspended particles and the surrounding fluid plays a vital role in both the particle kinematics and the interphase coupling. The slip velocity determines the exchanged forces and the local flow regime around a particle (Maxey & Riley Reference Maxey and Riley1983; Calzavarini et al. Reference Calzavarini, Kerscher, Lohse and Toschi2008; Bellani et al. Reference Bellani, Byron, Collignon, Meyer and Variano2012), and is critical for the momentum two-way coupling with the fluid (Sahu et al. Reference Sahu, Hardalupas and Taylor2016). It is also a key quantity for modelling sub-grid-scale dynamics in large-eddy simulations of particle-laden turbulence (Soldati & Marchioli Reference Soldati and Marchioli2009). Due to the abovementioned difficulties in measuring both phases, the slip velocity has been reported by only a few experimental studies, notably Kiger & Pan (Reference Kiger and Pan2002) and Khalitov & Longmire (Reference Khalitov and Longmire2003) in wall-bounded turbulent flows, and Yang & Shy (Reference Yang and Shy2005) and Sahu et al. (Reference Sahu, Hardalupas and Taylor2016) in homogeneous turbulence.

Here we evaluate the slip velocity as $\boldsymbol{u}_{slip}=\boldsymbol{v}_{\boldsymbol{p}}-\boldsymbol{u}_{\boldsymbol{f}}(\boldsymbol{x}_{\boldsymbol{p}})$ , where $\boldsymbol{v}_{\boldsymbol{p}}$ is the particle velocity vector. In the present measurements we only have access to the in-plane projection of those vectors. Here we only use data from the small and intermediate FOVs. We remark that the relevant definition of slip velocity (as it appears, e.g. in Stokes drag law) is the difference between the particle velocity and the undisturbed fluid velocity evaluated at the particle location. Because in the considered cases $d_{p}<\unicode[STIX]{x1D702}$ and $Re_{p,0}\leqslant 1$ , the region of perturbation from an individual particle (as estimated, e.g. by the Oseen solution) is expected to be a few particle diameters, i.e. typically smaller than the Kolmogorov scale. Under this assumption, the interpolated fluid velocity approximates the undisturbed velocity. However, as we shall see, the instantaneous particle Reynolds number $Re_{p}$ can reach much larger values than $Re_{p,0}$ .

Figure 24(a) presents PDFs of the slip velocity magnitude normalized by the fluid r.m.s. fluctuation. At the higher $St_{\unicode[STIX]{x1D702}}$ , $\boldsymbol{u}_{slip}$ can significantly exceed $u^{\prime }$ , suggesting that the potential energy transferred to the fluid by the fast falling particles is considerable with respect to the turbulent kinetic energy. We will return to this point in § 5.5. In figure 24(b) we plot the PDF of the vertical component of the slip velocity $W_{slip}$ , normalized by the still-air particle settling velocity $W_{0}$ . The mean vertical slip is approximately equal to $W_{0}$ , as one expects from (5.1). However, the standard deviation is large and, especially for the lower $St_{\unicode[STIX]{x1D702}}$ , greatly exceeds the mean value; indeed, the probability of an upward particle slip is significant. This behaviour is consistent with the simulations of particle-laden wall-bounded turbulence by Zhao, Marchioli & Andersson (Reference Zhao, Marchioli and Andersson2012), who found that the r.m.s. fluctuations of both streamwise and wall-normal slip velocities were several times greater than their mean. The r.m.s. slip velocities are presented in figure 24(c,d) separating horizontal and vertical components and normalizing by the r.m.s. fluid velocity fluctuations in their respective directions. The slip velocity fluctuations increase with $St$ , also consistent with the results of Zhao et al. (Reference Zhao, Marchioli and Andersson2012): the heavier particles increasingly move independently from the fluid, and their r.m.s. slip velocities reach levels comparable to the r.m.s fluid fluctuations. The vertical r.m.s. is consistently higher than the horizontal. This is likely a consequence of the vertical r.m.s. velocity of the particles being larger than the horizontal one, as discussed in the previous section.

Figure 24.

Probability density functions of slip velocity magnitude normalized by fluid r.m.s. velocity (a), and vertical slip velocity (b) normalized by the still-fluid settling. (c,d) Plots of horizontal and vertical r.m.s. particle slip velocity, respectively, normalized by the fluid r.m.s. fluctuations along the corresponding direction.

The magnitude of the in-plane slip velocity can be used to calculate the instantaneous particle Reynolds number. The PDFs of $Re_{p}$ for various cases are plotted in figure 25(a), with the vertical lines indicating the respective theoretical $Re_{p,0}$ . The large slip variance results in long tails of the distributions: while the wake-shedding regime is never achieved, there is a sizeable probability of $Re_{p}$ being an order of magnitude larger or smaller than the nominal value. The values are in fact somewhat underestimated, because the out-of-plane slip is not accounted for.

The slip velocity and Reynolds number can be used to evaluate the drag coefficient on the settling particles. In the presence of background turbulence, this is expected to differ from the standard drag in a steady uniform flow (Bagchi & Balachandar Reference Bagchi and Balachandar2003). The latter can be estimated following the Schiller & Naumann correction for finite Reynolds number:

(5.2) $$\begin{eqnarray}C_{d,SN}=\frac{24}{Re_{p,0}}(1+0.15Re_{p,0}^{0.687}).\end{eqnarray}$$

As suggested by Balachandar & Eaton (Reference Balachandar and Eaton2010), for particles settling in zero-mean-flow homogeneous turbulence, one can evaluate an effective drag coefficient:

(5.3) $$\begin{eqnarray}C_{d,eff}=\frac{24}{\langle Re_{v}\rangle }\frac{1}{\langle W_{p}\rangle }\left|\langle \boldsymbol{u}_{slip}\rangle +0.15\left(\frac{d_{p}}{\unicode[STIX]{x1D708}}\right)^{0.687}\langle \boldsymbol{u}_{slip}^{1.687}\rangle \right|,\end{eqnarray}$$

where $Re_{v}=d_{p}|\boldsymbol{v}_{\boldsymbol{p}}|/\unicode[STIX]{x1D708}$ . This may differ from the $C_{d,SN}$ owing to three main factors: the mean fluid velocity seen by the particles being non-zero (e.g. due to preferential sweeping); the nonlinear relation between drag and slip velocity; and possible two-way coupling. Most previous studies concerned with the influence of ambient turbulence on $C_{d}$ considered larger $Re_{p}$ compared to the present case (typically of order $10^{2}{-}10^{3}$ ; e.g. Wu & Faeth Reference Wu and Faeth1994; Warnica, Renksizbulut & Strong Reference Warnica, Renksizbulut and Strong1995; Bagchi & Balachandar Reference Bagchi and Balachandar2003). In these cases, the effect on the mean drag was found to be small. Bagchi & Balachandar (Reference Bagchi and Balachandar2003), however, recognized the potential importance of preferential flow sampling for particles falling freely through turbulence. Here we use (5.2) and (5.3) and plot the ratio $C_{d,eff}/C_{d,SN}$ in figure 25(b). There is a significant reduction in effective drag over the $St_{\unicode[STIX]{x1D702}}$ range displaying preferential sweeping. This is consistent with the early measurements of Rudoff & Bachalo (Reference Rudoff and Bachalo1988), who found substantially reduced drag for liquid droplets in turbulent air. We also observe a clear increasing trend with $St_{\unicode[STIX]{x1D702}}$ . This is likely a consequence of diminishing preferential sweeping effects, giving way to those of nonlinear drag and possibly loitering (Good et al. Reference Good, Ireland, Bewley and Bodenschatz2014). We remark, however, that the high- $St$ data points are affected by larger uncertainty due to two factors: non-negligible perturbation of the local fluid velocity by the particles (because of their relatively large $Re_{p}$ ), and limited spatial resolution (because those are imaged with the intermediate FOV). Thus, further analysis and measurements are warranted for those regimes.

Figure 25.

PDFs of instantaneous particle Reynolds number. Vertical lines represent average Reynolds particle number based on the Stokes still-fluid settling velocity for each respective grouping of experiments. Panel (b) shows the ratios between drag coefficients (5.2) and (5.3), showing reduced drag when the ‘effective drag’ model is used.

We then consider the spatial orientation of the particle slip velocity. Figure 26(a) shows the PDF of cos( $\unicode[STIX]{x1D703}$ ), $\unicode[STIX]{x1D703}$ being the angle between $\boldsymbol{v}_{\boldsymbol{p}}$ and $\boldsymbol{u}_{\boldsymbol{f}}(\boldsymbol{x}_{\boldsymbol{p}})$ . For moderate particle inertia ( $St_{\unicode[STIX]{x1D702}}<10$ ) the fluid and particle velocity vectors tend to be closely aligned, with $\unicode[STIX]{x1D703}<30^{\circ }$ in more than 90 % of the instances. As expected, the more inertial particles are more likely to display a velocity orientation substantially different from the fluid velocity. The probability of oppositely aligned vectors is also non-negligible, possibly due to weak upward gusts that are unable to reverse the particle settling motion. The probability distribution of slip orientation and magnitude is illustrated by the wind rose diagrams in figure 26(b) for selected cases. The length of each spoke represents the probability associated with its orientation; additionally, each spoke is broken down into segments representing the probability of a certain slip velocity magnitude (normalized by $W_{0}$ ). In all cases the slip velocity favours the downward direction, as expected. In the larger $Re_{\unicode[STIX]{x1D706}}$ cases, however, the strong fluid velocity fluctuations mitigate this tendency and the azimuthal distribution of the particle slip is less skewed. As discussed above, particles with $St_{\unicode[STIX]{x1D702}}\approx 1$ can have slip velocities several times larger than their still-air fall speed. On the other hand, $\boldsymbol{u}_{slip}$ for the more inertial cases is closer in magnitude to $W_{0}$ and mostly oriented downward (typically within $\pm 45^{\circ }$ from the vertical), confirming that these particles follow more ballistic downward trajectories.

Figure 26.

PDF of the cosine of the angle between the particle velocity and the fluid velocity at the particle location (a). (b) Three example cases showing the probability distribution of the slip velocity orientation and magnitude.

5.5 Turbulence modification by particles

We finally consider the scale-by-scale effect of the particles on the fluid turbulent fluctuations. We consider the second-order structure functions, which contain analogous information as the energy spectra used in most previous studies:

(5.4) $$\begin{eqnarray}D_{ii}(\boldsymbol{r})=\overline{[u_{i}(\boldsymbol{x}+\boldsymbol{r})-u_{i}(\boldsymbol{x})]^{2}}.\end{eqnarray}$$

Here $\boldsymbol{x}$ is the generic location on the imaging plane, $\boldsymbol{r}$ is the separation vector and $u_{i}$ is the $i$ th component of the velocity fluctuations. Thus, we denote as $D_{uu}$ and $D_{ww}$ the structure functions associated with the horizontal and vertical components. We focus on longitudinal structure functions (i.e. with separations parallel to the velocity components), the transverse structure functions yielding similar results. For this analysis we use only small and intermediate FOVs, in order to resolve a sufficient fraction of the fine-scale fluctuations.

Figure 27 displays horizontal and vertical structure functions for representative cases, comparing the laden and unladen measurements. At very low volume fractions ( $\unicode[STIX]{x1D719}_{v}=O(10^{-7})$ ), the particles marginally affect the energy distribution across scales, and the expected scaling $D_{ii}(r)\sim r^{2/3}$ is approached entering the inertial range according to Kolmogorov theory (figure 27 a). At significant loadings ( $\unicode[STIX]{x1D719}_{v}=O(10^{-5})$ ), we observe an increase of turbulent energy at small scales. For $St_{\unicode[STIX]{x1D702}}=O(1)$ (figure 27 b) the laden and unladen curves approach each other at larger scales, possibly leading to a cross-over which however is not captured within the FOV. At $St_{\unicode[STIX]{x1D702}}=O(10)$ , the cross-over happens at smaller scales as demonstrated in figure 27(d), in which particles appear to excite turbulent fluctuations at the small scales and modulate them at the large scales. This behaviour (sometimes termed pivoting) was reported by several numerical studies, as reviewed in detail by Poelma & Ooms (Reference Poelma and Ooms2006). Most of those (e.g. Sundaram & Collins Reference Sundaram and Collins1999; Ferrante & Elghobashi Reference Ferrante and Elghobashi2003) showed the cross-over wavenumber to increase with increasing particle response time, consistent with our observations. The only previous experimental study to clearly demonstrate the pivoting effect was, to our best knowledge, the one from Poelma et al. (Reference Poelma, Westerweel and Ooms2007). However, a direct comparison is difficult as these authors considered spatially decaying turbulence at maximum $Re_{\unicode[STIX]{x1D706}}\approx 29$ .

Figure 27.

Horizontal (a,b,d) and vertical (c,e) longitudinal structure functions for laden (colour coded by volume fraction) and unladen (black) multiphase flows. The shaded regions represent 95 % bootstrap confidence bounds.

The vertical velocity structure functions (figure 27 c,e) display the same behaviour as their horizontal counterpart (figure 27 b,d), suggesting that the turbulence responds to the presence of the particles in similar ways in both directions. It is indeed verified that the anisotropy ratio $u_{x}^{\prime }/u_{z}^{\prime }$ is not significantly affected by the particles. Poelma et al. (Reference Poelma, Westerweel and Ooms2007) and Frankel et al. (Reference Frankel, Pouransari, Coletti and Mani2016) found that particles enhanced the fluid velocity fluctuations aligned with gravity compared to the transverse ones, whereas in Hwang & Eaton (Reference Hwang and Eaton2006b ) the level of isotropy was marginally influenced even at particle mass loadings that heavily modulated the turbulence intensity. Among these studies, Hwang & Eaton considered a much higher Reynolds number ( $Re_{\unicode[STIX]{x1D706}}\approx 240$ ), comparable to the present study. One may conjecture that the wide range of scales allows for a redistribution of the energy among the different components; certainly, further studies that systematically investigate the role of $Re_{\unicode[STIX]{x1D706}}$ are needed.

Even for the higher volume fractions, the changes in the fluid turbulence intensities in both vertical and horizontal directions (not presented) are scattered within the measurement uncertainty, with no consistent trend in terms of $St_{\unicode[STIX]{x1D702}}$ , $Sv_{\unicode[STIX]{x1D702}}$ or $\unicode[STIX]{x1D719}_{v}$ . This suggests that the turbulence augmentation/attenuation, over different scales and from different mechanisms, largely compensate each other, at least in the considered range of parameters. On one hand, particles increase the overall inertia of the mixture and locally enhance the dissipation around them, modulating the turbulence; on the other hand, the drag force exerted by the falling particles perturbs the fluid an excites the turbulence (Balachandar & Eaton Reference Balachandar and Eaton2010). Both effects are amplified for larger and heavier particles (which provide higher loading, distort more the flow around them and shed more energetic wakes). For comparison, Hwang & Eaton (Reference Hwang and Eaton2006a ) investigated particles with $St_{\unicode[STIX]{x1D702}}\approx 50$ settling in homogeneous turbulence and measured a reduction of r.m.s. fluid velocity of approximately 10–15 % for a mass fraction of 0.1 (close to the highest loading considered in the present measurements). The balance between opposing effects was demonstrated by measurements obtained in the same facility in micro-gravity (Hwang & Eaton Reference Hwang and Eaton2006b ), where the turbulence modulation was found to be greater than in the fixed laboratory frame. Hwang & Eaton (Reference Hwang and Eaton2006b ) also compared their laboratory results against DNS studies which showed significantly smaller turbulence attenuation. While this suggested that the point-particle approach used in the simulations missed important physics, the DNS were also at significantly smaller $St_{\unicode[STIX]{x1D702}}$ , between 1 and 11 and comparable to the present cases.