Predicting the impact of particle-particle collisions on turbophoresis with a reduced number of computational particles

TL;DR

This paper presents a novel method to accurately simulate particle-particle collisions in turbulent flows with fewer computational particles by artificially enlarging collision radii, preserving collision dynamics across different regimes.

Contribution

It introduces a collision treatment that maintains accurate collision statistics with reduced particles by scaling collision radii, enabling efficient simulations of turbophoresis effects.

Findings

The method reproduces full simulation results with fewer particles.

It is effective across low and high Stokes number regimes.

The approach preserves deterministic collision modeling.

Abstract

A common feature of wall-bounded turbulent particle-laden flows is enhanced particle concentrations in a thin layer near the wall due to a phenomenon known as turbophoresis. Even at relatively low bulk volume fractions, particle-particle collisions regulate turbophoresis in a critical way, making simulations sensitive to collisional effects. Lagrangian tracking of every particle in the flow can become computationally expensive when the physical number of particles in the system is large. Artificially reducing the number of particles in the simulation can mitigate the computational cost. When particle-particle collisions are an important aspect determining the simulation outcome, as in the case when turbophoresis plays an active role, simply reducing the number of particles in the simulation significantly alters the computed particle statistics. This paper introduces a computational particle treatment for particle-particle collisions which reproduces the results of a full simulation with a reduced number of particles. This is accomplished by artificially enhancing the particle collision radius based on scaling laws for the collision rates. The proposed method retains the use of deterministic collision models and is applicable for both low and high Stokes number regimes.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}=-\frac{1}{{\rho }_f}\nabla p+{\nu }_f{\nabla }^2\mathbf{u}+\mathbf{f},\nabla \cdot \mathbf{u}=0,$$

$$\frac{d{\mathbf{x}}^{(i)}}{dt}={\mathbf{v}}^{(i)},\frac{d{\mathbf{v}}^{(i)}}{dt}=\mathbf{a}\left({\mathbf{x}}^{(i)},t,{\mathbf{v}}^{(i)}\right).$$

$$S{t}^{+}=\frac{1}{18}\frac{{\rho }_p}{{\rho }_f}{d_p^{+}}^2.$$

$$f\left(\mathbf{x},\mathbf{v};t\right)=⟨\sum _{i=1}^{N_p}\delta \left(\mathbf{x}-{\mathbf{x}}^{(i)}(t)\right)\delta \left(\mathbf{v}-{\mathbf{v}}^{(i)}(t)\right)⟩,$$

$$\frac{\partial f}{\partial t}+{\nabla }_{\mathbf{x}}\cdot \left(\mathbf{v}f\right)+{\nabla }_{\mathbf{v}}\cdot \left(⟨\mathbf{a}|\mathbf{x},\mathbf{v}⟩f\right)={\dot{f}}_{coll},$$

$$C(\mathbf{x};t)=\int f(\mathbf{x},\mathbf{v};t)d\mathbf{v}.$$

$$C_0=\frac{N_p}{V}=\frac{1}{V}\int C(\mathbf{x};t)d\mathbf{x}=\frac{1}{V}\iint f(\mathbf{x},\mathbf{v};t)d\mathbf{v}d\mathbf{x},$$

$${\Phi }_V=\frac{\pi }{6}{d}_p^3{C}_0,{\Phi }_M=\frac{{\rho }_p}{{\rho }_f}{\Phi }_V.$$

$${\varphi }_v=\frac{\pi }{6}{d}_p^{3}C,{\varphi }_m=\frac{{\rho }_p}{{\rho }_f}{\varphi }_v.$$

$$\frac{\partial C}{\partial t}+{\nabla }_{\mathbf{x}}\cdot \left(⟨\mathbf{v}|\mathbf{x}⟩C\right)=0.$$

$$\frac{d\left(⟨v_y|y⟩C\right)}{dy}=0,$$

$$\frac{\partial \left(⟨\mathbf{v}|\mathbf{x}⟩C\right)}{\partial t}+{\nabla }_{\mathbf{x}}\cdot \left(⟨\mathbf{v}\otimes \mathbf{v}|\mathbf{x}⟩C\right)-⟨\mathbf{a}|\mathbf{x}⟩C=0$$

$$\frac{d\left(⟨v_y^2|y⟩C\right)}{dy}-⟨a_y|y⟩C=0,$$

$$⟨v_y^2|y⟩\frac{dC}{dy}=\left(⟨a_y|y⟩-\frac{d⟨v_y^2|y⟩}{dy}\right)C.$$

$$\begin{gathered} C(y)=\mathcal{N}\exp \left({\int }_{}^y\frac{⟨a_y|\eta ⟩}{⟨v_y^2|y⟩}d\eta -{\int }_{}^y\frac{d⟨v_y^2|\eta ⟩}{d\eta }d\eta \right) \\ =\frac{\mathcal{N}}{⟨v_y^2|\eta ⟩}\exp \left({\int }_{}^y\frac{⟨a_y|\eta ⟩}{⟨v_y^2|y⟩}d\eta \right), \end{gathered}$$

$$I_{bias}={\int }_{}^y\frac{⟨a_y|\eta ⟩}{⟨v_y^2|y⟩}d\eta ,I_{turb}={\int }_{}^y\frac{d⟨v_y^2|\eta ⟩}{d\eta }d\eta .$$

$$g\left(\mathbf{x},\mathbf{v};t\right)=⟨\sum _{i=1}^{N_c}{W}^{(i)}\delta \left(\mathbf{x}-{\mathbf{x}}^{(i)}(t)\right)\delta \left(\mathbf{v}-{\mathbf{v}}^{(i)}(t)\right)⟩,$$

$$C(\mathbf{x};t)=\frac{1}{V}\int g(\mathbf{x},\mathbf{v};t)d\mathbf{v},$$

$$C_0(t)=\frac{N_p}{V}=\int C(\mathbf{x};t)d\mathbf{x}=\frac{N_{c}W}{V},$$

$$\frac{\partial g}{\partial t}+{\nabla }_{\mathbf{x}}\cdot \left(\mathbf{v}g\right)+{\nabla }_{\mathbf{v}}\cdot \left(⟨\mathbf{a}|\mathbf{x},\mathbf{v}⟩g\right)={\dot{g}}_{coll}.$$

$${\dot{f}}_{coll}={\iint }_0^{d_p}bdbd\theta \int d{\mathbf{v}}_1|{\mathbf{v}}_1-\mathbf{v}|\left[f({\mathbf{v}}^{\prime })f({\mathbf{v}}_1^{\prime })-f(\mathbf{v})f({\mathbf{v}}_1)\right].$$

