4.1 Cluster topology

In order to characterize the general features of the particle distribution in each cluster, we use singular value decomposition (SVD) to find the principal axes and corresponding singular values for the particle set in each coherent cluster. SVD is a natural method to achieve a low-order description of a set of points in multi-dimensional spaces, the physical three-dimensional space being a clear example. By definition, the first principal axis of a set lies along the direction of the greatest spread of particle distances from the cluster centroid; the second principal axis lies along the direction of greatest spread within a plane normal to the first axis; the third principal axis is orthogonal to the first two and lies in the direction of minimum spread. The corresponding singular values are a measure of particle spread along each principal axis; specifically, a singular value normalized by the number of particles in the cluster gives the variance of the particle distances from the cluster centroid projected along that axis. An example of the principal axes of a coherent cluster is shown in figure 7, in which the length of each vector arrow is proportional to the corresponding singular value. The complex spatial distribution of the clustered particles is apparent.

Figure 7.

Principal axes of a cluster found by singular value decomposition. The length of each vector arrow is proportional to the corresponding singular value.

The singular values $s_{1}\geqslant s_{2}\geqslant s_{3}$ are used here to define two aspect ratios for each cluster, $s_{2}/s_{1}$ and $s_{3}/s_{2}$ . These can be interpreted by considering the limiting ( $s_{2}/s_{1}$ , $s_{3}/s_{2}$ ) cases: (0, 0) corresponds to a thin ribbon; (1, 0) corresponds to a thin sheet; (0, 1) corresponds to a thin rod; and (1, 1) corresponds to a sphere. The joint PDF of $s_{2}/s_{1}$ and $s_{3}/s_{2}$ is shown figure 8. For all considered cases the most likely combination of aspect ratios is about (0.5, 0.5), which corresponds to proportions 4 : 2 : 1. The cases with gravity exhibit somewhat lower aspect ratios, which reflects the tendency to form elongated streaks. The local maximum at very low values of $s_{3}/s_{2}$ present in some of the cases is due to small clusters with relatively low particle counts, which tend to be thin in one dimension.

Figure 8.

Joint PDFs of aspect ratios $s_{2}/s_{1}$ and $s_{3}/s_{2}$ of the coherent clusters, based on singular value decomposition. Labels in the first panel indicate the topology in the extreme cases in which the aspect ratios are either 0 or 1.

This analysis is similar to the one carried out by Moisy & Jiménez (Reference Moisy and Jiménez2004) for turbulent flow structures characterized by high vorticity and high strain. However, these authors used a different set of length scales: for each structure of volume V they defined three lengths scales $r_{1}\leqslant r_{2}\leqslant r_{3}$ such that $V=r_{1}r_{2}r_{3}$ ; the outer scale $r_{3}$ is the size of the smallest embedding box that fully contains the cluster, the inner scale $r_{1}$ is the size of the largest box that can be fully embedded in the cluster, and the intermediate scale $r_{2}$ is computed as $V/(r_{1}r_{3})$ . The set of aspect ratios $r_{1}/r_{2}$ and $r_{2}/r_{3}$ has the same geometrical interpretation as $s_{3}/s_{2}$ and $s_{2}/s_{1}$ . We also applied this method, and the average values $\langle r_{1}/r_{2}\rangle$ and $\langle r_{2}/r_{3}\rangle$ are shown in figure 9. The inner scale $r_{1}$ is computed from a variation of the box-counting method, and is reliably obtained only for relatively large clusters. We therefore only consider clusters with $V_{C}>100\unicode[STIX]{x1D702}^{3}$ . There is not a large spread between cases, but compared to the SVD-based aspect ratios, $\langle r_{1}/r_{2}\rangle$ and $\langle r_{2}/r_{3}\rangle$ indicate more tubular and generally thinner shapes.

Figure 9.

Mean aspect ratios $\langle r_{1}/r_{2}\rangle$ and $\langle r_{2}/r_{3}\rangle$ of coherent clusters with a volume larger than $100\unicode[STIX]{x1D702}^{3}$ , based on embedding and embedded boxes. The topologies are indicated in the extreme cases in which the aspect ratios are either 0 or 1.

We remark that the method of Moisy & Jiménez (Reference Moisy and Jiménez2004) does not take into account the interior structure of the particle distribution within the cluster: as the characteristic lengths are quantified by embedding and embedded boxes, $r_{3}$ and $r_{1}$ represent estimates of the ‘maximum length’ and ‘minimum waist’ of the clusters, respectively. The SVD method on the other hand uses the entire particle set in each cluster. While both methods provide useful indications of the clusters’ aspect ratio, neither of them can provide a full description of their shape, because the complex fractal structure of the particle distribution is not accounted for. This is discussed in § 4.3.

4.2 Particle concentration and distribution

To determine how the particle concentration varies with cluster size, scatter plots of cluster volumes $V_{C}$ against particle counts $N_{PC}$ are displayed in figure 10. The concentration within a cluster is given by $C_{C}=N_{PC}/V_{C}$ , so the mean concentration is given by

