# Clustering of Inertial Particles in Turbulent Flow Through a Porous Unit
Cell

Sourabh V. Apte1<EMAIL_ADDRESS>Thibault Oujia2  Keigo
Matsuda3  Benjamin Kadoch4  Xiaoliang He5  Kai Schneider2  1School of
Mechanical, Industrial and Manufacturing Engineering, Oregon State University,
Corvallis, OR 97330, USA 2 Aix–Marseille Université, I2M–CNRS, Marseille,
France 3 Research Institute for Value-Added-Information Generation (VAiG),
Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokohama
236-0001, Japan 4 Aix–Marseille Université, IUSTI–CNRS, Marseille, France 5
Pacific Northwest National Laboratory, Richland, WA 99352, USA

###### Abstract

Direct numerical simulation is used to investigate effects of turbulent flow
in the confined geometry of a face-centered cubic porous unit cell on the
transport, clustering, and deposition of fine 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 collision of particles with
pore walls is modeled as perfectly elastic with specular reflection. Tools for
studying inertial particle dynamics and clustering developed for homogeneous
flows are adapted to take into account the embedded, curved geometry of the
pore walls. The pattern and dynamics of clustering are investigated using the
volume change of Voronoi tesselation in time to analyze the divergence and
convergence of the particles. Similar to the case of homogeneous, isotropic
turbulence, the cluster formation is present at large volumes, while cluster
destruction is prominent at small volumes and these effects are amplified with
Stokes number. However, unlike homogeneous, isotropic turbulence, formation of
large number of very small volumes was observed at all Stokes numbers and is
attributed to the collision of particles with the pore wall. Multiscale
wavelet analysis of the particle number density showed peak of clustering
shifts towards larger scales with increase in Stokes number. Scale-dependent
skewness and flatness quantify the intermittent void and cluster distribution,
with cluster formation observed at small scales for all Stokes numbers, and
void regions at large scales for large Stokes numbers.

##  1 Introduction

Several applications including catalysis in the chemical synthesis and process
industries (Dixon & Nijemeisland, 2001; Aris, 1999), high temperature nuclear
reactor cooling (Shams et al., 2013), hyporheic exchange of pollutants at the
sediment-water interface (Hester et al., 2017), 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 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., 1992). 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 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., 2004), and onset of
sanding triggered by rock failure (Yi et al., 2005). 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, 1974; Mahmud et al., 2019).
Particle rentention, bridge formation, jamming, and deposition has been
studied experimentally in simplified configurations such as microchannels
(Valdes & Santamarina, 2006; Dai & Grace, 2010; Agbangla et al., 2012) and
packed beds (Ramachandran & Fogler, 1999; Pandya et al., 1998) 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 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., 2019) explored, using direct
numerical simulation (DNS), how this change of flow structure impacts the
turbulence (i.e. the turbulent kinetic energy distribution and transport
mechanisms) in a face-centered 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 presence of three-dimensional (3D) 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 turbulent kinetic energy, 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 and wall-bounded channel flows (Maxey, 1987;
Eaton & Fessler, 1994; Toschi & Bodenschatz, 2009; Monchaux et al., 2012). 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 (Esmaily-Moghadam &
Mani, 2016; Bec et al., 2007). 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 et al. (2020) 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 analyze
clustering as it is directly related to the particle collision rates (Sundaram
& Collins, 1997; Wang et al., 2000). The PCF typically shows a power-law
behavior 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 (2006) in which particles
are swept by large-scale flow motion while sticking to stagnation points of
Lagrangian fluid acceleration (Coleman & Vassilicos, 2009). Probability
distribution function (PDF) of particle mass density and coarse graining
techniques have been used to investigate scale dependence of particle
distribution (Bec et al., 2007). Bassenne et al. (2017) 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 wave number contributions of
turbulent flow fields. The wavelet representation can be used to analyze
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 a review on wavelets and turbulence
refer to Farge (1992); Schneider & Vasilyev (2010) and more particularly on
wavelet-based statistics to Farge & Schneider (2015). Recently, Matsuda et al.
(2021) 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 effect of turbulent flow in the
confined geometry of a face-centered cubic 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 collision of particles with pore walls is
modeled 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 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., 2020; Matsuda et al., 2021). 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 Section 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 Section 4.

(a)

(b)

Figure 1: (a) Schematic of a face-centered 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 Simulation setup

Several different definitions of Reynolds numbers have been used in porous
media literature (He et al., 2019; Wood et al., 2020) based on different
length scales such as the particle diameter ($D_{B}$), the hydraulic diameter
($D_{H}$), or the permeability. In this work, Reynolds number ($Re_{H}$) based
on $D_{H}$ is used. The hydraulic diameter is related to the particle diameter
as,

$D_{H}={\frac{2}{3}}{\frac{\phi}{1-\phi}}D_{B}.$ (2.1)

Then, $Re_{H}$ is defined as (dropping the factor 2/3),

$Re_{H}=\frac{\overline{\langle u_{x}\rangle^{f}}D_{B}}{\nu_{f}}$ (2.2)

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 $x$-(streamwise)
direction. The spatial averaging operation is denoted by
$\langle~{}\cdot~{}\rangle$, the superscript or subscript $f$ indicates the
fluid phase, and $\overline{~{}\cdot~{}}$ is the time averaging operator. For
clarity, note that the intrinsic average velocity is defined by,

$\left.\langle u_{x}\rangle^{f}\right|_{t}=\frac{1}{V_{f}}\int_{{\bf
x}\in\Omega_{f}}u_{x}({\bf x},t)\,dV({\bf x}),$ (2.3)