$${\dot{f}}_{coll}\sim \delta v_p{d}_p^2{N}_p^2,$$

$${\dot{g}}_{coll}={\iint }_0^{d_c}bdbd\theta \int d{\mathbf{v}}_1|{\mathbf{v}}_1-\mathbf{v}|\left[g({\mathbf{v}}^{\prime })\frac{g({\mathbf{v}}_1^{\prime })}{W}-g(\mathbf{v})\frac{g({\mathbf{v}}_1)}{W}\right],$$

$${\dot{g}}_{coll}\sim \delta v_c{d}_c^2{N}_p^2/W.$$

$$v_c\sim d_c^h.$$

$$\underset{{\tau }_p\to 0}{\lim }h=1,\underset{{\tau }_p\to \infty }{\lim }h=0.$$

$${\dot{g}}_{coll}\sim d_c^{2+h}{N}_p^2/W,$$

$$d_c=d_p{W}^{1/(2+h)},$$

$$\underset{{\tau }_p\to 0}{\lim }\frac{d_c}{d_p}=\sqrt[3]{W},\underset{{\tau }_p\to \infty }{\lim }\frac{d_c}{d_p}=\sqrt{W}.$$

$$\delta v_c=\delta v_p{W}^{h/(2+h)}.$$

Full text

1 Introduction

Encounters with particle-laden flows are common in both engineering and natural sciences. Because of the wide range of length and time scales typically involved in fluid flows with dispersed particles, it is often infeasible to directly simulate these flows with the sub-particle scale resolution necessary to capture the particle-fluid interface and apply the relevant boundary or interface conditions (Burton and Eaton, 2005; Ardekani et al., 2016; Sierakowski and Prosperetti, 2016; Horne and Mahesh, 2019). On the other hand, treating the collection of particles as a continuum provides a vast reduction in computational cost, but such an approximation may not always be appropriate for accurate simulations (Druzhinin and Elghobashi, 1998; Fox et al., 2008). In between these two approaches, the popular Eulerian-Lagrangian method employs a ‘point-particle’ or ‘discrete element’ description of individual particles as they are tracked along their trajectories through the flow simulated on fixed grid (Squires and Eaton, 1990; Capecelatro and Desjardins, 2013; Subramaniam, 2013). Each particle is associated with a number of computational degrees of freedom describing its location, velocity, and other relevant properties. The interaction of a particle with the surrounding flow or other particles is not resolved (by the Eulerian grid), so closures are needed in the form of drag laws, collision models, and other dynamical relationships for integrating ordinary differential equations describing the evolution of the particle (Maxey and Riley, 1983; Gatignol, 1983). The computational cost grows (at least) linearly with the number of particles, and can grow faster in cases where identifying and computing particle-particle collisions is important.

Given that systems with billions, trillions, and more particles are not uncommon, Eulerian-Lagrangian approach to particle-laden flows can be inconveniently, or even prohibitively, expensive. A reduction in cost may be achieved by having a single computational particle ‘represent’ many particles in the physical system (Elghobashi and Truesdell, 1993; Elghobashi, 1994). Generally, the goal of such methods to have the statistics of the computational particles accurately reproduce the particle statistics of the case where the physical number of particles is simulated. Garg et al. (2009) described this as statistical equivalence and proscribed a method for dynamically assigning statistical weights to computational particles. In the absence of two-way coupling or particle-particle collisions, a minimal constraint is that a sufficient number of particles should be used to obtain enough samples for computing statistical quantities of interest (Elghobashi and Truesdell, 1992, 1993). Garg et al. (2009) used evolving statistical weights with annihilation and cloning of computational particles to reduce the non-uniformity in the computational particle number density, thus making statistical errors relatively uniform in space. This scheme also helps with load balancing for parallelization across many processors. The physics-based necessity of including two-way coupling at high mass fractions and four-way coupling at high volume fractions leads to additional complexities. The two-way coupling force for computational particles must be amplified in proportion to the number of physical particles represented (Garg et al., 2009). This can exacerbate the induced fluid velocity at the particle location and highlight the need for estimating an undisturbed flud velocity (Horwitz and Mani, 2016, 2018; Horwitz, 2019). Moin and Apte (2006) presented a scheme to prevent the number of computational particles from exceeding a specified threshold during liquid jet breakup, without considering the impact on collisions or coalescence. In this work, the treatment of collisional effects is isolated for study.

The present work is motivated by the impact of turbophoresis and non-uniform particle concentrations on radiation transport and absorption in internal wall-bounded turbulent flows (Frankel, 2017; Banko, 2018). Turbophoresis is a phenomenon in which small, heavy particles experience a net drift from regions of high turbulence intensity toward lower intensity regions in a turbulent shear flow (Caporaloni et al., 1975; Reeks, 1983; Guha, 2008). In wall-bounded flows, the fluid velocity vanishes at the wall due to the no-slip and no-penetration boundary conditions. Within a thin layer of flow adjacent to the wall, therefore, the turbulent fluctuations are less intense. Particles in this region near the wall tend to have lower velocity fluctuations themselves compared with the bulk of the flow away from the wall, leading to a differential in dispersion rates depending on distance from the boundary. Particles are rapidly dispersed throughout most of the turbulent flow but are ‘trapped’ near the wall for longer durations. This results in higher particle concentrations in a thin layer adjacent to the wall and interactions between particles and near wall coherent structures become important (Marchioli and Soldati, 2002; Sikovsky, 2014).

Particle-particle collisions have a strong impact on turbophoresis, signficantly reducing the non-uniformity of particle concentration profiles as the volume fraction is increased (Li et al., 2001). The increase in particle concentration levels leads to higher collision rates in the near wall region (Kuerten, 2016). Particles in low-speed and high-speed streaks experience elevated relative velocities. Oblique collisions redistribute particle streamwise velocity fluctuations into the wall-normal direction (Caraman et al., 2003). Thus, the wall-normal mixing of particles near the wall is greatly enhanced, flattening the concentration profiles even at relatively low bulk volume fractions $\sim 10^{-4}$ (Yamamoto et al., 2001; Kuerten and Vreman, 2015). The quadratic scaling of collision rate with particle number density creates the scenario in which the particle concentration profile, as a function of wall-normal distance, is very sensitive to the bulk volume fraction of the particle-laden flow (Johnson, 2018; Johnson et al., 2019). Any approach which relies on reducing computational cost by decreasing the number of simulated particles will necessarily struggle to provide an accurate concentration field without additional treatment of particle-particle collisions to account for this effect. Because these collisions occur in a very thin layer near the wall, where the flow and particle motions are highly anisotropic, and the obliqueness of collisions is important, stochastic approaches for collisions are less promising.