The exponent of the fitting power law (red line in figure 10) is close to one, indicating that the volume and particle count of clusters are nearly linearly related, and therefore the mean concentration is approximately constant throughout the range of cluster sizes. We note that, for the small clusters, there is significant scatter around the linear trend. This is partly due to the fact that, when the number of particles in a cluster is relatively small, the calculated volume (which is in fact the volume of the set of its Voronoï cells) is a less accurate approximation of the cluster size. The value of the mean concentration for each St, $Sv$ case is given in table 4 along with the mean volume of the coherent clusters. The concentration in the clusters is about one order of magnitude higher than the domain average, with the highest values found for St close to unity, as expected. For $St<1$ , gravitational settling somewhat reduces the in-cluster concentration (the change for $St=0.18$ being marginal); for $St>1$ , the opposite is true, and settling increases the in-cluster concentration. Such a trend was found in multiple recent studies (e.g. Bec, Homann & Ray Reference Bec, Homann and Ray2014; Gustavsson, Vajedi & Mehlig Reference Gustavsson, Vajedi and Mehlig2014; Ireland et al. Reference Ireland, Bragg and Collins2016; Matsuda, Onishi & Takahashi Reference Matsuda, Onishi and Takahashi2017). The consensus is that strongly inertial and fast-falling particles become uncorrelated from the turbulent flow structures, but the mechanism for enhanced clustering in these condition is still debated. The mean cluster volume generally increases with St and $Sv$ , confirming the trend visible in figures 3 and 5. Finding larger clusters at higher St is expected, because more inertial particles respond to larger eddies (Goto & Vassilicos Reference Goto and Vassilicos2006; Yoshimoto & Goto Reference Yoshimoto and Goto2007). On the other hand, the trend with $Sv$ suggests the presence of a clustering mechanism other than the classic centrifugal effect, such as the sweep-stick mechanism (Goto & Vassilicos Reference Goto and Vassilicos2008) and/or the path-history effect (Bragg & Collins Reference Bragg and Collins2014), as we discuss in § 5.

Figure 10.

Cluster volume normalized by the Kolmogorov length scale versus cluster particle counts normalized by the total number of particles in the domain.

Table 4.

Mean particle concentration within clusters normalized by the global concentration, and mean cluster volume normalized by the Kolmogorov length scale.

To characterize the nature and spatial extent of particle clustering, we can use the radial distribution function (RDF):

where $P_{i}$ is the number of pairs of particles separated by a distance $r_{i}\pm \unicode[STIX]{x0394}r/2$ , $v_{i}=4\unicode[STIX]{x03C0}[(r_{i}+\unicode[STIX]{x0394}r/2)^{3}-(r_{i}-\unicode[STIX]{x0394}r/2)^{3}]/3$ is the volume of the spherical shell of mean radius $r_{i}$ , and $P=N_{P}(N_{P}-1)/2$ is the total number of particle pairs in the simulation domain $v$ . In order to study anisotropic effects due to gravity, the angular distribution function (ADF, Gualtieri et al. Reference Gualtieri, Picano and Casciola2009) can be derived by binning the particle pairs using both the distance and the angle spanned by the inter-particle separation vector and the vertical:

Here $P_{i,j}$ is the number of particle pairs separated by distance $r_{i}\pm \unicode[STIX]{x0394}r/2$ and whose displacement vector has an orientation $\unicode[STIX]{x1D703}_{j}\pm \unicode[STIX]{x0394}\unicode[STIX]{x1D703}/2$ relative to the vertical, and $v_{i,j}$ is the volume of the part of the $i$ th spherical shell contained in the $j$ th angular bin. A power-law behaviour $g(r)\propto r^{-\unicode[STIX]{x1D6FC}}$ of the RDF or ADF indicates self-similar clustering. The value of $\unicode[STIX]{x1D6FC}$ is used to find the correlation dimension $D_{2}=3-\unicode[STIX]{x1D6FC}$ , which characterizes the clustering behaviour of a particle distribution (see for example Fessler, Kulick & Eaton Reference Fessler, Kulick and Eaton1994; Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Gualtieri et al. Reference Gualtieri, Picano and Casciola2009). A steeper slope (i.e. a lower value of $D_{2}$ ) implies a higher degree of local accumulation.

We first consider the RDF for the full set of P-I particles. Figure 11 presents the dependence of $D_{2}$ (calculated from a power-law fit over the range $0.2<r/\unicode[STIX]{x1D702}<1.5$ ) with Stokes number. The normalized standard deviation of the PDF of Voronoï cell size $\unicode[STIX]{x1D70E}_{V}/\langle V\rangle$ is also plotted. For the cases without gravity, $D_{2}$ shows a minimum and $\unicode[STIX]{x1D70E}_{V}/\langle V\rangle$ a maximum around $St\approx 1$ . This is expected, as both quantities are regarded as measures of clustering (Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2010). For the cases with gravitational settling, $D_{2}$ has a similar (although less pronounced) behaviour, while $\unicode[STIX]{x1D70E}_{V}/\langle V\rangle$ increases with St. This suggests that the standard deviation of the Voronoï volumes PDF should be interpreted with caution. In fact, besides the local concentration in the clusters, the latter is also a function of the size of the voids, which in turn depends in a non-trivial way on the number and spatial distribution of the clusters.

Figure 11.

Normalized standard deviation of the Voronoï cell size $\unicode[STIX]{x1D70E}_{V}/\langle V\rangle$ and correlation dimension $D_{2}$ for the full set of P-I particles versus Stokes number.