where $V_{f}$ is the volume of the fluid phase within the unit cell and
$\Omega_{f}$ denotes the fluid domain. The overbar represents traditional time
averaging of a quantity $q$ defined as,

$\overline{q}(t)=\frac{1}{T}\int\limits_{~{}~{}~{}~{}~{}~{}t^{*}=t-T/2}^{~{}~{}~{}~{}~{}~{}t^{*}=t+T/2}q(t^{*})dt^{*}.$
(2.4)

###  2.1 Porous geometry

A porous face-centered cubic (FCC) unit cell (Fig. 1) is used based on our
prior work on turbulence in porous media (He et al., 2018, 2019). It has a
half sphere entering at each face of the cube, and a half quarter sphere at
each corner. The face-centered cubic arrangement creates the lowest possible
porosity ($\phi$) to be 0.26 for the structured packings. 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.
Majority of the flow enters the cubic cell through the upstream open corners,
converges into the center 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 direct numerical simulation, 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.

###  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 (Ref. Apte et al. (2009)). 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 (Ref. Finn & Apte (2013)). The following
governing equations are solved over the entire domain, including the region
within the solid bed, and a rigidity constraint force, $\bf f$, is applied
that is non-zero only in the solid region.

The governing equations read as:

$\displaystyle\nabla\cdot{\bf u}$ $\displaystyle=0,$ (2.5)
$\displaystyle\rho_{f}\bigg{[}\frac{\partial{\bf u}}{\partial t}+\left({\bf
u}\cdot\nabla\right){\bf u}\bigg{]}$ $\displaystyle=-\nabla
p+\mu_{f}\nabla^{2}{\bf u}+{\mathbf{f}}\>,$ (2.6)