The rate at which collisions between two particles occur scales quadratically with the local number density of particles (Wang et al., 2000; Reade and Collins, 2000), so simply reducing the number of particles and assigning a statistical weight to each particle leads to under-prediction of collisional effects. Others approaches have considered stochastic models for particle-particle collisions, but this increases the difficulty and imprecision of the collision treatment (Schmidt and Rutland, 2000; O’Rourke et al., 2009). Sommerfeld (2001) introduced an approach based on the creation of fictitious collision partners with random obliqueness and relative velocity specified by observed trends in homogeneous turbulence. It is not clear how universal such a treatment would be, particularly in highly anisotropic and inhomogeneous regions near the wall in a turbulent flow.

In this paper, a novel method is introduced to address particle-particle collisions within the general computational particle framework of Garg et al. (2009). The deterministic treatment of collisions is retained. The method consists in artificially enhancing the collisional radius of the particles so as to maintain accurate collision rates as the number of computational particles is reduced. The collisional radius determining the instance and outcome of a collision between two particles. Scaling arguments for the collision rate can be used to determine how the altered collision radius depends on its statistical weight. The scaling may be directly determined in the opposite limits of low and high Stokes numbers, and a hybrid model is proposed for bridging between these limits. The resulting approach dynamically computes the collision radius for each individual collision based on the properties of the particle and fluid involved.

The rest of the paper is organized as follows. The sensitivity of the particle concentration field to particle-particle collisions in a canonical flow is considered in detail in §2. In §3, the theoretical basis for reducing the number of computational particles is reviewed and extended to include collisional effects by introducing scaling relations for an enhanced collisional radius treatment. Simulation results in §4 for particle-laden turbulent channel flow demonstrate that such scaling arguments work in the low and high Stokes number limits. Further, a hybrid model bridging these limits is constructed and demonstrated to work well in both regimes. Finally, §5 provides a summary and concluding discussion.

2 Turbophoresis and particle-particle collisions

This section demonstrates the difficulties associated with reducing the number of computational particles in wall-bounded turbulent flow simulations. The basic physics of turbophoresis are reviewed (§2.3) in a statistical framework (§2.2) and the impact of particle-particle collisions is demonstrated (§2.4), with particular attention to the particle concentration profile near the wall. To set the stage, a fully-developed particle-laden turbulent channel flow is introduced in §2.1 as the canonical flow for the current study. The physical model and numerical approach, which combines a Lagrangian description of the particle dynamics with direct numerical simulation (DNS) of the fluid on an Eulerian grid, are also described in §2.1.

2.1 Canonical flow: particle-laden turbulent channel

The turbulent channel flow is adopted to demonstrate the basic physical principles of turbophoresis and the impact of particle-particle collisions. The flow is bounded only in the $y$ (wall-normal) direction by two smooth, flat walls separated by a distance $2h$, where $h$ denotes the half-height of the channel. A pressure gradient, $-dp/dx=\rho_{f}u_{*}^{2}/h$, forces the mean flow in the $x$ (streamwise) direction. The statistics of the turbulent flow and the particles embedded in the flow are homogeneous with respect to the $x$ and $z$ (spanwise) directions. Here, $\rho_{f}$ is the fluid mass density and $u_{*}$ is the friction velocity which characterize the typical turbulent velocity fluctuation magnitude. The friction Reynolds number, $Re_{*}=u_{*}h/\nu_{f}$, where $\nu_{f}$ is the kinematic viscosity of the fluid, characterizes the range of active turbulent length scales between the geometry-imposed length scale, $h$, and viscous length scale, $\delta_{*}=\nu_{f}/u_{*}$. It also encapulates the range of active timescales in the turbulent fluctuations, from the large-eddy turnover time, $h/u_{*}$, to the viscous timescale, $\tau_{*}=\nu_{f}/u_{*}^{2}$.

The numerical simulations used in this work employ periodic boundary conditions with domain size $L_{x}=4\pi h$ and $L_{z}=2\pi h$. The flow is assumed incompressible and is described by a velocity field, $\mathbf{u}(\mathbf{x},t)$, and a pressure field, $p(\mathbf{x},t)$, which evolve according to,

[TABLE]

where and $f_{i}$ represents the mean pressure gradient forcing as well as any two-way coupling force from the particles. Numerical simulations for this paper solve Eq. (1) on a staggered Cartesian grid with uniform spacing $\Delta x^{+}=\Delta x/\delta_{*}\approx 11$, $\Delta z^{+}=\Delta z/\delta_{*}\approx 7$ in the streamwise and spanwise directions, respectively. The grid is non-uniform in the wall-normal direction with $\Delta y_{min}^{+}\approx 0.5$ adjacent to the wall and $\Delta y_{max}^{+}\approx 7$ at the center of the channel. The discretization of spatial derivatives uses second-order central differencing and a second-order Runge-Kutta scheme is used in a fractional step method for time advancement.

Each particle is described by its center of mass, $\mathbf{x}(t)$, and velocity, $\mathbf{v}(t)$, which evolve according to,

[TABLE]

where $\mathbf{a}$ represents acceleration due to fluid forces. The Stokes drag force, $\mathbf{a}=\mathbf{a}_{St}=\tau_{p}^{-1}\left(\mathbf{u}(\mathbf{x}^{(i)},t)-\mathbf{v}^{(i)}\right)$, applies when the particle Reynolds number is very small, $Re_{p}=|\mathbf{u}-\mathbf{v}|d_{p}/\nu_{f}\ll 1$. The particle adjusts to the background flow on the timescale, $\tau_{p}=\rho_{p}d_{p}^{2}/(18\mu_{f})$. A better model at finite Reynolds number is the Schiller-Naumann drag (Schiller and Naumann, 1933; Balachandar and Eaton, 2010), $\mathbf{a}=\mathbf{a}_{SN}=\mathbf{a}_{St}\left(1+0.15Re_{p}^{0.687}\right)$, which is used in all the simulations in this paper. Particles are characterized by the Stokes number, $St^{+}=\tau_{p}/\tau_{*}$, particle-to-fluid density ratio, $\rho_{p}/\rho_{f}$, and non-dimensional diameter, $d_{p}^{+}=d_{p}/\delta_{*}$. Only two of these dimensionless parameters are independent,

[TABLE]