To assess the effect of the cluster identification process, in figure 12 the RDFs of the four particle categories P-I through P-IV are shown for the $St=0.18$ , $Sv=0$ case. The values of $D_{2}$ follow a generally decreasing trend from P-I ( $D_{2}=2.73$ ) to P-IV ( $D_{2}=2.38$ ), reflecting the higher degree of spatial organization and accumulation of the clustered particles. The trend is the same for all St, $Sv$ cases. Focusing on the coherent clusters, figure 13(a) shows RDFs of the P-IV particle distributions. As expected, the steepest slope is found for St near unity, corresponding to a minimum value of $D_{2}$ , both in the cases with and without gravity. The range over which the RDF remains greater than one indicates the length scale over which the clustering occurs. This increases monotonically with St, consistently with the above-mentioned tendency of more inertial particles to form larger clusters. To highlight the spatial organization of the clusters in the presence of gravity, figure 13(b) displays the ADF calculated for the $Sv>0$ cases for angular bins $|\unicode[STIX]{x1D703}|=[0^{\circ },10^{\circ }]$ (vertical direction) and $|\unicode[STIX]{x1D703}|=[80^{\circ },100^{\circ }]$ (horizontal direction). For the high-St, high- $Sv$ cases, the vertical ADF decays over much larger distances than the horizontal ADF, a signature of the vertical particle streaks, which can have lengths greater than $100\unicode[STIX]{x1D702}$ (figure 5).

Figure 12.

RDFs of the four particle categories P-I–P-IV for the $St=0.18$ , $Sv=0$ case.

Figure 13.

(a) RDFs of the P-IV particle distributions for all cases, with filled symbols for cases without gravity and empty symbols for cases with gravity. (b) ADFs calculated for the $Sv>0$ cases for angular bins $|\unicode[STIX]{x1D703}|=[0^{\circ },10^{\circ }]$ (vertical direction, empty symbols) and $|\unicode[STIX]{x1D703}|=[80^{\circ },100^{\circ }]$ (horizontal direction, filled symbols). In both plots, the cases are as follows: $St=0.18$ (black squares), $St=0.9$ (green triangles), $St=3.6$ (blue diamonds), and $St=11.2$ (red circles).

Taken together, these results demonstrate how, besides exhibiting a power-law decay of their size distribution (which indicates statistical self-similarity among the various objects), the individual coherent clusters also display geometrically self-similar structures. Their fractal nature is also reflected in the highly convoluted boundaries, which may have important implications for, for example, the mixing between laden and unladen fluid. The general shape of the clusters (tubular rather than planar) also appears to depend on the object size. Overall, the results suggest a link between the geometry of the coherent clusters and that of the underlying turbulent motions. To corroborate these observations, it will be of interest to use larger simulations and evaluate the effect of the turbulence Reynolds number (and therefore of the range of scales) on the quantities calculated here.

4.3 Fractal dimensions of individual clusters

Being self-similar objects, the coherent clusters display a complex topology which can be described in a compact manner by their fractal dimension. As discussed below, different aspects of their features (global shape, small-scale structure, enveloping surface) can be described via different fractal dimension definitions; here we use some of the more common, often employed to characterize physical and biological systems (Falconer Reference Falconer2004). This quantitative characterization may also help in the development of low-order models of the cluster features and behaviour. The fractal structure of inertial particles fields in turbulence has been the subject of previous investigations (e.g. Bec Reference Bec2003; Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Calzavarini et al. Reference Calzavarini, Kerscher, Lohse and Toschi2008). Here we instead focus on the fractal dimension of individual clusters. The relation between the fractal description of individual structures and that of the global field in the domain is not straightforward, see Moisy & Jiménez (Reference Moisy and Jiménez2004).

Several definitions of fractal dimension are commonly used, and are in general not equivalent except for mathematically constructed fractal sets. Moreover, inertial particles in turbulence have been shown to resemble multifractal sets (Bec Reference Bec2005), i.e. to possess different fractal properties at different scales. Multifractal sets are characterized by a range of fractal dimensions, among which are the well-known box-counting dimension $D_{0}$ and correlation dimension $D_{2}$ (Grassberger & Procaccia Reference Grassberger and Procaccia1983). Here we calculate both dimensions for the sets of particles belonging to individual clusters, as well as the fractal dimension of the cluster surface (defined in the sense specified in § 3). Only the number $N_{C}^{\prime }$ of clusters larger than $100\unicode[STIX]{x1D702}^{3}$ are analysed; smaller clusters did not have sufficient separation of scales to obtain good fractal dimension estimates. Depending on the case, $N_{C}^{\prime }$ ranges from several hundred to more than one thousand.

To determine the cluster box-counting dimension $D_{0C}$ , the portion of the domain containing each cluster is divided into non-overlapping cubes of side length r, and the number of boxes $N(r)$ containing at least one particle belonging to the cluster is counted. For large enough clusters, the number of boxes N needed to contain the full particle set follows a power law $N\propto r^{-D_{0}}$ over some range of r. As an example, plots of $N(r)$ for the four largest clusters of the $St=0.9$ , $Sv=0$ case are shown in figure 14(a). A fit through the power-law region in the range $1<r/\unicode[STIX]{x1D702}<8$ is used to calculate $D_{0C}$ for each cluster. To determine the correlation dimension of individual clusters $D_{2C}$ (not to be confused with $D_{2}$ , calculated over the entire particle set), the RDF of the particles in each cluster is considered and a best fit is applied to the power-law region in the range $0.5<r/\unicode[STIX]{x1D702}<1.5$ . Figure 14(b) shows RDFs for the same four clusters as in figure 14(a). Finally, the surface area dimension $D_{S}$ is found by fitting a power law to the scatter plots of cluster surface area against cluster volume (see figure 4).

