3.2.1. Voronoi tessellation for particle clustering

Voronoi tessellation (see, e.g. (Aurenhammer Reference Aurenhammer1991)) is a technique to construct a decomposition of the fluid domain, into a finite number of Voronoi cells. If there are a finite number of points (particles) dispersed in space, a Voronoi cell is defined as a region of all points that are closer to a particle than any other particles (see figure 5b). The volume of the Voronoi cell is referred to as the Voronoi volume, $V_p$. The magnitude of this volume can be used to quantify particle clustering and void regions in a 3-D space, a smaller volume indicating pronounced clustering. Voronoi tessellation techniques have been first applied to inertial particles in turbulence by Monchaux, Bourgoin & Cartellier (Reference Monchaux, Bourgoin and Cartellier2010). Here, 3-D Voronoi tessellation is applied to the particle data obtained from the present DNS data using the Quickhull algorithm provided by the Qhull library in python (Barber, Dobkin & Huhdanpaa Reference Barber, Dobkin and Huhdanpaa1996), which has a computational complexity of ${{O}}(N_p {\rm log}N_p)$. Results are compared with those of Oujia et al. (Reference Oujia, Matsuda and Schneider2020) obtained in the case of HIT.

Figure 5.

Particle distribution for $St = 1$ in a slice of thickness $1/100$ with spheres and mirror particles (a). A magnified view with Voronoi tessellation (b).

To construct the Voronoi tessellation of particles in the presence of the embedded bead boundaries, a special treatment is needed for particles near the bead walls. Figure 5 shows how the geometry of the spherical beads is taken into account. All particles located at a given small distance from each sphere boundary is assigned a mirror particle. Introducing these ghost particles then accounts for the bead boundary, and the Voronoi tessellation can then be constructed making use of the mirror particles. To verify that this approach works, a random distribution of particles in the fluid domain is first considered. For randomly distributed particles, the p.d.f. of the normalized Voronoi volume can be approximated with a gamma distribution (Ferenc & Néda Reference Ferenc and Néda2007). For the 3-D case, the p.d.f. of the Voronoi volume is approximated by $\varGamma (5,1/5)$, where $\varGamma (k,\theta )$ corresponds to the gamma distribution,

(3.3)\begin{equation} f_{V_p}(x) = \varGamma(k)^{{-}1}\theta^{{-}k}x^{k-1}\exp({-}x/\theta), \end{equation}

where $k$ and $\theta$ are shape parameters of the p.d.f., respectively and $\varGamma (k)$ denotes the gamma function. Figure 6(a) shows the p.d.f. of the Voronoi volume $V_p$, normalized by the mean volume $\overline {V_p}$, for different Stokes numbers as well as for randomly distributed particles. The p.d.f. for random distribution follows closely with the gamma distribution, as expected. It should be noted that in the present fictitious domain method, the bead boundaries are smeared over a grid control volume owing to the interpolation function between the bead material points and the computational grid typical of penalization based techniques. Accordingly, to perfectly match the random particle p.d.f. to the gamma distribution, the sphere radius had to be increased by approximately 1.1 % ($0.35749$ instead of $0.353553$), only for the Voronoi analysis, which is comparable to the grid resolution used. This small modification in sphere geometry is applied to all particle distributions obtained from different Stokes numbers.

Figure 6.

Probability distribution function of volumes of Voronoi cells (a) in log–log representation normalized by the mean volume for different Stokes numbers as well as for particles distributed randomly following a Poisson distribution. The dotted vertical line represents $\eta ^3/\overline {V_p}$. Mean Voronoi volume of particle clusters as a function of the distance from the bead surface (b) normalized by the bead radius.