The particle phase is computed by advancing (2) with a second-order Runge-Kutta method, using trilinear interpolation to compute the fluid velocity at particle locations. Higher-order Lagrange interpolation schemes are tested and resulted in little change to the quantities of interest in this paper (Johnson et al., 2019). Particle-wall collisions are computed using a unity restitution coefficient. Particle-particle collisions were computed using a hard-sphere collision model with a specified restitution coefficient. Unless otherwise stated, particle-particle collisions also used a unity restitution coefficient. Particle and collision statistics shown in this paper are computed in wall-normal discretized bins of width $0.5$ viscous units, which equal to the minimum grid spacing of the Eulerian fluid solution.

2.2 Particle statistics

The (single-particle) statistics of monodisperse particles are described by the particle distribution function,

[TABLE]

where $\delta(\mathbf{x})$ is the Dirac delta function and $N_{p}$ is the number of particles in the domain. In this work, $N_{p}$ remains constant in time, and the numerical simulations reintroduce particles exiting the domain at the complementary location along the opposite boundary, consistent with the periodic boundary conditions. The evolution of the particle statistics is governed by,

[TABLE]

where $\dot{f}_{coll}=\dot{f}_{coll}(\mathbf{x},\mathbf{v};t)$ represents the change in particle velocities caused by particle-particle collisions.

The particle concentration is defined by

[TABLE]

The bulk concentration of particles is

[TABLE]

where $V$ is the volume of the domain. The volume fraction and mass fraction characterizing the particle ensemble are, respectively,

[TABLE]

Local volume and mass fractions can also be defined,

[TABLE]

Integrating (5) in velocity and dividing by volume gives the evolution equation for the concentration field,

[TABLE]

The collisional term vanishes because individual collisions do not add or remove particles (mass conserving). At steady state in the channel flow, this implies,

[TABLE]

which means that $\langle v_{y}|y\rangle=0$.

2.3 Turbophoresis

Likewise, an evolution equation expressing the conservation of momentum for the particle phase can be obtained by multiplying (5) by $\mathbf{v}$ before integrating in velocity,

[TABLE]

For the steady-state channel flow, this implies,

[TABLE]

or

[TABLE]

The two right side terms represent two causes of non-uniform concentration profiles in the turbulent channel flow. The first term is the average wall-normal force on all particles at a particular distance from the wall. Using Stokes drag to illustrate the basic idea, $\langle a_{y}|y\rangle=\tau_{p}^{-1}\langle u_{y}|y\rangle$ (recall that $\langle v_{y}|y\rangle=0$ at steady-state), where $\langle u_{y}|y\rangle>0$ as inertial particles tend to preferentially sample ejection events (Johnson, 2018; Johnson et al., 2019). This biased sampling of the flow leads to a net force on particles away from the wall. The biased sampling force counteracts the turbophoresis pseudo-force, the second term on the right, which pushes particles toward the wall, down the gradient of particle wall-normal velocity variance, $\tfrac{d}{dy}\langle v_{y}^{2}|y\rangle$. Equation (14) can be formally solved for the concentration,

[TABLE]

where $\mathcal{N}$ is an integration constant for normalizing the integral of $C/C_{0}$. The concentration profile can thus be written in terms of phoresis integrals for biased sampling and turbophoresis,

[TABLE]

These phoresis integrals allow a direct evaluation of turbophoresis and biased sampling in determining the concentration profile.

Figure 1(a) shows simulation results for the concentration profiles for a wide range of $St^{+}$ (holding $d_{p}^{+}=0.5$ constant and changing the density ratio). The simulations shown here neglect two-way coupling and particle-particle collisions, effectively representing the case of $\Phi_{V}=0$. The $St^{+}=0$ particles (Lagrangian particles) have a uniform distribution, though a slight deviation in the viscous sublayer illustrates the minor effects of discretization and interpolation errors. At finite but small $St^{+}=2$, turbophoresis begins to raise the concentration in the near-wall region up to $10$ times the bulk concentration. As $St^{+}$ is further increased, the near-wall concentration exceeds $100$ times the mean concentration.

The change in concentration profile with $St^{+}$ is analyzed in Figure 1(b) using phoresis integrals, Eq. (16). At $St^{+}=2$, the biased sampling force is comparable in magnitude to turbophoresis, attenuating the large near-wall concentrations seen at higher $St^{+}$. In fact, it can be shown analytically that the biased sampling term exactly cancels turbophoresis in the $St^{+}\rightarrow 0$ limit (Johnson, 2018; Johnson et al., 2019). As $St^{+}$ increases, both turbophoresis and biased sampling decrease in magnitude, though turbophoresis decreases much more slowly. The result is the predominance of turbophoresis over biased sampling at higher $St^{+}$, leading to $C(y)\sim 1/\langle v_{y}^{2}|y\rangle$, Eq. (15).

2.4 Particle-particle collisions

Particle-particle collisions are known to effect particle dynamics in wall-bounded turbulent flows as volume fractions become high enough (Li et al., 2001; Yamamoto et al., 2001; Caraman et al., 2003; Kuerten and Vreman, 2015, 2016; Johnson et al., 2019). In fact, the volume fraction at which collisions have a significant effect on the concentration profile can be quite ‘low’, as shown in Figure 2. Relative concentration profiles for three different $8\leq St^{+}\leq 128$ are shown in the four panels, each one containing a range of $\Phi_{V}$. As the $\Phi_{V}$ is increased, the near-wall peak in concentration is attenuated. At high enough $\Phi_{V}$, the concentration profile becomes almost flat. The volume fractions shown for each $St^{+}$ value in Figure 2 and 3 were chosen to best capture the transition away from the $\Phi_{V}=0$ concentration profile.

These changes in the relative concentration profiles with $\Phi_{V}$ reflect how particle-particle collisions alter turbophoresis. Indeed, the main effect of increasing $\Phi_{V}$ is to increase the rate at which a given particle collides with another particle. Each collision redistributes the participating particles in velocity-space, impacting $\langle v_{y}^{2}|y\rangle$ and hence changing how the turbophoresis forms the concentration profile. On average, particle-particle collisions near the wall transfer particle fluctuation energy from the streamwise to wall-normal components (return-to-isotropy for particle velocity fluctuations), enhancing the near-wall $\langle v_{y}^{2}|y\rangle$ and decreasing the turbophoretic drift (Caraman et al., 2003; Johnson et al., 2019). The resulting profiles of $\langle v_{y}^{2}|y\rangle$ are shown in Figure 3 for the same $St^{+}$ and $\Phi_{V}$ in Figure 2. At low $\Phi_{V}$, the particle velocities near the wall follow closely the asymptotic behavior of fluid velocities, i.e., $v_{y}\sim y^{2}$. This is particularly true at lower $St^{+}$. As $\Phi_{V}$ increases, the near-wall $\langle v_{y}^{2}|y\rangle$ increases dramatically due to particle-particle collisions, resulting in the reduction of the near-wall concentration peak seen in Figure 2.