Figure 14.

(a) Number of boxes versus box size normalized by the Kolmogorov length for the four largest clusters of the $St=0.9$ , $Sv=0$ case. The box-counting dimensions were found by fitting power laws in the range $1<r/\unicode[STIX]{x1D702}<8$ . (b) RDFs for the same four clusters. The correlation dimensions were found by fitting power laws in the range $0.5<r/\unicode[STIX]{x1D702}<1.5$ .

Figure 15.

$D_{0C}$ , $D_{2C}$ , and $D_{S}$ as a function of Stokes number, calculated as the average over the $N_{C}^{\prime }$ coherent clusters with volumes larger than $100\unicode[STIX]{x1D702}^{3}$ .

Figure 15 shows the values of $D_{0C}$ , $D_{2C}$ , and $D_{S}$ as a function of St, calculated as the average over the $N_{C}^{\prime }$ coherent clusters for each case. $D_{2C}$ has a similar trend to $D_{2}$ (see figure 11), but the values are lower due to the higher degree of spatial organization. Macroscopically, both the correlation dimension and the box-counting dimension are measures of the dimensionality of the particle distributions: values near 1 correspond to clusters resembling rod-like objects, values near 2 correspond to sheet-like objects, and values near 3 correspond to globular objects. However, both dimensions characterize different geometrical aspects of the small-scale distribution. In particular, $D_{2C}$ gives more weight to regions with high particle number densities (Gustavsson & Mehlig Reference Gustavsson and Mehlig2016). The difference between the two is thus not unexpected. Still, assuming the particles in each cluster follow the multifractal formalism, one would expect $D_{2C}\leqslant D_{0C}$ (Grassberger & Procaccia Reference Grassberger and Procaccia1983), in contrast to what is found here. However, one should consider that, as we are well outside of the limit $St\ll 1$ , the clusters may not exactly follow a multifractal behaviour. Moreover, we note that large uncertainties are likely associated to $D_{0C}$ , which is known to be very sensitive to discreteness effects and other sources of error that can significantly limit its robustness (Falconer Reference Falconer2004). Nevertheless, it is remarkable that the present values of $D_{0C}$ are close to the value of 1.7 obtained by Moisy & Jiménez (Reference Moisy and Jiménez2004) for individual high-shear structures in homogeneous isotropic turbulence. $D_{S}$ characterizes the surface of the clusters rather than their particle distribution, and is bounded by $2\leqslant D_{S}\leqslant 3$ . The limit $D_{S}=2$ is recovered for clusters enveloped by perfectly smooth surfaces, and the dimension increases as the cluster surface becomes more convoluted and therefore more space-filling. $D_{S}$ is found to be 2.85 or larger for all cases, confirming the visual impression that the borders of the clusters are extremely convoluted (see figure 6). As shown in figure 4, the values of $D_{S}$ are robust to possible subsampling bias.

Figure 16.

Correlation dimension of the coherent clusters conditionally averaged on their volume, for clusters with volumes larger than $100\unicode[STIX]{x1D702}^{3}$ . The $St=0.18$ cases are not shown because the number of large clusters for each data point is not sufficient for a statistically robust estimate.

Given that coherent clusters exist over a wide range of scales, it is of interest to investigate whether their fractal dimension changes with their size. Figure 16 shows $D_{2C}$ conditionally averaged on $V_{C}$ . The $St=0.18$ cases are not shown because the number of clusters for each data point is not sufficient for a statistically robust estimate. It is observed that the correlation dimension markedly increases with the cluster volume, and seems to plateau at values dependent on St and $Sv$ (although the limited range of scales cannot rule out further increase at larger volumes). For the cases with $St=O(1)$ , $D_{2C}$ for the large clusters tends to values close to the 2.6 value found by Calzavarini et al. (Reference Calzavarini, Kerscher, Lohse and Toschi2008) for the Kaplan–Yorke dimension in similar regimes. The dependence of $D_{2C}$ with the cluster size suggests that smaller clusters tend to be more tubular, and larger clusters tend to be more planar or globular. Under the assumption that clusters are a reflection of the eddy structures contributing to their formation, this trend appears consistent with the picture emerging from recent investigations of the geometry of eddy structures in homogeneous isotropic turbulence. In particular, Bermejo-Moreno, Pullin & Horiuti (Reference Bermejo-Moreno, Pullin and Horiuti2009) analysed isosurfaces of the curvelet transform of enstrophy and dissipation fields, and showed a transition from tubular to blob-like shapes with increasing scale. A similar conclusion was reached by Leung, Swaminathan & Davidson (Reference Leung, Swaminathan and Davidson2012) using a band-pass filtering approach at various scales. One should, however, be cautious when describing the topology, especially of large clusters, using $D_{2C}$ , as the latter is based on inter-particle correlation at small distances. Given the complex and multi-scale nature of the clusters, the use of multiple geometrical descriptors seems necessary to achieve a complete picture.

4.4 Preferential concentration