Figure 6(a) shows that the p.d.f.s of Voronoi volumes for different Stokes numbers intersect with the gamma distribution, i.e. the one for the random particles. An equivalent figure can be found for HIT in Oujia et al. (Reference Oujia, Matsuda and Schneider2020), where it is observed that as the Stokes number increases, the number of the large Voronoi cells increases and then stabilizes. The number of small Voronoi cells increases as the Stokes number increases, and then decreases for $St> 1$. In the present case of a porous cell, the observed behaviour is similar to that found for HIT, but with some key differences. The number of large Voronoi cells increases with increasing Stokes numbers and then stabilizes similar to the HIT case. It is also seen that as the Stokes number increases and gets closer to $1$, the number of small normalized volumes ($10^{-2}-5 \times 10^{-1}$) also increases similar to the HIT case. However, a large number of very small volumes ($10^{-4}$–$10^{-2}$) are observed for nearly all Stokes numbers, even for very small Stokes numbers. Such a behaviour was not observed in HIT, wherein the number of very small Voronoi volumes would decrease monotonically. This is attributed to the interaction of particles with the bead surfaces in the present case. Figure 6(b) shows the mean Voronoi volume of particle clusters as a function of the distance from the bead surfaces normalized by the bead radius. It illustrates that on average the volumes close to the beads are smaller than the volumes further away. Particles with finite, non-zero Stokes numbers interact with the bead surface and undergo specular reflection. In addition, as the particles approach the bead surfaces, their velocities are slowed down significantly as the fluid velocity itself is smaller owing to the no-slip condition. Thus, existence of a large number of very small volumes is mainly attributed to the collision of the particles with the bead surfaces. The small-scale behaviour of the particles has its limitation implied by the used point particle and the collision model. Considering solid particles in air with a density ratio of 1000 gives particle size of the order of $0.04 \eta$ for a Stokes number of $0.1$. In this case, the behaviour of the small volumes would change depending on the surface condition at the beads. Note that, fluid particles (or tracers) would not collide with the bead surfaces, except at the stagnation points and at large times. Thus, the distribution of fluid tracer particles will follow the one of random particles, i.e. the gamma distribution.

3.2.2. Voronoi-based divergence of particle velocity

To understand the clustering dynamics of inertial particles, the particle number density $n$, as a continuous function, is commonly used (Oujia et al. Reference Oujia, Matsuda and Schneider2020). It satisfies the conservation equation

(3.4)\begin{equation} \frac{{\rm D}n}{{\rm D}t} = \frac{\partial n}{\partial t} + {\boldsymbol u}_p\boldsymbol{\cdot} \boldsymbol{\nabla} n ={-} n \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol u}_p . \end{equation}

The divergence of the particle velocity ${\boldsymbol u}_p$ appears as a source term in the number density equation. However, finding the divergence of the particle velocity is not straightforward as the discrete particle distribution is not continuous, and the particle velocity is only known at the discrete particle locations, ${\boldsymbol u}_{p,j} = {\boldsymbol u}({\boldsymbol x}_{p,j})$, but not everywhere else. Moreover, it can be multi-valued. Oujia et al. (Reference Oujia, Matsuda and Schneider2020) proposed a method to compute the divergence of the particle velocity ${\mathcal {D}} = \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol u}_p$ in a discrete manner using a Lagrangian approach. From the conservation equation for the particle number density (3.4), one obtains ${\mathcal {D}} = -({1}/{n})({{\rm D}n}/{{\rm D}t})$.

To calculate the Lagrangian derivative of $n$, Oujia et al. (Reference Oujia, Matsuda and Schneider2020) defined the local number density $n_p$ as the number density averaged over a Voronoi cell, which is given by the inverse of the Voronoi volume $V_p$; i.e. $n_p =1/V_p$. Then it was shown that the divergence is obtained to the first-order approximation as

(3.5)\begin{equation} {\mathcal{D}}_p = \frac{2}{\Delta t} \, \frac{V^{k+1}_p-V^{k}_p}{V^{k+1}_p+V^{k}_p} + O(\Delta t) , \end{equation}

where $V^{k}_p$ denotes the Voronoi volume at time instant $t^k$. This shows that the divergence of the particle velocity can be estimated from subsequent Voronoi volumes, provided the time step is sufficiently small and the number of particles is sufficiently large. To obtain the subsequent Voronoi volumes, the particle positions were linearly advanced by ${\boldsymbol u}_p$; i.e. ${\boldsymbol x}_p^{k+1}= {\boldsymbol x}_p^{k}+{\boldsymbol u}_p\Delta t$. The time step was set to be the same as the flow solver time step, which is sufficiently small for the present DNS study. The influence of the step size has been checked and found to give the same result as in Oujia et al. (Reference Oujia, Matsuda and Schneider2020). For the different time steps, the same probability distribution of the divergence was obtained, except that the extreme values of the divergence are changing with $\Delta t$.