It can be seen from Figures 2 and 3 that collisions alter the relative concentration profile at quite low $\Phi_{V}$. While the trend with $\Phi_{V}$ is clearly demonstrated, a more detailed view of the relative importance of collisions can consider the local volume fraction, $\phi_{v}$. Figure 4 plots $\phi_{v}$ as function of wall-normal distance. For $St^{+}\geq 32$, the near-wall peak volume fraction appears to saturate even as the bulk volume fraction, $\Phi_{V}$, increases. This likely also occurs for lower $St^{+}$, but the simulations for this paper did not reach high enough $\Phi_{V}$ to observe it. As the volume (and mass) fraction continues to increase, at some point two-way coupling could alter this observation, but the two-way coupling has a negligible impact on the concentration profile for the simulation parameters shown in this paper (Johnson et al., 2019). The saturation of local volume fraction near the wall suggests that particle-particle collisions limit turbophoresis in such a way as to effectively set a maximum local volume fraction above which turbophoresis cannot sustain the near-wall peak. This maximum local volume fraction, $\phi_{v,max}$, evidently decreases with increasing $St^{+}$. A simple explanation for this trend is as follows.

Roughly speaking, $\phi_{v,max}$ may be related to the volume fraction as which collisional effects begin to dominate over drag in determining the trajectories and statistics of particles, e.g., $\langle v_{y}^{2}|y\rangle$. The relative impact of collisions and drag can be quantified by the dimensionless parameter, $c=\tau_{p}/\tau_{coll}$, where $\tau_{coll}$ is the typical time between collisions for a given particle. If $c$ is small, particles have plenty of time in between collisions to adjust back to the velocity of the nearby fluid. If $c$ is large, however, particles behave ballistically between collisions and fluid drag play a minimal role in establishing particle statistics such as $\langle v_{y}^{2}|y\rangle$. It follows, then, that the maximum (local) volume fraction that turbophoresis (alone) can support can be approximated by a critical value of this parameter, $c=c_{cr}\sim\mathcal{O}(1)$. For wall-bounded flows, $\tau_{coll}\sim\lambda_{mf}/u_{*}$, where $\lambda_{mf}$ is the mean-free-path of the particles and the friction velocity, $u_{*}$, establishes the typical relative velocity between a particle that has just collided and other nearby particles. If the swept volume between collisions, $\sim d_{p}^{2}\lambda_{mf}$, is set equal to the typical volume containing one particle on average, $\sim d_{p}^{3}/\phi_{v}$, then the mean-free-path may be estimated by $\lambda_{mf}\sim d_{p}/\phi_{v}$. Combining these approximations, $c\sim St^{+}\phi_{v}/d_{p}^{+}\sim d_{p}^{+}\phi_{m}$. In this way, it can be argued that the maximum volume fraction for significant collisional effects scales as $\phi_{v,max}\sim d_{p}^{+}/St^{+}$. This scaling is roughly observed in Figure 2, though certainly a more nuanced picture could be gained from more detailed analysis which is not pursued here. As a result of this trend, that $\phi_{v,max}$ decreases with $St^{+}$, collisional effects become noticeable in the relative concentration profiles at lower $\Phi_{V}$ as $St^{+}$ increases. As a result, the deviation from $\Phi_{V}=0$ behavior in Figures 2 and 3 begins at lower volume fractions as $St^{+}$ increases. This qualitative trend was recognized by Elghobashi (1994) when updating the regime diagram from Elghobashi (1991) to include four-way coupling at lower volume fractions as the particle Stokes number increases.

3 Reducing the number of computational particles: ‘super-particles’

In the previous section, it was shown how turbophoresis generates a near-wall peak in particle concentration. In particular, the impact of particle-particle collisions on this process was explored, showing that the relative concentration profile in a wall-bounded turbulent flow can depend strongly on the volume fraction. Therefore, reducing the number of computational particles in the simulation will lead to unphysical changes to the particle concentration field if the collision treatment is not modified. This section introduces a simple method for recovering collisional effects without resorting to stochastic collision models.

The statistical framework of §2.2 is modified in §3.1 to address the situation where a reduced number of computational particles are simulated, with each computational particle carrying a statistical weight signifying that it effectively represents multiple physical particles. Here, these computational particles are referred to as ‘super-particles’. In other contexts, they are called particle ‘parcels’ (Subramaniam, 2013). In §3.2, collisions between super-particles are considered, leading to proposed scaling laws for enhancing the collision radius of super-particles in §3.3 in the limit of low and high Stokes numbers.

3.1 Super-particle statistics

A super-particle denotes a computational particle which does not necessarily represent one physical particle, but represents more than one physical particle through the concept of a statistical weight assigned to each particle. In particle-laden flow simulations with Lagrangian tracking, Eq. (2), the computational cost of simulating the particle phase scales with $N_{p}$, the number of physical particles in the system of interest. The super-particle approach instead tracks $N_{c}$ computational particles, whose statistics are meant to be representative of the full physical system by assigning each super-particle a statistical weight of $W^{(i)}$. In some sense, the super-particle approach for the particle phase is akin to large-eddy simulations (LES) for the fluid phase, in that a description of the system with a reduced number of computational degrees of freedom is sought.

The statistics of the super-particles are described by,

[TABLE]

which reduces to (4) when $W^{(i)}=1$ and $N_{c}=N_{p}$. The present work follows and extends the theoretical approach of Garg et al. (2009) and Subramaniam (2013), which used ‘$h$’ to denote this super-particle distribution function – a notation which is not followed here for reasons that will become obvious later. For now, the statistical weight attached to each particle, $W^{(i)}$, will be kept constant. Further extension along the lines of Garg et al. (2009) to vary $W^{(i)}$ in order to achieve faster numerical convergence is possible but not pursued at the present. All single-particle statistics can be written in terms of $g$. For example, the concentration is,

[TABLE]

and

[TABLE]

which fixes the correct choice of $W=N_{p}/N_{c}$. Thus, $W$ represents the factor by which the number of computational particles is reduced. The volume and mass fractions can be defined using Eqs. (8) and (9).