To examine the preferential concentration of clustered particles, the local concentration is conditionally averaged on the scalar vorticity magnitude $\unicode[STIX]{x1D6FA}$ and the strain rate $S$ for P-I and P-IV particles. Here $\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D6FA}_{ij}\unicode[STIX]{x1D6FA}_{ij}^{1/2}$ and $S=S_{ij}S_{ij}^{1/2}$ , where $\unicode[STIX]{x1D6FA}$ and $S$ are the rotation and strain-rate tensors, respectively. Following Wang & Maxey (Reference Wang and Maxey1993), we count the number of particles in each computational cell (centred at each grid point) in the domain, and bin those based on the local values of $\unicode[STIX]{x1D6FA}$ and $S$ . The concentration is then calculated as the number of particles counted for each bin, divided by the sum of the cell volumes associated to that bin. The conditionally averaged concentration is finally normalized by the global particle concentration of P-I and P-IV particles as appropriate, while $\unicode[STIX]{x1D6FA}$ and $S$ are normalized by their r.m.s. values. Results are shown in figure 17.

Figure 17.

Local concentration conditionally averaged on the scalar vorticity magnitude $\unicode[STIX]{x1D6FA}$ (empty symbols) and the strain rate $S$ (filled symbols) for P-I (black triangles) and P-IV (red circles) particles. Local concentration is normalized by the global concentration, and $\unicode[STIX]{x1D6FA}$ and $S$ are normalized by their r.m.s. values.

For the $St=0.18$ and $St=0.9$ cases without gravitational effects, both P-I and P-IV particles are concentrated in areas of high strain and low vorticity. This is also true for the $St=0.18$ , $Sv=1.9$ case. Such effects (particularly the preferential sampling of high-strain regions) are highly intensified for particles in coherent clusters as compared to generic inertial particles. This is consistent with the classic picture of clusters being formed by preferential concentration of particles in high-strain, low-vorticity regions of the flow (e.g. Eaton & Fessler Reference Eaton and Fessler1994). At higher St, in the absence of gravity there is still a tendency for the particles to accumulate in low-vorticity regions, but no preference for high strain. Finally, for high-St settling particles there is no relation between local concentration and the local vorticity or strain: their large fall speed makes them uncorrelated from small-scale flow structures.

Figure 18.

PDFs of $\langle Q\rangle _{C}$ (black dots) for the cases without gravity using unfiltered and filtered velocity fields. PDFs of $Q$ sampled at all P-I particle locations (red dashed line) and in the entire domain (green dotted line) are also shown for comparison. Values of $Q$ are normalized by their r.m.s. values.

The high-St particles are expected to be affected by large-scale flow structures (Yoshimoto & Goto Reference Yoshimoto and Goto2007), at least for zero or moderate settling. Preferential concentration is therefore examined for the particles in coherent clusters using a low-pass-filtered velocity field. The velocity field is filtered using a Gaussian kernel at two different length scales, $\unicode[STIX]{x1D6E5}=30\unicode[STIX]{x1D702}\approx L/2$ and $\unicode[STIX]{x1D6E5}=60\unicode[STIX]{x1D702}\approx L$ . Then the second invariant of the velocity gradient tensor, $Q=-1/2\;\text{tr}(\unicode[STIX]{x1D64E}^{2}-\unicode[STIX]{x1D734}^{2})$ , which is a measure of excess of vorticity over strain (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988), is computed at each grid point in the flow field. A cluster-averaged value $\langle Q\rangle _{C}$ is obtained for each coherent cluster by averaging $Q$ over the grid points contained within the volume of that cluster. The PDF of $\langle Q\rangle _{C}$ is shown for the cases without gravity in figure 18 for unfiltered and filtered velocity fields. The PDFs of $Q$ sampled at all P-I particle locations and in the entire domain are also shown for comparison. All PDFs are normalized by $Q_{RMS}$ obtained from the flow field filtered with the corresponding kernel. In all cases, the PDF of $Q$ sampled at all particle locations is nearly identical to the PDF of $Q$ at all fluid grid points. In contrast, the PDF of $\langle Q\rangle _{C}$ shows that particles in coherent clusters oversample low- $Q$ regions, where strain is high compared to rotation. For the $St=0.18$ and $St=0.9$ cases this effect can be most clearly seen in the PDFs for the unfiltered velocity fields. For the $St=3.6$ and $St=11.2$ cases, the effect is visible in the PDFs of the velocity fields filtered at increasing scales, confirming the argument that particles cluster over scale with which they ‘resonate’ based on their response time (see, e.g., Yoshimoto & Goto Reference Yoshimoto and Goto2007). Figure 19 shows similar plots for the $St=0.18$ , $St=0.9$ , and $St=3.6$ cases with gravity, using unfiltered velocity fields. These confirm that, when $Sv\gtrsim 10$ , the clusters are uncorrelated from low- $Q$ regions due to the high settling velocity. This is true even considering the low-pass-filtered velocity fields at $\unicode[STIX]{x1D6E5}=30\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D6E5}=60\unicode[STIX]{x1D702}$ (not shown).

Figure 19.

PDFs of $\langle Q\rangle _{C}$ (black dots) for the $St=0.18$ , $St=0.9$ , and $St=3.6$ cases with gravity using unfiltered velocity fields. PDFs of $Q$ sampled at all P-I particle locations (red dashed line) and in the entire domain (green dotted line) are also shown for comparison. Values of $Q$ are normalized by their r.m.s. values.

Figure 20.

Local concentration for P-I (black triangles) and P-IV particles (red circles), conditionally averaged on the vertical fluid velocity for the cases with gravity. The local concentration is normalized by the global concentration, and velocities are normalized by the theoretical terminal velocity.

4.5 Particle settling