where $\bf u$ is the velocity vector (with components given by ${\bf
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 (Ref. Apte et al. (2009)). 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 (Ref. Finn & Apte (2013); He
et al. (2018, 2019)).

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 (2002), a constant
pressure gradient $\nabla P$ in the main direction of the flow ($x$
direction), proportional to the body force $F$, is used to drive the flow,

$\displaystyle\partial_{x}P$
$\displaystyle=\frac{18\>\mu_{f}\>c\>U_{int}}{D_{B}^{2}}F\>,$ (2.7)

where, $\mu_{f}$ is the dynamic viscosity of the fluid and $c$ the solid
volume fraction, defined as $\frac{2}{3}\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 (2002) and is given as,

$\displaystyle F$
$\displaystyle=462+9.85\>\left(\frac{1-\phi}{2}Re_{H}\right)\>~{}~{}~{}(Re_{H}\geq
216)\>.$ (2.8)

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., 2019).
A uniform, cubic grid is used with resolution chosen such that the first grid
point is at $y^{+}<1$ (i.e., in the viscous sub-layer) 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, 3D
DNS studies were performed at $Re_{H}=500$ in a unit cell of face-centered
cubic (FCC) spheres with 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 (turbulent kinetic energy, 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 to 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., 2019).

###  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 modeled using simple closure models. In the present work, the
particle motion is assumed to be governed by a simple, linear drag force. For
different Stokes numbers ($St$), the particle motion is given by,

${d{\mathbf{x}}_{p}\over dt}={\mathbf{u}}_{p};~{}~{}~{}{d{\mathbf{u}}_{p}\over
dt}={{\mathbf{u}}_{f,p}-{\mathbf{u}}_{p}\over St\cdot\tau_{\eta}},$ (2.9)

where ${\mathbf{x}}_{p}$ and ${\mathbf{u}}_{p}$ are particle position and
velocity, $\tau_{\eta}$ is Kolmogorov time scale, ${\mathbf{u}}_{f,p}$ is the
instantaneous fluid velocity interpolated to the particle location, i.e.,
${\bf u}_{f,p}={\bf u}({\bf 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 homogeneous isotropic turbulence (Matsuda et
al., 2021; Oujia et al., 2020).

Particles with five different Stokes numbers, $St=$ 0.01, 0.1, 0.5, 1, and 2
are simulated. 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. 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 center 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 modeled using Snell’s law of specular
reflection assuming perfectly elastic collisions. The direction of particle
reflection, ${\hat{\mathbf{s}}}_{r}$, is determined from the incident
direction, ${\hat{\mathbf{s}}}_{i}$, and the inward face normal at the
spherical bead surface, ${\hat{\mathbf{n}}}$ as,

${\hat{\mathbf{s}}}_{r}={\hat{\mathbf{s}}}_{i}+2|{\hat{\mathbf{s}}}_{i}\cdot{\hat{\mathbf{n}}}|{\hat{\mathbf{n}}}\,.$
(2.10)

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 center of the bead.

##  3 Results and discussion

In this section, the mean and turbulent flow structure is described briefly
using the the mean and rms (root mean square) velocity fields inside the pore,
and the integral length and time scales are estimated based on the Eulerian
and Lagrangian two-point auto-correlations, respectively. Next, the
probability distribution functions 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 setup, the flow through the triply
periodic domain is driven by a body force in the flow direction based on the
correlations given by Eq. (2.8). After an initial transient, a stationary
state is reached and the computation is continued for several flow through
times, $T_{f}=L/U_{int}$ 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 about $100T_{f}$ for $Re_{H}=500$.
After a stationary state has been established, the computations were performed
for additional $80T_{f}$ for each Reynolds number 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
turbulent kinetic energy 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 center slices are chosen as representative sections
for visualizations. One is the center $xy-$plane, where the mean flow is going
from left to right; the other the center $yz-$plane, where the mean flow in
going into the slice.

(a)

(b)

(c)

(d)

Figure 2: Visualization of the mean velocity magnitude (top) and turbulent
kinetic energy (bottom) at $Re_{H}=500$: (a–c) $xy-$plane, and (b–d)
$yz-$plane.

The mean flow enters the center pore from the left-hand side corners (regions
labeled by (i1) and (i2) in figure 2a) and accelerates to region (ii) due to
the geometric constrictions; then the mean flow at the center 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 (iv1) and (iv2)). 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 2b shows the distributions of the mean
velocity magnitude on the center $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 center pore.

The turbulent kinetic energy $k$, defined as
${\frac{1}{2}\left({\overline{{{{u^{\prime}_{x}}}^{2}}}+\overline{{{u^{\prime}_{y}}}^{2}}+\overline{{{u^{\prime}_{z}}}^{2}}}\right)}$,
is normalized by the square of the mean interstitial velocity
$({\overline{\langle u_{x}\rangle^{f}}})^{2}$. Here the temporal fluctuation
velocity $u_{x}^{\prime}$ is defined as $u_{x}-\overline{u_{x}}$. Figure 2(c)
and 2(d) illustrate the distribution of normalized TKE on the center $xy-$ and
$yz-$ planes, respectively. The overall magnitude of the TKE in the center
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 et al., 2018) wherein
the TKE reaches a peak value near the boundaries and then decreases in the
center region. The normalized TKE distributions on the center $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 particles aligned together seem to form a
locally duct-like flow pattern in this section. However, it disappears in
other sections away from the center $yz-$plane, owing to the three-dimensional
nature of the spherical particles.

(a)

(b)

(c)

Figure 3: Contours of turbulence intensity components normalized by the square
of interstitial velocity on the center $xy-$plane: (a) $x$-component, (b)
$y$-component, and (c) $z$-component. Mean flow is from left to right.

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 rms velocity,
$u^{\prime}_{x-rms}$ is computed by

$u^{\prime}_{x-rms}=\sqrt{\overline{u_{x}^{2}}-{\overline{u}_{x}}^{2}},$ (3.1)

likewise for $u^{\prime}_{y-rms}$ and $u^{\prime}_{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 3(a)), which is mostly
attributed to the impingement-like flow in this region. In this center
$xy-$plane, the $x$-component of turbulence intensity, corresponding to the
main flow direction, is weaker compared to the other two components. The
maximum value of $x$-component turbulence intensity is about 80% and 50% of
that for $y$\- and $z$-components, respectively. by Patil & Liburdy (2013). On
the center $xy-$plane, these two components are distributed in various
patterns. The $y$-component mostly concentrates near the left-hand-side
particle 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 center $xy-$plane are interchanged on the center $xz-$plane, due to the
homogeneous and symmetric geometry of the FCC packing; and the $x$-component
distribution remains the same.

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 rms velocities (normalized by the square of intrinsic
averaged velocity) was computed. Because of the periodicity and symmetry in
the computational domain, the volume-averaged rms velocity in $z$-direction is
similar to that in $y$-direction and hence not included here. The turbulence
intensity in both streamwise ($x$-) and spanwise ($y$-)directions are
increasing with Reynolds number, as expected. However, the rms velocity in the
spanwise direction ($y$\- or $z$-) has a larger magnitude (0.288) than that in
the $x$-direction (0.222). This behavior is quite different from the well-
studied turbulent channel or duct flows even with complex boundary shapes
(Orlandi et al., 2018), wherein the rms velocities in the wall normal
directions are much smaller than the streamwise component. This again is
caused by the three-dimensional effect of the complex configuration of the
packed spheres, and the tortuous mean flow patterns within the pore.
Qualitatively, a similar behavior that larger values of turbulence intensities
in non-streamwise direction has also been observed experimentally in randomly-
packed porous media (Patil & Liburdy, 2013).

Finally, the Eulerian and Lagrangian auto-correlations are used to compute the
integral length and time scales, respectively. To compute the Lagrangian auto-
correlations, fluid tracer particles are tracked to obtain the Lagrangian
trajectories. The Lagrangian auto-correlations are then computed according to
Eq. (3.2) (Ref. Monin & Yaglom (1965)),

$\displaystyle\begin{aligned}
\rho^{L}_{ij}(\tau)=\frac{{{\langle{v}^{\prime}_{i}({X_{0}},t)\>\>{v}^{\prime}_{j}({X_{0}},t+\tau)\rangle}}}{{{\Big{[}\langle{{v_{i}}^{\prime
2}}({X_{0}},t)\rangle\>\langle{{v_{j}}^{\prime
2}}({X_{0}},t+\tau)\rangle\Big{]}^{1/2}}}}\end{aligned}$ (3.2)

where $\rho^{L}_{ij}$ is the Lagrangian auto-correlation, $v^{\prime}_{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 auto-
correlation function.

The Eulerian integral length scale ($L_{11}^{E}$) normalized by the bead
diameter ($D_{B}$) is found to be $0.0884$, whereas the Lagrangian integral
time scale normalized by the flow time scale ($T_{11}^{L}$) based on the
interstitial velocity ($U_{int}$) and the bead diameter ($D_{B}$) is about
0.356 for the present flow conditions. The integral length scale is only about
10% of the sphere diameter, indicating that the coherent structures are
confined within the pore. Such observation supports the pore scale prevalence
hypothesis (PSPH) and the results reported in Jin et al. (2015). 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 case
$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 auto-correlation.
Assuming isotropic turbulence, which is not the case for the flow in porous
media, and using the same definition as in Matsuda et al. (2021), an even
smaller value of $Re_{\lambda}\sim 21$ is obtained for the Taylor microscale
Reynolds number.

###  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. 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 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.

(a)

(b)

(c)

(d)

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 (1991), is a technique to construct
a decomposition of the fluid domain, into a finite number of Voronoi cells. If
there are 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. 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 three-dimensional space, smaller
volume indicating pronounced clustering. Three-dimensional 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 et al., 1996), which has a computational complexity of
${\mathcal{O}}(N_{p}{\rm log}N_{p})$.

Figure 5:  Particle distribution for $St=1$ in a slice of thickness $1/100$
with spheres and mirror particles (left). A magnified view with Voronoi
tessellation (right).  Figure 6: PDF of volumes of Voronoi cells 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 PDF of randomly distributed particles is also compared with the gamma
distribution and shows perfect agreement.

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 PDF of the Voronoi volume
becomes a gamma distribution (Ferenc & Néda, 2007). For the 3D case, the PDF
of Voronoi volume is given by $\Gamma(5,1/5)$, where $\Gamma(k,\theta)$
corresponds to the gamma distribution,

$f_{V_{p}}(x)=\Gamma(k)^{-1}\theta^{-k}x^{k-1}{\rm exp}(-x/\theta),$ (3.3)

where $k$ and $\theta$ are shape parameters of the PDF, respectively. Figure 6
shows the PDF 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 PDF 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 PDF to the gamma
distribution, the sphere radius had to be increased by about 1.1% ($0.35749$
instead of $0.353553$) 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 shows that the PDFs of Voronoi volumes for different Stokes numbers
intersect with the gamma distribution, i.e. the one for the random particles.
The observed behavior is similar to that found for homogeneous isotropic
turbulence (HIT), see (Oujia et al., 2020), but with some key differences. The
number of large Voronoi cells increases with increasing Stokes numbers and
then stabilizes. This behavior is resembling 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\cdot 10^{-1}$) also increases similar
to the HIT case. However, 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 behavior 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. 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 large number of very small volumes is mainly attributed to the
collision of the particles with the bead surfaces. 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.,
2020). It satisfies the conservation equation,