Determining the divergence with the above Lagrangian approach requires two time instants of the Voronoi volumes. Above it was shown that the bead geometry, i.e. the presence of curved walls, necessitates some special treatment for computing the Voronoi tesselation. Mirror particles are introduced to account for the geometry. For computing the discrete divergence, the particles are then linearly advanced in time to obtain the volume at a subsequent time instant. This procedure implies that some particles could enter the beads and the flow geometry is not respected anymore. Consequently, the volume of adjacent cells are impacted and the computed divergence value is erroneous. To remove this artifact, particles having a distance less than $0.35782$ from each of the bead centres were not taken into account, for the Voronoi-based divergence analysis. This value has been determined by considering the statistics of the divergence as a function of the wall distance. While for the mean value, which is close to zero, no significant influence was found, for the higher order statistics (variance, skewness and flatness), a significant change with the distance, even by orders of magnitude, was observed. For wall distances larger than $0.35782$, the values were found to be stable and remained almost constant.

In Oujia et al. (Reference Oujia, Matsuda and Schneider2020) the p.d.f.s of the Voronoi-based divergence are almost symmetric, centred around zero with stretched exponential tails. The variance of the divergence increases as the Stokes number increases. For Stokes numbers less than or equal to 0.2, the divergence values are close to those of the obtained values for fluid particles. But the divergence of the fluid particles should be equal to zero because the fluid is incompressible. For this reason, Oujia et al. (Reference Oujia, Matsuda and Schneider2020) considered that for Stokes less than 0.2, the results are dominated by a geometrical error due to the Voronoi tessellation, so these values were not analysed in the rest of their paper. In the case of a porous cell, the p.d.f.s of the Voronoi-based divergence for different Stokes numbers are shown in figure 7. They are centred around 0 and confined between $-40$ and $40$. It can be observed that extrema correspond to ${\pm }2/\Delta t$ which is an upper/lower bound, even a rigorous bound when neglecting $O(\Delta t)$ terms, similar to what is found for HIT (Oujia et al. Reference Oujia, Matsuda and Schneider2020). The divergence should become closer to zero as the Stokes number decreases to zero because the fluid particles in an incompressible flow are divergence free. However, in figure 7 the divergence for $St=0.01$ is comparable to that of $St=0.1$. It can be deduced that the divergence for these Stokes numbers is mainly caused by a geometrical effect due to Voronoi tessellation, which is also discussed in Oujia et al. (Reference Oujia, Matsuda and Schneider2020). However, for larger Stokes numbers, physical effects do predominate. Table 1 assembles the variance and flatness of the divergence ${\mathcal {D}}_p$ as a function of the Stokes number. The flatness decreases as the Stokes number increases (with the exception of $St=0.1$), this implies that the tails of the p.d.f.s decay faster as the Stokes number increases. In contrast the variance has the opposite behaviour and, thus, the p.d.f.s are becoming wider and wider for increasing Stokes number, confirmed in figure 7. Thereafter the Voronoi analysis will be discussed only for Stokes numbers larger or equal to $0.5$, to avoid the influence of the geometrical effect.

Figure 7.

Probability distribution function of divergence of particle velocity for different Stokes numbers.

Table 1.

Variance and flatness of divergence ${\mathcal {D}}_p$ as a function of the Stokes number.

The mean of the divergence as a function of the volume is defined as

(3.6)\begin{equation} \langle {\mathcal{D}}_p \rangle_{Vp} = \frac{1}{P(V_p/\overline{V_p})}\int_{-\infty}^{+\infty} {\mathcal{D}}_p \, P({\mathcal{D}}_p,V_p/\overline{V_p}) \, {\rm d}{\mathcal{D}}_p. \end{equation}

The values are shown in figure 8 for three Stokes numbers, $St=0.5, 1$ and $2$. The insert is a zoom to focus on the negative values. Recall that negative values correspond to convergence of the particles and positive values to divergence. It can be seen that the mean of the divergence is positive for small volumes, negative for large volumes and, similarly to Oujia et al. (Reference Oujia, Matsuda and Schneider2020), the amplitude increases with the Stokes number. This implies the presence of cluster formation for large volumes, cluster destruction for small volumes and that these two behaviours are amplified when the Stokes number increases. The zero-crossing point is $V_p/\overline {V_p} \approx 0.06-0.2$, a value close to what is observed in Oujia et al. (Reference Oujia, Matsuda and Schneider2020). These observations are consistent with those found in Oujia et al. (Reference Oujia, Matsuda and Schneider2020) where equilibrium distribution is argued to be the main reason. For the statistically stationary distribution of particles, the value of the divergence must be positive for clusters of particles and negative in the case of void regions. Otherwise the set of particles would converge to a set of discrete points in space in a finite time.