To investigate the effect of turbulence on the settling of the particles in the coherent clusters, in figure 20 the local concentration is conditionally averaged on the vertical fluid velocity W (taken positive when directed upward), with the same method used for $\unicode[STIX]{x1D6FA}$ and $S$ . Velocities are normalized by the still-fluid terminal velocity $V_{t}$ , and both P-I and P-IV particles are considered for comparison. Preferential sweeping, i.e. the tendency of the particles to favour downward regions of the flow (Wang & Maxey Reference Wang and Maxey1993), is apparent for the $St=0.18$ , $Sv=1.9$ case and, to a lesser extent, the $St=0.9$ , $Sv=9.2$ case. In the former case, P-IV particles experience much stronger preferential sweeping than non-clustered particles, consistent with the stronger preferential concentration discussed in § 4.4. There is no evidence of preferential sweeping for the high-St cases; on the contrary, at $St=11.2$ , $Sv=28.7$ there is a clear tendency to preferentially sample regions of upward-moving fluid, especially for the P-IV particles. This behaviour, which is also visible (albeit weakly) for the clustered particles at $St=3.6$ , $Sv=36.7$ , is sometimes termed loitering (Nielsen Reference Nielsen1993; Good, Gerashchenko & Warhaft Reference Good, Gerashchenko and Warhaft2012) and is thought to be due to the fact that fast-falling particles spend more time to overcome upward-moving regions of a vortical flow field.

Figure 21.

Local concentration for the settling P-I (black triangles) and P-IV particles (red circles), conditionally averaged on the normalized vertical fluid velocity using the low-pass-filtered velocity field.

Since we have shown that the coherent clusters of highly inertial particles sample high-strain low-vorticity regions in the low-pass-filtered velocity field, it is interesting to explore how they sample downward or upward regions of the filtered field. This is illustrated in figure 21, which is obtained from the velocity filtered at $\unicode[STIX]{x1D6E5}=60\unicode[STIX]{x1D702}$ . The clustered particles at $St=0.18$ , $Sv=1.9$ show much weaker preferential sweeping by the coarse-grained fluid velocity. The opposite is true for $St=0.9$ , $Sv=9.2$ particles, which appear to ‘resonate’ better with the filtered velocity field and clearly favour downward flow regions. The clustered particles at $St=3.6$ , $Sv=36.7$ exhibit a remarkable non-monotonic trend: for moderate vertical velocities they preferentially sample downward regions of the filtered field, while for strong vertical velocity they are rather found in upward regions. Finally, $St=11.2$ , $Sv=28.7$ particles clearly favour upward flow regions. Filtering at $\unicode[STIX]{x1D6E5}=30\unicode[STIX]{x1D702}$ returns the same qualitative picture. Taken together, these trends suggest that the effect of turbulence on the settling speed results from the combination of multiple mechanisms. The particles are centrifuged out of vortices with which they ‘resonate’ based on their response time (and as allowed by the time they can spend in the vortex vicinity, which is limited by the settling). Their fall speed is in this case enhanced (preferential sweeping). However the heavier particles cannot respond to the more intense and fast events, and they instead fall through these quick eddy motions, spending more time in the upward regions (loitering). Which mechanism prevails depends thus not only on the particle inertia, but also on the specific range of scales within the spectrum of the turbulent flow events. This behaviour is evidenced by particles belonging to coherent clusters much more so than by generic inertial particles.

Figure 22 displays PDFs of normalized vertical velocity for P-I and P-IV particles, whose mean values are reported in table 5. Consistently with previous studies, a significant increase in settling velocity is found for St of order 1 and $Sv$ below 10 (Good et al. Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014; Rosa et al. Reference Rosa, Parishani, Ayala and Wang2016). For the $St=0.18$ , $Sv=1.9$ case the distribution is more positively skewed, and the settling velocity in the coherent clusters is more strongly enhanced compared to generic P-I particles. Indeed the mean settling velocity, which for the P-I particles is only 17 % higher than the theoretical still-fluid value $V_{t}$ , is almost twice as high as $V_{t}$ for the particles belonging to coherent clusters. Aliseda et al. (Reference Aliseda, Cartellier, Hainaux and Lasheras2002) also observed that inertial particles reached larger settling velocities when located in high-concentration regions, but given their volume fraction range they attributed this behaviour to the collective action of the clustered particles. In the present one-way coupled simulations no collective effect can ensue, therefore the large fall speed of the clusters must rather be a consequence of the more prominent preferential sweeping (figure 20).

Figure 22.

PDFs of normalized vertical particle velocity for P-I (black triangles) and P-IV particles (red circles). The velocity range in the different cases has been adjusted to better represent the shape of the distributions.

Table 5.

Mean normalized vertical particle velocity for P-I and P-IV particles.

The settling velocity distributions at $St=0.9$ , $Sv=9.2$ and $St=3.6$ , $Sv=36.7$ are fairly symmetric, and the mean values are negligibly different from $V_{t}$ for both P-I and P-IV particles, confirming that for these St and $Sv$ levels the particles become uncorrelated from the flow structures that cause preferential sweeping. At $St=11.2$ , $Sv=28.7$ the mean settling velocity for the P-IV particles is 2 % smaller than $V_{t}$ . More evidently, the PDFs show strong negative skewness. This indicates that the particles are more likely to sample upward fluid motions which reduce their fall speed, but occasionally encounter downward sweeps, leading to intermittent and very rapid settling events. We interpreted this as further evidence of loitering, which typically acts on particles falling too quickly through the flow to react to the eddies and spend more time on average in upward-moving regions of the flow (Nielsen Reference Nielsen1993; Good et al. Reference Good, Gerashchenko and Warhaft2012). In their DNS simulations, Good et al. (Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014) found that reduced settling only occurred when accounting for nonlinear drag. The present results suggest that this can in fact occur even when neglecting nonlinear drag (as done in the present simulations), but is only apparent when focusing on particles belonging to clusters.

