1. Introduction
Several applications including catalysis in the chemical synthesis and process industries (Aris Reference Aris1999; Dixon & Nijemeisland Reference Dixon and Nijemeisland2001), high temperature nuclear reactor cooling (Shams et al. Reference Shams, Roelofs, Komen and Baglietto2013), hyporheic exchange of pollutants at the sediment-water interface (Hester et al. Reference Hester, Cardenas, Haggerty and Apte2017), sand-migration in oil/gas wells, involve unsteady transitional and turbulent flows through confined spaces and porous structures. In these examples, the contribution of the inertial terms in the fluid flow equations can dramatically change the topology of the flow field leading to the formation of jets, vortices, dead zones, etc. within the pores. Such flow features can substantially alter the dispersion characteristics of pollutants and play critical roles in the transport of reactants and products to and from active reaction sites.
Dynamics of small inertial particles, such as sand particles, through porous media is of importance in oil and gas production (Carlson et al. Reference Carlson, Gurley, King, Price-Smith and Waters1992). Sand production is considered one of the most important issues facing hydrocarbon wells as it can erode hardware, plug tubulars, cease flow, create downhole cavities and needs separation and disposal at the surface. Variations in the reservoir pressure and completion permeability leads to high velocity in the perforation tunnels, significant inertial and turbulence effects and kinetic energy losses in the completion region (Cook et al. Reference Cook, Lee, DiGiovanni, Bronowski, Perkins and Williams2004), and onset of sanding triggered by rock failure (Yi, Valkó & Russell Reference Yi, Valkó and Russell2005). Minimizing pressure drop means that the gravel bed porosity should be as large as possible. However, to act as an effective filter for sand grains, the gravel also has to be small enough to restrain formation sand (Saucier Reference Saucier1974; Mahmud, Van Hong & Lestariono Reference Mahmud, Van Hong and Lestariono2019). Particle retention, bridge formation, jamming and deposition has been studied experimentally in simplified configurations such as microchannels (Valdes & Santamarina Reference Valdes and Santamarina2006; Dai & Grace Reference Dai and Grace2010; Agbangla, Climent & Bacchin Reference Agbangla, Climent and Bacchin2012) and packed beds (Pandya, Bhuniya & Khilar Reference Pandya, Bhuniya and Khilar1998; Ramachandran & Fogler Reference Ramachandran and Fogler1999) to show that rate of plugging by bridging has a nonlinear dependence on particle concentration. This bridging effect depends on flow velocity, particle diameter-to-pore size ratio and flow geometry.
The goal of the present work is to investigate and understand the clustering dynamics of inertial particles in a turbulent flow through confined geometries representative of a porous medium. Specifically, how does the interaction of particles with the solid walls affect the clustering and deposition and how this changes with particle Stokes number is of critical importance. In addition, with high-Reynolds-number flows through porous media, the pore-scale flow structure can change significantly owing to the inertial and unsteady flow features. Previous work (He et al. Reference He, Apte, Finn and Wood2019) explored, using direct numerical simulation (DNS), how this change of flow structure impacts the turbulence (i.e. the turbulent kinetic energy (TKE) distribution and transport mechanisms) in a face-centred cubic (FCC) lattice with very low porosity. The flow geometry gives rise to rapid acceleration and deceleration of the mean flow in different regions with the presence of three-dimensional (3-D) helical motions, weak wake-like structures behind the bead spheres, stagnation and jet-impingement-like flows together with merging and spreading jets in the main pore space. The jet-impingement and weak wake-like flow structures give rise to regions with negative production of total TKE, a feature observed in jet-impingement-like configurations. Unlike flows in complex shaped ducts, the turbulence intensity levels in the cross-stream directions were found to be larger than those in the streamwise direction. Furthermore, due to the compact nature and confined geometry of the FCC packing, the turbulent integral length scales were estimated to be less than 10 % of the bead diameter even for the lowest, transitional Reynolds number, indicating the absence of macroscale turbulence structures for this configuration. This finding suggests that even for the highly anisotropic flow within the pore, the upscaled flow statistics are captured well by the representative volumes defined by the unit cell. In the present work the previous analysis of turbulence in porous geometry is extended to investigate dynamics and transport of inertial particles and quantify the effect of geometric confinement on particle clustering.
Clustering of inertial particles in turbulent flows has been well studied in homogeneous isotropic turbulence (HIT) and wall-bounded channel flows (Maxey Reference Maxey1987; Eaton & Fessler Reference Eaton and Fessler1994; Toschi & Bodenschatz Reference Toschi and Bodenschatz2009; Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2012). The divergence of particle velocity, which differs from the divergence-free fluid velocity in an incompressible flow, plays a crucial role in clustering of inertial particles. The divergence of particle velocity has been shown to be proportional to the second invariant of flow velocity gradient tensor for sufficiently small Stokes numbers, which is defined as the ratio of the particle relaxation time, 
$\tau _p$, to the Kolmogorov time, 
$\tau _{\eta }$. Particles tend to concentrate in low vorticity and high strain regions in turbulence resulting in the preferential concentration. Divergence of particle velocity has been used to quantify clustering mechanisms (Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Esmaily-Moghadam & Mani Reference Esmaily-Moghadam and Mani2016). An inherent difficulty for determining the divergence of the particle velocity is its discrete nature, i.e. it is only defined at particle positions. Recently, Oujia, Matsuda & Schneider (Reference Oujia, Matsuda and Schneider2020) computed the particle velocity divergence from the position and velocity of a large number of particles, using a Voronoi tessellation technique. They proposed a model for quantifying the divergence using tessellation of the particle positions. The corresponding time change of the volume is shown to yield a measure of the particle velocity divergence.
Pair correlation function (PCF) has also been widely used to analyse clustering as it is directly related to the particle collision rates (Sundaram & Collins Reference Sundaram and Collins1997; Wang, Wexler & Zhou Reference Wang, Wexler and Zhou2000). The PCF typically shows a power-law behaviour at sub-Kolmogorov scales and the slope is dependent on the Stokes number. Scale similarity of particle distribution has been explained by the sweep-stick mechanism proposed by Goto & Vassilicos (Reference Goto and Vassilicos2006) in which particles are swept by large-scale flow motion while sticking to stagnation points of Lagrangian fluid acceleration (Coleman & Vassilicos Reference Coleman and Vassilicos2009). The probability distribution function (PDF) of particle mass density and coarse graining techniques has been used to investigate the scale dependence of particle distribution (Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007). Bassenne et al. (Reference Bassenne, Urzay, Schneider and Moin2017) proposed a wavelet-based method to extract coherent clusters of inertial particles in fully developed turbulence. Wavelets decompose turbulent flow, scalar and vector valued fields, in scale, position and possibly direction, complementary to Fourier techniques which yield insight into wavenumber contributions of turbulent flow fields. The wavelet representation can be used to analyse spatial intermittency and quantify spatial fluctuations at different scales. This is not easily possible with Fourier transform, owing to the global character of the basis functions. For reviews on wavelets and turbulence, we refer the reader to Farge (Reference Farge1992), Schneider & Vasilyev (Reference Schneider and Vasilyev2010) and more particularly on wavelet-based statistics to Farge & Schneider (Reference Farge and Schneider2015). Recently, Matsuda, Schneider & Yoshimatsu (Reference Matsuda, Schneider and Yoshimatsu2021) obtained scale-dependent statistics of the particle distribution and insights into the multiscale structure of clusters and voids in particle-laden, homogeneous isotropic turbulent flow using orthogonal wavelet decomposition of the Eulerian particle density field at high Reynolds numbers.
Although particle clustering has been observed and studied in several particle-laden turbulent flows, dynamics and clustering of small, inertial particles in complex and confined configurations of densely packed porous beads has not been explored. Whether the commonly observed heavy particle clustering mechanisms at unit particle Stokes numbers also appear in a confined geometry, and how collisional interactions of these particles with the bead surfaces affect such clustering has not been studied. To address these issues, DNS is used to investigate the effect of turbulent flow in the confined geometry of a FCC porous unit cell on the transport of fine, inertial particles at different Stokes numbers (
$St = 0.01, 0.1, 0.5, 1, 2$) and at a pore Reynolds number of 500. Particles are advanced using one-way coupling and the collision of particles with pore walls is modelled as perfectly elastic specular reflection. Detailed analysis of clustering and void formation statistics is then conducted based on (i) Voronoi tessellation, and (ii) wavelet analysis of the particle number density field in the porous geometry. The clustering statistics are compared and contrasted against those in a homogeneous, isotropic turbulence flow (Oujia et al. Reference Oujia, Matsuda and Schneider2020; Matsuda et al. Reference Matsuda, Schneider and Yoshimatsu2021). Specifically, the impact of geometric confinement on particle clustering and void formation at different Stokes numbers is evaluated.
The rest of the paper is arranged in the following way. Section 2 provides the details of the porous geometry, simulation parameters, the numerical approach, as well as Lagrangian tracking and motion of inertial particles. Mean velocity field, TKE distributions, integral scales are described in § 3. Analysis of particle clusters and voids is then conducted using Voronoi tessellation and wavelet-based multiscale, scale-dependent statistics of particle number density. Finally, summary and conclusions are given in § 4.
2. Simulation set-up
Several different definitions of Reynolds numbers have been used in porous media literature (He et al. Reference He, Apte, Finn and Wood2019; Wood, He & Apte Reference Wood, He and Apte2020) based on different length scales such as the particle diameter (
$D_B$), the hydraulic diameter (
$D_H$) or the permeability. In this work the Reynolds number (
$Re_H$) based on 
$D_H$ is used. The hydraulic diameter is related to the particle diameter as
\begin{equation} D_H = {\frac{2}{3}} {\frac{\phi}{1-\phi}} D_B. \end{equation}
Then, 
$Re_H$ is defined as (dropping the factor 2/3)
\begin{equation} Re_H = \frac{\overline{\langle u_x \rangle ^f} D_B}{\nu_f}, \end{equation}
where 
$\nu _f$ is the kinematic viscosity of fluid, 
$\phi$ is the porosity of the medium defined as the ratio of the void volume (which corresponds to the fluid volume) to the total volume, 
$\overline {\langle u_x \rangle ^f}$ is the time-averaged interstitial (intrinsic average) velocity of flow in porous media and 
$u_x$ is the instantaneous velocity component in the 
$x$-(streamwise) direction. The spatial averaging operation is denoted by 
$\langle \cdot \rangle$, the superscript or subscript 
$f$ indicates the fluid phase and 
$\bar {\cdot }$ is the time-averaging operator. For clarity, note that the intrinsic average velocity is defined by
\begin{equation} \left.\langle u_x \rangle^f\right|_t = \frac{1}{V_f}\int_{{\boldsymbol x}\in \varOmega_f} u_x({\boldsymbol x},t)\, {\rm d}V({\boldsymbol x}), \end{equation}
where 
$V_f$ is the volume of the fluid phase within the unit cell and 
$\varOmega _f$ denotes the fluid domain. The overbar represents traditional time averaging of a quantity 
$q$ defined as
\begin{equation} \bar{q}(t)=\frac{1}{T} \int_{t^*=t-T/2}^{t^*=t+T/2} q(t^*)\,{\rm d}t^*. \end{equation}2.1. Porous geometry
A porous FCC unit cell (figure 1) is used based on our prior work on turbulence in porous media (He et al. Reference He, Apte, Schneider and Kadoch2018, Reference He, Apte, Finn and Wood2019). It has a half-sphere entering at each face of the cube, and a half-quarter sphere at each corner. The FCC arrangement creates the lowest possible porosity (
$\phi$) to be 0.26 for the structured packings. The bead surfaces touch each other forming a point contact, and the exact radii of spheres are used in the fluid flow simulation. Due to this extreme compactness, the flow through the unit cell experiences rapid expansion and contraction. A pressure gradient is imposed to drive the flow through the bed and a triply periodic boundary condition is applied for the unit cell. The majority of the flow enters the cubic cell through the upstream open corners, converges into the centre pore resulting in strong accelerations and decelerations, and then leaving the unit cell along downstream corners. In this work the flow at one pore Reynolds numbers, 
$Re_H = 500$ is simulated using DNS, resulting in a turbulent flow within the pore. Emphasis is placed on dynamics of inertial particles at different Stokes numbers (0.01, 0.1, 0.5, 1, 2), the definition is given below.
Figure 1. (a) Schematic of a face-centred porous unit cubic cell showing the bead arrangement and surface geometry, (b) the isosurface of swirling strength at 
$\lambda = 0.25 \lambda _{max}$ for 
$Re_H = 500$.
2.2. Numerical approach and grid convergence
The numerical approach is based on a fictitious domain method to handle arbitrarily shaped immersed objects without requiring the need for body-fitted grids (Apte, Martin & Patankar Reference Apte, Martin and Patankar2009). Uniform Cartesian grids are used in the entire simulation domain, including both fluid and solid phases. An additional body force is imposed on the solid part to enforce the rigidity constraint and satisfy the no-slip boundary condition. The absence of highly skewed unstructured mesh at the bead surface has been shown to accelerate the convergence and lower the uncertainty (Finn & Apte Reference Finn and Apte2013). The following governing equations are solved over the entire domain, including the region within the solid bed, and a rigidity constraint force, 
$\boldsymbol {f}$, is applied that is non-zero only in the solid region.
The governing equations read as
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol u} = 0, \end{gather}
\begin{gather}\rho_f \left[\frac{\partial {\boldsymbol u}}{\partial t} + \left({\boldsymbol u}\boldsymbol{\cdot} \boldsymbol{\nabla}\right) {\boldsymbol u} \right] ={-}\boldsymbol{\nabla} p + \mu_f \nabla^2{\boldsymbol u} + {\boldsymbol f} , \end{gather}
where 
$\boldsymbol u$ is the velocity vector (with components given by 
${\boldsymbol {u}}=(u_x,u_y,u_z)$, 
$\rho _f$ the fluid density, 
$\mu _f$ the fluid dynamic viscosity and 
$p$ the pressure. A fully parallel, structured, collocated grid solver has been developed and thoroughly verified and validated for a range of test cases including flow over a cylinder and sphere for different Reynolds numbers, flow over touching spheres at different orientations, flow developed by an oscillating cylinder, among others. The details of the algorithm as well as very detailed verification and validation studies have been published elsewhere (Apte et al. Reference Apte, Martin and Patankar2009). In addition, the solver was also used to perform direct one-to-one comparison with a body-fitted solver with known second-order accuracy for steady inertial, unsteady inertial and turbulent flow through porous media (Finn & Apte Reference Finn and Apte2013; He et al. Reference He, Apte, Schneider and Kadoch2018, Reference He, Apte, Finn and Wood2019).
For the present studies, the flow is driven by a pressure drop as a body force in a triply periodic domain. According to Hill & Koch (Reference Hill and Koch2002), a constant pressure gradient 
$\boldsymbol {\nabla } P$ in the main direction of the flow (
$x$ direction), proportional to the body force 
$F$, is used to drive the flow,
\begin{equation} \partial_x P = \frac{18 \mu_f c \overline{\langle u_x \rangle ^f}}{ D_B^2} F , \end{equation}
where 
$\mu _f$ is the dynamic viscosity of the fluid and 
$c$ the solid volume fraction, defined as 
$\frac {2}{3} {\rm \pi}(D_B/L)^3$ (
$L$ is the length of the unit cube). The body force 
$F$ changes with the pore Reynolds number according to the linear fit obtained by Hill & Koch (Reference Hill and Koch2002) and is given as
\begin{equation} F = 462+9.85 \left(\frac{1-\phi}{2}Re_H\right)\ (Re_H \geq 216) . \end{equation} A posteriori calculation of body force needed to balance the shear stress on the sphere surfaces for different Reynolds numbers exhibits a good agreement with the above correlation and has been published elsewhere (He et al. Reference He, Apte, Finn and Wood2019). A uniform, cubic grid is used with resolution chosen such that the first grid point is at 
$y^{+}<1$ (i.e. in the viscous sublayer) to accurately capture the wall layers, where 
$y^+=yu_{\tau }/\nu$ indicates the normalized distance from the sphere surface, 
$u_\tau =\sqrt {\tau _\omega /\rho }\approx 0.5\|\overline {\langle u_x \rangle ^f}\|$ is the friction velocity and 
$\tau _\omega$ is the estimated wall shear stress. To obtain a more direct estimate on grid resolution requirements in the present DNS simulations, 3-D DNS studies were performed at 
$Re_H=500$ in a unit cell of FCC spheres with a systematic grid refinement study using 
$48$, 
$64$, 
$96$, 
$112$, 
$128$ and 
$144$ grid points per bead diameter 
$D_B$. A grid converged solution was obtained for first-order (mean flow) as well as second-order statistics (TKE) for 
$Re_H=500$ at a resolution of 
$D_B/96$. For 
$Re_H=500$, the mean flow converged at 
$D_B/96$, whereas, TKE showed small changes compared with coarser mesh indicating that a grid converged solution can be expected with a resolution of around 
$D_B/100$. However, in order to obtain a high resolution DNS study and provide sufficient resolution in the bead contact region, a refined grid based on 
$D_B/\delta = 250$ (
$\delta$ is the grid resolution in one direction) was used to resolve the pore-scale flow structures (He et al. Reference He, Apte, Finn and Wood2019).
2.3. Particle tracking algorithm
A point-particle approach is used to model the motion of small, heavy inertial particles in a Lagrangian frame. In this approach the particle size is assumed to be much smaller than the Kolmogorov scale, and the forces acting on them are modelled using simple closure models (Maxey Reference Maxey1987). The particle motion is assumed to be governed by a simple, linear drag force. One-way coupling, wherein the inertial particles do not affect the fluid flow, is used by assuming that the particle size and concentration is small. In addition, inter-particle interactions are also neglected. It should be noted that for sediment-laden beds consisting of low particle-to-fluid density ratios, the above assumptions may not hold, and the particle sizes may become of the order of the Kolmogorov scales with high volume loadings. Studies accounting for two-way coupling and fluid displacement effects (Apte, Mahesh & Lundgren Reference Apte, Mahesh and Lundgren2008; Finn, Li & Apte Reference Finn, Li and Apte2016) can impact the particle dynamics and clustering and will be considered in future studies.
For different Stokes numbers (
$St$), the particle motion is given by
\begin{equation} \frac{{\rm d} {\boldsymbol x}_p }{{\rm d}t} = {\boldsymbol u}_p;\quad \frac{{\rm d} {\boldsymbol u}_p }{{\rm d}t} = \frac{{\boldsymbol u}_{f,p} - {\boldsymbol u}_p }{St\cdot\tau_{\eta}}, \end{equation}
where 
${\boldsymbol x}_p$ and 
${\boldsymbol u}_p$ are particle position and velocity, 
$\tau _{\eta }$ is Kolmogorov time scale, 
${\boldsymbol u}_{f,p}$ is the instantaneous fluid velocity interpolated to the particle location, i.e. 
${\boldsymbol u}_{f,p}={\boldsymbol u}({\boldsymbol x}_p)$, and 
$\tau _p = St \tau _{\eta }$ is the particle relaxation time. The particles are advanced using the instantaneous fluid velocity interpolated to the particle location. The effect of inhomogeneity in the fluid velocity and the confinement of the bead walls on inertial particle motion can thus be quantified by qualitatively comparing results to inertial particles in HIT (Oujia et al. Reference Oujia, Matsuda and Schneider2020; Matsuda et al. Reference Matsuda, Schneider and Yoshimatsu2021).
Particles with five different Stokes numbers, 
$St=$ 0.01, 0.1, 0.5, 1 and 2 are simulated. After a stationary turbulent flow is achieved within the porous cell, inertial particles are injected into the fluid domain. One particle is added at the centre of each control volume in the fluid region, giving about 
$N_p \approx 10.4\times 10^6$ particles for each Stokes number. The fluid velocity interpolated to the particle location is obtained from tri-linear interpolation, and a fourth-order Runge–Kutta (RK4) scheme is implemented to advance the particle locations in time. Interactions of particles with the bead walls in the unit cell are modelled using Snell's law of specular reflection assuming perfectly elastic collisions. The direction of particle reflection, 
${\hat {\boldsymbol s}}_r$, is determined from the incident direction, 
${\hat {\boldsymbol s}}_i$, and the inward face normal at the spherical bead surface, 
${\hat {\boldsymbol n}}$, as
\begin{equation} {\hat{\boldsymbol s}}_r = {\hat{\boldsymbol s}}_i + 2 |{\hat{\boldsymbol s}}_i \boldsymbol{\cdot} {\hat{\boldsymbol n}}| {\hat{\boldsymbol n}} . \end{equation}The particle velocity after reflection is also modified according to the above equation. The inward normal to the sphere surface can be easily obtained by knowing the location on the spherical bead surface that the particle crossed and the centre of the bead. In order to ensure that the particles do not cross the bead surface boundary and enter the solid region, a small gap proportional to the grid size near the bead surface is used for implementing the specular reflection.
3. Results and discussion
In this section the mean and turbulent flow structure is described briefly using the mean and root mean square (r.m.s.) velocity fields inside the pore, and the integral length and time scales are estimated based on the Eulerian and Lagrangian two-point autocorrelations, respectively. Next, the PDFs of Voronoi tesselation volumes and the divergence of particle velocity field are evaluated and discussed. In addition, the characteristics of inertial particle clustering in the porous media at various Stokes numbers are discussed using the multiscale, wavelet analyses of the particle number density, wavelet spectra, higher-order statistics of flatness and skewness. The Voronoi tesselation and wavelet analysis results for inertial particle statistics in the FCC unit cell are contrasted with those from isotropic turbulence to identify effects of geometric confinement.
3.1. Turbulent flow statistics
As described earlier in the simulation set-up, the flow through the triply periodic domain is driven by a body force in the flow direction based on the correlations given by (2.8). After an initial transient, a stationary state is reached and the computation is continued for several flow through times, 
$T_f = L/\overline {\langle u_x \rangle ^f}$, where 
$L$ is the length of the unit cube. For each case, the flow was first computed for several flow through times to ensure that a stationary state has been reached. This is confirmed by monitoring the total kinetic energy in the domain, which starts out to be a large value and then decreases and remains more or less constant after an initial transient period. This initial transient period was approximately 
$100 T_f$ for 
$Re_H=500$. After a stationary state has been established, the computations were performed for an additional 
$80 T_f$ to collect flow statistics which was found to be large enough to obtain converged statistics.
In order to get a good understanding of the simulated flow topology, distributions of the mean velocity magnitude 
$\overline U_m$ and the TKE are presented first (see figure 2). The mean velocity 
$\overline U_m$ is calculated as 
$\sqrt {\overline {u_x}^2 + \overline {u_y}^2 +\overline {u_z}^2}$ and normalized by the mean interstitial velocity 
$\overline {\langle u_x \rangle ^{f}}$. Two centre slices are chosen as representative sections for visualizations. One is the centre 
$xy$-plane, where the mean flow is going from left to right; the other the centre 
$yz$-plane, where the mean flow in going into the slice.
Figure 2. Visualization of the mean velocity magnitude (a,b) and TKE (c,d) at 
$Re_H=500$: (a–c) 
$xy$-plane, and (b–d) 
$yz$-plane.
The mean flow enters the centre pore from the left-hand side corners (regions labelled by (i
$_1$) and (i
$_2$) in figure 2a) and accelerates to region (ii) due to the geometric constrictions; then the mean flow at the centre starts to decelerate as it encounters the particle at region (iii), similar to an impinging jet. The flow then accelerates in the spanwise directions and leaves the pore through the right-hand side corners (regions (iv
$_1$) and (iv
$_2$)). The pattern of mean velocity distribution indicates how the mean flow in the pore is affected by the geometry. High-Reynolds-number flow features such as wake and large-scale vortex shedding in an external flow behind a spherical particle, are not observed in this low porosity bed. A weak wake-like structure (small negative mean velocity just behind the left-hand side sphere near region (ii) in figure 2a) is present, but is confined to a small region. The closely packed solid beads in the low porosity FCC configuration tend to break down large-scale flow structures and prevent the generation of a significant wake region. Figure 2(b) shows the distributions of the mean velocity magnitude on the centre 
$yz$-plane. On this slice, the mean flow is moving perpendicularly into the page and there is a distinguished region with high velocity magnitude near each corner of the centre pore.
The TKE 
$k$, defined as 
${({\overline {{{{u'_x}}^2}}+\overline {{{u'_y}}^2}+\overline {{{u'_z}}^2}})}/2$, is normalized by the square of the mean interstitial velocity 
$({\overline {\langle u_x\rangle ^f}})^2$. Here the temporal fluctuation velocity 
$u_x'$ is defined as 
$u_x -\overline {u_x}$. Figures 2(c) and 2(d) illustrate the distribution of normalized TKE on the centre 
$xy$- and 
$yz$- planes, respectively. The overall magnitude of the TKE in the centre pore region remains substantially high. It also suggests that the flow through the pore, although bounded by curved particle walls, is different from simple channel/duct flows even with complex boundaries (Orlandi, Davide & Pirozzoli Reference Orlandi, Davide and Pirozzoli2018) wherein the TKE reaches a peak value near the boundaries and then decreases in the centre region. The normalized TKE distributions on the centre 
$yz$-plane presented in figure 2(d) shows some similarities with flow through a duct. Away from the wall, the TKE increases and reaches a peak value; and then decreases as the core region of the pore is approached. The homogeneous particle packing involving four bead particles aligned together seem to form a locally duct-like flow pattern in this section. However, it disappears in other sections away from the centre 
$yz$-plane, owing to the 3-D nature of the spherical particles. Prior work (He et al. Reference He, Apte, Finn and Wood2019) focused on exploring the turbulence characteristics within the FCC unit cell for three different Reynolds numbers (
$Re_H = 300$, 500 and 1000). The detailed TKE budgets, production-dissipation mechanisms within the porous bed, turbulence anisotropy characteristics and distribution were investigated in detail. Specifically, it was shown that for these porous geometries, there was an important region near the bead surfaces with net negative production of the TKE, typically observed in jet-impingement-like flows. In addition, it was also found that the TKE levels were high within the pore region. These flow features illustrated the complexity of the turbulent flow within the porous cell.
The turbulence intensity components in streamwise and spanwise directions are estimated to illustrate the anisotropy caused by the confined geometry, they are shown in the 
$xy$-plane in figure 3. The 
$x$-component r.m.s. velocity, 
$u'_{x-rms}$ is computed by
\begin{equation} u'_{x-rms} = \sqrt{\overline{u_x^2}-{\bar{u}_x}^2}, \end{equation}
likewise for 
$u'_{y-rms}$ and 
$u'_{z-rms}$. The 
$x$-component turbulence intensity is only heavily pronounced in the region close to the right-hand side particle (region (i) in figure 3a), which is mostly attributed to the impingement-like flow in this region. In this centre 
$xy$-plane the 
$x$-component of turbulence intensity, corresponding to the main flow direction, is weaker compared with the other two components. The maximum value of 
$x$-component turbulence intensity is approximately 80 % and 50 % of that for 
$y$- and 
$z$-components, respectively, by Patil & Liburdy (Reference Patil and Liburdy2013). On the centre 
$xy$-plane, these two components are distributed in various patterns. The 
$y$-component mostly concentrates near the left-hand side sphere in figure 3(b). However, for the 
$z$-component, the distribution is substantially different. There are two distinct regions with high magnitude of 
$z$-component turbulence intensity in figure 3(c) near the left-hand side of the pore. More importantly, the shape of such regions is almost identical to the high TKE regions illustrated in figures 2(c), indicating that the 
$z$-component turbulence intensity has the most important contribution to total TKE at this particular plane. Although not shown here, it is noteworthy that the distributions of 
$y$- and 
$z$-components of turbulence intensity on the centre 
$xy$-plane are interchanged on the centre 
$xz$-plane, due to the homogeneous and symmetric geometry of the FCC packing; and the 
$x$-component distribution remains the same.
Figure 3. Contours of turbulence intensity components normalized by the square of interstitial velocity on the centre 
$xy$-plane: (a) 
$x$-component, (b) 
$y$-component and (c) 
$z$-component. Mean flow is from left to right.
To investigate the overall characteristics of both the streamwise and spanwise components of turbulence intensity and for purposes of comparison, intrinsic spatial average of these r.m.s. velocities (normalized by the square of intrinsic averaged velocity) was computed. Because of the periodicity and symmetry in the computational domain, the volume-averaged r.m.s. velocity in the 
$z$-direction is similar to that in the 
$y$-direction and, hence, not included here. It is observed that the r.m.s. velocity in the spanwise direction (
$y$- or 
$z$-direction) has a larger magnitude (0.288) than that in the 
$x$-direction (0.222). This behaviour is quite different from the well-studied turbulent channel or duct flows even with complex boundary shapes (Orlandi et al. Reference Orlandi, Davide and Pirozzoli2018), wherein the r.m.s. velocities in the wall-normal directions are much smaller than the streamwise component. This again is caused by the 3-D effect of the complex configuration of the packed spheres, and the tortuous mean flow patterns within the pore. Qualitatively, a similar behaviour that larger values of turbulence intensities in non-streamwise direction has also been observed experimentally in randomly packed porous media (Patil & Liburdy Reference Patil and Liburdy2013).
Finally, the Eulerian and Lagrangian autocorrelations are used to compute the integral length and time scales, respectively. To compute the Lagrangian autocorrelations, fluid tracer particles are tracked to obtain the Lagrangian trajectories. The Lagrangian autocorrelations are then computed according to (3.2) (Monin & Yaglom Reference Monin and Yaglom1965)
\begin{equation} \rho^L_{ij} (\tau) = \frac{{ {\langle {v}'_i({ X_0}, t) \, {v}'_j({ X_0} ,t+\tau) \rangle} }}{{ {\left[ \langle {{v_i}'^2}({ X_0}, t)\rangle\,\langle {{v_j}'^2}({ X_0}, t+\tau)\rangle \right]^{1/2}} }}, \end{equation}
where 
$\rho ^L_{ij}$ is the Lagrangian autocorrelation, 
$v'_i$ the 
$i$th component of the particle fluctuation velocity and 
$\langle \cdot \rangle$ represents ensemble averaging. The Lagrangian integral time scale, 
$T_{11}^{L}$, is simply given by the integral of the autocorrelation function.
The longitudinal Eulerian integral length scale (
$L_{11}^E$) is computed by integrating the Eulerian autocorrelations over the abscissa. Then it is normalized by the bead diameter (
$D_B$) and found to be 
$0.0884$. The flow time scale, defined as 
${D_B}/{\overline {\langle u_x \rangle ^f}}$, is used to normalize the Lagrangian integral time scale 
$T_{11}^L$, which is found to be approximately 0.356 for the present flow conditions. The integral length scale is only approximately 10 % of the sphere diameter, indicating that the coherent structures are confined within the pore. Such observation supports the pore-scale prevalence hypothesis and the results reported in Jin et al. (Reference Jin, Uth, Kuznetsov and Herwig2015). This implies that, the turbulence in the pore scale is strongly affected by the porosity, and restrained by the pore size. The integral time scale is also smaller than the flow time scale, suggesting that the Lagrangian coherent structures are restricted by the pore size. As a result, the single periodic unit cell considered in the present work is sufficient. The overall dissipation rate 
$\langle \epsilon \rangle$ was estimated from the TKE budget first. Then the Kolmogorov time scale is computed as 
$\tau _\eta =\sqrt {{\nu }/{\langle \epsilon \rangle }}$.
For the present case of 
$Re_H=500$, 
$Re_\lambda$ can be estimated to be around 32, which is computed as the drape of the fitted parabola to the Eulerian autocorrelation. Assuming isotropic turbulence, which is not the case for the flow in porous media, and using the same definition as in Matsuda et al. (Reference Matsuda, Schneider and Yoshimatsu2021), an even smaller value of 
$Re_\lambda \sim 21$ is obtained for the Taylor microscale Reynolds number. Particle clustering dynamics in HIT for a wide range of 
$Re_{\lambda }$ has been studied by Sumbekova et al. (Reference Sumbekova, Cartellier, Aliseda and Bourgoin2017) showing significant impact on the cluster and void size. The flow through a porous unit cell is highly anisotropic; however, as is shown later, the qualitative behaviour of clustering characteristics observed in the present work are similar to HIT, e.g. at higher Stokes numbers, the energy density is higher and the peak shifts to larger scales. For higher 
$Re_{\lambda }$, qualitatively similar trends are expected. The main focus of the present work is on investigating the effect of geometric confinement and interaction of inertial particles with the bead surfaces on the clustering dynamics. Effect of 
$Re_{\lambda }$, volume loading and porous bead arrangements will be part of future studies.
3.2. Inertial particle dynamics and clustering
Figure 4 shows the instantaneous distribution of particles for three different Stokes numbers in the 
$xy$-plane after a stationary state is reached. A uniform random distribution of particles in the porous geometry, representative of fluid particles in the limit of 
$St\rightarrow 0$, is also shown for comparison. Significant particle clustering is observed for 
$St=1$ as expected, whereas for 
$St=0.01$, particles are more uniformly distributed with only a few pockets of voids and clusters near the bead boundaries. The clustering of inertial particles and effect of the bead walls is evaluated by conducting a multiscale analysis of the particle number density. Figure 4 shows a small gap around the bead surface where inertial particles are not observed. Although the fluid flow computations were performed using exact radii of the bead surfaces and without any gap, to ensure that the inertial particles are not trapped inside the solid bead region, a small gap corresponding to the grid resolution is used to enforce particle collisions with the bead surface.
Figure 4. Instantaneous distribution of particles in the 
$xy$-plane: (a) 
$St=1$, (b) 
$St=0.1$, (c) 
$St=0.01$ and (d) uniform random distribution.
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,
\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
\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
\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
\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. Eulerian field: particle number density
The number density of the discrete particle positions is obtained using a histogram method by binning particles onto an equidistant cubic grid of 
$N_g^3$ grid points (Matsuda et al. Reference Matsuda, Schneider and Yoshimatsu2021). Irrespective of the grid resolution used in the DNS, 
$N_g=2^8$ was used to calculate the number density as
\begin{align} n({\boldsymbol x},t) = \sum_{i_1,i_2,i_3=0}^{N_g-1}\left\{\int_{\varOmega}K_h({\boldsymbol x}_{i_1,i_2,i_3}-{\boldsymbol x}^{\prime}) \frac{1 }{n_0}\sum_{m=1}^{N_p}\delta({\boldsymbol x}^{\prime} - {\boldsymbol x}_{p,m}(t))\,{\rm d}{\boldsymbol x}^{\prime}\right\}h^3 K_h({\boldsymbol x}-{\boldsymbol x}_{i_1,i_2,i_3}), \end{align}
where 
${\boldsymbol x}_{i_1,i_2,i_3}=h(i_1+1/2,i_2+1/2,i_3+1/2)$ is the box position, and 
$K_h({\boldsymbol x}) = 1/h^3$ for 
$-h/2\leq x_i \leq h/2$ (
$i=1,2,3)$, while 
$K_h({\boldsymbol x})=0$ otherwise, is a piecewise constant function, 
$h=L/N_g$, and 
$n_0 = N_p/L^3$ is the mean dimensional number density, where 
$N_p$ is the total number of particles and 
$L$ is the side length of the cubic domain. With non-dimensionalization by the mean number density, the above equation satisfies 
$\left \langle n\right \rangle =1$.
It should be noted that the entire computational domain volume is used to compute the number density, even though some part of the volume is a solid bead region to simplify the number density computation. The porosity of the unit cell can be used to relate this to the number density calculated based on the fluid volume only. The consequence of the presence of solid beads within the computational domain is that a uniform distribution of inertial particles in the fluid domain results in non-uniform number density variations across the bead surface. These gradients in number density, even for a uniform inertial particle distribution, can result in non-zero wavelet decomposition. To analyse the true clustering of inertial particles, this effect of pseudo variations in number density due to the bead geometry needs to be removed before performing the wavelet decomposition. Therefore, the number density gap from 
$n({\boldsymbol x})$ is subtracted as
\begin{equation} n^\prime({\boldsymbol x}) = n({\boldsymbol x}) - \phi^{{-}1} \langle n \rangle~\left( 1- \chi({\boldsymbol x}) \right), \end{equation}
with 
$\chi (\boldsymbol {x})$ the mask function,
\begin{equation} \chi(\boldsymbol{x}) = \left\{ \begin{array}{@{}ll} 1 & \text{if }\boldsymbol{x} \in \varOmega_s,\\ 0 & \text{if }\boldsymbol{x} \in \varOmega_f,\end{array} \right. \end{equation}
where 
$\varOmega _s$ is the solid domain, 
$\varOmega _f$ the fluid domain and 
$\varOmega =\varOmega _f \cup \varOmega _s$ the computational domain. For the sake of clarity, 
$n^\prime ({\boldsymbol x})$ is denoted by 
$n({\boldsymbol x})$.
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,
\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
\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
\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
\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
\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
\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
\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:
\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
\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
\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
\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.
4. Summary and conclusion
Fully resolved DNS of a turbulent flow through the confined geometry of a triply periodic, FCC porous unit cell was performed using a Cartesian grid-based fictitious domain method. The low porosity of the flow geometry gives rise to rapid acceleration and deceleration of the mean flow with the presence of 3-D helical motions, weak wake-like structures behind the bead spheres, stagnation and jet-impingement-like flows together with merging and spreading jets in the main pore. Details of turbulence characteristics in this confined geometry of the porous cell for a range of Reynolds numbers spanning unsteady inertial, transitional and turbulent flow were characterized in detail in prior work (He et al. Reference He, Apte, Finn and Wood2019). In this work emphasis is placed on clustering dynamics of inertial particles in a turbulent flow through the porous cell. Specifically, how particle–wall interactions affect clustering and deposition mechanisms at different particle Stokes number were studied. To this end, point particles were advanced for five different Stokes numbers 
$St$ using one-way coupling and assuming perfectly elastic wall collisions for a single pore Reynolds number of 
$Re_H = 500$, corresponding to inertial, turbulent flow. About ten million particles for each Stokes number were introduced into the flow and tracked over several flow through times to provide meaningful statistics on inertial particle dynamics of clustering, void formation and transport inside the confined geometry of the porous medium.
Tools for studying inertial particle dynamics and clustering, previously developed for homogeneous flows, have been adapted by taking into account the curved flow geometry of the bead walls in the porous cell. Mirror particles were used in the Voronoi analysis and adjusted for tesselation in the presence of bead walls. The probability distribution of the Voronoi volumes for the different Stokes numbers quantified the departure from the gamma distribution and allowed us to assess the influence of both the geometry and the flow, in comparison to isotropic turbulence. It was found that the geometry only impacts the small volumes below 
$10^{-2}$ to the mean Voronoi volume. Cluster formation and destruction was quantified by analysing the time change of the Voronoi volumes for the different Stokes numbers. This Lagrangian approach yields a time discrete measure of the spatial divergence of the particle velocity. It was found that the p.d.f.s of the divergence, which are symmetric, are becoming wider for increasing Stokes number, as the variance increases. In contrast, the tails of the p.d.f.s are becoming shallower with Stokes number. The conditional average of the divergence as a function of volume is found to be positive for small volumes and negative for large volumes. This explains that cluster formation is present for large volumes, while cluster destruction is more prominent for small volumes. Moreover, these effects are amplified with the Stokes number. These findings are similar to what has been found for isotropic turbulence (Oujia et al. Reference Oujia, Matsuda and Schneider2020). Wavelet-based multiscale statistical analyses were applied to particle number density fields in the flow through the porous geometry. By decomposing the number density fields into orthogonal wavelets, scale-dependent statistics of the number density distribution were computed. The wavelet energy spectra showed that the peak of clustering gradually shifts towards a larger scale as the Stokes number increases. To reduce the influence of bead geometry, the difference of each energy spectrum for inertial particles from that for randomly distributed particles only in a fluid domain was also computed. This allowed us to identify the peaks in the spectra for 
$St \le 0.1$. Scale-dependent skewness and flatness of the particle density quantified the intermittent void and cluster distribution statistically. The positive skewness values at smaller scales found for all Stokes numbers confirm the observed small-scale prominent clusters. Negative skewness values for 
$St >0.5$ quantify the presence of prominent void regions. The flatness values which are increasing with decreasing scale and that the values become even larger for larger Stokes numbers confirm the strongly intermittent cluster distribution.
In conclusion both static and dynamic analyses of particle clustering in a porous cell have been performed. With scale-dependent analyses of snapshots of particle density distributions (static analyses), voids and clusters were quantified statistically. Focusing on intermittency, a signature of void and clusters was observed in higher-order statistics. With the Lagrangian analysis of Voronoi tesselations (dynamic analysis), the convergence and divergence of particle velocity were computed thus providing an explanation for cluster formation and destruction. A comparison with results for HIT, showed many similarities and also pointed out differences due to the flow geometry, in particular for small volumes.
Voronoi and wavelet analysis of the particle density are somewhat connected. It was shown that Voronoi tesselation of the particle positions contains length scale information, reflected by small and large cells. Wavelet analysis of the corresponding particle density field, which is obtained by averaging the particle positions (using histograms), can also detect the presence of small scales and large scales, reflected in strong wavelets coefficients measuring variation (or gradients) in the particle density. Thus thresholding of the wavelet coefficients would yield an adaptive grid, as shown in Bassenne et al. (Reference Bassenne, Urzay, Schneider and Moin2017). Such locally refined grids, are expected to be present in regions of small Voronoi cells and coarse grids are expected in regions of large Voronoi cells. A qualitative agreement can thus be expected. In contrast to the wavelet representation, in the Voronoi analysis, the scales are not separated, all scales are simultaneously given by the tesselation. Combining the multiscale statistical analyses with the Lagrangian formulation of the Voronoi tesselation constitutes an interesting perspective of this work for quantifying scale-dependent divergence and convergence of the particle velocity and the related inter-scale transfer. The approach developed in this work should be directly applicable to cases involving polydispersed, and randomly packed beads with periodic conditions. Effect of flow Reynolds number, volume loading, particle-size relative to the Kolmogorov scale and different types of particle–bead interactions can affect the clustering dynamics in these confined configurations and will be the focus of future work. The analysis presented can also be extended in the future to vorticity fields to obtain 3-D directional information at different scales.
Acknowledgements
This work was initiated during S.V.A.'s visiting scientist/Professor position at Aix-Marseille Université, and CNRS-I2M, France. S.V.A. acknowledges kind hospitality during his visit and stay in Marseille. Simulations were performed at the Texas Advanced Computing Center's (TACC) Stampede2 and Frontera systems. B.K., T.O. and K.S. thankfully acknowledge Centre de Calcul Intensif d'Aix-Marseille for granting access to its high performance computing resources.
Funding
S.V.A. gratefully acknowledges CNRS, France and Aix-Marseille Université for the visiting scientist/Professor position at Aix-Marseille Université. Partial funding from US Depart of Energy, Office of Basic Energy Sciences (Geosciences) under award number DE-SC0021626 as well as US National Science Foundation award #205324 is acknowledged. X.H. thanks Dr T. Scheibe from Pacific Northwest National Laboratory (PNNL) for supporting funding from the DOE Office of Biological and Environmental Research, Subsurface Biogeochemical Research program, through the PNNL Subsurface Science Scientific Focus Area project (http://sbrsfa.pnnl.gov/). T.O. and K.S. acknowledge financial support from Agence Nationale de la Recherche, project ANR-20-CE46-0010-01.
Declaration of interests
The authors report no conflict of interest.
