On the settling and clustering behaviour of polydisperse gas–solid flows

Abstract

Sedimenting flows occur in a range of society-critical systems, such as circulating fluidised bed reactors and pyroclastic density currents (PDCs), the most hazardous volcanic process. In these systems, mass loading is sufficiently high ($\gg \mathcal {O}(1)$) and momentum coupling between the phases gives rise to mesoscale behaviour, such as formation of coherent structures capable of generating and sustaining turbulence in the carrier phase and directly impacting large-scale quantities of interest, such as settling time. While contemporary work has explored the physical processes underpinning these multiphase phenomena for monodispersed particles, polydispersed behaviour has been largely understudied. Since all real-world flows are polydisperse, understanding the role of polydispersity in gas–solid systems is critical for informing closures that are accurate and robust. This work characterises the sedimentation behaviour of two polydispersed gas–solid flows, with properties of the particles sampled from historical PDC ejecta. Highly resolved data at two volume fractions (1 % and 10 %) are collected using an EulerLagrange framework and is compared with monodisperse configurations of particles with diameters equivalent to the arithmetic mean of the polydisperse configurations. From these data, we find that polydispersity has an important impact on cluster formation and structure and that this is most pronounced for dilute flows. At higher volume fraction, the effect of polydispersity is reduced. We also propose a new metric for predicting the degree of clustering, termed ‘surface loading’, and a model for the coefficient of drag that accurately captures the settling velocity observed in the high-fidelity data.

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 2007; Mettler, Vlachos & Dauenhauer 2012), transport of contaminants within ground water (Miller et al. 1998) or the atmosphere (Ravishankara 1997; Arganat & Perminov 2020) and several other environmentally relevant systems, such as volcanic processes (Lube et al. 2020) and avalanches (Sovilla, McElwaine & Köhler 2018; Cicoira et al. 2022; Gardezi et al. 2022; Zhuang et al. 2023). 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 2017; Breard, Dufek & Lube 2018; Sovilla et al. 2018) and computationally (Capecelatro, Desjardins & Fox 2014; Lube et al. 2020), 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 (1987a , 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 1992; Bosse, Kleiser & Meiburg 2006; Breard & Lube 2017; Breard et al. 2018). 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 2015).

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 2019), alter thermal exchange efficiency (Guo & Capecelatro 2019) and entrance lengths (Beetham & Capecelatro 2023) and enhance the mean particle settling velocity (Wang & Maxey 1993; Aliseda et al. 2002; Ferrante & Elghobashi 2003; Yang & Shy 2003 Bosse et al. 2006).

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 2020). 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. 2015). 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. 2014; Beetham, Fox & Capecelatro 2021). 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. 2021).

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 1989; Sundaram & Collins 1994; Cao & Ahmadi 1995; Dasgupta, Jackson & Sundaresan 1998; Zeng & Zhou 2006; Rao et al. 2012), 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 1941) and underscores the shortcomings of extending single-phase turbulence models to multiphase flows. To this end, work by Fox (2014) 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. 2021) 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 2010; Subramaniam 2020) 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 2024), however, several of these studies often were limited to bi-disperse assemblies of particles (Municchi & Radl 2017), 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. 2019).

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. (2023) 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. 2023). 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. 2018).

In recent years, several one-dimensional (1-D) (Degruyter & Bonadonna 2012; Dufek, Esposti Ongaro & Roche 2015; Folch, Costa & Macedonio 2016; Pouget et al. 2016), 2-D (Pfeiffer, Costa & Macedonio 2005; Barsotti, Neri & Scire 2008; Williams, Stinton & Sheridan 2008) and 3-D multiphase flow models (Neri et al. 2007; Dufek & Bergantz 2007b ; Costa 2016; Lube et al. 2020; Cao et al. 2021) 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 1982; Basson, Berres & Bürger 2009; Shankar, Pandey & Shukla 2021). 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. 2020; Breard et al. 2023; de’Michieli Vitturi et al. 2023). To this end, de’Michieli Vitturi et al. (2023) 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 1994) 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. 2016; Breard & Lube 2017; Lube et al. 2020).

Particle density is constant across all configurations and is consistent with the density typically found in the intermediate layer of PDCs (Lube et al. 2015; Breard et al. 2016). Particle diameter distributions are chosen following the data presented in (Lube et al. 2014). 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.

Table 1. 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 2007a ; Lube et al. 2015, 2020; Breard et al. 2016; Breard & Lube 2017; Breard et al. 2018), 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 2002). 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 2016) 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 1967). 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. 2024) 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

(2.1) \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

(2.2) \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

(2.3) \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

(2.4) \begin{align} \frac {\text {d}\boldsymbol {x}_p^{(i)}}{\text {d}t}&=\boldsymbol {u}_p^{(i)} \end{align}

and

(2.5) \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,

(2.6) \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 (2011) that takes into account local volume fraction and Reynolds number effects and is given by

(2.7) \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

(2.8) \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

(2.9) \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

(2.10) \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 1979) 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 (2013). 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

(2.11) \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

(2.12) \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. 2008), 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 2013; Capecelatro et al. 2014, 2015). 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 (2001). 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 (2013).

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. 2016) 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. (2015) and Beetham et al. (2021).

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

(3.1) \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. 2015; Beetham et al. 2021) 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

(3.2) \begin{align} k_p &= \frac {1}{2}\langle \boldsymbol {u}_p^{\prime \prime } \cdot \boldsymbol {u}_p^{\prime \prime } \rangle _p \end{align}

and

(3.3) \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 (2013).

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

(3.4) \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

(3.5) \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

(3.6) \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.

Figure 1. 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 (1994) and defined as

(4.1) \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.

Table 2. 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

(4.2) \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

(4.3) \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. 2015; Beetham et al. 2021), 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

(4.4) \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.

Figure 2. 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

(4.5) \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

(4.6) \begin{align} \mathbb {E}\left (\langle \alpha _p \rangle \right ) &= -8.2\langle \alpha _p \rangle + 0.9 ,\end{align}

(4.7) \begin{align} \mathbb {F}\left (\langle \alpha _p \rangle \right ) &=76.0\langle \alpha _p \rangle - 0.8 ,\end{align}

(4.8) \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. 2021) 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

(4.9) \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.

Table 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,

(4.10) \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.

Figure 3. 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:

(4.11) \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

(4.12) \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

(4.13) \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 (1994). 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.

Table 4. 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.

1. (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.

2. (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.

3. (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.

Figure 4. 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.

Figure 5. 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

(4.14) \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. 2011). 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.

Figure 6. 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.

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

Figure 8. 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$ .

Figure 9. 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. (2023), 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 1994) was proposed. This model is given as

(4.15) \begin{align} \mathcal {V}_s &= \sqrt { \frac {4}{3C_D}g {d}_p \left (\frac {\rho _p-\rho _f}{\rho _f}\right )} ,\end{align}

(4.16) \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.

Table 5. Summary of the mean settling velocity for each configuration compared with Stokes law and the settling law of de’Michieli Vitturi et al. (2023) 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. (2023) (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$

(4.17) \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

(4.18) \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.

Table 6. 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.

Table 7. 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$ .

Figure 10. 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. 2023) 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.

Figure 11. 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 (1994). 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 2001) 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. (2022). The relative contribution to the drag coefficient for each of these terms is demonstrated for the model corresponding to distribution $A$ in figure 11.

Figure 12. Summary of the numerical implementation of the Wiener process (Higham 2001). 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. (2023) 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. (2023).

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.