Recent attempts have been made to further our understanding of particle–turbulence interaction under gravity effects, in particular by analysing the topology of the flow experienced by settling particles. Bec et al. (Reference Bec, Homann and Ray2014) provided an elegant analytical argument to show that, for $St\ll 1$ , the concentration of falling particles along horizontal planes is correlated with downward vertical fluid velocity, as already demonstrated numerically by Wang & Maxey (Reference Wang and Maxey1993) in homogeneous isotropic turbulence. They also showed that fast-falling particles with $St\geqslant 1$ experience more intense clustering as a consequence of an effectively two-dimensional dynamics along horizontal planes. Focusing on a similar range of parameters, Park & Lee (Reference Park and Lee2014) argued that fast-falling particles are more exposed to regions of converging fluid motion along horizontal planes, identified by negative in-plane divergence, i.e. $(\unicode[STIX]{x2202}U/\unicode[STIX]{x2202}x+\unicode[STIX]{x2202}U/\unicode[STIX]{x2202}y)=-\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}z<0$ . In figure 23(a) we explore the correlation between particle concentration and the vertical strain rate $\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}z$ (normalized by the mean strain rate $\langle S\rangle$ ). For the $St=0.18$ , $Sv=1.9$ case, P-I particles show weak correlation with the in-plane divergence, while P-IV particles are more abundant in regions where $\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}z<0$ , i.e. in regions of compressive vertical strain (i.e. positive horizontal-plane divergence). No significant correlation between concentration and $\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}z$ is found for particles at $St=0.9$ , $Sv=9.2$ or higher (not shown). The correlation between concentration and compressive vertical strain contrasts the proposition of Park & Lee (Reference Park and Lee2014), but is indeed in agreement with the classic view of preferential sweeping. Considering a simple cellular flow (similar to the one studied first by Maxey & Corrsin (Reference Maxey and Corrsin1986) and Maxey (Reference Maxey1987), and recently revisited by Bergougnoux et al. (Reference Bergougnoux, Bouchet, Lopez and Guazzelli2014)), inertial particles that tend to favour the downward sweeps also tend to cluster in regions of $\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}z<0$ , as illustrated in figure 23(b). While turbulent velocity fields are much more complex, the same qualitative picture for preferential sweeping was deduced by Wang & Maxey (Reference Wang and Maxey1993) by analysing DNS simulations. This behaviour is hardly visible when considering the totality of the inertial particles, but becomes apparent when considering the particles belonging to coherent clusters.

Figure 23.

Local concentration of P-I (black triangles) and P-IV particles (red circles) conditioned on $\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}z$ normalized by the global average of $S$ (a) for the $St=0.18$ , $Sv=1.9$ case with unfiltered velocity, the $St=0.9$ , $Sv=9.2$ case with unfiltered velocity, and the $St=0.9$ , $Sv=9.2$ case with low-pass-filtered velocity. The mechanism is illustrated in a sketch of a simple cellular flow (b), where inertial particles will tend to accumulate in regions of compressive vertical strain (marked with a star) rather than in regions of extensive vertical strain (marked with a diamond).

In figure 21 we have shown that clustered particles with $St=0.9$ , $Sv=9.2$ , which fall too fast to experience substantial preferential sweeping, still display a significant tendency to favour downward flow regions of a low-pass-filtered velocity field. In a similar manner, figure 23(a) shows that particles at $St=0.9$ , $Sv=9.2$ concentrate in regions where $\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}z<0$ when the velocity field is low-pass-filtered at $\unicode[STIX]{x1D6E5}=30\unicode[STIX]{x1D702}$ . Once again, such effect is only visible when considering P-IV particles. Even heavier particles are uncorrelated with $\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}z$ , independently of the filter size. Besides agreeing with the phenomenology of Maxey and co-workers, the correlation of the clustered particles with regions of compressive vertical strain may help explain the increase in cluster size with gravity for $St=0.18$ and $St=0.9$ (see table 4), as it could provide a mechanism for the particles to ‘pile up’ as they fall through these flow regions. It does not, however, explain the larger and vertically elongated clusters generated at $St=3.6$ and $St=11.2$ , since in these cases no significant correlation with $\unicode[STIX]{x2202}W/\unicode[STIX]{x2202}z$ is found. The cluster size and orientation in presence of gravity is further discussed in the next section.

4.6 Cluster orientation

The identification of coherent clusters allows an investigation of their individual spatial orientation. In particular, to determine the effect of gravity on their orientation, we calculate the angle $\unicode[STIX]{x1D703}_{g}$ between the vertical axis and the primary cluster axis determined by the SVD method. The PDF of the cosine of $\unicode[STIX]{x1D703}_{g}$ is shown in figure 24 for each case. Heavier particles have a stronger tendency to form vertically oriented clusters, in agreement with the visual evidence that fast-falling particles form long, vertical streaks (figure 5) and with the ADF analysis (figure 13). As expected, the degree of this effect appears to be determined by $Sv$ rather than by St.