Figure 8.

Mean divergence $\langle {\mathcal {D}}_p \rangle _{Vp}$ as a function of the Voronoi volume for different Stokes numbers.

In summary, the above findings for the Voronoi analysis confirm the results found in Oujia et al. (Reference Oujia, Matsuda and Schneider2020) for HIT, in particular the $St$ dependence of cluster formation and destruction. Differences are found for small Voronoi volumes below 0.01, they behave differently due to the influence of the bead geometry. The zero-crossing point of the mean value of the divergence is likewise shifted towards smaller values.

3.3.1. Scale-dependent wavelet analysis of number density

Inertial particle clustering and related multiscale statistics are quantified using the orthogonal wavelet decomposition (Daubechies Reference Daubechies1993; Mallat Reference Mallat2009) of the particle number density. Consider the particle number density, $n({\boldsymbol x},t)$, at a given instant $t$, within the computational domain of a triply periodic $L^3$ cubic box. This scalar field is decomposed into a 3-D orthogonal wavelet series to unfold into scale, positions and seven directions ($\mu = 1, 2,\ldots, 7$). The 3-D mother wavelet, $\psi _{\mu }({\boldsymbol x})$, is hereby based on a tensor product construction and a family of wavelets, $\psi _{\mu,{\boldsymbol \lambda }}({\boldsymbol x})$ can be generated by dilatation and translation. The multi-index ${\boldsymbol \lambda }= (j,i_1, i_2, i_3)$ denotes the scale $2^{-j}$ and position $L\times 2^{-j}{\boldsymbol i} = L\times 2^{-j}(i_1,i_2,i_3)$ of the wavelets for each direction, where $i_{\ell } = 0,\ldots,2^{j-1}\ (\ell =1,2,3)$. The wavelets are well localized in space around position $L\times 2^{-j}(i_1,i_2,i_3)$, and scale $2^{-j}$, oscillating, and smooth. A periodization technique (Mallat Reference Mallat2009) is applied to the wavelets. The spatial average of $\psi _{\mu,{\boldsymbol \lambda }}({\boldsymbol x})$, defined by $\langle \psi _{\mu,{\boldsymbol \lambda }} \rangle = L^{-3}\int _{{\mathbb {T}}^3} \psi _{\mu,{\boldsymbol \lambda }}({\boldsymbol x})\,{\rm d}{\boldsymbol x}$, vanishes for each index, which is a necessary condition for being a wavelet.

Similar to Matsuda et al. (Reference Matsuda, Schneider and Yoshimatsu2021), the number density field $n({\boldsymbol x},t)$ sampled on $N_g^3 = 2^{3J}$ equidistant grid points can be developed into an orthogonal wavelet series,

(3.10)\begin{equation} n({\boldsymbol x}) = {\bar{n}}_{000}({\boldsymbol x}) + \sum_{{ j}=0}^{J-1} n_{j} ({\boldsymbol x}), \end{equation}

where $n_{j}({\boldsymbol x})$ is the contribution of $n({\boldsymbol x})$ at scale $2^{-j}$ defined as

(3.11)\begin{equation} n_{j}({\boldsymbol x}) = \sum_{{ \mu}=1}^7 \sum_{{ i_1=0}}^{2^{ j}-1} \sum_{{i_2=0}}^{2^{j}-1} \sum_{{ i_3=0}}^{2^{ j}-1} {\tilde n}_{ { j,} {i_1, i_2, i_3} }^{ \mu} {\psi^{\mu}_{j, i_1, i_2, i_3}} ({\boldsymbol x}), \end{equation}

with

(3.12)\begin{equation} {\tilde n}_{ {j} {i_1 i_2 i_3} }^{ \mu} = \left\langle n({\boldsymbol x}), {\psi}_{ {j,} {i_1, i_2, i_3} }^{ \mu}({\boldsymbol x}) \right\rangle, \end{equation}