The trajectory of each individual super-particle is described by Eq. (2). The evolution of the super-particle distribution function is,

[TABLE]

In order to achieve accurate statistics from the super-particle approach, the distribution functions $f$ and $g$ should have the same evolution for the same initial and boundary conditions. This is the principle of statistical equivalence (Subramaniam, 2013). Another way to state statistical equivalence is to say that the evolution of $g$ should be invariant with respect to changes in $W$. It is immediately clear, then, that the drag law used for super-particles, $\mathbf{a}$, must be kept invariant with respect to $W$. That is, the super-particles must use the same drag law as physical particles.

It is evident from Figures 2 and 3, along with the discussion of the underlying physics, that the particle-particle collisions represented by $\dot{g}_{coll}$ break the invariance of (20) with respect to $W$. That is, if a reduced number of computational particles are simulated with the same drag law and same collision model, the result is as if a lower volume fraction is simulated. In this way, the variation of relative concentration, $C(y)/C_{0}$, and particle wall-normal velocity variance, $\langle v_{y}^{2}|y\rangle$, shown in Figures 2 and 3 respectively may be interpreted as errors of a naive super-particle approach. A closer look at $\dot{g}_{coll}$ and how its breaks invariance with respect to $W$ is warranted.

3.2 Super-particle collisions

The analysis pursued below treats only binary collisions. Collisions involving three or more particles are assumed to be rare enough so as to be negligible in their impact on the quantities of interest. A sketch of a representative binary collision is shown in Figure 5. A reference particle with pre-collision velocity $\mathbf{v}$ is impacted by a colliding particle with pre-collision velocity $\mathbf{v}_{1}$. The particles collide with an offset $b<d_{p}$, which determines the obliqueness of the collision. After the collision, the particles have velocities $\mathbf{v}^{\prime}$ and $\mathbf{v}_{1}^{\prime}$ respectively. The post-collision velocities depend on the restitution of coefficient. The collisional source term in Eq. (5) may be modeled as,

[TABLE]

This form is a bit simplified, for example, in assuming that the velocities of the colliding particles are statistically independent. Removing this assumption requires knowledge of two-particle statistics, a richer description than provided by Eq. (4) or (17). Nonetheless, Eq. (21) illustrates the basic scaling of $\dot{f}_{coll}$, which is all that is used for the modeling in this paper. In particular,

[TABLE]

where $\delta v_{p}$ is the the relative velocity at collision, i.e., relative velocity when the particles center of masses are separated by a distance $d_{p}$. More generally, it can be seen that $\dot{f}_{coll}$ scales with the collision rate, i.e., the number of collisions per unit time. Indeed, previous estimates of the collision rate (Sundaram and Collins, 1997; Wang et al., 2000) lead to Eq. (22) as well, though may differ from each other in whether the relative velocity magnitude or relative radial velocity is used for $\delta v_{p}$.

In the same way, the collisional source term in Eq. (20) may be modeled as,

[TABLE]

or

[TABLE]

In other words, if the collisional diameter of the super-particles, $d_{c}$, is kept equal to the physical diameter of the particles, $d_{p}$, and the relative velocity statistics are assumed unchanged, the collision rate will be under-estimated proportional to $W^{-1}$. The underlying reason is as follows. The single-particle distribution, Eq. (4), scales linearly with the number of particles, so that particle statistics can be represented by a reduced number of particles using the (linear) correction factor $W=N_{p}/N_{c}$. The binary collision rate, and hence $\dot{f}_{coll}$, scales quadratically with the number of particles. Therefore, the linear correction factor $W$ in Eq. (17) only corrects for one of these two factors.

Said another way, each super-particle statistically represents $W$ physical particles, and hence when it does collide with another particle, this effectively represents the effect of $W$ particle-particle collisions. However, all other features held constant (relative velocity statistics, preferential concentration, etc…), the collision rate falls by a factor of $W^{-2}$ in the super-particle system. Therefore, the net effect of collisions on the resulting simulation statistics falls by a factor of $W^{-1}$, Eq. (24). Without further treatment, the effect of collisions will not be represented accurately in the super-particle approach leading to errors in, for instance, concentration (Figure 2) and velocity (Figure 3).

3.3 Enhanced collision radius model

A simple approach to modeling collisional effects in a super-particle approach is to try to restore invariance of $\dot{g}_{coll}$ with respect to $W$. One way this can be done is by enhancing the collisional diameter of the particles to be larger than the physical diameter of the particles, $d_{c}>d_{p}$. Because the velocity at collision, $\delta v_{c}$, may depend on the collisional diameter, a Hölder exponent-like relation is assumed for the velocity difference between two particles separated by a distance $d_{c}$,

[TABLE]

In the limit of very low Stokes number, the perturbation solution of Maxey (1987) may be used to approximate the relative velocity of the particles in the vicinity of collision. Here the effects of two-way coupling, particle screening, and lubrication are not explicitly considered. In this approximation, the particle velocities are represented as a smooth, single-valued field, so that $h=1$ in Eq. (25). In the opposite limit of high Stokes number, the colliding particles do not adjust at all to the local fluid velocity field as they approach. As a result, the collision velocity is independent of their separation at collision and $h=1$ for (25). These two limits may be summarized as,

[TABLE]

With this approximation for the collision velocity, the super-particle collisional source term scales as,

[TABLE]

and therefore the invariance of $\dot{g}_{coll}$ with respect to $W$ may be upheld by the choice,

[TABLE]

which in the two limits discussed above results in,

[TABLE]

As an example, consider a super-particle simulation with $64$ times fewer particles than the physical system under investigation. If the particles are at very low Stokes numbers, then the analysis here suggests scaling the collisional diameter by a factor of $4$ and proceeding with a deterministic collision algorithm. For very high Stokes number particles, this analysis suggests that more accurate results could be obtained by scaling the collision diameter by a factor of $8$.

The scaling for the collisional diameter, Eq. (28) along with Eq. (25) implies a scaling for the collisional velocity,

[TABLE]

Because the change of momentum of particles due to collision is proportional to the relative (radial) velocity at collision, this ratio can be factored into the coefficient of restitution used for the collision of super-particles,

[TABLE]

where $e_{p}$ is the physical coefficient of restitution. The proposed model, then, consists of computing collisions based on an enhance collisional diameter according to Eq. (28) with a modified restitution coefficient from Eq. (31). The model has a single parameter $h$ which must be chosen based on the expected behavior of the relative velocities of collisions (smooth vs. ballistic). While the collision radius is increased for particle-particle collisions, it is important to note that the treatment of particle-wall collisions continues to use the physical diameter of the particle, $d_{p}$, because particle-wall collision rates scale linearly with $N_{c}$ and naturally do not break invariance with $W$.