Figure 24.

PDFs of the cosine of cluster orientation angle $\unicode[STIX]{x1D703}_{g}$ for coherent clusters.

The presence of large, vertically elongated clusters of heavy, fast-falling particles has been noted in several previous studies also focused on one-way coupled point-particle simulations (Woittiez et al. Reference Woittiez, Jonker and Portela2009; Dejoan & Monchaux Reference Dejoan and Monchaux2013; Bec et al. Reference Bec, Homann and Ray2014; Park & Lee Reference Park and Lee2014; Ireland et al. Reference Ireland, Bragg and Collins2016). Clearly, given the nature of the modelling approach, this behaviour cannot be explained by long-range multi-particle hydrodynamic interactions, as in columnar settling of suspensions of particles at relatively large Reynolds numbers (e.g. Uhlmann & Doychev Reference Uhlmann and Doychev2014). Previously the vertical alignment of the clusters has been related to a two-dimensionalization of the dynamics in the horizontal plane (Bec et al. Reference Bec, Homann and Ray2014; Park & Lee Reference Park and Lee2014), and to the influence of settling on the balance between inward drift and outward diffusion of particle pairs (Ireland et al. Reference Ireland, Bragg and Collins2016). However, a conclusive explanation has not yet emerged, especially to justify the large extension of these structures. The skewed distribution of vertical velocities for the $St=11.2$ , $Sv=28.7$ case in figure 22(d) suggests an alternative mechanism for the formation of the vertically elongated clusters: while high-St particles spend more time in upward-moving fluid regions (hence experiencing loitering), they sometimes encounter a strong downward gust of turbulence which, adding to their already high fall speed, projects them ballistically through the flow. These correspond to the extreme events of highly negative velocity that make up the left tail of the PDF in figure 22(d). For the duration of such gust, all particles that pass through the corresponding flow region are channelled vertically downward, creating columnar structures. This reasoning can also be used to explain the vertically elongated clusters at lower Stokes number, although it is more effective for more inertial particles. The mechanism would cease to play a role for particles so heavy as to be completely unaffected by the flow, which would remain randomly distributed. While this view (which shares conceptual similarities with the path-history effects discussed in Bragg & Collins (Reference Bragg and Collins2014) and Ireland et al. (Reference Ireland, Bragg and Collins2016)) is consistent with the observations, it awaits to be confirmed by dedicated investigations on the cluster formation and evolution.

Figure 25.

PDFs of the cosine of the angles $\unicode[STIX]{x1D703}_{i}$ between the cluster primary axis and the cluster-averaged vectors $\boldsymbol{e}_{1}$ (black squares), $\boldsymbol{e}_{2}$ (green triangles), $\boldsymbol{e}_{3}$ (blue diamonds), and $\unicode[STIX]{x1D6FA}$ (red circles) for coherent clusters.

In order to determine how cluster orientation is affected by the local velocity field, we calculate the angle $\unicode[STIX]{x1D703}_{i}$ between the primary axis of each coherent cluster and four vectors describing the local flow topology: the extensive, intermediate, and compressive eigenvectors of the strain-rate tensor ( $\boldsymbol{e}_{1}$ , $\boldsymbol{e}_{2}$ , and $\boldsymbol{e}_{3}$ , respectively), and the local vorticity vector, $\unicode[STIX]{x1D6FA}$ . For each cluster, these quantities are obtained from the average over the grid points contained within the volume of each cluster. The PDFs of the cosine of the four angles are shown in figure 25. For zero or moderate $Sv$ , the clusters tend to be strongly aligned with the local vorticity vector. It is well known that the vorticity tends to be indifferently aligned with $\boldsymbol{e}_{1}$ , strongly aligned with $\boldsymbol{e}_{2}$ , and strongly anti-aligned with $\boldsymbol{e}_{3}$ (Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987; Vincent & Meneguzzi Reference Vincent and Meneguzzi1991; Tsinober, Kit & Dracos Reference Tsinober, Kit and Dracos1992; Hamlington, Schumacher & Dahm Reference Hamlington, Schumacher and Dahm2008; Buxton & Ganapathisubramani Reference Buxton and Ganapathisubramani2010). The PDFs show that the clusters follow the same alignment with the strain-rate eigenvectors as the vorticity vector, confirming that the alignment with the latter is an inherent feature of coherent clusters. This result is consistent with (but not implied by) the classic picture of the particles being centrifuged outside of vortex cores and collecting around vortices (Eaton & Fessler Reference Eaton and Fessler1994; Calzavarini et al. Reference Calzavarini, Kerscher, Lohse and Toschi2008). Moreover, Ni, Ouellette & Voth (Reference Ni, Ouellette and Voth2014) recently showed that the vorticity vector tends to align with the direction of maximum Lagrangian stretching. It is therefore possible that, over times comparable with the clusters’ lifetime, groups of particles with St of order unity (and marginally or moderately affected by gravity) become elongated as they are carried by the fluid that undergoes vortex stretching, and therefore end up forming a cluster aligned with the vorticity vector. To confirm this proposition, a detailed investigation of the Lagrangian cluster properties is required, which is beyond the scope of the present study.

Another observation prompted by figure 25 is that the high- $Sv$ clusters, while showing no preferential alignment with the vorticity vector, have a noteworthy tendency to align with the compressive strain-rate eigenvector, $\boldsymbol{e}_{3}$ . The dynamics behind this result are not clear and deserve further investigation.