where $\langle \cdot\ , \cdot \rangle$ denotes an inner product. At scale $2^{-j}$, there are $7\times 2^{3j}$ wavelet coefficients for $n({\boldsymbol x})$. Thus, in total there are $N_g^3$ coefficients for each component of the vector field corresponding to the $N_g^3-1$ wavelet coefficients and the non-vanishing mean value. These coefficients are efficiently computed for $N_g^3-1$ grid points for $n({\boldsymbol x})$ using fast wavelet transform, which has linear computational complexity. The scale from the wavelet transform and wavenumber, $k_j$, from the Fourier transform are related as

(3.13)\begin{equation} k_j = \frac{2{\rm \pi}}{L} k_{\psi} 2^{{-}j}, \end{equation}

where $k_{\psi } =0.77$ is the centroid wavenumber of the chosen Coiflet 12 wavelet.

Figure 9(a,b) shows instantaneous, normalized mean number density contours in the $xy$-plane for two Stokes numbers, $St=1$ and $0.01$. A large number density, indicative of highly clustered regions, is clearly observed for $St=1$, but is absent for $St=0.01$. The instantaneous number density field is decomposed using wavelet transform and the scale-dependent fields ($n_j$) normalized by its variance ($\sigma [n_j]=\sqrt {M_2[n_j]}$) are shown for different scales in figure 9(c–f). It can be seen that clusters are prominent as scales become smaller (larger $j$), whereas large void regions (negative $n_j$) of size comparable to the pore size are seen at larger scales (smaller $j$). Prominent cluster regions are also seen near the bead surfaces, even for a low Stokes number ($St=0.01$). For intermediate scales, both clusters and voids are distributed more intermittently in space. The multiscale nature of clusters and voids is qualitatively clear from these figures. Scale-dependent statistics are computed to quantify these differences at various scales.

Figure 9.

Instantaneous distribution of normalized particle number density (a,b) and scale contributions ($n_j/\sigma [n_j]$) at $j=2$ (c,d), $j=4$ (e,f) and $j=6$ (g,h) in the $xy$-plane for $St=1$ (a,c,e,g) and $St=0.01$ (b,d,f,h). Results are shown for (a) $n^\prime$, $St=1$; (b) $n^\prime$, $St=0.01$; (c) $j=2$, $St=1$; (d) $j=2$, $St=0.01$; (e) $j=4$, $St=1$; (f) $j=4$, $St=0.01$; (g) $j=6$, $St=1$; (h) $j=6$, $St=0.01$.

3.3.2. Wavelet spectra, scale-dependent skewness and flatness

Wavelet-based statistics of the particle number density can be computed as described below. The $q$th moment of $n_j({\boldsymbol x})$ is defined as

(3.14)\begin{equation} M_q[n_j] = \left\langle (n_j)^q\right\rangle, \end{equation}

where the mean values $\left \langle n_j \right \rangle =0$, by giving a central moment and are related to the $q$th order structure functions (Schneider, Farge & Kevlahan Reference Schneider, Farge and Kevlahan2004).

The wavelet energy spectrum of $n_j({\boldsymbol x})$ can be defined using the second-order moment, $M_2[n_j]$, as

(3.15)\begin{equation} E[n_j] = \frac{1 }{\Delta k_j} M_2[n_j], \end{equation}

where $\Delta k_j = (k_{j+1}-k_j){\rm ln}2$ as given by Meneveau (Reference Meneveau1991). The energy spectrum obtained from the above equation has the particle number dependence due to the Poisson noise. Matsuda et al. (Reference Matsuda, Schneider and Yoshimatsu2021) analytically obtained the effect of the Poisson noise on $M_2[n_j]$ and succeeded in removing the Poisson noise from the wavelet energy spectrum. The second-order moment for randomly distributed particles in the cubic domain $\varOmega = \varOmega _f \cup \varOmega _s$ is $M_{2,{random},\varOmega }[n_j] = (7\times 2^{3j}/N_{p,\varOmega })\langle n \rangle _{\varOmega }^2$. For the present case, particles exist only in the fluid domain $\varOmega _f$, and $N_p = \phi N_{p,\varOmega }$ and $\langle n \rangle = \phi \langle n \rangle _{\varOmega }$. The energy of the Poisson noise is also reduced by a factor of $\phi$, i.e. $M_{2,{random}}[n_j] = \phi M_{2,{random},\varOmega }[n_j]$, yielding

(3.16)\begin{equation} M_{2,{{random}}} [ n_j ] = \frac{7\times 2^{3j}}{N_p}. \end{equation}

