1. Introduction
Sedimenting multiphase flows occur in a wide range of society-critical applications, such as the upgrading of biomass into biofuels in circulating fluidised bed reactors (Cui & Grace Reference Cui and Grace2007; Mettler, Vlachos & Dauenhauer Reference Mettler, Vlachos and Dauenhauer2012), transport of contaminants within ground water (Miller et al. Reference Miller, Christakos, Imhoff, McBride, Pedit and Trangenstein1998) or the atmosphere (Ravishankara Reference Ravishankara1997; Arganat & Perminov Reference Arganat and Perminov2020) and several other environmentally relevant systems, such as volcanic processes (Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020) and avalanches (Sovilla, McElwaine & Köhler Reference Sovilla, McElwaine and Köhler2018; Cicoira et al. Reference Cicoira, Blatny, Li, Trottet and Gaume2022; Gardezi et al. Reference Gardezi, Xing, Bilal, Zhuang, Muhammad and Janjua2022; Zhuang et al. Reference Zhuang, Bilal, Xing, Li, He and Zhang2023). All of these inherently multiphase systems represent problems of great societal interest that are also historically challenging to study, both computationally and experimentally. It is these challenges that, to date, have hindered our ability to fully understand and accurately predict their behaviour, thus impeding more optimal engineering designs and systems, as well as more effective procedures and standards for safety and risk mitigation.
As is often the case in multiphase systems of interest, large-scale heterogeneity arises in the form of turbulence in the carrier phase and cohesive structures such as clusters in the dispersed phase. Clustering, in particular, has been extensively examined both experimentally (Breard & Lube Reference Breard and Lube2017; Breard, Dufek & Lube Reference Breard, Dufek and Lube2018; Sovilla et al. Reference Sovilla, McElwaine and Köhler2018) and computationally (Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2014; Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020), and has been shown to arise through mechanisms that are related to the Stokes number (the ratio of particle to fluid time scales) and the volume fraction,
$\langle \alpha _p \rangle$
, or mass loading,
$\varphi = (\rho _p \alpha _p)/(\rho _f \alpha _f)$
. Four traditionally recognised regimes emerge depending upon the values of the volume fraction and the Stokes number. Namely, when volume fraction is sufficiently low (
$\lt \mathcal {O}(10^{-6})$
) particles are considered to be ‘one-way coupled’ and particles have a negligible effect on the fluid, but the fluid dictate the motion of particles. In this regime, clustering can emerge through the mechanism of preferential concentration, as originally described by Maxey (Reference Maxey1987a
, Reference Maxey
b
). It is of course worth noting that, in order to form cohesive structures through the mechanism of preferential concentration, there must be a background turbulence or the existence of regions with increased shear. As the volume fraction is increased, there is an increasing interplay between particles and the carrier fluid. For moderate volume fractions (
$ \langle \alpha _p \rangle \in \lbrack \mathcal {O}(10^{-6}), \mathcal {O}(10^{-3})\rbrack$
), suspensions are considered to be ‘two-way coupled’. For this regime, there exist two sub-regimes based upon the Stokes number. For particles with Stokes number less than unity, particles tend to contribute to turbulent energy dissipation. For particles with Stokes number greater than unity, particles have an opposite effect and tend to generate turbulent kinetic energy. Finally, when the volume fraction surpasses
$\mathcal {O}(10^{-3})$
, the system is considered to be strongly coupled (Elghobashi & Truesdell Reference Elghobashi and Truesdell1992; Bosse, Kleiser & Meiburg Reference Bosse, Kleiser and Meiburg2006; Breard & Lube Reference Breard and Lube2017; Breard et al. Reference Breard, Dufek and Lube2018). Here, not only do both phases effect on one another’s behaviour, but collisions between particles also becomes important. In this regime, there is less of an effect on flow behaviour due to the Stokes number and clustering primarily arises due to the reduced drag felt by particles (Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2015).
As one may expect, the emergence of heterogeneity has an important effect on nearly all quantities of interest for a multiphase system, such as settling velocity and degree of mixing in the flow. Several recent works have highlighted this effect and shown that heterogeneity can impede chemical conversion (Beetham & Capecelatro Reference Beetham and Capecelatro2019), alter thermal exchange efficiency (Guo & Capecelatro Reference Guo and Capecelatro2019) and entrance lengths (Beetham & Capecelatro Reference Beetham and Capecelatro2023) and enhance the mean particle settling velocity (Wang & Maxey Reference Wang and Maxey1993; Aliseda et al. Reference Aliseda, Cartellier, Hainaux and Lasheras2002; Ferrante & Elghobashi Reference Ferrante and Elghobashi2003; Yang & Shy Reference Yang and Shy2003 Bosse et al. Reference Bosse, Kleiser and Meiburg2006).
In addition to the complexity generated by heterogeneity in sedimenting flows and the coupling between the phases, another key challenge is the breadth of length and time scales at which physical processes occur. At the particle scale, which is often of the order of micrometres, small wakes and boundary layer-induced instabilities are formed at the surface of particles and in systems that exhibit heat and mass transfer or are chemically reactive, these physical processes take place at very small time scales. These multi-physics processes drive the need for high computational resolution both in space and time in order to fully resolve relevant physics. Computational studies that aim to capture this level of detail often require grids fine enough that there are tens of grid points across the particle diameter (Subramaniam Reference Subramaniam2020). Juxtaposed with this, nearly all systems of societal interest span length scales on orders relevant to a reactor or volcanic topography (metres or kilometres) and time scales of the order of minutes or hours. In light of this, it is obvious that, even with modern computing capabilities, fully resolved computations of full-scale systems of interest are intractable. Because of this, several computational strategies exist for simulating multiphase flows at a range of scales. Particle resolved direct numerical simulation (PR-DNS), while model free, is useful for studying behaviour at the microscale (i.e. systems of the order of hundreds of particles) and for developing the models required by coarser-grained methodologies, such as highly resolved Euler–Lagrange simulations, where the fluid is treated in the Eulerian sense, particles are individually tracked as Lagrangian bodies and the phases are coupled through interphase exchange terms (Capecelatro et al. Reference Capecelatro, Desjardins and Fox2015). In these coarser-grained simulations, the grid is of the order of 2–3 times larger than the particles, and models (informed by PR-DNS) are required to capture sub-grid behaviour. Importantly, these methodologies allow simulations at the mesoscale with hundreds of millions of particles (Capecelatro et al. Reference Capecelatro, Desjardins and Fox2014; Beetham, Fox & Capecelatro Reference Beetham, Fox and Capecelatro2021). Systems at this size are sufficiently large to observe the mesoscale heterogeneity that microscale PR-DNS simulations cannot. While this is a dramatic improvement in the scale of systems that can be considered, mesoscale systems are also computationally limited, making them incapable of tractably considering systems at scales of interest (i.e. full reactor scale or the several kilometres required to capture volcano or avalanche runout). Thus, in the same way that PR-DNS data are useful for informing model closures for mesoscale (EulerEuler or Euler–Lagrange) simulations, the resulting data can then in turn be used to formulate closures for even coarser-grained strategies, such as the multiphase Reynolds-averaged Navier–Stokes (multiphase RANS) equations (Beetham et al. Reference Beetham, Fox and Capecelatro2021).
It is worth emphasising that this incremental and careful approach – moving strategically from PR-DNS to highly resolved Euler–Lagrange to even coarser-grained multiphase RANS simulations – is a necessity for retaining relevant mesoscale physics. While early work in multiphase flow modelling attempted to augment existing single-phase turbulence models (Sinclair & Jackson Reference Sinclair and Jackson1989; Sundaram & Collins Reference Sundaram and Collins1994; Cao & Ahmadi Reference Cao and Ahmadi1995; Dasgupta, Jackson & Sundaresan Reference Dasgupta, Jackson and Sundaresan1998; Zeng & Zhou Reference Zeng and Zhou2006; Rao et al. Reference Rao, Curtis, Hancock and Wassgren2012), this approach fails beyond very dilute suspensions of particles. This primarily owes to the fact that in multiphase systems, turbulent energy is generated at the scale of the particles and cascades up to generate large-scale turbulence. This phenomenon is in direct opposition to classical turbulence theory (Kolmogorov Reference Kolmogorov1941) and underscores the shortcomings of extending single-phase turbulence models to multiphase flows. To this end, work by Fox (Reference Fox2014) pointed out that deriving the multiphase RANS equations directly from the microscale equations (Navier–Stokes) is incapable of capturing the momentum exchange between the phases. Instead, it is necessary to formulate the multiphase RANS equations by averaging the mesoscale equations. This strategy, while exact, leads to several unclosed terms that require modelling. While recent work (Beetham et al. Reference Beetham, Fox and Capecelatro2021) has proposed closures for these terms, these closures are limited to unbounded flows of spherical, monodisperse particles of a single density ratio and the development of accurate, widely applicable closures is still very much in its infancy.
While a large body of work in the last decades (Balachandar & Eaton Reference Balachandar and Eaton2010; Subramaniam Reference Subramaniam2020) has shed light on the physics underlying multiphase processes as well as formulations for guiding reduced-order modelling of their behaviour, most contemporary studies consider particles that are monodisperse. In contrast to this, all real-world flows of interest are comprised of polydispersed particles, where particle parameters often span a wide range of diameters, shapes, densities, etc. We note that some recent studies have considered particulate phases that were not monodisperse (see, e.g., Fox Reference Fox2024), however, several of these studies often were limited to bi-disperse assemblies of particles (Municchi & Radl Reference Municchi and Radl2017), quantities of particles that were too few to capture mesoscale behaviour or were limited to particles too small to be considered strongly coupled (Islam et al. Reference Islam, Saha, Gemci, Yang, Sauret, Ristovski and Gu2019).
In addition to computational work, recent experimental studies have also studied the settling behaviour of gas–solid flows that exhibit mesoscale behaviour, such as clustering. Penlou et al. (Reference Penlou, Roche, Manga and van der Wildenberg2023) considered how the Stokes number impacted settling. In the cases of similarly sized small particles (78
$\unicode{x03BC} \textrm{m}$
), particles were observed to experience enhanced settling behaviour due to cluster formation, whereas larger particles (467
$\unicode{x03BC}\textrm{m}$
) demonstrated hindered settling (Penlou et al. Reference Penlou, Roche, Manga and van der Wildenberg2023). Other work, which also considered similarly sized suspensions of particles over a range of particle diameters observed that the maximum particle volume fraction due to clustering behaviour scaled with the Reynolds number as Re
$^{1/5}$
or Re
$^{1/2}$
depending on the region of the flow considered (Weit et al. Reference Weit, Roche, Dubois and Manga2018).
In recent years, several one-dimensional (1-D) (Degruyter & Bonadonna Reference Degruyter and Bonadonna2012; Dufek, Esposti Ongaro & Roche Reference Dufek, Esposti Ongaro and Roche2015; Folch, Costa & Macedonio Reference Folch, Costa and Macedonio2016; Pouget et al. Reference Pouget, Bursik, Singla and Singh2016), 2-D (Pfeiffer, Costa & Macedonio Reference Pfeiffer, Costa and Macedonio2005; Barsotti, Neri & Scire Reference Barsotti, Neri and Scire2008; Williams, Stinton & Sheridan Reference Williams, Stinton and Sheridan2008) and 3-D multiphase flow models (Neri et al. Reference Neri, Esposti Ongaro, Menconi, De’Michieli Vitturi, Cavazzoni, Erbacci and Baxter2007; Dufek & Bergantz Reference Dufek and Bergantz2007b ; Costa Reference Costa2016; Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020; Cao et al. Reference Cao, Bursik, Yang and Patra2021) have been used to quantify key physical processes in polydispersed, multiphase flows (such as pyroclastic density currents, or PDCs) that are often inaccessible by large-scale experimentation. In addition to quantifying the physics of such flows, large-scale simulations such as these have also been used to develop reduced-order forecasting models. That said, while the current state-of-the-art models mentioned previously can reveal some aspects of the internal structure of PDCs, to date they have fallen short of capturing the heterogeneous effects of clustering. This represents a critical missing element for accurate models because in order for a reduced-order model of such flow configurations to be successful, they must be representative of polydisperse heterogeneity and settling.
Due to the lack of more comprehensive models, the Stokes settling velocity,
$\mathcal {V}_0 = \tau _p g$
, with
$\tau _p$
, the particle response time defined as
$\rho _p {\textrm{d}}_p^2/(18\unicode{x03BC} _f)$
, has been traditionally used in many works to predict settling behaviour in multiphase flows (Dietrich Reference Dietrich1982; Basson, Berres & Bürger Reference Basson, Berres and Bürger2009; Shankar, Pandey & Shukla Reference Shankar, Pandey and Shukla2021). Here,
$g$
represents the magnitude of gravity,
$d_p$
and
$\rho_p$
are the particle diameter and density, respectively and
$\mu_f$
is the fluid dynamic viscosity. As an improvement to a Stokes prediction, which makes gross oversimplifications, contemporary work has focused on extending depth-averaged equations that were originally used in the context of shallow-water systems to gas–solid flows (Esposti Ongaro et al. Reference Esposti Ongaro, Cerminara, Charbonnier, Lube and Valentine2020; Breard et al. Reference Breard, Dufek, Charbonnier, Gueugneau, Giachetti and Walsh2023; de’Michieli Vitturi et al. Reference de’Michieli Vitturi, Esposti Ongaro and Engwell2023). To this end, de’Michieli Vitturi et al. (Reference de’Michieli Vitturi, Esposti Ongaro and Engwell2023) generalised the depth-averaged equations to include a dispersed phase and employed an Eulerian–Eulerian approach to simulate a wide range of volcanic phenomena. As a result of this work, a depth-averaged settling velocity, denoted as
$\mathcal {V}_s$
, with a model for the drag coefficient originally developed using kinetic theory arguments (Gidaspow Reference Gidaspow1994) was proposed.
In this work, we present an analysis of the clustering and settling behaviour observed in high-fidelity Euler–Lagrange simulations for two polydisperse assemblies of particles, each at two volume fractions. In addition, comparisons are drawn against the behaviour observed in analogous monodisperse configurations in an effort to isolate the effect of polydispersity on heterogeneous multiphase behaviour. This work represents, to the authors’ knowledge, the most highly resolved study of polydisperse clustering and settling behaviour to date. Further, we propose a new model for the mean settling velocity of strongly coupled gas–solid flows based on particle size, volume fraction and degree of polydispersity. This model represents an initial step toward improved reduced-order models.
2. System description
While polydisperse settling has broad applications, as previously discussed, this work is motivated by (PDCs – the gravity current resulting from the collapse of a volcanic column. The interested reader is referred to Appendix A for a brief summary of this particular class of flows. Given this motivation, we leverage historical ejecta information to inform the parameters for the highly resolved simulations under study.
In particular, we aim to quantify the effect of polydispersity on clustering and settling behaviour. To this end, we consider four configurations containing a polydisperse assembly of particles with parameters sampled from historical PDC ejecta data. In each of these configurations, the particles are assumed to be rigid spheres with diameter,
${\textrm{d}}_p$
, and density,
$\rho _{p}$
. The gas has a constant density,
$\rho _{f}$
, and dynamic viscosity,
$\nu _{f}$
, which are specified based on the properties of air at 400
$^{\circ}$
C, the average temperature in the intermediate region (Breard et al. Reference Breard, Lube, Jones, Dufek, Cronin, Valentine and Moebis2016; Breard & Lube Reference Breard and Lube2017; Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020).
Particle density is constant across all configurations and is consistent with the density typically found in the intermediate layer of PDCs (Lube et al. Reference Lube, Breard, Cronin and Jones2015; Breard et al. Reference Breard, Lube, Jones, Dufek, Cronin, Valentine and Moebis2016). Particle diameter distributions are chosen following the data presented in (Lube et al. Reference Lube, Breard, Cronin, Procter, Brenna, Moebis, Pardo, Stewart, Jolly and Fournier2014). In this study, we consider two log-normal distributions corresponding to
$(\phi , \phi _{{sorting}}) = (0, 1)$
and
$(1.5, 1)$
(for a description of
$\phi$
and
$\phi _{{sorting}}$
, refer to Appendix B), denoted as
$A$
and
$B$
throughout the manuscript. While the details of both distributions are summarised in table 1, it is important to note that the particles for each configuration are selected from the described distribution between the values of
$d_{{min}}$
and
$d_{{max}}$
. This results in the mean and standard deviation of the particle assemblies in the simulations that deviate from the mean and standard deviation for a log-normal distribution described by
$\unicode{x03BC}$
and
$\sigma$
with infinite support. Here,
$\mu$
and
$\sigma$
are the distribution parameters associated with log-normal data having the probability distribution function
$\frac{1}{x \sigma \sqrt{2 \pi}}\exp{\left(- \frac{(\ln{x}-\mu)^2}{2\sigma^2} \right)} $
, for some random variable,
$x$
, which in this case is the particle diameter.
Summary of simulation parameters under consideration. The PDFs shown above the table represent the distribution from which particles were sampled and are specified according to log-normal parameters
$\mu$
and
$\sigma$
. The dashed vertical lines represent the diameter for the corresponding monodisperse simulations.
For each configuration
$A$
and
$B$
, we draw comparisons with a monodisperse analogue, denoted as
$A_0$
and
$B_0$
, where the diameter of the particles in each of these is equivalent to the mean diameter of distributions
$A$
and
$B$
, respectively. Here, these diameters are specified by the expected diameter of the distributions
$A$
and
$B$
with full support.
In addition to studying the effect of two different degrees of polydispersity and two analogous monodisperse configurations, we also consider each configuration at two volume fractions. For consistency with the conditions typical of the intermediate region of PDCs (Dufek & Bergantz Reference Dufek and Bergantz2007a
; Lube et al. Reference Lube, Breard, Cronin and Jones2015, Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020; Breard et al. Reference Breard, Lube, Jones, Dufek, Cronin, Valentine and Moebis2016; Breard & Lube Reference Breard and Lube2017; Breard et al. Reference Breard, Dufek and Lube2018), we consider two global particle-phase volume fractions:
$\langle \alpha _p \rangle = 0.01$
and
$0.1$
. Here, angled brackets denote a domain average and the particle volume fraction,
$\alpha _p$
, is defined as the ratio of the volume occupied by particles to the volume of the total domain.
To isolate the effect of polydispersity on clustering and settling behaviour, we consider the eight assemblies of particles described above in an unbounded domain and subjected to gravity. Due to the relative simplicity of such a configuration, these flows can be characterised by a small set relevant non-dimensional quantities: the particle Reynolds number, Re
$_{\! p}$
, the Froude number, Fr, the Stokes number, St, and the Archimedes number. The particle Reynolds number, Re
$_{\! p}$
, is defined as
$\mathcal {V}_0 {d}_p/\nu _f$
, where
$\mathcal {V}_0$
is the Stokes settling velocity for a given particle diameter, defined as
$\mathcal {V}_0 = \tau _p g$
, where
$\tau _p = \rho _p {d}_p ^2/(18 \unicode{x03BC} _f)$
is the Stokes settling time. The Froude number, which quantifies the balance between gravitational and inertial forces, is defined as Fr
$=\tau _p^2 g/{d}_p$
. Finally, the Archimedes number is defined as Ar
$= g {d}_p^3 \rho _f(\rho _p-\rho _f)/(\unicode{x03BC} _f^2)$
.
The particulates commonly observed in PDCs have sizes in the range
$\pm 5 \phi$
units (corresponding to a range from 0.032 to 32 mm) in diameter and densities spanning 900 to 3300 kg m−
$^3$
(Taddeucci & Palladino Reference Taddeucci and Palladino2002). Temperature of the fluidising gas is considered to be dry air at approximately 400
$^{\text {o}}$
C. These ranges of parameters result in typical settling velocities of the order of
$\mathcal {O}(10^{-2})$
to
$\mathcal {O}(10^{4})$
and Archimedes numbers of the order of
$\mathcal {O}(10^{-1})$
to
$\mathcal {O}(10^{8})$
. To this end, we have adjusted the value of gravity and the upper and lower end cutoffs for diameter such that our simulations fall within this range while also reducing
$\mathcal {L}$
and maintaining a grid spacing such that a less extensive domain size is both statistically relevant and computationally tractable.
It is important to make note here of the computational challenge of highly resolved simulations of polydisperse gas–solid particles on a uniform grid – a necessity considering that an adaptive mesh refinement strategy based on particles in a two-phase flows does not yet exist. Thus, for a uniform grid, the grid spacing must be sufficiently fine such that the smallest particles are resolved (
$\delta x \approx 1.75-3.0 {d}_{p,min}$
) and the extent of the domain must be large enough that large-scale heterogeneity is also resolved. Prior work (Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2016) has shown that for smaller, monodisperse particles, an overall domain size should scale with the characteristic cluster length,
$\mathcal {L} = \tau _p^2 g$
. However, it is not clear that this scaling extends to larger or polydisperse particle suspensions and this study is reserved for future work.
A list of these non-dimensional numbers and other simulation parameters are summarised in table 1. Note that quantities appended with a subscript of ‘10’ corresponds to that quantity evaluated with a particle diameter equal to
$D_{32}$
, the Sauter mean particle diameter. A complete description of the definitions of statistical diameters can be found in Appendix C.
In this work, we utilise an Euler–Lagrange framework to capture high fidelity behaviour of each of the particle configurations under study. The mathematical formulation for this framework is summarised in the subsequent sections.
2.1. Eulerian fluid-phase description
To account for the presence of particles in the fluid phase without resolving the boundary layers around individual particles, we consider the volume-filtered, incompressible Navier–Stokes equations (Anderson & Jackson Reference Anderson and Jackson1967). This procedure replaces the point microscale variables (fluid velocity, pressure, etc.) with smooth, locally filtered fields and requires models to be used for sub-grid-scale behaviour, such as particle drag (see, e.g., Xian et al. Reference Xian, Lin, Guo and Yu2024) and particle–particle collision forces. Vectors and tensors in equations are represented in bold text. The volume-filtered conservation of mass and momentum equations are given by
\begin{equation} \frac {\partial \left ( \alpha _f \rho _f \right )}{\partial t} + \boldsymbol{\nabla} \cdot (\alpha _f \rho _f \boldsymbol {u}_f) = 0, \end{equation}
and
\begin{equation} \frac {\partial \left ( \alpha _f \rho _f \boldsymbol {u}_f \right )}{\partial t}+ \boldsymbol{\nabla} \cdot (\alpha _f \rho _f \boldsymbol {u}_f \otimes \boldsymbol {u}_f)= \boldsymbol{\nabla} \cdot \boldsymbol {\tau }_f + \alpha _f \rho _f \boldsymbol {g} + \boldsymbol {\mathcal {F}} + \boldsymbol {F}_{{mfr}}, \end{equation}
where
$\boldsymbol{\mathcal {F}}$
accounts for the two-way coupling between the fluid and particle phases and is defined in § 2.3. Here,
$\boldsymbol {F}_{{mfr}}$
is a source term imposed to ensure the system maintains a net constant mass flow rate in order to achieve a statistically stationary state. The fluid-phase viscous stress tensor is defined as
\begin{equation} \boldsymbol {\tau }_f= -p_f \mathbb {I} + \unicode{x03BC} _f [\boldsymbol{\nabla} \boldsymbol {u}_f + (\boldsymbol{\nabla} \boldsymbol {u}_f)^T - \frac {2}{3} \boldsymbol{\nabla} \cdot \boldsymbol {u}_f \mathbb {I}], \end{equation}
where
$\mathbb {I}$
is the identity matrix.
2.2. Lagrangian particle-phase equations
The dispersed phase is solved using Lagrangian mechanics. The position and velocity of the particles are advanced according to Newton’s second law
\begin{align} \frac {\text {d}\boldsymbol {x}_p^{(i)}}{\text {d}t}&=\boldsymbol {u}_p^{(i)} \end{align}
and
\begin{align} m_p^{(i)}\frac {\text {d}\boldsymbol {u}_p^{(i)}}{\text {d}t} &= \boldsymbol {F}^{(i)}_{{inter}}+\boldsymbol {F}^{(i)}_{{col}} + m_p^{(i)} \boldsymbol {g} , \end{align}
where the superscript
$(i)$
denotes the
$i$
th particle. Throughout the rest of this section, square brackets indicate a fluid quantity evaluated at the centre of the
$i$
th particle’s location.
As shown in (2.5), particles are subject to three forces: the force due to interphase momentum exchange,
$\boldsymbol {F}_{{inter}}$
, the force due to inter-particle collisions,
$\boldsymbol {F}_{{col}}$
, and the body force. Here, the interphase exchange term is given by,
\begin{equation} \boldsymbol {F}_{{inter}}^{(i)}= V_p^{(i)} \boldsymbol{\nabla} \cdot \boldsymbol {\tau }_f\lbrack \boldsymbol {x}_p^{(i)}\rbrack + \frac {m_p^{(i)} \alpha _f\lbrack \boldsymbol {x}_p^{(i)}\rbrack F_d\left (\alpha _f, \textit{Re}^{(i)}_p\right )}{\tau _p^{(i)}} (\boldsymbol {u}_f\lbrack \boldsymbol {x}_p^{(i)}\rbrack - \boldsymbol {u}_p^{(i)}), \end{equation}
where the rightmost term is the drag force felt by the
$i$
th particle. Here,
${F}_d$
is the non-dimensional drag correction of Tenneti, Garg & Subramaniam (Reference Tenneti, Garg and Subramaniam2011) that takes into account local volume fraction and Reynolds number effects and is given by
\begin{equation} F_d(\alpha _f , Re_p) =\frac {1 +0.15Re_p^{0.687}}{\alpha _f^2} + \alpha _fF_1(\alpha _f) + \alpha _fF_2(\alpha _f , Re_p). \end{equation}
Here, the local particle Reynolds number is defined as
\begin{equation} \textit{Re}_p^{(i)} = \frac {\alpha _f\lbrack \boldsymbol {x}_p^{(i)}\rbrack \left \vert \boldsymbol {u}_f\lbrack \boldsymbol {x}^{(i)} \rbrack - \boldsymbol {v}_p^{(i)}\right \vert d_p}{\nu _f} ,\end{equation}
and the remaining two terms in (2.7) are given by
\begin{equation} F_1(\alpha _f)= \frac {5.81 \alpha _p}{\alpha _f^3} + \frac {0.48 \alpha _p^{1/3}}{\alpha _f^4}, \end{equation}
and
\begin{equation} F_2(\alpha _f, \textit{Re}_p)= \alpha _p^{3} \textit{Re}_p \left (0.95 + \frac {0.61 \alpha _p^3}{\alpha _f^2}\right ). \end{equation}
The force of collisions is accounted for using a soft-sphere collision model (Cundall & Strack Reference Cundall and Strack1979) and particles are treated as inelastic and frictional with a coefficient of restitution of 0.85 and coefficient of friction of 0.1.
2.3. Two-way coupling
The fluid-phase equations introduced in § 2.1 contains two interphase exchange terms: the particle volume fraction,
$\alpha _p$
, and the momentum exchange term,
$\boldsymbol {\mathcal {F}}$
. Each of these terms requires projection of the Lagrangian particle information to the Eulerian mesh. We make use of the two-step filtering approach proposed by Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013). In this approach, the particle-localised data are extrapolated to the nearest grid points and is implicitly smoothed using a Gaussian filter kernel (denoted
$\mathcal {G}$
in (2.11) and (2.12)) with a width equal to 7 times the mean particle diameter.
Given this, the particle-phase volume fraction is defined as
\begin{equation} \alpha _f = 1- \sum _{i=1}^{N_p} \mathcal {G}\left (\vert \boldsymbol {x} - \boldsymbol {x}^{(i)}_p \vert \right ) V_p^{(i)} , \end{equation}
where
$N_p$
is the total number of particles and
$V^{(i)}_p$
is the volume of the
$i$
th particle.
Similarly, the interphase momentum exchange is defined as
\begin{equation} \boldsymbol {\mathcal {F}} = - \sum _{i=1}^{N_p}\mathcal {G}\left (\vert \boldsymbol {x} - \boldsymbol {x}^{(i)}_p \vert \right ) \boldsymbol {F}_{{inter}}. \end{equation}
2.4. Computational methodology
The equations presented in §§ 2.1–2.3 are evolved using an in-house code, NGA (Desjardins et al. Reference Desjardins, Blanquart, Balarac and Pitsch2008), a fully conservative, low-Mach-number finite volume solver that has been used extensively for strongly coupled gassolid flows and has been shown to produce physically relevant results (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013; Capecelatro et al. Reference Capecelatro, Desjardins and Fox2014, Reference Capecelatro, Desjardins and Fox2015). A pressure Poisson equation is solved to enforce continuity via fast Fourier transforms in all three periodic directions. The fluid equations are solved on a staggered grid with second-order spatial accuracy and advanced in time with second-order accuracy using the semi-implicit Crank–Nicolson scheme of Pierce & Moin (Reference Pierce and Moin2001). Lagrangian particles are integrated using a second-order Runge–Kutta method. Fluid quantities appearing in § 2.2 are evaluated at the position of each particle via trilinear interpolation and particle data are projected onto the Eulerian mesh using the two-step filtering process described in Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013).
Gravity is oriented in the negative
$x$
direction and is specified such that flow statistics are properly resolved. The domain size and grid spacing is specified such that clustering statistics are properly resolved. Namely, the grid spacing is set to be no greater than 1.75
${d}_{p,{min}}$
and the domain length in the gravity-aligned direction is specified to be several times larger than the expected length of clusters. These domain parameters are summarised in table 1.
It is important to note that, while prior work (Capecelatro et al. Reference Capecelatro, Desjardins and Fox2016) has suggested the domain size for gravity-driven gas–solid flows should be at least 32
$\tau _p^2 g$
long in the streamwise direction in order to properly resolve flow statistics, this guidance originated from monodisperse particles of a smaller diameter (90
$\unicode{x03BC}\textrm{m}$
) than the mean particles considered in this work. Thus, we have considered a domain length large enough to contain a sufficient number of particles to observe clustering behaviour.
The domain for all the configurations under study are triply periodic. Particles are initially randomly distributed in the domain and the fluid phase is initially quiescent. After a transient period, the flow reaches a statistically stationary state, which is assessed by monitoring the mean particle settling velocity and the mean variance in particle volume fraction.
3. Phased-averaged quantities of interest
As described in § 1, coarse-grained models often make use of averaging. Thus, phase-averaged terms of high-fidelity data are useful for constructing improved models. In this section, we present several flow equations and quantities to aid in the quantification of the configurations under study and lay the foundation for improved models. This section serves as a brief summary, however, more details regarding the quantities presented here can be found in Capecelatro et al. (Reference Capecelatro, Desjardins and Fox2015) and Beetham et al. (Reference Beetham, Fox and Capecelatro2021).
For multiphase flows, it is convenient to introduce a phase average after taking the Reynolds average of the volume-filtered questions. These phase averages are analogous to Favré averages for variable-density flows and are denoted by
$\langle (\cdot ) \rangle _p = \langle \alpha _p (\cdot ) \rangle / \langle \alpha _p \rangle$
and
$\langle (\cdot ) \rangle _f = \langle \alpha _f (\cdot ) \rangle / \langle \alpha _f \rangle$
. Phase averages are denoted with angled brackets and a subscript
$f$
or
$p$
to indicate the phase with respect to which the average has been taken. Angled brackets without a subscript indicate a standard Reynolds average, which consists of an average overall spatial dimensions in this work. Fluctuations about particle phase-averaged quantities are denoted with two primes, e.g.
$\boldsymbol {u}_p^{\prime \prime }=\boldsymbol {u}_p- \langle \boldsymbol {u}_p \rangle _p$
, with
$\langle \boldsymbol {u}_p^{\prime \prime } \rangle _p=0$
. Fluctuations from fluid phase-averaged quantities are denoted with three primes, e.g.
$\boldsymbol {u}_f^{\prime \prime \prime }=\boldsymbol {u}_f- \langle \boldsymbol {u}_f \rangle _f$
.
For fully developed gravity-driven flow, phase-averaged variables are statistically stationary. As a consequence of this, the phase-averaged continuity equation implies
$\alpha _f$
is constant. Further, the fluid-phase momentum equation reduces to
$\langle u_f \rangle _f=0$
. Taking the phase average of (2.5) projected to the Eulerian mesh results in the particle phase-averaged momentum equation. The
$x$
-direction component of this expression, the only non-zero component, is given as
\begin{equation} \frac {\partial \langle u_p \rangle _p}{\partial t} = \frac {1}{\tau _p^*}(\langle {u_f}\rangle _{p}-\langle {u_p}\rangle _{p}) + \frac {1}{\rho _p} \left (\left \langle \frac {\partial \sigma _{f,xi}}{\partial x_i}\right \rangle _p - \left \langle \frac {\partial p_f}{\partial x} \right \rangle _p \right ) +g. \end{equation}
Here, we note that the transport of the phase-averaged particle velocity results in a balance between the forces exchanged by the fluid (i.e. drag model and resolved surface stresses) and gravity. Prior work (Capecelatro et al. Reference Capecelatro, Desjardins and Fox2015; Beetham et al. Reference Beetham, Fox and Capecelatro2021) has shown that the fluid stress and fluid pressure gradient terms are small enough for gas–solid flows to be reasonably neglected. This implies that at steady state,
$\langle u_p \rangle _p \approx \langle u_f \rangle _p + \tau _p^{\ast} g$
. Here,
$\langle u_f\rangle _p$
is often thought of as the fluid velocity seen by the particles. It is also notable that we incorporate the nonlinearities associated with the drag in
$\tau _p^*=\langle \tau _p\rangle / \langle F_d \rangle _p$
, where
$\langle F_d \rangle _p(\langle \alpha _f \rangle , \langle Re_p \rangle )$
and
$\langle \tau _p \rangle (\langle {d}_p \rangle )$
are the nonlinear drag correction and particle response time with averaged flow quantities used as argument to (2.7) and the definition of
$\tau _p$
. In the following sections, we will make use of this relationship to propose a model for settling velocity that follows the formulation that the mean settling velocity is equal to the Stokes settling velocity plus an unclosed term.
Phase averaging also gives rise to two Reynolds stress tensors, one in each of the particle and fluid phases,
$\langle \boldsymbol {u}_p^{\prime \prime } \boldsymbol {u}_p^{\prime \prime } \rangle _p$
and
$\langle \boldsymbol {u}_f^{\prime \prime \prime } \boldsymbol {u}_f^{\prime \prime \prime } \rangle _f$
, respectively. Analogous to single-phase turbulence, the trace of each of these tensors yields the turbulent kinetic energy in the particle and fluid phases, respectively. These are denoted as
$k_p$
and
$k_f$
and defined as
\begin{align} k_p &= \frac {1}{2}\langle \boldsymbol {u}_p^{\prime \prime } \cdot \boldsymbol {u}_p^{\prime \prime } \rangle _p \end{align}
and
\begin{align} k_f &= \frac {1}{2}\langle \boldsymbol {u}_f^{\prime \prime \prime } \cdot \boldsymbol {u}_f^{\prime \prime \prime } \rangle _f. \end{align}
It is notable that in the case of the particle-phase averages, Lagrangian particle quantities, such as
$\boldsymbol {u}_p$
, must be evaluated on the Eulerian mesh. This is done by trilinear interpolation and the application of a Gaussian filter kernel, as described in § 2.3 and Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013).
In addition to the phase-wise Reynolds stresses, it is also important to note that the Reynolds average of the inner product of the particle velocity fluctuations gives way to the total energy fluctuations in the particle phase, or the total granular energy,
$\kappa _p$
, given by
\begin{equation} \kappa _p = \frac {1}{2}\langle \boldsymbol {u}_p^{\prime } \cdot \boldsymbol {u}_p^{\prime } \rangle . \end{equation}
Here, angled brackets denote a Reynolds decomposition according to
$\boldsymbol {u}_p = \langle \boldsymbol {u}_p \rangle + \boldsymbol {u}_p^{\prime }$
, with a single prime denoting a fluctuation from the Reynolds-averaged quantity.
When comparing the particle turbulent kinetic energy,
$k_p$
, and the total granular energy,
$\kappa _p$
, it is observed that there is an additional term that represents the uncorrelated, random component to the total particle energy that exists at the particle scale. This is the information that is lost when filtering Lagrangian particle data to the Eulerian mesh to evaluate
$k_p$
and is termed ‘granular temperature’,
$\Theta$
. This yields an expression for the total granular energy
\begin{equation} \kappa _p = k_p + \frac {3}{2} \langle \Theta \rangle _p, \end{equation}
that is equal to the sum of the correlated and uncorrelated contributions,
$k_p$
and
$3\langle \Theta \rangle _p/2$
.
We compute
$\Theta _p$
directly, by computing the volume-filtered particle volume fraction and velocity. Here, the particle volume, the product of velocity and particle volume and the product of velocity squared and particle volume are extrapolated to the Eulerian grid via trilinear interpolation. These fields are then divided by the Eulerian cell volume and filtered using the Gaussian kernel described in § 2.3, thus yielding the smoothed Eulerian quantities for particle volume fraction,
$\tilde {\alpha }_p$
, particle velocity,
$\tilde {\boldsymbol {u}}_p$
, and particle velocity squared,
$\widetilde {\boldsymbol {u}^2}_p$
. Additional details regarding this formulation for granular temperature can be found in ?. From this, the granular temperature is computed according to
\begin{equation} \Theta _p = \frac {1}{3}\left \lbrack \frac {\text {tr} \left (\widetilde {\boldsymbol {u}^2}_p \right )}{\tilde {\alpha }_p} - \frac {\text {tr}\left (\tilde {\boldsymbol {u}}_p^2\right )}{\tilde {\alpha }_p} \right \rbrack . \end{equation}
These quantities will be leveraged in the following discussion of our findings and utilised to propose models that capture polydisperse effects for both clustering and settling behaviour.
Representative snapshots from the statistically stationary period of each configuration under study. Each slice is an
$x$
$y$
plane at
$z = 0$
. Fluid-phase velocity (grey) is normalised with the polydisperse Stokes velocity,
$\mathcal {V}_{0,10}$
and particles are coloured by diameter (from blue (small) to yellow (large) for distribution A and from pink (small) to red (large) for distribution B).
4. Results
In this work, simulations are evolved from their initial, uncorrelated condition to a statistically stationary state, which is determined based on the temporal variation of the volume mean settling velocity and volume mean variance of volume fraction. After a stationary state is maintained for several characteristic time scales,
$\tau _p$
, representative snapshots are used for analysis. Representative
$x{-}y$
planes of these snapshots are shown in figure 1.
In our analysis, we examine the effect of both polydispersity and volume fraction on clustering behaviour and how this in turn implicates overall settling behaviour. In addressing both clustering and setting phenomena, we draw comparisons between observations in the polydisperse configurations and analogous monodisperse configurations as well as between the two volume fractions.
First, we begin this discussion with averaged behaviour on the scale of the full domain. In doing so, we leverage several mean quantities that are often useful in characterising heterogeneity. We present these quantities first, to paint a broad picture of the clustering and settling behaviour of the configurations studied, but note that examining mean flow quantities in isolation paints an incomplete picture. Then, to present a more complete analysis of flow phenomena, we present a statistical analysis on a particle-wise basis to illuminate the physics underlying domain-scale behaviour.
4.1. Characterisation of mean flow quantities
Mean degree of clustering is often characterised using the parameter,
$\mathcal {D}$
, introduced by Eaton & Fessler (Reference Eaton and Fessler1994) and defined as
\begin{equation} \mathcal {D} = \frac {\left (\sqrt {\langle \alpha _p^{\prime 2}\rangle } -\sqrt {\langle \alpha _p^{\prime 2}\rangle }\Big \vert _0 \right )}{\langle \alpha _p \rangle } \approx \frac {\sqrt {\langle \alpha _p^{\prime 2}\rangle } }{\langle \alpha _p \rangle }, \end{equation}
where
$\langle \alpha _p^{\prime 2}\rangle$
is the variance of the particle volume fraction of the fully developed configuration, also denoted as
$\text {var}(\alpha _p)$
in this work, and
$\langle \alpha _p^{\prime 2}\rangle \Big \vert _{0}$
is the variance of an uncorrelated assembly of particles, assumed here to be null. Given these definitions, the numerator of
$\mathcal {D}$
represents the deviation of the standard deviation of particle volume fraction from a random distribution of particles. As a collection of particles becomes increasingly correlated and heterogeneity in the flow develops, this quantity similarly increases.
Summary of surface loading and statistics on volume fraction for all configurations under study.
While
$\mathcal {D}$
has been traditionally used to characterise clustering behaviour in gas–solid flows, this description is statistically incomplete when the solid phase is polydisperse. The third and fourth moment means of the filtered particle-phase volume fraction are also relevant for characterising the extent of clustering. These quantities are tabulated in table 2 for each of the configurations and are defined as
\begin{align} \text {skew}(\alpha _p) &= \frac {\sum \limits _{i=1}^{N_{{{cells}}}} \left ({\alpha }_p^{(i)} - \langle {\alpha }_p\rangle \right )^3}{(N_{{cells}}-1)\left (\text {var}(\alpha _p)\right )^{3/2}} \end{align}
and
\begin{align} \text {kurt}(\alpha _p) &= \frac {\sum \limits _{i=1}^{N_{{cells}}} \left ({\alpha }_p^{(i)} - \langle {\alpha }_p\rangle \right )^4}{(N_{{cells}}-1)\left (\text {var}(\alpha _p)\right )^{2}}. \end{align}
We note that skewness represents the asymmetry of the distribution of volume fraction, with
$\text {skew}(\alpha _p)=0$
denoting perfect symmetry in the data. Values for skewness less than and greater than null signifying asymmetry skewed toward volume fractions more dilute and more concentrated than the global mean, respectively. In other words, the value of skewness quantifies the propensity for a given distribution of particles to produce correlated regions that are either very dense (positive skewness) or very dilute (negative skewness). Kurtosis is a measure of the ‘tailedness’ of the distribution and indicates how the volume fraction is distributed between the mean and tails. Null represents a normal distribution, while positive (leptokurtic) and negative (platykurtic) values are representative of distributions that are more tightly and loosely spread about the mean, respectively.
In computing the variance, skewness and kurtosis of the particle volume fraction for all the configurations under study, we make several observations. First, we find that the variance in volume fraction is substantially higher for dilute polydisperse configurations relative to their monodisperse counterparts. In particular, the ratios of the standard deviation for distribution
$A/A_0$
and
$B/B_0$
are 2.07 and 2.73, respectively. Interestingly, the standard deviation in particle volume fraction is nearly equivalent between the polydisperse configurations and their monodisperse counterparts at high volume fraction. The difference between
$A$
and
$A_0$
is 1.4 % and the difference between
$B$
and
$B_0$
is 2.5 %. A similar trend is observed when considering the skewness and the kurtosis in the particle volume fraction. Again, the dilute configurations demonstrate much higher deviations in the polydisperse configurations relative to their monodisperse counterparts, and this difference is diminished at higher volume fraction. As previously described, reduced drag is a primary mechanism of clustering. The marked difference in mean clustering behaviour between polydispersed and monodispersed assemblies of particles at dilute configurations and relatively little difference at higher concentrations owes to this. Specifically, at higher volume fractions, there are a sufficient number of particles to consistently disturb the flow and produce regions of reduced drag, thereby initiating clustering events. This, then, reduces the importance of larger particles in the flow. In contrast, for more dilute suspensions, when particles are polydispersed, the presence of larger particles induces larger regions of reduced drag compared with a monodispersed analogue. This then translates to more regions of heterogeneity and clustering.
As previously described,
$\mathcal {D}$
is routinely used as an a posteriori estimator of clustering. We also note that mass loading,
$\varphi = \rho _p \langle \alpha _p \rangle /(\rho _f \langle \alpha _f \rangle )$
, has historically been used as an a priori estimate for predicting clustering. This metric has been shown to be related to
$\mathcal {D}$
in the case of monodispersed assemblies of particles (e.g. Capecelatro et al. Reference Capecelatro, Desjardins and Fox2015; Beetham et al. Reference Beetham, Fox and Capecelatro2021), however, it is agnostic to polydispersity in the particle phase. To underscore this, the mass loading of all the configurations
$A$
,
$A_0$
,
$B$
and
$B_0$
is 50.5 at
$\langle \alpha _p \rangle = 0.01$
and 555.5 for all configurations at
$\langle \alpha _p \rangle = 0.10$
. This, combined with the wide variation in
$\mathcal {D}$
across each of these configurations, highlights the notion that mass loading is incapable of providing differentiation in the clustering behaviour between assemblies of particles that have the same mass loading but have the particle-phase mass partitioned differently. Thus, while mass loading still may serve as an initial indicator of whether or not clustering will occur, our data motivate the need for an alternative a priori quantity that is capable of giving a more complete prediction of clustering behaviour. It is also prudent to note that for a given flow condition, within a statistically stationary period, the clustering behaviour continuously changes, making any a priori predictor a many-to-one mapping to clustering. However, in this work, we emphasise the prediction of the mean clustering behaviour from this stationary regime, thus making it more tractable to reduce this many-to-one mapping. Further, we note that as with any model, the extent of its applicability is reliant on the data used to inform it. Thus, the improved models we develop herein are limited to the flow configurations described in this work. Expanding these notions to different flows, such as those with background turbulence, is reserved for future work.
To this end, we note that a key difference between configurations that contain differing particle distributions at equivalent mass loading is the number of particles present, or alternatively, the amount of particulate surface area present in the flow. We also mention that total particle surface area is likely to be an important predictor of clustering, due to this quantity’s role in particle drag. Based on the data collected as a part of this study, we propose an alternative a priori predictor for the degree of clustering that can be expected for a given gas–solid flow, based on the degree of surface area contact between the phases. We term this ‘surface loading’ and define it as
\begin{align} \mathcal {S} &= \left (\frac {1}{\langle \alpha _f \rangle A_{{cross}}}\right )\left (\frac {\rho _p}{\rho _f}\right )\frac {\pi }{4}\underbrace {\frac {1}{N_p }{\sum \limits _{i=1}^{N_p} \left ({d}_p^{(i)}\right )^2}}_{\left (D_{30}\right )^3/D_{32}} ,\end{align}
where
$A_{\text {cross}}$
is the area of the cross-stream plane and the term in brackets is the square of the surface mean diameter (see Appendix C for more details). It is important to note here, that when
$\mathcal {S}$
tends to zero, this is indicative of extremely fine, dilute particles and when
$\mathcal {S}$
approaches infinity, this represents very concentrated, very large particles. Due to this, we expect that
$\mathcal {D}$
should tend to zero in both limits. In addition, there should exist some critical surface loading for which the degree of heterogeneity is a maximum.
Shown in figure 2, we observe that for each void fraction there is a critical surface loading for which the degree of clustering,
$\mathcal {D}$
, attains a maximum with respect to the surface loading,
$\mathcal {S}$
. For the configurations studied here, this maximum occurs at lower surface loading for higher mean volume fraction, but the degree of clustering overall is higher for the dilute suspensions. Additionally, at higher volume fraction we find that
$\mathcal {D}$
remains nearly constant, further underscoring that clustering behaviour is less sensitive to changes in
$\mathcal {S}$
(e.g. polydispersity) at higher volume fraction.
Deviation of normalised particle-phase volume fraction as a function of surface loading. Circles represent data for
$\langle \alpha _p \rangle = 0.01$
and squares represent data for
$\langle \alpha _p \rangle = 0.10$
. The loosely and densely dashed lines represent the model described in (4.5) for
$\langle \alpha _p \rangle = 0.01$
and
$0.10$
, respectively.
In light of our data, we propose a model for
$\mathcal {D}$
as a function of the surface loading, shown as dashed lines in figure 2, and defined as
\begin{align} \frac {\sqrt {\langle \alpha _p^{\prime ^2} \rangle }}{\langle \alpha _p \rangle } &= \frac {1}{\mathbb {E}\; \mathcal {S}} \exp {\left ( \frac {-\left (\ln {\left (\mathcal {S}\right )}-\mathbb {F} \right )^2}{\mathbb {G}}\right )} ,\end{align}
where the coefficients
$\mathbb {E}$
,
$\mathbb {F}$
and
$\mathbb {G}$
, depend upon the mean volume fraction as
\begin{align} \mathbb {E}\left (\langle \alpha _p \rangle \right ) &= -8.2\langle \alpha _p \rangle + 0.9 ,\end{align}
\begin{align} \mathbb {F}\left (\langle \alpha _p \rangle \right ) &=76.0\langle \alpha _p \rangle - 0.8 ,\end{align}
\begin{align} \mathbb {G}\left (\langle \alpha _p \rangle \right ) &=164.0\langle \alpha _p \rangle - 0.9. \end{align}
It is important to note, however, that since this study considered only two volume fractions, the dependence of the coefficients
$\mathbb {E}$
,
$\mathbb {F}$
and
$\mathbb {G}$
can at most only be described as linear. A more comprehensive model that considers a wider volume fraction is reserved for future work. Further, it is important to note here that prior work (Beetham et al. Reference Beetham, Fox and Capecelatro2021) has noted that there is also a (likely nonlinear) relationship between
$\mathcal {D}$
and the Archimedes (or Froude) number, which has not been captured through
$\mathcal {S}$
. This dependence is also reserved for future work.
In addition to the global clustering behaviour, we also consider global settling behaviour for each configuration studied. As was previously discussed in § 3, the mean settling velocity can also be described as
\begin{equation} \langle u_p \rangle = \langle u_f \rangle _p + \mathcal {V}_0^{\ast}, \end{equation}
where
$\mathcal {V}_0^{\ast}$
is defined as
$\tau _p^{\ast}g$
and
$\tau _p^{\ast}$
is the
$\tau _p/F$
with mean flow quantities,
$\langle d_p \rangle$
and
$\langle \alpha _p \rangle$
, used as argument to the previously defined expressions for the settling time and drag correction factor. The way settling velocity is partitioned into each of these terms is summarised in table 3 and plotted in figure 3.
Summary of the domain mean settling velocities and contributions from
$\langle u_f \rangle _p$
and
$\mathcal {V}_0^{\ast}$
along with the domain Stokes and fluid-phase integral Reynolds numbers.
In view of the data reported in table 3, we observe that the mean settling velocity is greater for each of the polydisperse configurations, relative to their analogous monodisperse configuration for distribution
$B$
, but that this behaviour is inverted for distribution
$A$
, in which the monodisperse configurations have a greater settling velocity compared with their polydisperse counterparts. This effect is more pronounced in the configurations where
$\langle \alpha _p \rangle = 0.01$
. Additionally, we observe that the dilute configurations demonstrate a greater correlation between the degree of clustering,
$\mathcal {D}$
, and the normalised mean settling velocity, as shown in figure 3. In contrast, the configurations with
$\langle \alpha _p \rangle = 0.10$
span a wider range of settling velocities for a much narrow range of
$\mathcal {D}$
. When comparing the normalised mean settling velocity against the surface loading parameter,
$\mathcal {S}$
, however, we note that our data collapse onto a single curve, characterised by,
\begin{equation} \frac {\langle u_p \rangle }{\mathcal {V}_{0,10}} = \frac {2.5}{\left (\mathbb {I}\mathcal {S}\sqrt {2\pi }\right )}\exp \left (-\frac {(\ln (\mathcal {S})-\mathbb {H})^2}{2\mathbb {I}^2}\right ) ,\end{equation}
where
$\mathbb {H} = 0.15$
and
$\mathbb {I}=0.8$
. Notably, we observe that our data collapse for both volume fractions studied in this work, resulting in parameters
$\mathbb {H}$
and
$\mathbb {I}$
that are constant coefficients.
Summary of the normalised mean settling velocity with respect to (a) degree of clustering,
$\mathcal {D}$
and (b) surface loading
$\mathcal {S}$
. The normalised contribution of
$\langle u_f\rangle _p$
relative to
$\mathcal {S}$
is shown in (c). Configurations with
$\langle \alpha _p \rangle = 0.01$
are denoted with solid circles and configurations with
$\langle \alpha _p \rangle = 0.10$
are denoted with solid diamonds. The dashed lines represents a one-to-one correlation in (a) and the model prescribed by (4.10) in (b) and (4.11) in (c).
Noting that we may also consider the mean settling velocity as the sum of
$\mathcal {V}_0^{\ast}$
and
$\langle u_f \rangle _p$
, it is apparent that
$\mathcal {V}_0^{\ast}$
is closed but a model is required for
$\langle u_f \rangle _p$
. Here, we note that the relative contribution to the unclosed portion of mean settling velocity is 2.4 times larger for the higher volume fraction configurations. In light of this, we introduce an alternative formulation for the mean settling velocity, which uses the closed definition for
$\mathcal {V}_0^{\ast}$
and the following model for the fluid velocity seen by the particles:
\begin{equation} \frac {\langle u_f \rangle _p}{\mathcal {V}^{\ast}_{0}} = \frac {\mathbb {J}}{\left (0.6\mathcal {S}\sqrt {2\pi }\right )}\exp \left (-\frac {(\ln (\mathcal {S}))^2}{2\mathbb {I}^2}\right ) + 1, \end{equation}
where the coefficient
$\mathbb {J}$
is 1.65 and 3.9 for
$\langle \alpha _p \rangle = 0.01$
and
$0.10$
, respectively. As noted previously, additional studies at a wider range of volume fractions is needed to develop an accurate model for the dependence of
$\mathbb {J}$
on
$\langle \alpha _p \rangle$
.
In addition the way settling velocity is partitioned into
$\mathcal {V}_0^{\ast} = \tau _p^{\ast} g$
and
$\langle u_f\rangle _p$
, we also report the Stokes number, St, and the Taylor Reynolds number, Re
$_{\lambda }$
, in table 3. Here, the Stokes number is defined as St
$=\tau _p/\tau _{\eta }$
, where the Kolmogorov time scale,
$\tau _{\eta }$
is defined as
$\sqrt {\nu _f/\varepsilon _f}$
and the fluid-phase viscous dissipation is defined as one half the trace of the viscous dissipation tensor given as
\begin{equation} \varepsilon _{f,ij} = \frac {1}{\rho _f} \left \langle \sigma _{f,ik} \frac {\partial u^{\prime \prime \prime }_{f,j}}{\partial x_k} \right \rangle, \end{equation}
with the viscous stress tensor defined as
\begin{equation} \sigma _{f,ij} = \unicode{x03BC} _f \left \lbrack \frac {\partial u_{f,j}}{\partial x_i} + \frac {\partial u_{f,i}}{\partial x_j} - \frac {2}{3}\frac {\partial u_{f,k}}{\partial x_k} \delta _{ij} \right \rbrack . \end{equation}
It is important to note here that, while the Stokes number for the flows considered in this work are quite small, because of the values of the volume fraction considered, we still expect and ultimately observe clustering as a result of the flow being four-way coupled as predicted by Elghobashi (Reference Elghobashi1994). Finally, Re
$_{\lambda } = u_{{rms}}\lambda /\nu _f$
also makes use of the fluid viscous dissipation to compute the Taylor micro-scale,
$\lambda = \sqrt {15 \nu _f/\varepsilon _f}u_{{rms}}$
.
In addition to mean settling velocity, the contributions to turbulent energy in both phases is also helpful in characterising the flow. In addition to the relative streamwise and cross-stream components of the fluid- and particle-phase Reynolds stress tensors, we also present the contributions to total granular energy,
$\kappa$
, from both the particle-phase turbulent kinetic energy (TKE) (
$k_p$
) as well as from the granular temperature. As shown in table 4, we observe that the dilute suspensions exhibit greater contributions from granular temperature to the overall granular energy as compared with the denser suspensions. This is indicative of denser suspensions containing a larger proportion of their particles in clusters, in which their granular energy is strongly correlated with bulk settling. In contrast, dilute suspensions have more particles that are not correlated with clusters. This phenomenon is also evident in the higher ratio of fluid to particle-phase turbulent kinetic energy in the denser suspensions. This is attributable to the fact that larger clusters are able to generate and sustain greater degrees of fluid-phase turbulence.
Summary of the mean turbulent kinetic energy in both phases as well as the relative contributions to total granular energy from particle-phase TKE and granular temperature as well as the relative magnitude of fluid-phase TKE relative to granular energy.
4.2. A brief portrait of clustering
In the following section, we present a brief overview of the implications of polydispersity and volume fraction on cluster formation and composition. The following is based upon a more thorough discussion, which can be found in Appendix E, which draws upon statistical analysis of the Eulerian filtered particle volume fraction and the use of isosurfaces. Because what constitutes a cluster is not strictly defined, we make use of the following regions of clustering and consider a particle volume fraction
$\geqslant 50$
% larger than the global mean to be ‘clustered’:
-
Region A: most dilute regime,
$\alpha _p \leqslant 1.5 \langle \alpha _p \rangle$
. -
Region B: loosely clustered regime,
$1.5 \langle \alpha _p \rangle \lt \alpha _p \leqslant 2.25 \langle \alpha _p \rangle$
. -
Region C: moderately clustered regime,
$2.25 \langle \alpha _p \rangle \lt \alpha _p \leqslant 3.0 \langle \alpha _p \rangle$
. -
Region D: densely clustered regime,
$\alpha _p \gt 3.0 \langle \alpha _p \rangle$
.
We employ these partitions in the remainder of the manuscript and each are designated with colours increasing from white (most dilute) to dark grey (most dense). In addition, it is important to identify an appropriate reference volume fraction with respect to which we can normalise particle-phase concentrations. Here, we utilise the random close-packing efficiency corresponding to distribution
$B$
as it represents the maximum of the particle distributions considered in this work. A justification for this choice is discussed in detail in Appendix D. As previously mentioned, we leave a more extensive discussion of clustering for the interested reader in Appendix E and note here the following key findings relevant to polydisperse clustering.
-
(i) Polydisperse collections of particles form clusters that attain higher packing efficiencies in their cores, as smaller particles can fill voids that would remain empty in a monodisperse assembly of particles with diameters equal to the statistical mean of the polydispersed case.
-
(ii) Large particles are statistically more likely to be found in the densest regions of clusters for all polydisperse configurations considered, however, this effect was more pronounced in the dilute suspensions as compared with the denser suspensions.
-
(iii) We believe that because lower volume fraction flows have fewer particles, the disturbance to the flow caused by large particles is more substantial, as compared with higher volume fraction flows. In the latter case, the clustering patterns are more similar between polydisperse and monodisperse particles due to the amount of disturbance in the fluid phase due to the existence of a greater number of particles.
Summary of the particle size distributions based on local volume fraction for the polydisperse distributions
$A$
(teal) and
$B$
(mustard) at
$\langle \alpha _p \rangle = 0.01$
(left two columns) and
$\langle \alpha _p \rangle = 0.10$
(right two columns). The probability distribution functions (PDFs) are computed based on local volume fraction cutoffs corresponding to regions A, B, C and D, as noted. All plots show the normalised PDF (shaded bars) of the particle diameters in each region of the flow, with the normalised distribution of all the particles shown as a solid black line.
These findings are substantiated by the profiles for both particle size (figure 4) and granular temperature (figure 5) in each of the regions outlined above. Specifically, as the amount of clustering increases (i.e. as one traverses from region A to D), we observe an increase in the likelihood of observing large particles, particularly in the dilute configurations. Additionally, we observe that has the local volume fraction increases, there is a marked decline in granular temperature, an expected result and a validation for the metric in which clustering was assessed.
Summary of the granular temperature based on local volume fraction for the polydisperse distributions
$A$
(teal) and
$B$
(mustard) at
$\langle \alpha _p \rangle = 0.01$
(left two columns) and
$\langle \alpha _p \rangle = 0.10$
(right two columns). Distributions are computed based on local volume fraction cutoffs corresponding to regions A, B, C and D, as noted. All plots show the normalised distribution (shaded bars) of the particle diameters in each region of the flow, with the normalised distribution of all the particles shown as a solid black line.
4.3. A statistical description of particle settling
In this section, we examine the settling behaviour of each of the configurations under study, with a particular focus on how clustering behaviour implicates settling. To this end, we consider the distributions of particle velocity in regions A–D, in direct analogy with the discussion in the previous section. Figure 6 shows these distributions, where the solid lines denote the distribution for each component of velocity for the full domain and the shaded distributions denote the distribution for the given region of the flow.
Beginning with a comparison of the streamwise velocity,
$u_p$
, for the dilute cases (left two columns of figure 6), we find that particles in the most dilute region of the flow (region A) tend to have statistically smaller settling velocities. This is due to two primary factors: (i) the particles are generally lone and would be expected to have smaller velocities for this reason and (ii) particles may be swept up in jet bypassing, which causes much smaller velocities, and occasionally velocities that oppose gravity. As cluster density increases, particles in these regions are increasingly more likely to have larger velocities, due to their correlation with other particles and entrapment in coherent structures. Naturally, the particles that are found in the densest region (region D), attain the largest settling velocities.
In contrast with the dilute cases, the denser cases (left two columns of figure 6) exhibit distributions similar to the global distribution for regions A and B, and only moderately depart from this trend for regions C and D. Again, in these denser regions, particles are statistically more likely to attain greater settling velocities, but the extent for this preference toward larger settling velocities is less prominent than for the dilute cases.
In the cross-stream directions, particles in region A are more likely to attain cross-stream velocities (
$v_p$
and
$w_p$
) further from null, indicating that they are more susceptible to being entrained in the underlying turbulence of the flow. This observation is consistent across all four polydispersed configurations. As the particle density increases, particles increasingly have velocities very near zero. This is indicative of granular temperature that is higher on the outer regions of clusters, as compared with the interior, consistent with what was reported in the previous section. This overall behaviour is mirrored across distributions
$A$
and
$B$
and at both mean volume fractions, however, distribution
$A$
has a narrower band of cross-stream velocities in region D compared with distribution
$B$
, indicating that distribution
$A$
clusters are more ‘rigidly’ packed and particles are less able to move translationally.
While this discussion illuminates the behaviour of particles in each of the regions of the flow based on local volume fraction, it is also instructive to connect settling behaviour with particle size. In § 4.2, we have already demonstrated that particle size and clustering potential are correlated, thus connections between particle size and settling velocity are related through a particle’s likelihood to be involved in clustering.
When carrying out this analysis, we leverage the fact that at steady state, all the forces acting on the particles are in equilibrium. In particular, we draw attention to the drag force felt by the particles. As described in § 2.2, the drag force on a particle is given as
\begin{equation} \boldsymbol {F}^{(i)}_{{drag}} = \frac {m_p^{(i)} \alpha _f\lbrack \boldsymbol {x}_p^{(i)}\rbrack F_D}{\tau _p} \left (u_f[x_p^{(i)}] - u_p^{(i)}\right ), \end{equation}
where
$F_D$
is the correction based on particle Reynolds number and local volume fraction (Tenneti et al. Reference Tenneti, Garg and Subramaniam2011). Since settling velocity is related to the balance between drag (which opposes settling) and mass (see (2.6)), when the drag force increases, settling velocity is hindered. Conversely, when drag is decreased, settling is enhanced. Since drag depends upon not only the particle size, but also the slip velocity and local volume fraction, we note that complex behaviour occurs as particles become correlated.
To illustrate this behaviour, figures 7 and 8 show the relationship between individual particle settling velocities and drag forces against local volume fraction for six bins of particle size ranging from very fine to coarse. By examining particle settling velocity in this way, we note that the finest particles attain the widest range of settling velocities for all configurations and that these particles are located primarily in dilute regions of the flow for lower volume fraction configurations, and exist across all regions – from dilute to dense – in the higher volume fraction configurations. For the denser configurations, these small particles are also exclusively involved in entertainment in strong jet by-passing, as they are the only particles that achieve positive velocities.
Distributions of particle velocities (gravity direction in lighter shading and representative cross-stream velocity in darker shading) for all 4 polydisperse configurations considered. Distributions are shown for dilute to dense regions of the flow (top to bottom). The dark lines represent the full domain distribution of velocities.
This effect is more pronounced in the higher mean volume fraction cases, due to the fact that the clusters formed are larger, thereby producing larger turbulent wakes. The fine particles that exist at increasingly large volume fraction have an increasingly narrow range of attained velocities, indicating that they have become increasingly correlated with surrounding particles in a cluster.
At higher volume fraction, the largest particles exist over a substantially broader range of volume fractions in distribution
$A$
as compared with distribution
$B$
. This suggests that particles are more ‘mixed’ in distribution
$A$
than
$B$
, where large particles exist only in higher volume fraction areas. This indicates that the clusters formed in distribution
$B$
tend to contain large particles only in the cluster ‘cores’, whereas in distribution
$A$
, large particles are located throughout clusters. For the dilute configurations, we note that there exists a similar, but more pronounced effect of what we observe in their denser counterparts. In particular, particles of increasing size have substantially more overlap as volume fraction increases for distribution
$A$
, however, the larger particle sizes are extremely stratified in distribution
$B$
. This is suggestive of the existence of clusters formed by particles of like size in the case of larger particles, and clusters formed by the smaller to moderately sized particles.
Settling velocity of individual particles plotted against normalised local particle volume fraction for all four polydisperse cases under consideration. Colours correspond to six increasing diameters
$d_p = (0.3, 0.8, 1.3, 0.9, 2.4, 2.9)$
(mm) and the colour maps for distributions
$A$
and
$B$
used throughout. The data plotted represent the particles within the diameters listed
$\pm 0.03$
mm.
Drag force on individual particles plotted against normalised local particle volume fraction for all four polydisperse cases under consideration. Colours correspond to six increasing diameters
$d_p = (0.3, 0.8, 1.3, 0.9, 2.4, 2.9)$
and the colour maps for distributions
$A$
and
$B$
are used throughout. The data plotted represent the particles within the diameters listed
$\pm 0.03$
mm.
Figure 8 shows a similar analysis for the drag felt by each particle. Here, we find that the force of drag increases with both volume fraction and particle diameter for the dilute cases. This is due to the fact that clusters are smaller and composed of fewer particles, making these particles more susceptible to experiencing increased drag due to a higher slip velocity. In contrast, at higher mean volume fraction, when the clusters are larger, only the outer particles in the cluster – those particles that compose the outer ‘shell’ of the cluster – will observe this effect, while the inner core will experience reduced drag due to the number of particles surrounding them and the resulting reduced slip velocity.
Finally, we consider the settling velocity as a function of particle size, directly. To this end, figure 9 plots the settling velocity against diameter for each particle in the system. A solid dark line represents the moving average of settling velocity as a function of diameter. Here, we note that, on average, the settling velocity of particles increases with increasing particle size before eventually approaching an asymptotic limit. This is observed for all the configurations considered in this work. While this is the trend on average, we also note that the smallest particles experience the widest range of velocities, with some particles having positive velocities (indicative of being swept up in ‘jet bypassing’ events – the upward-moving gas that results from a large cluster moving downward) and others having high downward velocities, due to entrainment in clusters. This phenomenon is observed across all four polydisperse configurations, however, we note that the settling velocities are more widespread for distribution
$A$
compared with
$B$
. We also note that the spread of particle velocities for very fine particles is more widespread at
$\langle \alpha _p \rangle = 0.10$
.
Settling velocity of particles with respect to particle diameter (dots) for all four polydisperse configurations (top figures). The PDF of particle diameters are shown in light shading in the background of each panel. The horizontal dashed lines represent the predicted settling velocity for particles corresponding to
$D_{10}$
,
$D_{32}$
and
$D_{43}$
using (4.15). The heavy solid line is the mean velocity as a function of particle diameter. The vertical line delineates the diameter at which settling is enhanced for smaller particles and hindered for larger particles when using (4.15) as a predictor (shown as a densely dotted line). The loosely dashed line represents the prediction of Stokes velocity as a function of particle size and the red solid is the prediction of settling velocity according to the proposed model shown in (4.17) and (4.18). The bottom of each figure shows the distribution of settling velocity corresponding to the coarse particle bins shown beneath.
As previously summarised, we point out that all contemporary models for polydispersed settling to our knowledge are agnostic to clustering in the particulate phase. One such example of a commonly used model is that of de’Michieli Vitturi et al. (Reference de’Michieli Vitturi, Esposti Ongaro and Engwell2023), who generalised the depth-averaged shallow-water equations to include a dispersed phase and employed an Eulerian–Eulerian approach to simulate a wide range of volcanic phenomena. As a result of this work, a depth-averaged settling velocity, denoted as
$\mathcal {V}_s$
, with a model for the drag coefficient originally developed using kinetic theory arguments (Gidaspow Reference Gidaspow1994) was proposed. This model is given as
\begin{align} \mathcal {V}_s &= \sqrt { \frac {4}{3C_D}g {d}_p \left (\frac {\rho _p-\rho _f}{\rho _f}\right )} ,\end{align}
\begin{align} C_D &= \begin{cases} \dfrac {24}{Re}(1 + 0.15 {Re}^{0.687}) & {Re}\leqslant 1000, \\ 0.44 & {Re} \gt 1000, \end{cases} \end{align}
where we note that Re
$= \mathcal {V}_s d_p/\nu _f$
, which requires (4.15) to be solved numerically.
When using this expression to predict mean settling behaviour, it is necessary to choose a mean diameter as argument. Several statistical mean diameters can be chosen for polydisperse assemblies of particles (see Appendix C), and the predictions of both of these models with three different mean diameters as argument (
$D_{10}$
,
$D_{20}$
and
$D_{32}$
) are summarised in table 5. In examining these predictions, we note that the
$\mathcal {V}_{s}$
prediction is a more reliable predictor for the monodispersed configurations, however, still not consistently accurate. For the polydispersed configurations, using the surface mean diameter,
$D_{20}$
, and the model for
$\mathcal {V}_s$
provides the best approximation for mean settling behaviour. The use of the surface mean diameter for more accurate results is perhaps intuitive, as clustering and drag are intimately connected and drag is related to the degree of surface area contact between fluid and particles.
Summary of the mean settling velocity for each configuration compared with Stokes law and the settling law of de’Michieli Vitturi et al. (Reference de’Michieli Vitturi, Esposti Ongaro and Engwell2023) with several mean diameters used as input. Note that for monodisperse assemblies, the single value for each settling law is reported under the columns corresponding to the result using the
$d_{10}$
diameter as argument. All velocities are shown in metres per second.
While
$\mathcal {V}_{s,20}$
is reasonably predictive for global mean settling behaviour, it is important to note that both
$\mathcal {V}_0$
and
$\mathcal {V}_s$
fall short when predicting the settling behaviour as a function of particle diameter. To illustrate this, figure 9 shows the predictions of both Stokes (loosely dashed line) and the model of de’Michieli Vitturi et al. (Reference de’Michieli Vitturi, Esposti Ongaro and Engwell2023) (densely dashed line) as continuous functions of particle diameter.
While the model described by (4.15) demonstrates a clear improvement over a Stokes assumption, neither model captures the appropriate mean settling behaviour of particles over the range of diameters considered. Additionally, the concavity of these predictions is in direct opposition to the Euler–Lagrange data. The existing models suggest that the settling velocity should continue to increase with increasing diameter, whereas our data suggest that settling velocity approaches an asymptotic limit, due to cluster formation. Interestingly, we notice that there exists a critical particle diameter such that particles smaller than
${d}_p^{{crit}}$
exhibit enhanced settling, and particles larger than
${d}_p^{{crit}}$
exhibit hindered settling. This behaviour is directly connected to clustering behaviour. Smaller particles on average experience reduced drag due to their proximity to or presence within clusters, thus resulting in enhanced settling. Conversely, larger particles experience hindered settling because of their entrainment in coherent structures much larger in size than the particles themselves. This then results in particles that feel increased drag compared with the drag they would feel as a lone particle thus yielding hindered settling velocities.
For both distributions
$A$
and
$B$
, the critical diameter that delineates between enhanced and hindered settling is larger at
$\langle \alpha _p \rangle = 0.01$
than for
$\langle \alpha _p \rangle = 0.10$
. Interestingly, these values correspond approximately to
${d}_{20}$
and
${d}_{30}$
for distributions
$A$
and
$B$
at
$\langle \alpha _p \rangle = 0.10$
, respectively, perhaps indicating that for a distribution that favours fine particles, the surface area contact between the phases is more important than for distributions that contain more moderately sized particles, as in distribution
$A$
. At
$\langle \alpha _p \rangle = 0.01$
,
${d}_p^{{crit}}$
corresponds to values between
${d}_{30}$
and
${d}_{32}$
for both distributions. This suggests that for more dilute suspensions, higher-order statistical means are more relevant for predicting settling behaviour. We also note that the maximum settling velocity is slightly higher for the cases at lower volume fraction for both distributions
$A$
and
$B$
. We postulate that this is due to the nature of the cluster structures. In the denser suspensions, the clusters have a larger cross-sectional area which in turn results in a larger drag on the cluster itself.
In light of these data, we propose an improved model in which the expression for the settling velocity,
$\mathbb {V}_s$
, retains the same functional expression as for
$\mathcal {V}_s$
\begin{align} \mathbb {V}^{(i)}_s &= \sqrt { \frac {4}{3C_D}g {d}_p \left (\frac {\rho _p-\rho _f}{\rho _f}\right )} ,\end{align}
but the expression for
$C_D$
is modified as
\begin{align} \mathbb {C}_D &= \mathbb {M}\left \lbrack \frac {24}{\mathbb {N} {Re}_p}\left (0.2 + 0.01 \left ( \mathbb {N} {Re}_p\right )^{0.9}\right ) +0.35 \left (\mathbb {N} {Re}_p\right ) \right \rbrack \nonumber\\ &\quad - \frac {2}{\left (\mathbb {N} {Re}_p\right )^2 + 0.09} + \frac {\mathbb {P}\;\mathcal {W}}{1+ {Re}_p^2} ,\end{align}
where
$\mathbb {M}$
,
$\mathbb {N}$
and
$\mathbb {P}$
are constant coefficients specified using the highly resolved data in this study and summarised in table 7.
Predictions of the proposed model summarised alongside the mean settling velocity from the Euler–Lagrange studies and the Stokes prediction for the monodisperse assemblies. All velocities are shown in metres per second.
Summary of the model coefficients for the eight polydisperse configurations studied.
To develop this improved model, we begin by computing the coefficient of drag,
$C_D$
, required to result in the settling velocity observed in the highly resolved simulations, as prescribed by (4.17). These results are shown in figure 10 for all the particles in the configuration. The dark shaded line represents the mean of the data as a function of particle size. In this study, we did not observe Reynolds numbers greater than
$1000$
, and as such have left the upper limit equal to the historical model of
$C_D = 0.44$
.
Coefficient of drag required to ensure
$\mathbb {V}_s^{(i)} = u_p^{(i)}$
, using the expression in (4.17). Values for individual particles are shown as shaded dots, the mean of these data as a function of particle diameter is shown as a dark black line and the models for
$C_D$
are shown as dashed lines (black dashed is the baseline model of (de’Michieli Vitturi et al. Reference de’Michieli Vitturi, Esposti Ongaro and Engwell2023) and the red dashed is the proposed new model shown in (4.17).
Here, we note that the model for
$C_D$
closely follows the mean behaviour of the highly resolved data, with the exception of deviations at very small Reynolds numbers for all configurations studied and at high Reynolds numbers in the denser suspensions. Deviations at low Reynolds numbers were required to ensure accurate predictions of the settling velocity. This stems from the fact that the way in which
$C_D$
was computed does not take into account differences in the mean slip between the phases. This effect is particularly important for the small particles. Additionally, we observe that the behaviour for
$C_D$
at higher Reynolds number is not well captured by the proposed model for dense suspensions, however, the settling velocity for larger particles is not sensitive to these deviations in the drag coefficient.
Exemplary case (distribution
$A$
,
$\langle \alpha _p \rangle = 0.01$
) demonstrating the relative contribution of each of the terms in the proposed coefficient of drag model in (4.17).
In the proposed model for
$C_D$
(4.18), the first term is a modification of the original definition of
$C_D$
from Gidaspow (Reference Gidaspow1994). The second, linear term captures hindered settling for larger particles and the third term captures enhanced settling for smaller particles. Finally, the fourth term introduces stochasticity into the model through
$\mathcal {W}$
, a Wiener process (Higham Reference Higham2001) that depends upon the particle Reynolds number and is implemented numerically as shown in figure 12. Stochasticity has been similarly introduced into drag models in the literature, for example in the work of Lattanzi et al. (Reference Lattanzi, Tavanashad, Subramaniam and Capecelatro2022). The relative contribution to the drag coefficient for each of these terms is demonstrated for the model corresponding to distribution
$A$
in figure 11.
Summary of the numerical implementation of the Wiener process (Higham Reference Higham2001). This code snippet is written in the style of Matlab, where ‘randn’ represents a normally distributed random number bounded by [0, 1].
The settling velocity predictions resulting from this model are shown in figure 9 as solid red lines for the polydisperse configurations and summarised in table 5 for the monodispersed assemblies of particles.
Since this work considered only two volume fractions and two distributions of particle diameter, we leave the values for the model coefficients summarised in table 7 for each configuration, rather than proposing closures for each based on flow parameters. A study of additional points in the parameter space is required before a robust functional dependence of the coefficients on polydisperse and volume fraction can be developed. In particular, we postulate that
$\mathbb {M} = f(\sigma , \unicode{x03BC} , {d}_{10}, {d}_{20}, {d}_{30}, {d}_{32}, {d}_{43}, \ldots )$
and
$\mathbb {N} = f(\langle \alpha _p \rangle )$
, however, this is left for future work.
5. Conclusions
The work discussed here represents, to the authors’ knowledge, the most extensive study of clustering and settling behaviour of strongly coupled gas–solid flows at the mesoscale. In particular, we investigate the effects of polydispersity on particle clustering and settling behaviour using highly resolved Euler–Lagrangian simulations of two polydispersed distributions of particles and two analogous monodispersed distributions of particles. These four particle distributions were studied at two volume fractions (
$\langle \alpha _p \rangle = 0.01$
and
$\langle \alpha _p \rangle = 0.1$
) and all simulation parameters were sampled to align with values typical of PDCs.
Due to the strong coupling between the phases and the presence of a gravitational body force, coherent structures in the form of clusters spontaneously emerge, thereby altering settling behaviour as compared with the uncorrelated initial dispersion of particles. To date, the extent to which polydispersity is implicated in clustering structure and consequently on settling behaviour has been largely unquantified. To this end, this work demonstrated that mass loading – which has historically been used as a a priori estimate for predicting clustering – is insufficient for predicting the degree of clustering, particularly for polydispersed particles. This owes to the fact that mass loading is by definition agnostic to how mass is partitioned throughout a volume. Based on the data collected in this study, we propose an alternative a priori predictor for the degree of clustering expected in gas–solid flows, which we refer to as surface loading. This quantity leverages the mean surface diameter of particles and is shown to predict both the degree of clustering and mean settling velocity through two new models proposed in this work.
Additionally, this work identified that the surface mean diameter,
$D_{20}$
, is the best mean diameter to choose for use with existing models, such as the settling velocity model of de’Michieli Vitturi et al. (Reference de’Michieli Vitturi, Esposti Ongaro and Engwell2023) for the prediction of mean flow behaviour. This is perhaps intuitive due to the intimate connection between granular surface area and drag, and the relationship between drag and clustering. While using this diameter produces improved predictions for global mean settling behaviour, however, we demonstrate that existing models fail to capture the settling behaviour across a distribution of particle diameters. Importantly, we observe that fine-grained particles experience enhanced settling and coarse particles experience hindered compared with existing models. Further, existing models predict that settling velocity continues to increase with particle size, which is in contrast with our observation that settling velocity increases initially with increasing particle diameter, but quickly slows and approaches an asymptotic limit. This is attributed entirely to the existence of heterogeneity in the flow and how large particles initiate and become correlated through clustering. To this end, we propose a new model for the coefficient of drag that produces accurate settling velocity predictions as a function of particle diameter. Further, our proposed model can be immediately implemented in existing geophysical flow solvers, such as in de’Michieli Vitturi et al. (Reference de’Michieli Vitturi, Esposti Ongaro and Engwell2023).
Although this work represents the most highly resolved study of polydisperse clustering and settling behaviour to date and an initial step toward improved reduced-order models, future efforts are needed to build robustness into the models proposed herein. In particular, additional volume fractions and polydispersed assemblies of particles are required to formulate more comprehensive closures for the model coefficients proposed in this work.
Acknowledgements
The computing resources and assistance provided by Oakland University are greatly appreciated.
Funding
This work was supported by funding provided by the National Science Foundation (award number 2346972) and the National Aeronautics and Space Administration (NASA), under award number 80NSSC20M0124, Michigan Space Grant Consortium (MSGC). E.C.P.B. was supported by a NERC Independent Research Fellowship (NE/V014242/1).
Declaration of interests
The authors report no conflict of interest.
Appendix A. A brief background on pyroclastic density currents
One society-critical example of a polydisperse gravity-driven, sedimenting flow, and the motivation for this work, is DCs – the fast-moving, gravity-driven flow of particulate matter resulting from the collapse of an ejected volcanic column, collapse of a volcanic lava (dome) or proximal material perched on steep slopes (Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020). Pyroclastic density currents are the most destructive volcanic process and can cause extensive damage to human settlements, infrastructure and ecosystems (Breard & Lube Reference Breard and Lube2017; Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020; Breard et al. Reference Breard, Dufek, Charbonnier, Gueugneau, Giachetti and Walsh2023). Understanding the physics behind PDCs and accurately predicting their behaviour are essential due to their impact on society and the environment. However, as discussed previously, their inherent complexity, owing to their multiphase nature, makes accurate prediction and the formulation of robust models challenging.
In general, PDCs can be partitioned into three flow regions: An upper dilute buoyant (and loosely coupled) region, a dilute to intermediate density current region and a concentrated basal underflow (granular flow dominant) region (see figure 13). These regions are coupled and, as such, evolve together. In the uppermost dilute region, the flow is buoyant, forming a turbulent thermal cloud where particles are dispersed at very low mass loading, resulting in one-way coupled behaviour with the fluid phase (i.e. the particles tend to act primarily as tracers (Elghobashi Reference Elghobashi1991) and are unlikely to become correlated with one another, outside the mechanism of preferential concentrations due to underlying turbulent structures). In contrast, the bottom-most concentrated region is comprised of very high concentrations of particles (
$\alpha _p \gtrapprox 0.3$
) and is dominated by particle–particle interactions. Recent work (Breard et al. Reference Breard, Lube, Jones, Dufek, Cronin, Valentine and Moebis2016; Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020) showed that within PDCs, the non-turbulent underflow and fully turbulent ash-cloud areas were connected by an intermediate zone characterised by cluster-induced turbulence (CIT). In this region, the dispersed and continuous phases are two-way coupled due to sufficiently high concentrations of particles and mass loadings higher than unity throughout the layer. This strong coupling between the fluid and particles leads to the formation of mesoscale structures in the form of clusters (Ferrante & Elghobashi Reference Ferrante and Elghobashi2003), thereby resulting in altered particle settling when compared with a homogeneous particulate phase (Breard et al. Reference Breard, Lube, Jones, Dufek, Cronin, Valentine and Moebis2016; Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020). The mechanism for this development of heterogeneity is primarily due to the reduced drag particles are subjected to due to the presence of other nearby particles. When drag is reduced, these particles can attain higher settling velocities and become correlated with their neighbours, thereby forming coherent structures, such as clusters. Because this intermediate region is critical in determining the sedimentation rate into the concentrated basal underflow (CBU) region, it serves as the motivation and context for this work.
(a) Pyroclastic density current on 13 September 2012 at Fuego volcano (Guatemala, photo courtesy of V. Bejarano). (b) Sketch of the anatomy of a PDC, with CBU fed by settling of particulates from the dilute ash cloud (DAC) where abundant mesoscale clusters with CIT occurs. The upper part of the PDC is made of the co-PDC ash cloud, which forms thermals that rise buoyantly and feeds co-PDC plumes that can reach heights up to tens of kilometres in the atmosphere.
Historically, the study of PDCs has spanned field measurements, analogue experiments, and computational methods. Direct observation of PDCs is typically done post-eruption, providing an average of intrinsically transient behaviour. However, recent large-scale experiments (Dellino et al. Reference Dellino, Zimanowski, Büttner, La Volpe, Mele and Sulpizio2007; Lube et al. Reference Lube, Breard, Cronin and Jones2015; Breard & Lube Reference Breard and Lube2017; Breard et al. Reference Breard, Dufek and Lube2018, Reference Breard, Jones, Fullard, Lube, Davies and Dufek2019; Brosch et al. Reference Brosch, Lube, Esposti-Ongaro, Cerminara, Breard and Meiburg2022) have been conducted to understand the behaviour of PDCs and propose models for how PDC properties relate to quantities of interest. These experiments are capable of providing real-time measurements but the challenges associated with experimentally probing gas–solid flows remain and results are still more coarsely grained than computational methods.
To this end, recent 1-D (Bursik & Woods Reference Bursik and Woods1996), 2.5-D (Keim & de’Michieli Vitturi 2025) and 3-D multiphase flow models (Tonoyama et al. 2025) have been used to quantify the important physical processes of PDCs that are often inaccessible by large-scale experimentation and develop improved reduced-order forecasting models. To be successful, these models require accurate predictions of polydisperse settling behaviour. While current state-of-the-art models can reveal some aspects of the internal structure of PDCs, they have not captured the heterogeneous effects of clustering to date. This is largely because highly resolved studies of the settling behaviour and underlying physical processes present in polydisperse assemblies of particles has been, to date, largely understudied.
Appendix B. A brief note on distribution terminology
While the applications of polydisperse gas–solid flows are far reaching, ranging from natural to industrial flows, the parameters under study in this work are drawn from PDCs. In the geoscience community, it is common to quantify the distribution of sedimentary particles using a convention referred to as the ‘
$\phi$
scale’. Since this convention is less common in other areas of science and the applications of this study are far-reaching outside the geoscience community, we present here a very brief discussion connecting the
$\phi$
scale with a standard log-normal distribution defined in terms of millimetres.
As such, the
$\phi$
scale was developed on the notion that sediment behaviour is a function of particle diameter squared. Thus,
$\phi \equiv -\log _2({\textrm{d}}_p)$
is used as the basis for this scale (Krumbein Reference Krumbein1936). It is important to note that in the definition of
$\phi$
the particle diameter,
${d}_p$
, has units of millimetres. Since the basis of the definition of
$\phi$
is log
$_2$
, then the parameter
$\phi _{{sorting}}$
describes the normal distribution of the particles in terms of
$\phi$
. This also implies that the distribution of particle diameter in units of mm is log-normal.
To traverse these two definitions for particle distribution, one can make use of the following relationships:
\begin{align} \langle d_p \rangle &= 2^{-\phi } \nonumber \\ \unicode{x03BC} &= \ln (\langle {d}_p\rangle ) = \ln (2^{-\phi }) \nonumber \\ \sigma &= \phi _{{sorting}}^2 \ln (2)^2 ,\end{align}
where
$\langle d_p \rangle$
represents the mean particle diameter in millimetres and
$\unicode{x03BC}$
and
$\sigma$
represent the standard log-normal parameters, such that the probability distribution function for particle diameter is defined as
\begin{equation} f_{d_p} = \frac {1}{{d}_p \sigma \sqrt {2\pi }}\exp {\left ( -\frac {(\ln {d}_p - \unicode{x03BC} )^2}{2\sigma ^2}\right )} .\end{equation}
Appendix C. Statistical diameters
When considering assemblies of particles of varying size, as is done in this work, it is often informative to quantify how particles are distributed by using the following statistical diameters:
-
(i)
$D_{10}$
represents the mean diameter in the usual, arithmetic sense as(C1)This expression can be equivalently expressed in terms of mass fraction, where
\begin{equation} D_{10} = \frac {\sum \limits _{i=1}^{N_p} {d}_p^{(i)}}{N_p} .\end{equation}
$x_j$
and
$n_j$
, are the mass fraction and number of particles of size
${d}_j$
and
$J$
represents the total number of bins the particles are divided into(C2)
\begin{equation} D_{10} = \frac {\sum \limits _{j=1}^{J} {d}_j}{\sum \limits _{j=1}^{J} {3 x_j}/{(4 \rho _p {d}_j^3})}. \end{equation}
-
(ii)
$D_{20}$
represents the surface mean diameter. In other words, this is diameter of a monodisperse assembly of
$N_p$
particles with the same total surface area of the polydisperse assembly. It is defined as(C3)In terms of mass fraction, this can also be defined as
\begin{equation} D_{20} = \sqrt {\frac {\sum \limits _{i=1}^{N_p} \left ({d}_p^{(i)}\right )^2}{N_p}}. \end{equation}
(C4)
\begin{equation} D_{20} = \sqrt {\frac {\sum \limits _{j=1}^{J} x_j/{d}_j}{\sum \limits _{j=1}^{J} x_j/({d}_j^3)}}. \end{equation}
-
(iii)
$D_{30}$
represents the volume mean diameter. Similar to the surface mean diameter, this value represents the diameter of a monodisperse assembly of
$N_p$
particles with the same total volume as the polydisperse assembly. It is given as(C5)Equivalently in terms of mass fraction
\begin{equation} D_{30} = \Bigg (\frac {\sum \limits _{i=1}^{N_p} \left ({d}_p^{(i)}\right )^3}{N_p}\Bigg )^{1/3}. \end{equation}
(C6)
\begin{equation} D_{30} = \Bigg (\frac {1}{\sum \limits _{j=1}^{J}\left ( x_j/{d}_j^3\right )}\Bigg )^{1/3}.\end{equation}
-
(iv)
$D_{32}$
denotes the surface moment mean diameter, commonly referred to as the Sauter mean diameter. This is the diameter required for a monodisperse assembly of
$N_p$
particles to have the same ratio of volume to surface area as the polydisperse assembly(C7)In terms of mass fraction, this is defined as
\begin{equation} D_{32} = \frac {\sum \limits _{i=1}^{N_p} \left ({d}_p^{(i)}\right )^3}{\sum \limits _{i=1}^{N_p} \left ({d}_p^{(i)}\right )^2} = \frac {D_{30}^3}{D_{20}^2} .\end{equation}
(C8)
\begin{equation} D_{32} = \left (\sum \limits _{j=1}^{J} (x_j/{d}_j)\right )^{-1}. \end{equation}
-
(v)
$D_{43}$
denotes the volume moment mean diameter. This quantity is an indicator of which particle sizes contain a majority of the particle volume(C9)Similarly, the definition based on mass fraction is given as
\begin{equation} D_{43} = \frac {\sum \limits _{i=1}^{N_p} \left ({d}_p^{(i)}\right )^4}{\sum \limits _{i=1}^{N_p} \left ({d}_p^{(i)}\right )^3} .\end{equation}
(C10)
\begin{equation} D_{43} = \sum \limits _{j=1}^{J}\left ( x_j {d}_j\right ) .\end{equation}
The values corresponding to each of these statistical mean diameters for distributions
$A$
and
$B$
are summarised in table 8.
Summary of the arithmetic mean, surface area, volume, Sauter and volume moment mean diameters, distribution width,
$\tilde {\sigma }$
, and maximum random close-packing efficiency resulting from (D2) for the particles studied in each configuration. A complete description of the definitions of the statistical diameters can be found in Appendix C.
Appendix D. A brief discussion on the close-packing potential of polydispersed particles
An important comparison we draw when describing clustering is that of the local volume fraction compared with the theoretical close-packing volume fraction, denoted
$\alpha _{{rcp}}$
. As one might anticipate, the packing efficiency observed in clusters, along with their shape, is intimately connected with settling behaviour.
It has been previously established that the random close-packed (RCP) volume fraction of monodisperse spheres is
$\alpha _{{rcp}} = 0.64$
(Kansal, Torquato & Stillinger Reference Kansal, Torquato and Stillinger2002; Farr & Groot Reference Farr and Groot2009; Farr Reference Farr2013; Desmond & Weeks Reference Desmond and Weeks2014). This value represents the maximum achievable volume fraction for a randomly packed arrangement of spheres. While close-packing efficiencies for polydispersed assemblies of spheres are less established in comparison, several theoretical and computational studies have quantified the RCP efficiency for log-normal distributions of spheres (Kansal et al. Reference Kansal, Torquato and Stillinger2002; Farr & Groot Reference Farr and Groot2009; Farr Reference Farr2013; Desmond & Weeks Reference Desmond and Weeks2014). In Farr & Groot (Reference Farr and Groot2009), an algorithm that maps 3-D configurations of spheres onto a 1-D system of rods has shown success in accurately and efficiently estimating the RCP packing efficiency of log-normally distributed spheres. This approach, known as the rod packing (RP) algorithm (Farr & Groot Reference Farr and Groot2009), was subsequently validated (Farr Reference Farr2013) and also resulted in a closed form expression (D2). In this work, Farr (Reference Farr2013) observed that the RCP packing efficiency depended only upon a measure of the distribution ‘width’, denoted
$\tilde {\sigma }$
. This parameter is computed using the ratio of the volume moment diameter to the surface moment diameter, as
\begin{equation} \tilde {\sigma }=\sqrt {\ln \left (\frac {{d}_{4,3}}{{d}_{3,2}}\right )}. \end{equation}
Here, we make note that
${d}_{4,3}$
is sensitive to large particles and
${d}_{3,2}$
is sensitive to small particles (for more details, see Appendix C). When the ratio to these quantities is increasingly larger than 1 this implies a higher proportion of very small particles to larger particles. Naturally, the RCP efficiency for distributions containing a greater number of very small particles (i.e. distributions with increasingly large
$\tilde {\sigma }$
), will be higher. The model relating
$\alpha _{{rcp}}$
to
$\tilde {\sigma }$
based on results of the RP algorithm and developed by Farr (Reference Farr2013) is given as
\begin{align} \alpha _{\text {rcp}}(\tilde {\sigma }) &= 1-0.57\textrm{e}^{-\tilde {\sigma }} + 0.2135\textrm{e}^{\frac {-0.57\tilde {\sigma }}{0.2135}} \nonumber\\ &\quad + 0.0019\big(\cos \big(2\pi \big(1-\textrm{e}^{-0.75\tilde {\sigma }^{0.7} - 0.025\tilde {\sigma }^4}\big)\big)-1\big). \end{align}
We use this expression to infer the theoretical maximum RCP packing efficiency based on
$\tilde {\sigma }$
of the log-normal distributions under study. The statistical values relating to each of the polydisperse configurations studied are listed in table 8. For completeness, we note that while both volume fractions each for distributions
$A$
and
$B$
were sampled from the same log-normal distribution, since the number of discrete particles in the systems differ (fewer for
$\langle \alpha _p \rangle = 0.01$
compared with
$\langle \alpha _p \rangle = 0.10$
), there are minimal differences in the statistical descriptions between the two volume fractions.
As one might anticipate, the log-normally distributed particles exhibit a higher maximum RCP volume fraction, as small particles in the distribution occupy voids that would otherwise remain in a monodisperse packing arrangement of particles with a diameter equal to the mean particle diameter. Due to this, we use the RCP packing efficiency of distribution B,
$\alpha _{{rcp}}^{B} = 0.7038$
as the normalisation factor for volume fraction for all the configurations studied. Visualisations for the RCP packing efficiencies for distributions
$A$
,
$B$
,
$A_0$
and
$B_0$
are shown, along with their associated
$\alpha _{{rcp}}$
in figure 14.
It is notable that while the values reported in table 8 are the random close-packing efficiency limits for the entire collection of particles. However, when these particles spontaneously cluster, the subsets of particles contained within clusters do not necessarily mirror the overall distribution of particles. In other words, if larger particles are more likely to be found in clusters than smaller particles, then the close-packing efficiency of the subset of particles involved in a cluster will differ from what would be expected for an assembly representative of the full domain.
Random close-packed configurations for the configurations under study. Here,
$\alpha _{{rcp}}$
, denotes the RCP volume fraction from the RP algorithm. Visualisations generated by the Kansal Torquatoa Stillinger algorithm (Kansal et al. Reference Kansal, Torquato and Stillinger2002) and RP software (Farr Reference Farr2013).
Appendix E. An extensive discussion on clustering formation and behaviour
Overview of the clustering patterns observed in distributions
$A$
(left figures) and
$A_0$
(right figures) at
$\langle \alpha _p \rangle = 0.01$
(top row) and
$0.10$
(bottom row).
Following the four region convention for delineating regions of increasing correlation in the flow and normalisation with respect to the maximum random close-packing efficiency for distribution B as described above, we plot the local volume fraction and corresponding contours for four representative
$x{-}z$
planes of each configuration studied (see figures 15–17). Qualitatively, we observe that the contours defining clusters in all uniform distributions,
$A_0$
and
$B_0$
, are smoother and achieve lower packing efficiency in their cores. This owes to the increased potential packing efficiency for polydispersed assemblies of particles.
In comparing these uniform distributions of particles against each other, we observe that the configurations with larger particle diameter (distribution
$A_0$
) contain larger regions of correlated particles (i.e. larger clusters), as compared with a smaller particle diameter (distribution
$B_0$
) which results in a greater number of smaller clusters. Further, clusters tend to be denser in the case of smaller particles, however, this stems from the difference in the total number of particles present in the domain since monodispersed particles of any size have the same maximum packing efficiency.
In comparing the polydisperse configurations with each of their analogous uniform distributions (i.e. distributions
$A$
and
$A_0$
), we observe that at low particle volume fraction, the polydisperse distributions tend to have smaller, more numerous clustered regions with higher packing of particles, compared with their monodisperse counterparts (
$A_0$
and
$B_0$
), which exhibit less densely packed clusters that are more sparsely arranged throughout the domain. In the denser suspensions, cluster shape is qualitatively similar between the polydisperse and monodisperse configurations, however, the polydisperse clusters achieve a higher density of particles than their monodisperse analogues. These observations are illustrated in figures 15 and 16, where,
$x{-}z$
planes of the particle volume fraction are shown with the three contours discussed previously.
In addition to examining streamwise planes of data, we also consider the cross-stream constitution of clusters. To this end, for each
$y{-}z$
plane along the gravity-aligned direction,
$x$
, we sum the number of Eulerian cells containing a local particle volume fraction greater than
$1.5\langle \alpha _p \rangle$
. Multiplying this value by the
$y{-}z$
area of the cells thus yields a measure of total cross-stream cluster area. While this does not delineate whether this cross-sectional area is comprised of one or several clusters, it aids in quantifying the relative structural differences in clustering. In comparing the cluster cross-sectional area curves for
$A$
and
$A_0$
at
$\langle \alpha _p \rangle = 0.01$
, we notice that the maximum cluster cross-sectional area is substantially greater for configuration
$A_0$
, substantiated by the observation that clusters are qualitatively larger in size but fewer in quantity. We also observe that the cross-sectional cluster area for distribution
$A$
is more uniformly distributed throughout the streamwise direction, indicating that clusters are more evenly spread throughout the domain as compared with distribution
$A_0$
. In contrast to these observations at low volume fraction, we observe that the regions of maximum cross-sectional cluster area are comparable between distributions
$A$
and
$A_0$
at
$\langle \alpha _p \rangle = 0.10$
and that for both configurations, these peaks are observed at similar intervals.
Some consistency is seen in the distribution
$B$
and
$B_0$
configurations, with a few exceptions. First, at
$\langle \alpha _p \rangle = 0.01$
, we note that distribution
$B_0$
contains relatively larger peak values for cluster cross-sectional area, however, these peaks are similarly distributed as compared with distribution
$A$
. At higher volume fraction,
$\langle \alpha _p \rangle = 0.10$
, the polydisperse assembly, distribution
$B$
, achieves slightly higher cross-sectional cluster area values as compared with distribution
$B_0$
, which is consistent with the qualitative observation that a few of the polydisperse clusters are aligned more with the cross-stream, than the streamwise direction.
Finally, as detailed in figure 17, we note that the structures of the polydisperse clusters are more fragmented than their monodisperse counterparts. This is partly due to the existence of particles that range from very small to very large particles in the domain but also may point to the formation of larger clusters when two moderately sized clusters merge.
In addition to the qualitative analysis presented above, we also consider an alternative, quantitative method for parsing particle volume fraction information. Here, we bin Eulerian cells by their volume fraction and generate a domain-wise distribution based on this binning. These distributions are shown in figure 18 (for reference, the deviation, skewness and kurtosis for the full assemblies of particles have been previously summarised in table 2). In these figures, the shaded regions are consistent with regions A–D as described previously.
Overview of the clustering patterns observed in distributions
$B$
(left figures) and
$B_0$
(right figures) at
$\langle \alpha _p \rangle = 0.01$
(top row) and
$0.10$
(bottom row).
Cross-sections of particle volume fraction in the
$y{-}z$
plane with contours denoting volume fractions of
$(1.5, 2.25, 3.0)\langle \alpha _p \rangle$
. The colour map represents particle volume fraction and ranges from 0 (white) to
$\alpha _{\text {rcp}}$
(black). Cross-sections are the same as those detailed in figures 15 and 16, but are shown together here to aid in comparisons across the configurations studied at
$\langle \alpha _p \rangle = 0.10$
.
In examining clustering in this way, we note that for
$\langle \alpha _p \rangle = 0.01$
, the monodisperse configurations have distributions of particle volume fraction that are nearly normally distributed (recall skew
$(\alpha _p)$
is 0.66 and 0.82 for distributions
$A_0$
and
$B_0$
), while the skewness for the polydisperse simulations is much greater in comparison. This implies that the polydisperse configurations attain higher local volume fractions for more regions of the flow compared with their monodisperse analogues. In comparing distribution
$A$
and distribution
$B$
we note that distribution
$B$
contains more nearly void regions compared with
$A$
along with a longer tail in the concentrated regime. This implies that Dist
$B$
achieves a more stratified version of clustering, where clusters tend to be denser and more regions of the flow are absent of particles.
Interestingly, at
$\langle \alpha _p \rangle = 0.10$
, the polydispersed and monodispersed assemblies deviate significantly from a normal distribution. Further, the distributions of particle volume fraction are very similar, with the greatest differences occurring between
$A$
and
$A_0$
. This underscores the finding that polydispersity has a much greater effect on clustering behaviour at lower volume fraction. We postulate that when the overall number of particles is higher, as is the case for higher volume fraction flows, this allows for a more uniform introduction of flow disturbances, thereby dampening the effect of larger particles compared with smaller particles.
Distributions of
$\langle \alpha _p \rangle$
for all configurations at the dilute global volume fraction,
$\langle \alpha _p \rangle = 0.01.$
The shaded regions indicate regions of the flow for specific ranges of volume fraction:
$1.5\geqslant \alpha _p/\langle \alpha _p \rangle \lt 2.25$
(lightest grey),
$2.25\geqslant \alpha _p/\langle \alpha _p \rangle \lt 3$
(grey) and
$\alpha _p/\langle \alpha _p \rangle \geqslant 3$
(darkest grey).
Another important distinction we note in our qualitative observations of clustering is the number of clustered regions in the flow. A quantitative way to illustrate this clustering behaviour is to consider the connectivity of the particle volume fraction field. To this end, each Eulerian cell is assigned a binary value corresponding to whether the filtered particle volume fraction is above or below a threshold (corresponding to regions A–D). Cells that have two or more ‘connected’ cells (i.e. either a face, edge or vertex is shared) are considered to be connected and therefore part of a cluster. Connectivity mappings were carried out for all of the configurations at thresholds corresponding to
$\alpha _p = (1.5, 2.25, 3.0) \langle \alpha _p\rangle$
(regions B, C and D). These values are summarised in table 9 and illustrated in figure 19.
In keeping with our qualitative assessment of each configuration, we observe that the monodisperse distributions at
$\langle \alpha _p \rangle = 0.01$
contain fewer connected regions as the volume fraction threshold increases, whereas the polydisperse cases have an increasing number of connected regions. This implies that polydisperse clusters are comprised of a greater number of fragmented dense regions, whereas monodispersed clusters generally have a minority of regions in clusters that attain a dense packing efficiency. In contrast, at
$\langle \alpha _p \rangle =0.10$
, no clear trend exists, however, this is likely due to the fact that both
$A_0$
and
$B_0$
contain relatively very few connected regions in general.
Summary of the number of connected regions containing volume fractions at thresholds of
$(1.5, 2.25, 3.0)\langle \alpha _p \rangle$
. Eulerian cells containing volume fractions above the prescribed threshold with two or connected cells are considered a ‘connected region.’.
Summary of the trend in number of connected regions as a function of the particle volume fraction threshold. To aid in the comparison between configurations, the number of connected regions for each case is normalised by the maximum number of connectivities. Distributions at
$\langle \alpha _p \rangle = 0.01$
are shown on the left and
$\langle \alpha _p \rangle = 0.10$
are shown on the right. Distributions are denoted as:
$A$
(open circles),
$B$
(open diamonds),
$A_0$
(filled circles) and
$B_0$
(filled diamonds).
Finally, we consider the distribution of particle sizes (figure 4) and granular temperature (figure 5) within each of the regions A–D in the flow. Here, each particle is binned into regions according to the local Eulerian particle volume fraction interpolated to its centre. Then, particles belonging to each of the four regions are binned by diameter resulting in a distribution for each group. In each of the panels of figures 4 and 5, the distribution of all the particles in the domain are represented with a solid black line, which is consistent across the plots for each of regions A–D. The normalised distribution of the particles belonging to each of the four regions is denoted with shaded bars.
Considering first the distribution of particle diameters by region, we make several observations. First, smaller particles are preferentially found in the most dilute regime (i.e. region A) and this effect is enhanced for
$\langle \alpha _p \rangle = 0.01$
. Interestingly, for the dilute cases, all of the particles found in region A are less than
$D_{30}$
, though this does not extend to the higher volume fraction cases. In regions B and C, we note similar behaviour for both distributions
$A$
and
$B$
at
$\langle \alpha _p \rangle = 0.01$
. In these regions, there is a dramatic increase in the probability of mid-sized particles. As the density of the cluster increases, there is a much higher likelihood of encountering increasingly large particles and less likelihood of encountering very small particles. This is particularly true for distribution
$A$
, where in region D, almost all the particles are moderately large to very large. In distribution
$B$
we observe that in this region, there is a slight increase in the number of finer particles. This is likely attributable to the fact that distribution
$B$
has a great number of fine particles, and therefore more particles overall. These small particles can be more efficiently packed into a cluster and are also more likely to be entrained in clusters due to the greater number of small particles in the overall distribution.
At higher volume fraction, distribution
$A$
exhibits a smoother increase in preference toward moderate to large-sized particles, particularly in regions B and C. In fact, at high volume fraction, region C is the most likely region to encounter very large particles. In region D, there is a higher likelihood of finding very large particles, however, most notable is the reduction of mid-sized particles. At
$\langle \alpha _p \rangle = 0.10$
, distribution
$B$
also smoothly increases its preference toward larger particles as the local volume fraction increases, however, this increase persists through region D.
Importantly, we believe this preference for larger particles in region D indicates that for polydisperse configurations of particles, it is the larger particles in the domain that give rise to cluster formation, very likely due to the relatively large regions of reduced drag that these particles provide for smaller particles. As previously mentioned, this strong preference for large particles at the centre of clusters with smaller particles haloed around them which is more pronounced for lower volume fraction than at higher volume fraction, is likely observed because, at higher volume fraction, the regions of reduced drag stemming from large particles is obscured by the regions of reduced drag from the greater number of other nearby particles.
Open access