$\frac{D}{Dt}(n)=\frac{\partial n}{\partial t}+{\mathbf{u}}_{p}\cdot\nabla
n=-n\,\nabla\cdot{\mathbf{u}}_{p}\,.$ (3.4)

The divergence of the particle velocity ${\mathbf{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, ${\mathbf{u}}_{p,j}={\mathbf{u}}({\mathbf{x}}_{p,j})$, but
not everywhere else. Moreover it can be multi-valued. Oujia et al. (2020)
proposed a method to compute the divergence of the particle velocity
${\mathcal{D}}=\nabla\cdot{\bf u}_{p}$ in a discrete manner using a Lagrangian
approach. From the conservation equation for the particle number density (eq.
3.4), one obtains, ${\mathcal{D}}=-\frac{1}{n}\frac{Dn}{Dt}$.

To calculate the Lagrangian derivative of $n$, Oujia et al. (2020) 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,

${\cal D}_{p}=\frac{2}{\Delta
t}\,\frac{V^{k+1}_{p}-V^{k}_{p}}{V^{k+1}_{p}+V^{k}_{p}}+O(\Delta t)\;,$ (3.5)

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
${\mathbf{u}}_{p}$; i.e,
${\mathbf{x}}_{p}^{k+1}={\mathbf{x}}_{p}^{k}+{\mathbf{u}}_{p}\Delta t$. The
time-step was set to be 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 same result as in Oujia et al. (2020). 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$.

Figure 7: PDF of divergence of particle velocity for different Stokes numbers.

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.35749$ from each of the beads centers were not taken into account. 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.35749$, the values were found to be stable and remained almost constant.

The PDFs of the Voronoi-based divergence for different Stokes numbers are
shown in figure 7. They are centered around 0 and confined between $-40$ and
$40$. It can be observe 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., 2020). The divergence
should become closer to zero as the Stokes number decreases to zero because
the fluid particles in an incompressible flow are divrgence 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 Stoks numbers is mainly caused by
a geometrical effect due to Voronoi tessellation, which is also discussed in
Oujia et al. (2020). However for larger Stokes numbers, physical effects do
predominate. Table 1 assembles the variance and flatness of the divergence
${\cal 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 PDFs decay faster as the Stokes number increases. In
contrast the variance has the opposite behavior and thus the PDFs 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.

St | 0.01 | 0.1 | 0.5 | 1 | 2
---|---|---|---|---|---
Variance | 0.0083 | 0.0093 | 0.0952 | 0.2796 | 0.5423
Flatness | 5514.9 | 734.50 | 832.00 | 291.92 | 105.70
Table 1: Variance and flatness of divergence ${\cal D}_{p}$ as a function of
the Stokes number.  Figure 8: Mean divergence $\langle{\cal
D}_{p}\rangle_{Vp}$ as a function of the Voronoi volume for different Stokes
numbers.

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

$\langle{\cal
D}_{p}\rangle_{Vp}=\frac{1}{P(V_{p}/\overline{V_{p}})}\int_{-\infty}^{+\infty}{\cal
D}_{p}\,P({\cal D}_{p},V_{p}/\overline{V_{p}})\,d{\cal D}_{p},$ (3.6)

The values are shown in figure 8 for 3 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. (2020) the amplitude
increases with the Stokes number. This implies presence of cluster formation
for large volumes, cluster destruction for small volumes and that these two
behaviors 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. (2020).

In summary, the above findings for the Voronoi analysis confirm the results
found in Oujia et al. (2020) for homogeneous isotropic turbulence, 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., 2021). Irrespective of the grid
resolution used in the direct numerical simulation, $N_{g}=2^{8}$ was used to
calculate the number density as,

$n({\mathbf{x}},t)=\sum_{i_{1},i_{2},i_{3}=0}^{N_{g}-1}\left\\{\int_{\Omega}K_{h}({\mathbf{x}}_{i_{1},i_{2},i_{3}}-{\mathbf{x}}^{\prime}){1\over
n_{0}}\sum_{m=1}^{N_{p}}\delta({\mathbf{x}}^{\prime}-{\mathbf{x}}_{p,m}(t))d{\mathbf{x}}^{\prime}\right\\}h^{3}K_{h}({\mathbf{x}}-{\mathbf{x}}_{i_{1},i_{2},i_{3}}),$
(3.7)

where ${\mathbf{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}({\mathbf{x}})=1/h^{3}$ for $-h/2\leq x_{i}\leq
h/2$ ($i=1,2,3)$, while $K_{h}({\mathbf{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<n\right>=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 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
analyze the true clustering of inertial particles, this effect of pseudo
variations in number density due to the bead geometry, need to be removed
before performing the wavelet decomposition. Therefore the number density gap
from $n({\mathbf{x}})$ is subtracted as

$n^{\prime}({\mathbf{x}})=n({\mathbf{x}})-\phi^{-1}\langle
n\rangle~{}\left(1-\chi({\mathbf{x}})\right)$ (3.8)

with $\chi(\mathbf{x})$ the mask function :

$\chi(\mathbf{x})=\bigg{\\{}\begin{array}[]{l l}1&\quad\text{if
$\mathbf{x}\in\Omega_{s}$ }\\\ 0&\quad\text{if $\mathbf{x}\in\Omega_{f}$ }\\\
\end{array}$ (3.9)

where $\Omega_{s}$ is the solid domain, $\Omega_{f}$ the fluid domain and
$\Omega=\Omega_{f}\cup\Omega_{s}$ the computational domain. For sake of
clarity, $n^{\prime}({\mathbf{x}})$ is denoted by $n({\mathbf{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 (Mallat, 2009; Daubechies, 1993) of
the particle number density. Consider the particle number density,
$n({\mathbf{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 3D orthogonal wavelet series to unfold into scale, positions and seven
directions ($\mu=1,2,....,7$). The 3D mother wavelet,
$\psi_{\mu}({\mathbf{x}})$, is hereby based on a tensor product construction
and a family of wavelets, $\psi_{\mu,{\bm{\lambda}}}({\mathbf{x}})$ can be
generated by dilation and translation. an orthogonal basis of
$L^{2}({\mathbb{R}}^{3})$. The multi-index
${\bm{\lambda}}=(j,i_{1},i_{2},i_{3})$ denotes the scale $2^{-j}$ and position
$L\times 2^{-j}{\mathbf{i}}=L\times 2^{-j}(i_{1},i_{2},i_{3})$ of the wavelets
for each direction, where $i_{\ell}=0,...,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, 2009) is applied to the wavelets. The spatial
average of $\psi_{\mu,{\bm{\lambda}}}({\mathbf{x}})$, defined by,
$\left<\psi_{\mu,{\bm{\lambda}}}\right>=L^{-3}\int_{{\mathbb{T}}^{3}}\psi_{\mu,{\bm{\lambda}}}({\mathbf{x}})d{\mathbf{x}}$
vanishes for each index, which is a necessary condition for being a wavelet.

Similar to Matsuda et al. (2021) the number density field $n({\mathbf{x}},t)$
sampled on $N_{g}^{3}=2^{3J}$ equidistant grid points, can be developed into
an orthogonal wavelet series:

$n({\mathbf{x}})={\overline{n}}_{000}({\mathbf{x}})+\sum_{{j}=0}^{J-1}n_{j}({\mathbf{x}}),$
(3.10)

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

(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$

(a) $j=6$, $St=1$

(b) $j=6$, $St=0.01$

Figure 10: 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$ (e,f) in the $xy-$plane for $St=1$ (left) and $St=0.01$
(right).

$n_{j}({\mathbf{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}{\widetilde{n}}_{{j,}{i_{1},i_{2},i_{3}}}^{\mu}{\psi^{\mu}_{j,i_{1},i_{2},i_{3}}}({\mathbf{x}}),$
(3.11)

with,

${\widetilde{n}}_{{j}{i_{1}i_{2}i_{3}}}^{\mu}=\left<n({\mathbf{x}}),{\psi}_{{j,}{i_{1},i_{2},i_{3}}}^{\mu}({\mathbf{x}})\right>,$
(3.12)

where $\left<,\right>$ denotes an inner product. At scale $2^{-j}$, there are
$7\times 2^{3j}$ wavelet coefficients for $n({\mathbf{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({\mathbf{x}})$ using fast wavelet transform, which has linear
computational complexity. The scale from the wavelet transform and wave
number, $k_{j}$, from the Fourier transform are related as,

$k_{j}=\frac{2\pi}{L}k_{\psi}2^{-j},$ (3.13)

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

Figure 10a,b shows instantaneous, normalized mean number density contours in
the $xy-$plane for two Stokes numbers, $St=1$ and $0.01$. Large number
density, indicative of highly clustered regions, are clearly observed for
$St=1$, but are 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 10c–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 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.

####  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}({\mathbf{x}})$ is defined as,

$M_{q}[n_{j}]=\left<(n_{j})^{q}\right>,$ (3.14)

where the mean values $\left<n_{j}\right>=0$, by giving a central moment and
are related to the $q^{th}$ order structure functions (Schneider et al.,
2004).

The wavelet energy spectrum of $n_{j}({\mathbf{x}})$ can be defined used the
second-order moment, $M_{2}[n_{j}]$, as

$E[n_{j}]={1\over\Delta k_{j}}M_{2}[n_{j}],$ (3.15)

where $\Delta k_{j}=(k_{j+1}-k_{j}){\rm ln}2$ as given by Meneveau (1991). The
energy spectrum obtained from the above equation has the particle number
dependence due to the Poisson noise. Matsuda et al. (2021) 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
$\Omega=\Omega_{f}\cup\Omega_{s}$ is
$M_{2,\mathrm{random},\Omega}[n_{j}]=(7\cdot 2^{3j}/N_{p,\Omega})\langle
n\rangle_{\Omega}^{2}$. For the present case, particles exist only in the
fluid domain $\Omega_{f}$, and $N_{p}=\phi N_{p,\Omega}$ and $\langle
n\rangle=\phi\langle n\rangle_{\Omega}$. The energy of the Poisson noise is
also reduced by a factor of $\phi$, i.e., $M_{2,\mathrm{random}}[n_{j}]=\phi
M_{2,\mathrm{random},\Omega}[n_{j}]$, yielding

$M_{2,{\mathrm{random}}}[n_{j}]=\frac{7\cdot 2^{3j}}{N_{p}}.$ (3.16)

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

${E}[n_{j}]=\frac{1}{\Delta k_{j}}\left\\{M_{2}[n_{j}]-\frac{7\cdot
2^{3j}}{N_{p}}\right\\},$ (3.17)

where the influence of the Poisson noise has been removed. Note that the
analytical estimate of Eq. (3.16) does not contain 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 $\Omega_{f}$,
the combined effect of the Poisson noise and geometrical confinement in
$\Omega_{f}$ can be observed.

The asymmetry of the PDF of $n_{j}({\mathbf{x}})$ is quantified by the
skewness defined as,

$S[n_{j}]={M_{3}[n_{j}]\over(M_{2}[n_{j}])^{3/2}}.$ (3.18)

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

$F[n_{j}]={M_{4}[n_{j}]\over(M_{2}[n_{j}])^{2}},$ (3.19)

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. (2021) 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,\mathrm{random}}[n_{j}]$. They confirmed that the effect of the
Poisson noise on the statistics is negligibly small when 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 about 15 equally spaced snapshots of number
density over this time frame.

(a) $E\left[n_{j}\right]$

(b) $E\left[n_{j}\right]-(E\left[n_{j}\right])_{\rm Random,~{}FCC}$

Figure 11: 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).

Figure 11a 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}/\left<\epsilon\right>)^{1/4}$. Also shown is the spectrum for
random particle positions with uniform probability, where the PDF 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 behavior 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\leq 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 (Matsuda
et al., 2021). 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
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
11b, allows 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., 2014). 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.

(a) Skewness

(b) Flatness

Figure 12: Scale-dependent skewness $S[n_{j}]$ (left) and flatness $F[n_{j}]$
(right) 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}\geq 10$, and dotted lines are used otherwise.

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 12a,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 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 10e–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 to larger Stokes
numbers. Note that for lower Stokes numbers, particle clustering is typically
not significant (figure 10b). Thus, this increase in flatness at smaller
scales may be a result of collision of particles with bead surfaces.

For inertial particles at all Stokes numbers, skewness values also increase
with decreasing scales. As shown by Matsuda et al. (2021), positive skewness
of number density fluctuations indicate high probablity of large positive
values of $n_{j}({\mathbf{x}})$, that is prominent clusters, whereas negative
skewness corresponds to large excursions of negative values of
$n_{j}({\mathbf{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 10.

##  4 Summary and conclusion

Fully resolved direct numerical simulation of a turbulent flow through the
confined geometry of a triply periodic, face-centered cubic (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 presence of three-dimensional 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., 2019). 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 four 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 being taken 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
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 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 analyzing 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 PDFs of the divergence, which are
symmetric, are becoming wider for increasing Stokes number, as the variance
increases. In contrast the tails of the PDFs 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., 2020). 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 number
density distribution were computed. The wavelet energy spectra showed that the
peak of clustering gradually shifts towards 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 fluid domain was also computed. This allowed to identify the
peaks in the spectra for $St\leq 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
number 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 providing thus an explanation for cluster
formation and destruction. A comparison with results for homogeneous isotropic
turbulence, showed many similarities and also pointed out differences due to
the flow geometry, in particular for small volumes. 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 analysis presented can be
extended in the future to vorticity fields to obtain three-dimensional
directional information at different scales.

###### Acknowledgements.

This work was initiated during SVA’s visiting scientist/Professor position at
Aix-Marseille University. Funding for this position was provided by CNRS,
France and Aix-Marseille University. SVA acknowledges kind hospitality during
his visit and stay in Marseille. Funding from NSF award#2053248 is also
gratefully acknowledged. Simulations were performed at the Texas Advanced
Computing Center’s (TACC) Stampede2 and Frontera systems. BK, TO and KS
thankfully acknowledge Centre de Calcul Intensif d’Aix-Marseille for granting
access to its high performance computing resources. TO and KS thankfully
acknowledge financial support from Agence Nationale de la Recherche, project
ANR-20-CE46-0010-01. XH thanks Dr. Timothy 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/).

## References

  * Agbangla et al. (2012) Agbangla, Gbedo Constant, Climent, Éric & Bacchin, Patrice 2012 Experimental investigation of pore clogging by microparticles: Evidence for a critical flux density of particle yielding arches and deposits. Separation and purification technology 101, 42–48.
  * Apte et al. (2009) Apte, S. V., Martin, M. & Patankar, N. A. 2009 A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows. Journal of Computational Physics 228 (8), 2712–2738.
  * Aris (1999) Aris, Rutherford 1999 Elementary chemical reactor analysis. Courier Corporation.
  * Aurenhammer (1991) Aurenhammer, Franz 1991 Voronoi diagrams—a survey of a fundamental geometric data structure. ACM Computing Surveys (CSUR) 23 (3), 345–405.
  * Barber et al. (1996) Barber, C Bradford, Dobkin, David P & Huhdanpaa, Hannu 1996 The quickhull algorithm for convex hulls. ACM Transactions on Mathematical Software (TOMS) 22 (4), 469–483.
  * Bassenne et al. (2017) Bassenne, Maxime, Urzay, Javier, Schneider, Kai & Moin, Parviz 2017 Extraction of coherent clusters and grid adaptation in particle-laden turbulence using wavelet filters. Physical Review Fluids 2 (5), 054301.
  * Bec et al. (2007) Bec, Jeremie, Biferale, Luca, Cencini, Massimo, Lanotte, Alessandra, Musacchio, Stefano & Toschi, Federico 2007 Heavy particle concentration in turbulence at dissipative and inertial scales. Physical Review Letters 98 (8), 084502.
  * Carlson et al. (1992) Carlson, Jon, Gurley, Derrel, King, George, Price-Smith, Colin & Waters, Frank 1992 Sand control: Why and how? Oilfield Review 4 (4), 41–53.
  * Coleman & Vassilicos (2009) Coleman, S.W. & Vassilicos, J.C. 2009 A unified sweep-stick mechanism to explain particle clustering in two-and three-dimensional homogeneous, isotropic turbulence. Physics of Fluids 21 (11), 113301.
  * Cook et al. (2004) Cook, BK, Lee, MY, DiGiovanni, AA, Bronowski, DR, Perkins, ED & Williams, JR 2004 Discrete element modeling applied to laboratory simulation of near-wellbore mechanics. International Journal of Geomechanics 4 (1), 19–27.
  * Dai & Grace (2010) Dai, Jianjun & Grace, John R 2010 Blockage of constrictions by particles in fluid–solid transport. International Journal of Multiphase Flow 36 (1), 78–87.
  * Daubechies (1993) Daubechies, Ingrid 1993 Ten lectures on wavelets. Society of Industrial and Applied Mathematics.
  * Dixon & Nijemeisland (2001) Dixon, A. G. & Nijemeisland, M. 2001 CFD as a design tool for fixed-bed reactors. Industrial & Engineering Chemistry Research 40 (23), 5246–5254.
  * Eaton & Fessler (1994) Eaton, John K & Fessler, JR 1994 Preferential concentration of particles by turbulence. International Journal of Multiphase Flow 20, 169–209.
  * Esmaily-Moghadam & Mani (2016) Esmaily-Moghadam, Mahdi & Mani, Ali 2016 Analysis of the clustering of inertial particles in turbulent flows. Physical Review Fluids 1 (8), 084202.
  * Farge (1992) Farge, M. 1992 Wavelet transforms and their applications to turbulence. Annual Review of Fluid Mechanics 24, 395–457.
  * Farge & Schneider (2015) Farge, M. & Schneider, K. 2015 Wavelet transforms and their applications to MHD and plasma turbulence: a review. Journal of Plasma Physics 81 (6), 435810602.
  * Ferenc & Néda (2007) Ferenc, Járai-Szabó & Néda, Zoltán 2007 On the size distribution of poisson voronoi cells. Physica A: Statistical Mechanics and its Applications 385 (2), 518–526.
  * Finn & Apte (2013) Finn, J. & Apte, S. V. 2013 Relative performance of body fitted and fictitious domain simulations of flow through fixed packed beds of spheres. International Journal of Multiphase Flow 56, 54–71.
  * Goto & Vassilicos (2006) Goto, Susumu & Vassilicos, J.C. 2006 Self-similar clustering of inertial particles and zero-acceleration points in fully developed two-dimensional turbulence. Physics of Fluids 18 (11), 115103.
  * He et al. (2018) He, Xiaoliang, Apte, Sourabh, Schneider, Kai & Kadoch, Benjamin 2018 Angular multiscale statistics of turbulence in a porous bed. Physical Review Fluids 3 (8), 084501.
  * He et al. (2019) He, Xiaoliang, Apte, Sourabh V, Finn, Justin R & Wood, Brian D 2019 Characteristics of turbulence in a face-centred cubic porous unit cell. Journal of Fluid Mechanics 873, 608–645.
  * Hester et al. (2017) Hester, Erich T, Cardenas, M Bayani, Haggerty, Roy & Apte, Sourabh V 2017 The importance and challenge of hyporheic mixing. Water Resources Research 53 (5), 3565–3575.
  * Hill & Koch (2002) Hill, R. J. & Koch, D. L. 2002 The transition from steady to weakly turbulent flow in a close-packed ordered array of spheres. Journal of Fluid Mechanics 465, 59–97.
  * Jin et al. (2015) Jin, Y., Uth, M. F., Kuznetsov, A. V. & Herwig, H. 2015 Numerical investigation of the possibility of macroscopic turbulence in porous media: a direct numerical simulation study. Journal of Fluid Mechanics 766, 76–103.
  * Mahmud et al. (2019) Mahmud, Hisham Ben, Van Hong, Leong & Lestariono, Yuli 2019 Sand production: A smart control framework for risk mitigation. Petroleum 6 (1), 1–13.
  * Mallat (2009) Mallat, Stéphane 2009 A wavelet tour of signal processing, 3rd edn. Boston: Academic Press.
  * Matsuda et al. (2014) Matsuda, Keigo, Onishi, Ryo, Hirahara, Masaaki, Kurose, Ryoichi, Takahashi, Keiko & Komori, Satoru 2014 Influence of microscale turbulent droplet clustering on radar cloud observations. Journal of the Atmospheric Sciences 71 (10), 3569–3582.
  * Matsuda et al. (2021) Matsuda, Keigo, Schneider, Kai & Yoshimatsu, Katsunori 2021 Scale-dependent statistics of inertial particle distribution in high reynolds number turbulence. Physical Review Fluids 6, 064304\.
  * Maxey (1987) Maxey, Martin R 1987 The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields. Journal of Fluid Mechanics 174, 441–465.
  * Meneveau (1991) Meneveau, Charles 1991 Analysis of turbulence in the orthonormal wavelet representation. Journal of Fluid Mechanics 232, 469–520.
  * Monchaux et al. (2012) Monchaux, Romain, Bourgoin, Mickael & Cartellier, Alain 2012 Analyzing preferential concentration and clustering of inertial particles in turbulence. International Journal of Multiphase Flow 40, 1–18.
  * Monin & Yaglom (1965) Monin, A. S. & Yaglom, A. M. 1965 Statistical fluid mechanics: Mechanics of Turbulence. The MIT Press.
  * Orlandi et al. (2018) Orlandi, Paolo, Davide, Modesti & Pirozzoli, Sergio 2018 DNS of turbulent flows in ducts with complex shape. Flow, Turbulence and Combustion 100 (4), 1063–1079.
  * Oujia et al. (2020) Oujia, Thibault, Matsuda, Keigo & Schneider, Kai 2020 Divergence and convergence of inertial particles in high Reynolds number turbulence. Journal of Fluid Mechics 905, A14.
  * Pandya et al. (1998) Pandya, Vishal B, Bhuniya, S & Khilar, Kartic C 1998 Existence of a critical particle concentration in plugging of a packed bed. American Institute of Chemical Engineers. AIChE Journal 44 (4), 978.
  * Patil & Liburdy (2013) Patil, V. A. & Liburdy, J. A. 2013 Flow structures and their contribution to turbulent dispersion in a randomly packed porous bed based on particle image velocimetry measurements. Physics of Fluids 25 (11), 24.
  * Ramachandran & Fogler (1999) Ramachandran, Venkatachalam & Fogler, H Scott 1999 Plugging by hydrodynamic bridging during flow of stable colloidal particles within cylindrical pores. Journal of Fluid Mechanics 385, 129–156.
  * Saucier (1974) Saucier, RJ 1974 Considerations in gravel pack design. Journal of Petroleum Technology 26 (02), 205–212.
  * Schneider et al. (2004) Schneider, Kai, Farge, Marie & Kevlahan, Nicholas 2004 Spatial intermittency in two-dimensional turbulence: a wavelet approach. In Woods Hole Mathematics: Perspectives in Mathematics and Physics, pp. 302–328. World Scientific.
  * Schneider & Vasilyev (2010) Schneider, Kai & Vasilyev, Oleg V 2010 Wavelet methods in computational fluid dynamics. Annual Review of Fluid Mechanics 42, 473–503.
  * Shams et al. (2013) Shams, A, Roelofs, F, Komen, EMJ & Baglietto, E 2013 Quasi-direct numerical simulation of a pebble bed configuration. part i: Flow (velocity) field analysis. Nuclear Engineering and Design 263, 473–489.
  * Sundaram & Collins (1997) Sundaram, Shivshankar & Collins, Lance R 1997 Collision statistics in an isotropic particle-laden turbulent suspension. part 1. direct numerical simulations. Journal of Fluid Mechanics 335, 75–109.
  * Toschi & Bodenschatz (2009) Toschi, Federico & Bodenschatz, Eberhard 2009 Lagrangian properties of particles in turbulence. Annual Review of Fluid Mechanics 41, 375–404.
  * Valdes & Santamarina (2006) Valdes, Julio R & Santamarina, J Carlos 2006 Particle clogging in radial flow: Microscale mechanisms. SPE Journal 11 (02), 193–198.
  * Wang et al. (2000) Wang, Lian-Ping, Wexler, Anthony S & Zhou, Yong 2000 Statistical mechanical description and modelling of turbulent collision of inertial particles. Journal of Fluid Mechanics 415, 117–153.
  * Wood et al. (2020) Wood, Brian D., He, Xiaoliang & Apte, Sourabh V. 2020 Modeling turbulent flows in porous media. Annual Review of Fluid Mechanics 52 (1), 171–203, arXiv: https://doi.org/10.1146/annurev-fluid-010719-060317.
  * Yi et al. (2005) Yi, Xi, Valkó, PP & Russell, J.E. 2005 Effect of rock strength criterion on the predicted onset of sand production. International Journal of Geomechanics 5 (1), 66–73.