Hence, the following definition for the wavelet energy spectrum is used:

(3.17)\begin{equation} {E} [n_j] = \frac{1}{\Delta k_j} \left\{ M_2 [ n_j ] - \frac{7\times 2^{3j}}{N_p} \right\}. \end{equation}

Here the influence of the Poisson noise has been removed. Note that the analytical estimate of (3.16) does not contain a contribution of the beads geometry. It was found that the energy due to Poisson noise is orders of magnitude smaller for all scales and does not significantly affect the energy spectra of inertial particles. If $M_2[n_j]$ is computed from a realization of random particle distribution in the fluid domain $\varOmega _f$, the combined effect of the Poisson noise and geometrical confinement in $\varOmega _f$ can be observed.

The asymmetry of the p.d.f. of $n_j({\boldsymbol x})$ is quantified by the skewness defined as

(3.18)\begin{equation} S[n_j] = \frac{M_3[n_j] }{(M_2[n_j])^{3/2}}. \end{equation}

The scale-dependent flatness, which measures the intermittency at scale $2^{-j}$, is given as

(3.19)\begin{equation} F[n_j]= \frac{M_4[n_j] }{(M_2[n_j])^{2}}, \end{equation}

and is equal to three at all scales for a Gaussian distribution. It should be noted that the influence of the Poisson noise on $M_3[n_j]$ and $M_4[n_j]$ cannot be removed by subtracting the moments for randomly distributed particles. Matsuda et al. (Reference Matsuda, Schneider and Yoshimatsu2021) introduced the signal-to-noise ratio (SNR) defined as the ratio of the energy spectrum for inertial particles to that for randomly distributed particles, i.e. $\mathrm {SNR}=E[n_j]\Delta k_j/M_{2,{random}}[n_j]$. They confirmed that the effect of the Poisson noise on the statistics is negligibly small when the SNR is larger than 10.

Scale-dependent wavelet spectra, flatness and skewness of number density fluctuation are obtained by collecting data over a time span of $15 \, T_{11}^L \approx 30 \, \tau _{\eta }$, where $T_{11}^L$ is the Lagrangian time scale and $\tau _{\eta }$ is the Kolmogorov time scale. The scale-dependent statistics are averaged using approximately 15 equally spaced snapshots of number density over this time frame.Figure 10(a) presents the wavelet spectra of the number density fluctuations, $E[n_j]$, for different Stokes numbers as a function of the wavenumber, $k_j$, normalized by the Kolmogorov scale, $\eta =(\nu ^3/\langle \epsilon \rangle )^{1/4}$. Also shown is the spectrum for random particle positions with uniform probability, where the p.d.f. satisfies the Poisson distribution, resulting in $E[n_j]\propto k_j^2$. Moreover, random particles in the porous region of the FCC geometry are considered for comparison, and the spectrum shows a behaviour similar to the $St=0.01$ case. This is expected as random particles correspond to the case $St=0$, representative of fluid tracer particles. Comparing this with random particles in the cubic cell without the beads, i.e. the $k^2$ scaling, quantifies the influence of the geometrical confinement. At small scales increasing energy is observed and a similar magnitude as for $St \le 0.1$. For $St=$ 0.5, 1 and 2, it is seen that the spectrum first increases with a increasing $k_j \eta$ (large scales) and then gradually decreases for large $k_j\eta$ (small scales). The peak in the spectrum is found to gradually shift towards larger scales (smaller $k_j \eta$). This observation is consistent with the clustering of inertial particles in homogeneous, isotropic turbulence (HIT; Matsuda et al. Reference Matsuda, Schneider and Yoshimatsu2021). In addition, for higher Stokes numbers, the energy $E[n_j]$ is generally higher for all $k_j \eta$. For $St<0.5$, the peak observed in the spectrum at higher $St$ is not seen and $E[n_j]$ increases monotonically with $k_j \eta$. Subtracting the noise in the porous region, shown in figure 10(b), allows us to recover the peaks for small $St$. The peak values become higher for similar $k_j \eta$, which indicates that the void scale is almost constant. This is consistent with results for homogeneous, isotropic turbulence (Matsuda et al. Reference Matsuda, Onishi, Hirahara, Kurose, Takahashi and Komori2014). For isotropic turbulence, with low Stokes numbers, the spectra show a steeper slope close to $k^{-1}$ at small scales ($k_j \eta \gtrsim 0.4$). This difference in slopes at low $St$ in the presence of beads is attributed to the effect of inertial particle collisions with the bead surfaces as well as the geometric confinement effect of the flow.