4 Results: super-particles and turbophoresis

This section demonstrates the performance of the proposed enhanced collision radius approach. In particular, the low and high Stokes number limits of Eq. (26) are verified in §4.1. Following that, a hybrid model bridging the two limits is introduced and validated in §4.2. The advantage of the hybrid model is that it computes a unique value of $h$ for each binary collision based on the known properties of the flow and particles participating. Therefore, the hybrid model performs well for low, intermediate, and high Stokes numbers. As motivated in the introduction and §2, the concentration of particles as a function of wall-normal distance is treated as the main quantity of interest, though other features are shown when of particular interest.

4.1 Limiting behaviors

The enhanced collision radius model in §3.3 can be directly employed by assuming a value for $h$. This direct application makes sense for $h=1$ (for low $St^{+}$) or $h=0$ (for high $St^{+}$) by essentially assuming that all particle-particle collisions in a flow occurs either at low or high inertia relative to the surrounding flow. As a first demonstration, the direct application of $h=0$ and $h=1$ is used for monodisperse particle-laden channel flows with three different values of $St^{+}$.

Figures 6-8 show concentration profiles obtained from simulations with a reduced number of computational particles, indicated by increasing statistical weight, $W=N_{p}/N_{c}$, up to a factor of $16$. The order is from large to small Stokes number, beginning with $St^{+}=128$ in Figure 6 and decreasing to $St^{+}=8$ in Figure 8. The Stokes number is changed by altering the particle-to-fluid density ratio. The bulk volume fraction, $\Phi_{V}$, is increased as the Stokes number is decreased to keep the mass fraction from changing dramatically, keeping the simulations in a range where two-way coupling does not strongly influence the concentration profiles.

The results for $St^{+}=128$ are shown in Figure 6. In panel (a), the baseline result is obtained by making no correction to account for changes in particle-particle collisions. In that case, the sensitivity of the concentration profile to collision-driven changes in turbophoresis, as earlier demonstrated in Figure 2, is apparent. Simply decreasing the number of computational particles without correction leads to a lower number of collisions per particle. As a result, the turbophoretic drift is stronger and larger concentrations are observed near the wall. In fact, with a $16$x factor reduction in number of computational particles, the near wall concentration is over-predicted by an order of magnitude compared to the $N_{c}=N_{p}$ case.

Panels (b) and (c) of Figure 6 show results using the enhanced collision radius model with $h=1$ and $h=0$, respectively. When $h=1$, the discrepancies in predicted concentration profile with increasing $W$ are certainly less than in the case of no correction. Nevertheless, these discrepancies grow with increasing reduction in particle count, indicating the lack of statistical invariance with respect to changing $W$. In contrast, the concentration profiles using $h=0$ approximately collapse even as $W$ is increased. The reason $h=0$ succeeds in collapsing the concentration profiles is tied to the relatively large Stokes number. At $St^{+}=128$, the particles collide in a ballistic manner and therefore a larger correction of $d_{c}=d_{p}\sqrt{W}$ is necessary to maintain the same number of collisions per particle as the number of particles is reduced. The smaller correction of $d_{c}=d_{p}\sqrt[3]{W}$, being appropriate for smaller Stokes numbers, still under-predicts the number of collisions with reduced number of computational particles.

Figure 7 shows the results for $St^{+}=32$. Again, with no correction, panel (a) shows that the concentration profile changes dramatically as the number of computational particles is reduced. The low Stokes number correction, $h=1$, in panel (b) still shows significant increasing discrepancy with increasing $W$. In panel (c), the ballistic scale correction, $h=0$, demonstrates excellent performance for $St^{+}=32$, as was also true for $St^{+}=128$. Again, the concentration profiles collapse for all values of $W$, signaling that invariance with respect to $W$ as been achieved. Thus, $St^{+}=32$ and $St^{+}=128$ both appear to be safely in the high Stokes number regime anticipated in §3.3.

However, continuing to decrease the Stokes number leads to different results. Figure 8 shows concentration profiles for simulations at $St^{+}=8$ for decreasing number of particle concentrations. Results similar to the higher Stokes numbers are obtained when no correction is performed, panel (a). However, judging from panels (b) and (c) of Figure 8, the low Stokes number correction, $h=1$, appears to lead to better results. This qualitative trend is exactly as expected from the model development in §3.3. In the case of lower Stokes numbers, an increase in collision radius also increases the relative velocity with which the particles collide. As a result, the collision radius needs only to be scaled up as $d_{c}=d_{p}\sqrt[3]{W}$ to maintain accurate collisional activity with decreasing number of computational particles. This is reflected by the invariance of the concentration profile with respect to $W$ in Figure 8(b). Panel (c) shows that the $h=0$ scaling $d_{c}=d_{p}\sqrt{W}$ derived from high Stokes number assumptions actually over-corrects for collisional effects as $W$ is increased. In that case, the collision radius is made too large and the number of collisions actually increases with decreasing number of computational particles.

Taken together, the results shown in Figures 6-8 verify the expectations outlined for the large and small Stokes number limits, Eq. (26). A simple scaling correction applied universally across all particles and all collisions works surprisingly well in generating collisional effects to maintain accurate concentration profiles even as the number of particles is dramatically reduces. Before, leaving this section, a few notes are in order. First, results at $St^{+}=512$ (not shown) indicate that the $h=0$ correction continues to succeed at $St^{+}$ is further increased beyond what it shown here. At $St^{+}=2$, the decreasing slip velocity between the particle and surrounding fluid leads to a decrease in the turbophoretic drift, see Figure 1. As a result, concentration profiles at $St^{+}=2$ are not as sensitive to collisions in the volume fractions studied here. It follows, then, that even without a correction for the collisional radius, accurate concentration profiles can be obtained in the super-particle approach.

4.2 Hybrid model

While the choice $h=0$ seems to be effective for $St^{+}\geq 32$ but $h=1$ works better for $St^{+}\leq 8$, general simulation scenarios are not necessarily characterized by a single Stokes number due to complex flow inhomogeneities and multidisperse particle size distributions encountered in real flows. It is more desirable to specify a unique $h$ value for each collision based on the properties of the collision known to the simulation. The basic difference between $h=0$ and $h=1$ is whether particles colliding at a separation $d_{c}$ would have a significantly different (lower) collision velocity if allowed to approach to a separation of $d_{p}$ before colliding. Therefore, it is straightforward to propose a dimensionless parameter,