Figure 10.

Wavelet energy spectra of particle number density fluctuation $E[n_j]$ for different Stokes numbers. Also shown are energy spectra for Poisson noise ($k^2$ line) in a cubic box (Random, no beads) and uniform random data in the porous region of the FCC geometry (Random, FCC). Results are shown for (a) $E[n_j]$; (b) $E[n_j]-(E[n_j])_{Random,FCC}$.

The peak of the spectrum of the particle density shown in figure 10(a) yields a wavenumber which corresponds to a clustering length scale given as

(3.20)\begin{equation} \frac{L_c}{\eta} = \frac{\sum_{j=0}^{J-1} (k_j \eta)^{{-}1} E[n_j] \Delta k_j }{\sum_{j=0}^{J-1} E[n_j] \Delta k_j}. \end{equation}

Variation of this length scale in comparison to the Kolmogorov scale as a function of the Stokes number is reported in table 2. It is observed that as the Stokes number increases, the wavenumber where the peak occurs becomes smaller indicating a larger length scale of the clustering. This integral length scale approaches the value for random distribution of the particles in the FCC geometry as the Stokes number decreases. Compared with the integral length scale of the flow, $L_{11} = 25.58\eta$, the clustering length scale is an order of magnitude smaller for all Stokes numbers, suggesting that the clustering is mainly impacted by smaller scales.

Table 2.

Clustering integral length scale ($L_c$) normalized by the Kolmogorov scale ($\eta$) for different Stokes numbers and randomly distributed particles with and without the FCC bead geometry.

Flatness and skewness statistics of number density fluctuations are computed for various Stokes numbers, as well as randomly distributed particles for the cubic box and the porous region, and are shown in figure 11(a,b). The statistics of randomly distributed particles correspond to those of fluid tracer particles ($St=0$). The number density of the fluid particles is uniform due to the volume preserving nature of the incompressible flow and particle clustering is absent resulting in zero skewness and flatness close to 3 as expected. For random data in the porous region of the FCC geometry, similar results were obtained with only minor differences. For all inertial particles with different Stokes numbers, the flatness values increase with decreasing scales, showing that intermittency of particle clustering is high at smaller scales. This is also seen qualitatively in the contour plots (figure 9e–h) of number density fluctuations at smaller scales. For $St>0.5$, the flatness values for each scale ($k_j \eta$) decreases with decreasing Stokes number. However, for small Stokes numbers, there is an increase in intermittency at intermediate and smaller scales ($k_j \eta > 0.2$) compared with larger Stokes numbers. Note that for lower Stokes numbers, particle clustering is typically not significant (figure 9b). Thus, this increase in flatness at smaller scales may be a result of the collision of particles with bead surfaces.

Figure 11.

Scale-dependent skewness $S[n_j]$ (a) and flatness $F[n_j]$ (b) for different Stokes numbers and random particles in a cubic box (Random, no beads) and uniform random data in the porous region of the FCC geometry (Random, FCC). Solid lines connect the symbols at scales for which $\mathrm {SNR} \ge 10$, and dotted lines are used otherwise. (a) Skewness. (b) Flatness.

For inertial particles at all Stokes numbers, skewness values also increase with decreasing scales. As shown by Matsuda et al. (Reference Matsuda, Schneider and Yoshimatsu2021), positive skewness of number density fluctuations indicate a high probability of large positive values of $n_j({\boldsymbol x})$, that is prominent clusters, whereas negative skewness corresponds to large excursions of negative values of $n_j({\boldsymbol x})$, that is void-pronounced structures. Increase in positive skewness at smaller scales for all Stokes numbers implies more prominence of clusters. At $St=0.01$ and $0.1$, the skewness remains close to zero at intermediate and large scales, similar to the random particle statistics. However, at smaller scales there is an increase in positive skewness for these Stokes numbers. This again is conjectured to be clustering of these particles owing to collisions with the bead surfaces. Finally, small negative values of skewness observed for $St>0.5$ at large scales are indicative of large regions of voids. This is also qualitatively confirmed in the visualizations shown in figure 9.