[TABLE]

which encapsulates this behavior. The numerator is the stopping distance of a particle under Stokes drag, which indicates the relative distance the particles can travel before relaxing significantly toward local relative flow velocities. The denominator is simply a measure of the additional distance the particles must travel relative to each other before a collision based on their physical diameter. When $b$ is large, the particles respond very slowly to local relative flow near collision, indicating a mostly ballistic interaction ($h=0$). In the opposite limit of small $b$, the particles adjust nearly instantaneously to the local flow and thus experience a collision in the ‘smooth’ regime ($h=1$). A simple model containing these asymptotics is

[TABLE]

where $b_{tr}=32$ is set empirically. This leads to an implicit formula for $h$, namely,

[TABLE]

but this is observed to converge quickly under iteration. Since $d_{c}$ is not known until after collision properties are computed, collision detection during each time step is done assuming the maximum $d_{c}$, namely, assuming $h=0$ at first. Once $h$ is determined, the diameter and restitution coefficient of the collision are adjusted according to Eqs. (28) and (31).

The results of this super-particle collision model for the channel flow at $Re_{*}=150$ are shown in Figure 9 for three different Stokes numbers. It is apparent that, for each Stokes number, relatively good invariance is obtained with respect to $W$. This means that accurate concentration profiles are obtained even when the number of computational particles is dramatically reduced. The hybrid enhanced collision radius model, therefore, provides a realistic approach for the super-particle approach in situations where collisions are key to the accuracy of the simulation. To further support this assertion, additional quantities are shown in Figure 10 and 11. Note that wall-normal distances only up to $y^{+}=60$ are shown in these two figures due to the rarity of collisions, and hence noisy data further from the wall.

Figure 10 shows how collision frequency varies with distance from the wall for increasing statistical weight, $W$. The curves largely overlap for all three Stokes numbers, which is further evidence that particle statistics are roughly invariant with changes in $W$. Some discrepancy is noticeable within $y^{+}\leq 1$, and this is reflected in Figure 9 by some slight discrepancies in the near wall concentration. The discrepancy in collision rate, and hence concentration, in the very near wall region is due to the larger finite diameter relative to important flow features in the the viscous sublayer caused by the enhanced collision radius effect ($d_{p}^{+}=0.5$, so $d_{c}^{+}$ increases up to $2.0$ in the worst case shown). Although this effect remains minor for the cases shown here, this illustrates the main limitation of the enhanced collisional radius model developing in this paper. While the physical particle size may satisfy $d_{p}\lesssim\eta$ in a general turbulent flow, the distance at which super-particles collide may exceed this limitation if $W$ is taken too large. The collision radius for particle-wall collisions is unchanged in the super-particle approach, still computed based on the physical particle diameter. Tri-linear interpolation is used to obtain the fluid velocity at each particle location, even below the first grid point. Tests with higher-order interpolation showed no significant changes in the concentration profile (Johnson et al., 2019).

The average magnitude of relative velocity at collision is show in Figure 11. The noise in these plots serves as a reminder by illustrating that statistical convergence requires longer averaging times or more snapshots in the super-particle approach simply because less samples are available. Nonetheless, it can be seen that the super-particle approach with enhanced collision radius scaling performs reasonably well. The discrepancies, particularly for $y^{+}\leq 10$ are larger than for concentration or collision frequency in the $St^{+}=8$ and $St^{+}=32$ cases. At higher $St^{+}=128$, the collapse with varying $W$ is more impressive. The lower Stokes number particles interact more nontrivially with near-wall coherent structures of the turbulent flow, which leads to less accuracy for the finer details in a super-particle approach.

It should be appreciated that the collision treatment here is not stochastic in nature. Each collision is computed in the same deterministic way as in the simulation with the full number of particles. While a hard sphere collision model was used in this study, the model is directly applicable to other approaches such as the soft sphere collision model. In principle, it could also be extended to a resolved particle approach as well, though some of the benefits of resolving the flow around individual particles would necessarily be forfeited, such as the capturing of the lubrication interaction near collision or accurately simulating particles for $d_{p}\gtrsim\eta$.

5 Summary & Conclusions

The simulation of particles subject to turbophoresis in a turbulent wall-bounded flows motivated the present work. At very low volume fractions, turbophoresis leads to dramatically enhanced concentrations in a thin layer of fluid adjacent to the wall. As volume fraction increases, particle-particle collisions play a key role in limiting the turbophoretic drift and reducing the near-wall peak in concentration. It was demonstrated here that this trend occurs in such a way that the collisions effectively create a maximum near-wall concentration that can be sustained by turbophoresis. Therefore, simply using a reduced number of computational particles can lead to large errors in near-wall concentration levels.

In this paper, an enhanced collision radius model was introduced to yield simulation results approximately invariant with respect to the number of computational particles used in the simulation. This treatment of simulated collisions for super-particles extends of the method of Garg et al. (2009), see also Subramaniam (2013) by artificially enhancing the collision radius of the particles, i.e. the separation distance at which two particles collide. The primary challenge in prescribing such a treatment is to specify how quickly the collision radius must grow as the number of computational particles decreases so as to keep collision rates more or less constant.

The specification of the collision radius can be done using scaling arguments for collision rates to approximate the conditions for statistical equivalence between the physical and computational particle ensembles, which is based on assigning a statistical weight to each computational particle. The proposed scalings of the collision radius with statistical weight were based on analytically established limits for low and high Stokes numbers. The method was first tested in these limits and showed promising results. Further, a hybrid model was designed to bridge these limits by relating the scaling of the collision radius to the properties of each individual collision. The hybrid model locally detects the extent to which each particle-particle collision is of a ballistic nature. As a result, the proposed model results in particle statistics that produce very similar particle statistics even as the number of computational particles is reduced.

Acknowledgements

This investigation was funded by the Advanced Simulation and Computing program of the US Department of Energy’s National Nuclear Security Administration via the PSAAP-II Center at Stanford, Grant No. DE-NA0002373. Jeremy Horwitz is thanked for a detailed review of this manuscript and for insightful discussions about the super-particle framework and his complementary research into two-way coupling effects. This work also benefited from interaction with Parviz Moin, Mahdi Esmaily, Ali Mani, Maxime Bassenne, Javier Urzay, and Lluis Jofre-Cruanyes. I would also like to thank Cristian Marchioli (University of Udine, Italy) for sharing benchmark numerical results that were useful in verifying our simulation code.