2.1 Kinetic energy equations

In the context of multiphase flow, the conservation equations are derived from statistical approaches. One such statistical procedure used in gas–solid flow is the random field (RF) approach (Pai & Subramaniam Reference Pai and Subramaniam2009; Tenneti & Subramaniam Reference Tenneti and Subramaniam2014; Mehrabadi & Subramaniam Reference Mehrabadi and Subramaniam2017). In this approach the transport equations of mass, momentum and energy are multiplied by a phasic indicator function $I^{(\unicode[STIX]{x1D6FD})}$ with $\unicode[STIX]{x1D6FD}\in (f,p)$ which is 1 in phase $\unicode[STIX]{x1D6FD}$ and 0 otherwise. Here, $(f)$ and $(p)$ denote the fluid and particle phases, respectively. The equations are then ensemble averaged over all possible realizations in the event space. If the system is statistically homogeneous and ergodic, then ensemble averaging and volume averaging are equivalent. Therefore, the mean value of property $Q$ in phase $\unicode[STIX]{x1D6FD}$ is given by

(2.1) $$\begin{eqnarray}\langle Q^{(\unicode[STIX]{x1D6FD})}\rangle =\frac{\langle I^{(\unicode[STIX]{x1D6FD})}Q\rangle }{\langle I^{(\unicode[STIX]{x1D6FD})}\rangle }=\frac{\displaystyle \int I^{(\unicode[STIX]{x1D6FD})}Q\,\text{d}V}{\displaystyle \int I^{(\unicode[STIX]{x1D6FD})}\,\text{d}V},\end{eqnarray}$$

and subsequently the fluctuating component is defined as $Q^{\prime (\unicode[STIX]{x1D6FD})}=Q-\langle Q^{(\unicode[STIX]{x1D6FD})}\rangle$ . It is worth mentioning that in the above equation the denominator represents the volume fraction of phase $\unicode[STIX]{x1D6FD}$ .

Accordingly, the level of kinetic energy in phase $\unicode[STIX]{x1D6FD}$ is defined as $k^{(\unicode[STIX]{x1D6FD})}=\langle \boldsymbol{u}^{\prime (\unicode[STIX]{x1D6FD})}\boldsymbol{\cdot }\boldsymbol{u}^{\prime (\unicode[STIX]{x1D6FD})}\rangle /2$ with the velocity fluctuation being $\boldsymbol{u}^{\prime (\unicode[STIX]{x1D6FD})}=\boldsymbol{u}^{(\unicode[STIX]{x1D6FD})}-\langle \boldsymbol{u}^{(\unicode[STIX]{x1D6FD})}\rangle$ (in the present study, $\langle \boldsymbol{u}^{(\unicode[STIX]{x1D6FD})}\rangle$ is zero). The evolution equations of gas-phase and solid-phase kinetic energy for a homogeneous gas–solid suspension are, respectively, given by (Pai & Subramaniam Reference Pai and Subramaniam2009; Subramaniam et al. Reference Subramaniam, Mehrabadi, Horwitz and Mani2014; Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015)

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}^{(f)}(1-\unicode[STIX]{x1D719})\frac{\text{d}k^{(f)}}{\text{d}t}=\unicode[STIX]{x1D6F1}^{(f)}-\unicode[STIX]{x1D700}^{(f)}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}^{(p)}\unicode[STIX]{x1D719}\frac{\text{d}k^{(p)}}{\text{d}t}=\unicode[STIX]{x1D6F1}^{(p)}. & \displaystyle\end{eqnarray}$$

In the above equations, it is assumed that the rotational energy of particles is negligible compared to the translational kinetic energies. Also, the particles are assumed to be elastic, so that collisions make no net contribution to energy exchanges. In the above equations, $\unicode[STIX]{x1D70C}^{(f)}$ and $\unicode[STIX]{x1D70C}^{(p)}$ are the gas-phase and solid-phase mass densities, $\unicode[STIX]{x1D719}$ is the solid-phase volume fraction, and $\unicode[STIX]{x1D700}^{(f)}=2\unicode[STIX]{x1D707}^{(f)}\langle I^{(f)}\unicode[STIX]{x1D634}_{ij}\unicode[STIX]{x1D634}_{ij}\rangle$ ( $\unicode[STIX]{x1D668}$ being the strain-rate tensor) is the viscous dissipation (Subramaniam et al. Reference Subramaniam, Mehrabadi, Horwitz and Mani2014; Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015). In addition, $\unicode[STIX]{x1D6F1}^{(f)}$ and $\unicode[STIX]{x1D6F1}^{(p)}$ are respectively the fluid-phase and solid-phase interphase TKE transfer terms. It has been shown (Xu & Subramaniam Reference Xu and Subramaniam2007; Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015; Mehrabadi & Subramaniam Reference Mehrabadi and Subramaniam2017) that the interphase TKE transfer terms, arising from the velocity fluctuation–drag force covariance, play key roles in transfer of the kinetic energy due to interphase interactions. By definition, these transfer terms have the following forms:

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}^{(f)}=-\langle \boldsymbol{u}^{(f)}\boldsymbol{\cdot }\boldsymbol{f}\rangle , & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}^{(p)}=\langle \boldsymbol{u}^{(p)}\boldsymbol{\cdot }\boldsymbol{f}\rangle , & \displaystyle\end{eqnarray}$$

where $\boldsymbol{f}=-\unicode[STIX]{x1D748}\boldsymbol{\cdot }\unicode[STIX]{x1D735}I^{(p)}$ is the interfacial force per unit volume, and $\unicode[STIX]{x1D748}$ is the fluid stress tensor. It is worth mentioning that the gradient of the indicator function can be rewritten as $-\boldsymbol{n}^{(p)}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}^{(I)})$ (Drew & Passman Reference Drew and Passman1998), where $\boldsymbol{n}^{(p)}$ is the normal vector at each particle surface, and $\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}^{(I)})$ is a generalized delta function which is non-zero only at the fluid–solid interface. The alternative form of the drag force $\boldsymbol{f}=\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}^{(p)}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}^{(I)})$ is particularly useful since the computed drag force (integral of the fluid stress over the particle surface) in PR-DNS is reconcilable with this definition. In addition, it has been shown that this definition in conjunction with the kinematic condition at the interface leads to the conservation of the interphase TKE transfer principle which states that, for the case of zero mean-slip velocity, the interphase transfer terms are equal in magnitude and opposite in sign, i.e. $\unicode[STIX]{x1D6F1}^{(f)}+\unicode[STIX]{x1D6F1}^{(p)}=0$ (Xu & Subramaniam Reference Xu and Subramaniam2007; Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015; Mehrabadi & Subramaniam Reference Mehrabadi and Subramaniam2017). This implies that the rate of energy lost by one phase owing to interphase interactions at mutual fluid–solid interfaces is equal in magnitude to the rate of energy gained by the other phase. As a result, the mixture kinetic energy $e^{(m)}$ defined as $\unicode[STIX]{x1D70C}^{(f)}(1-\unicode[STIX]{x1D719})k^{(f)}+\unicode[STIX]{x1D70C}^{(p)}\unicode[STIX]{x1D719}k^{(p)}$ , which is the mass-weighted kinetic energy in the system governed by:

(2.6) $$\begin{eqnarray}\frac{\text{d}e^{(m)}}{\text{d}t}=\unicode[STIX]{x1D6F1}^{(f)}+\unicode[STIX]{x1D6F1}^{(p)}-\unicode[STIX]{x1D700}^{(f)},\end{eqnarray}$$

simplifies to

(2.7) $$\begin{eqnarray}\frac{\text{d}e^{(m)}}{\text{d}t}=-\unicode[STIX]{x1D700}^{(f)}.\end{eqnarray}$$

Therefore, the mixture energy decays monotonically due to the viscous dissipation in the absence of any driving force. This equation is also reconcilable with the PR-DNS formulation in which the drag force on each particle surface is directly recovered. Nevertheless, this conservation principle is violated by PP-DNS (or any other numerical approach) in which a particle drag model is used in place of direct enforcement of no-slip and no-penetration boundary conditions at particle surfaces.

2.2 Point-particle model kinetic energy equations

In the PP-DNS implementation, the drag force on particle $i$ , i.e. $\boldsymbol{F}_{i,pp}$ is provided by an assumed drag model. To project this drag force to the Eulerian grid, the following expression is used:

(2.8) $$\begin{eqnarray}\boldsymbol{f}_{pp}=\frac{1}{V_{cell}}\mathop{\sum }_{i=1}^{N_{p}}\boldsymbol{F}_{i,pp}{\mathcal{P}}(\unicode[STIX]{x1D6FF}(\boldsymbol{X}_{i}-\boldsymbol{x})).\end{eqnarray}$$

In the above equation, the subscript $pp$ denotes point-particle quantities. In addition, the operator ${\mathcal{P}}(\unicode[STIX]{x1D6FF}(\boldsymbol{X}_{i}-\boldsymbol{x}))$ is the numerical projection operator that is non-zero on the grid points of the corresponding Eulerian cell with the cell volume $V_{cell}$ . Considering the above definition for the point-particle drag force, the evolution equations for implied fluid-phase and solid-phase velocity fluctuations are (Sundaram & Collins Reference Sundaram and Collins1999; Subramaniam et al. Reference Subramaniam, Mehrabadi, Horwitz and Mani2014)

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}^{(f)}(1-\unicode[STIX]{x1D719})\frac{\text{d}k_{pp}^{(f)}}{\text{d}t}=\unicode[STIX]{x1D6F1}_{pp}^{(f)}-\unicode[STIX]{x1D700}_{pp}^{(f)}, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}^{(p)}\unicode[STIX]{x1D719}\frac{\text{d}k_{pp}^{(p)}}{\text{d}t}=\unicode[STIX]{x1D6F1}_{pp}^{(p)}. & \displaystyle\end{eqnarray}$$

In (2.9), $-\unicode[STIX]{x1D700}_{pp}^{(f)}\equiv (1-\unicode[STIX]{x1D719})(1/V)\int \unicode[STIX]{x1D707}^{(f)}\boldsymbol{u}^{(f)}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{(f)}\,\text{d}V$ , and the integration is over the whole fluid volume $V$ . Using the modelled drag force in (2.8), the implied interphase TKE transfer terms in the PP-DNS formulation now have the following forms:

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}_{pp}^{(f)}=(1-\unicode[STIX]{x1D719})\frac{1}{V}\mathop{\sum }_{i=1}^{N_{p}}\boldsymbol{u}^{(f)}(\boldsymbol{X}_{i})\boldsymbol{\cdot }\boldsymbol{F}_{i,pp}, & \displaystyle\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}_{pp}^{(p)}=\frac{1}{V}\mathop{\sum }_{i=1}^{N_{p}}\boldsymbol{V}_{i}\boldsymbol{\cdot }\boldsymbol{F}_{i,pp}. & \displaystyle\end{eqnarray}$$

In (2.11) and (2.12), $\boldsymbol{u}^{(f)}(\boldsymbol{X}_{i})$ represents the fluid velocity at the $i\text{th}$ particle location, and $\boldsymbol{V}_{i}$ is the particle velocity. The mixture energy for point particles, $e_{pp}^{(m)}=\unicode[STIX]{x1D70C}^{(f)}(1-\unicode[STIX]{x1D719})k_{pp}^{(f)}+\unicode[STIX]{x1D70C}^{(p)}\unicode[STIX]{x1D719}k_{pp}^{(p)}$ , can then be expressed as:

(2.13) $$\begin{eqnarray}\frac{\text{d}e_{pp}^{(m)}}{\text{d}t}=\unicode[STIX]{x1D6F1}_{pp}^{(f)}+\unicode[STIX]{x1D6F1}_{pp}^{(p)}-\unicode[STIX]{x1D700}_{pp}^{(f)}.\end{eqnarray}$$

The summation of the interphase TKE transfer terms in the above equations reveals that

(2.14) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{pp}^{(f)}+\unicode[STIX]{x1D6F1}_{pp}^{(p)}=\frac{1}{V}\mathop{\sum }_{i=1}^{N_{p}}(\boldsymbol{V}_{i}-\boldsymbol{u}^{(f)}(\boldsymbol{X}_{i}))\boldsymbol{\cdot }\boldsymbol{F}_{i,pp}+O(\unicode[STIX]{x1D719}^{2}).\end{eqnarray}$$

In (2.14), the additional term which is $O(\unicode[STIX]{x1D719}^{2})$ (see appendix A) is an error consistent with the neglect of volume fraction in the point-particle equations. The kinematic constraint of zero relative velocity at particle surfaces which is enforced by the PR-DNS method, does not hold for the PP-DNS method, i.e. $\boldsymbol{V}_{i}-\boldsymbol{u}^{(f)}(\boldsymbol{X}_{i})\neq 0$ . As a result, the summation of $\unicode[STIX]{x1D6F1}_{pp}^{(f)}$ and $\unicode[STIX]{x1D6F1}_{pp}^{(p)}$ in (2.14) is non-zero which means that the conservation of interphase TKE transfer principle does not hold. In the PP-DNS formulation, equation (2.14) is interpreted as an additional dissipation near particle surfaces, and is denoted as $\unicode[STIX]{x1D700}_{pp}^{\ast }$ (Sundaram & Collins Reference Sundaram and Collins1999; Xu & Subramaniam Reference Xu and Subramaniam2007). The additional dissipation can be interpreted as the model-form dissipation arising from the fact that not all of the true fluid dissipation is resolved by the term, $\unicode[STIX]{x1D700}_{pp}^{(f)}$ . Therefore, the use of a particle acceleration model for fluid–solid interactions in the PP-DNS approach leads to a different form of the mixture energy equation, that is

(2.15) $$\begin{eqnarray}\frac{\text{d}e^{(m)}}{\text{d}t}=-\unicode[STIX]{x1D700}_{pp}^{\ast }-\unicode[STIX]{x1D700}_{pp}^{(f)},\end{eqnarray}$$

which is now governed by the under-resolved viscous dissipation $\unicode[STIX]{x1D700}_{pp}^{(f)}$ and a model for the additional dissipation near particle surfaces $\unicode[STIX]{x1D700}_{pp}^{\ast }$ .

The formulations for the TKE equations presented here show that no direct comparison is possible between viscous dissipation obtained from PR-DNS and PP-DNS approaches (cf. equations (2.7) and (2.15)) unless the effect of the additional dissipation near particle surfaces $\unicode[STIX]{x1D700}_{pp}^{\ast }$ is accounted for (Sundaram & Collins Reference Sundaram and Collins1996, Reference Sundaram and Collins1999; Subramaniam et al. Reference Subramaniam, Mehrabadi, Horwitz and Mani2014; Horwitz & Mani Reference Horwitz and Mani2016). In this study it is of interest to investigate the effect of the gas–particle drag force in the PP-DNS approach on the implied statistics of turbulence when compared with those of the PR-DNS. Further details of the point-particle energetics derivation can be found in appendix A.

3.1 Particle-resolved direct numerical simulation

PR–DNS of freely evolving gas–solid suspensions are performed using the particle-resolved uncontaminated-fluid reconcilable immersed boundary method (PUReIBM) that has been developed to simulate flow past fixed particle assemblies (Garg et al. Reference Garg, Tenneti, Mohd-Yusof, Subramaniam, Pannala, Syamlal and O’Brien2011; Tenneti et al. Reference Tenneti, Sun, Garg and Subramaniam2013; Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015; Sun, Tenneti & Subramaniam Reference Sun, Tenneti and Subramaniam2015; Mehrabadi, Murphy & Subramaniam Reference Mehrabadi, Murphy and Subramaniam2016a ; Sun et al. Reference Sun, Tenneti, Subramaniam and Koch2016) and freely evolving suspensions (Tenneti et al. Reference Tenneti, Mehrabadi and Subramaniam2016; Mehrabadi, Tenneti & Subramaniam Reference Mehrabadi, Tenneti and Subramaniam2016b ). In PUReIBM, the entire physical domain is discretized using a uniform Cartesian grid and the governing equations of fluid flow are solved on all the grid points (including those lying inside the particles). The governing equations that are solved in PUReIBM are the continuity equation

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$

and the Navier–Stokes equations

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D70C}^{(f)}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D70C}^{(f)}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{u})=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}^{(f)}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{f}_{IB},\end{eqnarray}$$

where $\boldsymbol{u}$ and $p$ are respectively the instantaneous velocity and pressure fields. The third term on the right-hand side of (3.2) is the immersed boundary (IB) force $\boldsymbol{f}_{IB}$ that accounts for the presence of solid particles by enforcing the no-slip and no-penetration boundary conditions at the surface of the particles. Complete details of the computation of the IB force can be found in the work of Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011). The surface of the spherical particle is represented by a discrete number of Lagrangian points called boundary points that are parametrized in spherical coordinates. Two additional sets of points, termed as exterior and interior points are generated by projecting the boundary points onto spheres of radii $r+\unicode[STIX]{x0394}r$ and $r-\unicode[STIX]{x0394}r$ respectively, with $\unicode[STIX]{x0394}r$ chosen to be equal to the grid spacing. In PUReIBM, the IB force is computed only at the interior points so that a desired velocity $\boldsymbol{u}^{(d,k)}$ is obtained at the $k\text{th}$ interior point with respect to the $n\text{th}$ time step. Following the direct forcing method proposed by Mohd-Yusof (Reference Mohd-Yusof1996) the IB force $\boldsymbol{f}_{IB}^{(k)}$ at the $k\text{th}$ interior point is computed as:

(3.3) $$\begin{eqnarray}\boldsymbol{f}_{IB}^{(k)}=\unicode[STIX]{x1D70C}^{(f)}\frac{\boldsymbol{u}^{(k,d)}-\boldsymbol{u}^{(k,n)}}{\unicode[STIX]{x0394}t}+\unicode[STIX]{x1D70C}^{(f)}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{u})^{(k,n)}+\unicode[STIX]{x1D735}p^{(k,n)}-\unicode[STIX]{x1D707}^{(f)}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{(k,n)}.\end{eqnarray}$$

The desired velocity $\boldsymbol{u}^{(k,d)}$ depends on the velocity of the particle. For instance, for a fixed particle the desired velocity at the interior point is equal in magnitude but opposite in direction of the fluid velocity at the corresponding exterior point so that the velocity at the boundary point is zero. The IB force computed at all the interior points is distributed among the neighbouring Cartesian grid nodes, but limited only to those grid points in solid particles.

It has been shown by Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011) that volume-averaged equations in a periodic domain can be reconciled with the ensemble-averaged equations of statistically homogeneous gas–solid suspensions. Hence the governing equations (3.1) and (3.2) are solved numerically in a cubic domain with periodic boundaries, since the primary goal is to simulate statistically homogeneous suspensions. Similar to DNS of single-phase turbulence, it is simpler to solve for appropriately defined fluctuating variables. The velocity field is decomposed into a spatially uniform mean flow $\langle \boldsymbol{u}\rangle _{{\mathcal{V}}}(t)$ and a fluctuating velocity field $\boldsymbol{u}^{\prime }(\boldsymbol{x},t)$ , i.e.

(3.4) $$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\langle \boldsymbol{u}\rangle _{{\mathcal{V}}}(t)+\boldsymbol{u}^{\prime }(\boldsymbol{x},t),\end{eqnarray}$$

where the volumetric mean velocity is obtained by averaging the velocity field over the entire computational domain. Similar decompositions are written for the convective term, the pressure gradient $\unicode[STIX]{x1D735}p$ and IB force $\boldsymbol{f}$ . It should be noted that the mean values of the velocity, pressure gradient and immersed boundary forcing are zero in an isotropic turbulent flow. Substituting the above decompositions in (3.1) and (3.2), followed by averaging over the entire computational domain yields the volume-averaged mass and momentum conservation equations (detailed equations are given by Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011)). Evolution equations for the fluctuating variables are derived by subtracting the volume-averaged equations from their instantaneous counterparts. The resulting equations for the fluctuating variables are solved using a pseudo-spectral method, with the Crank–Nicolson scheme for the viscous terms, and an Adams–Bashforth scheme for the convective terms. A fractional time-stepping method based on the work of Kim & Moin (Reference Kim and Moin1985) is used to advance the fluctuating velocity fields in time.

In freely evolving gas–solid suspensions, the particles move under the influence of hydrodynamic and collisional forces. In PUReIBM, the particles are represented in a Lagrangian frame of reference. The position and translational velocity of the $i\text{th}$ particle at time $t$ are denoted by $\boldsymbol{X}_{i}(t)$ and $\boldsymbol{V}_{i}(t)$ respectively, and they evolve in time as:

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{X}_{i}(t)}{\text{d}t}=\boldsymbol{V}_{i}(t), & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle m_{p}\frac{\text{d}\boldsymbol{V}_{i}(t)}{\text{d}t}=\boldsymbol{F}_{i}(t)+\mathop{\sum }_{\substack{ j=1 \\ j\neq i}}^{N_{p}}\unicode[STIX]{x1D60A}_{ij}(t), & \displaystyle\end{eqnarray}$$

where $m_{p}=\unicode[STIX]{x1D70C}_{p}\unicode[STIX]{x03C0}d_{p}^{3}/6$ is the particle mass, and $N_{p}$ is the total number of particles used in the simulation. Here $\boldsymbol{F}_{i}$ is the hydrodynamic force calculated from the velocity and pressure fields at the particle surface. Also $\unicode[STIX]{x1D60A}_{ij}$ is the contact force on the $i\text{th}$ particle as a result of collision with the $j\text{th}$ particle, and is treated using a soft-sphere model originally proposed by Cundall & Strack (Reference Cundall and Strack1979). For the present study, the collisions are modelled as elastic. Nevertheless, since the suspension is very dilute in the cases studied here, particle collisions are infrequent. Some more discussion on collision frequency and the effect of collisions on point-particle results is presented in appendix B.

3.2 Point-particle direct numerical simulation

The point-particle algorithm is based on the code originally developed by Pouransari, Mortazavi & Mani (Reference Pouransari, Mortazavi and Mani2015). The fluid-phase governing equations are the Navier–Stokes equations, (equations (3.1) and (3.2)) but written in conservative form. However, the immersed boundary force $\boldsymbol{f}_{IB}$ in (3.2) is replaced by the negative of the gas–particle drag evaluated on the Eulerian grid by (2.8), that is

(3.7) $$\begin{eqnarray}-\boldsymbol{f}_{IB}\equiv \boldsymbol{f}_{pp}=\frac{1}{V_{cell}}\mathop{\sum }_{i=1}^{N_{p}}\boldsymbol{F}_{i,pp}{\mathcal{P}}(\unicode[STIX]{x1D6FF}(\boldsymbol{X}_{i}-\boldsymbol{x})).\end{eqnarray}$$

In addition, the particles are tracked in a Lagrangian frame using the same Newtonian laws of motion given in (3.5) and (3.6). The Navier–Stokes equations are solved on a staggered grid using second-order finite differences. The resulting Poisson equation for pressure is solved directly using the fast Fourier transform. Time advancement of both fluid and particle equations is accomplished by explicit fourth-order Runge–Kutta.

To complete the point-particle algorithm, it is necessary to specify the projection operator ${\mathcal{P}}\left(\boldsymbol{\cdot }\right)$ , and drag force $\boldsymbol{F}_{i,pp}$ . The particles explored in this study will have a size equal to the initial Kolmogorov scale of the fluid, that is $d_{p}/\unicode[STIX]{x1D702}_{0}=1$ . From a scaling perspective, this particle size ratio is borderline in terms of the need to model unsteady effects in the particle equation of motion (Ling, Parmar & Balachandar Reference Ling, Parmar and Balachandar2013). As we will show, some of these unsteady effects can be alleviated by seeding the particles with the local fluid velocity rather than at rest, as was done in Subramaniam et al. (Reference Subramaniam, Mehrabadi, Horwitz and Mani2014). Faxén corrections are ignored for two reasons. Faxén corrections are important when particles are large enough to see the curvature in the underlying fluid velocity field. Their importance can be shown to scale with the ratio of the particle size to Taylor microscale $(d_{p}/\unicode[STIX]{x1D706})^{2}$ (Calzavarini et al. Reference Calzavarini, Volk, Bourgoin, Leveque, Pinton and Toschi2009), which is much smaller than unity for the present problem. However, it is also important to assess whether the curvature induced by the disturbance fields of neighbour particles will be important. Since the volume fraction in this study is relatively low ( $10^{-3}$ ), the contribution to the Faxén correction owing to disturbance flows created by neighbour particles is small, which provides an additional reason to neglect Faxén corrections. Collisions are also neglected in the point-particle simulation owing to diluteness. In this study we only consider steady drag (Stokes, and Reynolds corrected) and neglect history, added mass, and fluid acceleration forces, to assess whether the simplest version of a point-particle method can find some agreement with a particle-resolved simulation, in the limit where the point-particle method is believed to be valid $d_{p}/\unicode[STIX]{x1D702}\leqslant 1$ . Our goal with this study is validation of the point-particle method, so we were not a priori resigned to only use steady drag. Rather we sought to add levels of complexity until statistics from the two approaches agreed within acceptable tolerance. What we will show in the results section, is that an appropriate numerical implementation of steady drag was found to yield good agreement between the two methods, for this particular non-dimensional problem. It is not a universal statement that the drag forces we consider here, and their numerical implementation, will be sufficient for validation problems in other regimes. In this study we consider Stokes (linear) drag, and Schiller–Naumann (nonlinear). Their formulas are:

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{i}^{St}=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}^{(f)}d_{p}(\tilde{\boldsymbol{u}}_{i}^{(f)}-\boldsymbol{V}_{i}), & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{i}^{SN}=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}^{(f)}d_{p}(\tilde{\boldsymbol{u}}_{i}^{(f)}-\boldsymbol{V}_{i})\;(1+0.15\times {Re_{p}}^{0.687}), & \displaystyle\end{eqnarray}$$

with the particle Reynolds number being defined as

(3.10) $$\begin{eqnarray}Re_{p}=\frac{d_{p}\unicode[STIX]{x1D70C}^{(f)}|\tilde{\boldsymbol{u}}_{i}^{(f)}-\boldsymbol{V}_{i}|}{\unicode[STIX]{x1D707}^{(f)}}.\end{eqnarray}$$

Note that in the above formulation, $\tilde{\boldsymbol{u}}_{i}^{(f)}$ is the undisturbed fluid velocity at the $i\text{th}$ particle location, the importance of which is discussed in the following section.

3.3 Consistent numerical implementation of the interphase momentum transfer

The drag formulas in (3.8) and (3.9) depend on the difference between the undisturbed fluid velocity $\tilde{\boldsymbol{u}}_{i}^{(f)}$ evaluated at the particle location, and the particle’s velocity. The undisturbed fluid velocity is the fluid velocity the particle would see in the absence of that particle. This means the undisturbed fluid velocity evaluated at particle $j$ contains a contribution owing to the summation of velocity disturbances created by all other particles $i\neq j$ , but not the disturbance created by particle $j$ . This is a modelling problem because the undisturbed fluid velocity is not readily accessible in the simulation. Usually the fluid velocity interpolated from the surrounding fluid grid points to the particle location is assumed to be approximately equal to the undisturbed fluid velocity. However, when the particle size relative to the grid spacing $\unicode[STIX]{x1D6EC}=d_{p}/\text{d}x$ is significant, then the fluid velocity interpolated to particle $j$ will contain a significant contribution owing to the disturbance velocity field created by particle $j$ . In other words, if this issue is left untreated, this particle will see its own disturbance field, and the drag force which will be used to update its position (and couple momentum and energy to the fluid) will be incorrectly calculated.

This observation was recognized by Horwitz & Mani (Reference Horwitz and Mani2016) who developed a procedure to estimate the undisturbed fluid velocity given a measured disturbed fluid velocity at the location of the particle. The idea is that the undisturbed fluid velocity as used in the above drag formulas for the $j\text{th}$ particle should not include the disturbance created by that particle, but may contain the disturbance information from the other particles. The procedure attempts to remove the self-disturbance portion from the measured disturbed fluid velocity which thereby provides an estimate for the undisturbed fluid velocity seen by that particle. The estimation procedure from that work is:

(3.11) $$\begin{eqnarray}\tilde{\boldsymbol{u}}_{p}^{(f)}=\boldsymbol{u}_{p}^{(f)}+C_{p}\,\text{d}x^{2}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{(f)}.\end{eqnarray}$$

The subscript $p$ in (3.11) denotes quantities evaluated at the location of the particle. Here, $C_{p}$ is an $O(1)$ correction coefficient reported by Horwitz & Mani (Reference Horwitz and Mani2016) and $\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{(f)}$ is the Laplacian of the fluid velocity interpolated to the particle position. Essentially, this correction estimates the undisturbed fluid velocity by adding to the measured fluid velocity at the location of the particle, $\boldsymbol{u}_{p}^{(f)}$ , an estimate for the disturbance created by the particle. The estimate is the Laplacian term. Such a correction is expected to work well provided that the curvature introduced into the fluid velocity field by the particle velocity disturbance is large compared with the underlying curvature already present in the flow. The Laplacian can be decomposed as

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{(f)}\approx \unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{P}+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{F},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{P}$ is the curvature in the fluid velocity created by the particle velocity disturbance and $\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{F}$ is the underlying curvature in the fluid velocity field in the absence of the disturbance flow. Horwitz & Mani (Reference Horwitz and Mani2016) have shown that the disturbance velocity scales as the characteristic fluid velocity and scales with the particle size relative to the grid $\unicode[STIX]{x1D6EC}$ , so that $|C_{p}\,\text{d}x^{2}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{P}|\approx \unicode[STIX]{x1D6EC}u^{\prime }$ , where $u^{\prime }$ here is the r.m.s. velocity of the turbulence. The ensemble average magnitude of the curvature in the fluid velocity field, in the absence of the particle may be estimated as: $|\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{F}|\sim u^{\prime }/\unicode[STIX]{x1D706}^{2}$ , where $\unicode[STIX]{x1D706}$ is the Taylor microscale. Therefore, we estimate the ratio of the inherent curvature in the flow to that produced by a point particle as:

(3.13) $$\begin{eqnarray}|C_{p}\,\text{d}x^{2}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{F}|/|C_{p}\,\text{d}x^{2}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{p}^{P}|\sim \frac{C_{p}\,\text{d}x^{2}u^{\prime }/\unicode[STIX]{x1D706}^{2}}{\unicode[STIX]{x1D6EC}u^{\prime }}\sim \frac{d_{p}^{2}}{\unicode[STIX]{x1D706}^{2}}\ll 1.\end{eqnarray}$$

In (3.13) we have made use of the fact that $C_{p}=O(1)$ , and for the present study $\unicode[STIX]{x1D6EC}=1$ . The particle size being equal to the initial Kolmogorov scale $(d_{p}=\unicode[STIX]{x1D702}_{0})$ means that the ratio in (3.13) is much less than unity ( $d_{p}^{2}/\unicode[STIX]{x1D706}_{0}^{2}\approx 10^{-2}$ for the present cases). In other words, the curvature created by the particles owing to their disturbance fields will be much larger than the inherent curvature in the fluid velocity field. An attempt to estimate the undisturbed fluid velocity using a curvature-based procedure is therefore non-dimensionally justified for this problem. Though we have used scaling analyses above to justify the employed correction procedure for the undisturbed fluid velocity (Horwitz & Mani Reference Horwitz and Mani2016), we have tested the sensitivity of this scheme in a variety of laminar (Horwitz & Mani Reference Horwitz and Mani2016, Reference Horwitz and Mani2018), and turbulent environments (Horwitz & Mani Reference Horwitz and Mani2015, Reference Horwitz and Mani2018), and have found the results to be fairly insensitive to non-dimensional particle size or Stokes number. Though not specifically tested in the context of finite Reynolds number corrections to steady drag, we expect the procedures developed by Gualtieri et al. (Reference Gualtieri, Picano, Sardina and Casciola2015) and Ireland & Desjardins (Reference Ireland and Desjardins2017) based on analytical solutions would also be applicable in the present configuration since the near field of the regularized point-force solution should vary slowly with Reynolds number (Batchelor Reference Batchelor1967). Nevertheless the analytical-based methods developed by Gualtieri et al. (Reference Gualtieri, Picano, Sardina and Casciola2015) and Ireland & Desjardins (Reference Ireland and Desjardins2017) rely on sufficiently large filter widths to remove errors associated with discretization of the Gaussian kernel, while the procedure developed in Horwitz & Mani (Reference Horwitz and Mani2016) was developed specifically taking into account the effect of discretization.

The procedure outlined above is formally accurate in the limit of small volume fraction and Reynolds number, and high Stokes number. It was found that the above formula (3.11), with no modifications works well for order-unity Reynolds number in the Schiller–Naumann equation (Horwitz & Mani Reference Horwitz and Mani2018), the error being of the same order as that which is present in the Schiller–Naumann correlation. In addition, this correction scheme was found to perform reasonably well over a range of Stokes numbers covering the range explored in this study. Therefore, the previous studies demonstrated this correction procedure is verifiable. For simple model problems, our previous work (Horwitz & Mani Reference Horwitz and Mani2016) suggests that particles will follow the correct trajectories and dissipate the correct rate of mechanical energy implied by those force models. In the next section we explore the validation question: can a verifiable numerical implementation of a drag model for point particles, which is valid for the non-dimensional parameters encountered in the present problem, reproduce the particle and fluid statistics implied by the PR-DNS?

Lastly, to complete the point-particle algorithm, it is necessary to specify the projection operator. The correction scheme developed by Horwitz & Mani (Reference Horwitz and Mani2016) to estimate the undisturbed fluid velocity is consistent with projecting the calculated drag force to the Eulerian grid using trilinear weights. For comparison, we will also perform the point-particle simulations assuming the disturbed fluid velocity is equal to the undisturbed fluid velocity (no correction). For these simulations, we will use trilinear weights both for interpolation of the fluid velocity to the particle location and for projection of the drag force back to the fluid grid. Trilinear interpolation/projection is chosen because it was found to be more accurate for two-way coupled problems than higher-order interpolation (Horwitz & Mani Reference Horwitz and Mani2016; Horwitz et al. Reference Horwitz, Rahmani, Geraci, Banko and Mani2016).

Before presenting the results, it is worth comparing the implied point-particle dissipation term with the form presented in Sundaram & Collins (Reference Sundaram and Collins1996). As noted by Sundaram & Collins (Reference Sundaram and Collins1996), when Stokes drag (3.8) is assumed for the point-force model, the sum of the inter-phase energy exchange takes a quadratic form viz:

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{pp}^{(f)}+\unicode[STIX]{x1D6F1}_{pp}^{(p)}=\unicode[STIX]{x1D70C}^{(p)}\unicode[STIX]{x1D719}\left\{\frac{(\boldsymbol{V}-\boldsymbol{u}_{p}^{(f)})\boldsymbol{\cdot }(\tilde{\boldsymbol{u}}_{p}^{(f)}-\boldsymbol{V})}{\unicode[STIX]{x1D70F}_{p}}\right\}_{N_{p}}.\end{eqnarray}$$

In (3.14), $\{\cdot \}_{N_{p}}$ denotes an average over all particles. The first slip velocity appearing in (3.14), ( $\boldsymbol{V}-\boldsymbol{u}_{p}^{(f)}$ ), accounts for the difference between the particle and fluid velocity at the location of the particle while the second slip velocity, ( $\tilde{\boldsymbol{u}}_{p}^{(f)}-\boldsymbol{V}$ ), involves the undisturbed fluid velocity arising due to the assumption of the Stokes drag model. Note that in the limit where the particle size to grid spacing is much smaller than unity $\unicode[STIX]{x1D6EC}\ll 1$ , then the undisturbed fluid velocity $\tilde{\boldsymbol{u}}_{p}^{(f)}$ evaluated at the particle location is approximately equal to the disturbed fluid velocity $\boldsymbol{u}_{p}^{(f)}$ which would be found by interpolating to the particle location from the surrounding grid points. In this limit, equation (3.14) is a quadratic form in the continuous sense, but the discrete version of this source term will only be numerically consistent if the operator used for interpolating the gas velocity to the particle has the same weights as those used to project the resulting force back to the fluid grid (Sundaram & Collins Reference Sundaram and Collins1996). Since this observation, symmetric interpolation–projection has become common practice in the literature. However, if the undisturbed fluid velocity is taken into account, as has been argued elsewhere (Horwitz & Mani Reference Horwitz and Mani2016, Reference Horwitz and Mani2018; Horwitz et al. Reference Horwitz, Rahmani, Geraci, Banko and Mani2016), the symmetry in the quadratic form disappears, so that there is no longer a reason to assume the same weights for interpolation as projection. Rather, in the authors’ opinion, it is best to compute the assumed force model as accurately as possible (which requires a good estimate of the undisturbed fluid velocity), and to then explore the consequences of that procedure. In other words, how good is the agreement between the point-particle statistics where a drag model is assumed (but computed accurately), and a particle-resolved simulation, where no force model is assumed? This is the question we explore in the next section. We conclude this section by noting that the additional particle dissipation term for an arbitrary force (acceleration) model will take the form:

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{pp}^{(f)}+\unicode[STIX]{x1D6F1}_{pp}^{(p)}=\unicode[STIX]{x1D70C}^{(p)}\unicode[STIX]{x1D719}\{(\boldsymbol{V}-\boldsymbol{u}_{p}^{(f)})\boldsymbol{\cdot }\boldsymbol{A}_{model}\}_{N_{p}}.\end{eqnarray}$$

4.1 Grid resolution

Following the recipe of Pope (Reference Pope2000), the grid resolution requirement for PR-DNS can be expressed as

(4.1) $$\begin{eqnarray}N^{3}=\left(\frac{{\mathcal{L}}}{\unicode[STIX]{x0394}x}\right)^{3}=\left(\frac{{\mathcal{L}}}{L_{11}}\frac{L_{11}}{L}\frac{L}{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x1D702}}{d_{p}}\frac{d_{p}}{\unicode[STIX]{x0394}x}\right)^{3},\end{eqnarray}$$

where ${\mathcal{L}}$ is the domain size, and $\unicode[STIX]{x0394}x$ is grid spacing. In addition, $L_{11}$ and $L$ are, respectively, the turbulent integral and large-eddy length scales (Xu & Subramaniam (Reference Xu and Subramaniam2010) reported a similar expression for the grid resolution requirement). Assuming that the length scale ratios $L_{11}/L$ and ${\mathcal{L}}/L_{11}$ are constant, equation (4.1) can be simplified by exploiting the turbulence scaling expressions $L/\unicode[STIX]{x1D702}=Re_{L}^{3/4}$ and $Re_{\unicode[STIX]{x1D706}}=\sqrt{20Re_{L}/3}$ (Pope Reference Pope2000), that is

(4.2) $$\begin{eqnarray}N^{3}\sim Re_{\unicode[STIX]{x1D706}}^{9/2}\left(\frac{\unicode[STIX]{x1D702}}{d_{p}}\right)^{3}\left(\frac{d_{p}}{\unicode[STIX]{x0394}x}\right)^{3}.\end{eqnarray}$$

The above expression reveals that unlike single-phase turbulence, the grid resolution requirement is strongly affected by the particle size and particle grid resolution. The computational cost increases as $d_{p}/\unicode[STIX]{x1D702}$ decreases. In addition, the overall grid resolution is proportional to the cube of the grid spacing across each particle $d_{p}/\unicode[STIX]{x0394}x$ . The immense computational demand imposed by the grid resolution requirement limits applicability of the PR-DNS of Kolmogorov-size particles to low-to-moderate Taylor microscale Reynolds numbers. Therefore, in the present study we choose $Re_{\unicode[STIX]{x1D706}}\approx 27$ , $d_{p}/\unicode[STIX]{x1D702}=1$ and $d_{p}/\unicode[STIX]{x0394}x=12$ for the PR-DNS. Based on these parameters, we estimated the size of the domain to be ${\mathcal{L}}/d_{p}=96$ which leads to $N^{3}=1152^{3}$ grid points for the PR-DNS. For the PP-DNS, the same parameters are considered, except that the grid spacing is $\unicode[STIX]{x0394}x=\unicode[STIX]{x1D702}$ which gives rise to $N^{3}=96^{3}$ grid points. It is worth mentioning that $\unicode[STIX]{x1D705}_{max}\unicode[STIX]{x1D702}_{0}$ ( $\unicode[STIX]{x1D705}_{max}$ being the highest wavenumber in the domain) for PR-DNS and PP-DNS are, respectively, 37.7 and 3.14. One last note concerns the PR-DNS grid resolution per particle diameter $d_{p}/\unicode[STIX]{x0394}x=12$ , which has been verified to be adequate as shown in appendix C. We also find this grid resolution to be consistent with similar particle-resolved studies (Schneiders et al. Reference Schneiders, Meinke and Schroder2016; Ozel et al. Reference Ozel, de Mottaa, Abbas, Fedea, Masbernat, Vincent, Estivalezes and Simonin2017).

4.2 Flow initialization

Owing to resolution requirements for the particle-resolved simulations, the initial Reynolds number is relatively low, $Re_{\unicode[STIX]{x1D706}}\approx 27$ . The initial condition for each simulation is a divergence-free random field based on the method of (Rogallo Reference Rogallo1981) whose energy spectrum obeys Pope’s model spectrum (Pope Reference Pope2000). The energy spectrum function has the following form:

(4.3) $$\begin{eqnarray}E(\unicode[STIX]{x1D705})=C\unicode[STIX]{x1D700}^{2/3}\unicode[STIX]{x1D705}^{-5/3}f_{L}(\unicode[STIX]{x1D705}L)f_{\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D705}\unicode[STIX]{x1D702}),\end{eqnarray}$$

where $C$ is a model constant, $\unicode[STIX]{x1D705}$ is the wavenumber, $\unicode[STIX]{x1D700}$ is the dissipation and $L$ and $\unicode[STIX]{x1D702}$ are large eddy and Kolmogorov length scales, respectively. The functions $f_{L}$ and $f_{\unicode[STIX]{x1D702}}$ determine the shape of energy-containing and dissipative range of the energy spectrum and have the form

(4.4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle f_{L}(\unicode[STIX]{x1D705}L)=\left(\frac{\unicode[STIX]{x1D705}L}{[(\unicode[STIX]{x1D705}L)^{2}+c_{L}]^{1/2}}\right)^{5/3+p_{0}}\\ \displaystyle f_{\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D705}\unicode[STIX]{x1D702})=\text{exp}\{-\unicode[STIX]{x1D6FD}\{[(\unicode[STIX]{x1D705}\unicode[STIX]{x1D702})^{4}+c_{\unicode[STIX]{x1D702}}^{4}]^{1/4}-c_{\unicode[STIX]{x1D702}}\}\}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The model constants $p_{0}$ and $\unicode[STIX]{x1D6FD}$ are the same as those suggested by Pope (Reference Pope2000), while $c_{L}$ and $c_{\unicode[STIX]{x1D702}}$ are determined by the multi-variable Newton–Raphson method such that the energy and dissipation rate from the spectrum matches those required for the chosen $Re_{\unicode[STIX]{x1D706}}$ with $d_{p}=\unicode[STIX]{x1D702}_{0}$ . The different resolutions and grid configurations (collocated for the particle-resolved and staggered for the point-particle algorithms) means the respective simulations are using different initial conditions, which implies each initial condition should be interpreted as different realizations of the same turbulence. When the parameters in the model spectrum were chosen to match integral length scales for the two initial conditions, both codes resulted in nearly identical evolution of fluid kinetic energy and dissipation rate in unladen simulations. Based on the aforementioned observations, we do not believe there is a sensitivity in the initial conditions that can account for any differences in the results of the particle-laden simulations. The particles are seeded randomly with uniform probability in the box with the initial particle velocity taken as the local fluid velocity at the particle location. The system is then allowed to evolve via fluid and particle equations. Since the sample paths taken by each particle will be different between the two simulation methodologies, the PR-DNS and PP-DNS studies will be said to be statistically equivalent in a weak sense if moments of the particle and fluid distributions can be matched between the two simulation approaches. The number of particles tracked limits the sampling error. While the initial condition does not correspond to a physical turbulent field, the skewness of the longitudinal velocity derivative reaches a physical value of $-0.5$ within $O(0.1)$ eddy turnover times so that all salient features in the decay are physical.

In this study, we examine particles whose diameter is equal to the initial Kolmogorov scale, $d_{p}/\unicode[STIX]{x1D702}_{0}=1$ . We also examine two Stokes numbers $St_{0}=1$ and $St_{0}=100$ where the particle Stokes number is defined as $St_{0}=(1/18)(\unicode[STIX]{x1D70C}^{(p)}/\unicode[STIX]{x1D70C}^{(f)})(d_{p}/\unicode[STIX]{x1D702}_{0})^{2}$ . These Stokes numbers are defined based on the initial Kolmogorov time scale. The particle size is held constant between the simulations, so that the increase in Stokes number is accomplished by increasing the particle to fluid density ratio from 18 to 1800 between the two simulation cases. The volume fraction is kept small $(\unicode[STIX]{x1D719}=0.001)$ to ensure particle–particle interactions are a second-order effect. With ${\mathcal{L}}/d_{p}=96$ and $\unicode[STIX]{x1D719}=0.001$ , the total number of particles in the domain is:

(4.5) $$\begin{eqnarray}N_{p}=\unicode[STIX]{x1D719}\frac{6}{\unicode[STIX]{x03C0}}\left\langle \frac{{\mathcal{L}}}{d_{p}}\right\rangle ^{3}=1689.\end{eqnarray}$$

In addition, the two simulation cases correspond to two mass loading ratios, $\unicode[STIX]{x1D6F7}_{m}=(\unicode[STIX]{x1D70C}^{(p)}/\unicode[STIX]{x1D70C}^{(f)})\unicode[STIX]{x1D719}=0.018$ and 1.8.

Figure 1.

Comparison of (a) fluid kinetic energy, (b) fluid dissipation rate and (c) particle kinetic energy versus time from the particle-resolved simulation with different point-particle models for the case $St_{0}=1$ . The fluid dissipation rate for the particle-resolved simulation, the right-hand side of (2.7), is compared against the sum of all terms on the right-hand side of (2.13).

4.3 Fluid and particle statistics for $St_{0}=1$

The evolution of the fluid energy, dissipation rate, and particle kinetic energy are shown in figure 1. Two point-particle models have been employed; in each, the particles are assumed to obey Stokes drag, however, one of the models accounts for the undisturbed fluid velocity (St corrected) and the other does not (St uncorrected). In comparing statistics, it is clear that both point-particle models are in excellent agreement with the particle-resolved simulation. The agreement between the point-particle simulations and the particle-resolved simulation with regard to fluid energy and dissipation rate is not surprising given the low mass loading. This means the particle–fluid system is effectively one-way coupled. More surprising however is the agreement in particle kinetic energy predicted by the point-particle models compared with the particle-resolved simulation. In the present study, since the particles are initialized with the local fluid velocity, and because the Stokes number is relatively small, the particles never acquire a significant enough slip velocity for finite Reynolds number corrections to steady drag to become important (as shall be contrasted with the results presented in the next section). This observation also explains why not accounting for the undisturbed fluid velocity has little effect on prediction of particle kinetic energy; at this Stokes number, particles act as tracers to large eddies which dominate contributions to the kinetic energy. This observation also helps explain some of the results presented in Subramaniam et al. (Reference Subramaniam, Mehrabadi, Horwitz and Mani2014) which compared point particles and resolved particles at $St_{0}=1$ in decaying isotropic turbulence where the particles were initialized at rest. Though fluid energy was matched well between the two simulation methodologies, the fluid dissipation rate and particle kinetic energies showed disagreement. Therefore it seems the manner in which particles are initialized plays an important role in determining whether inertial and/or unsteady contributions to drag will be significant, as well as how these forces ultimately correlate to other particle and fluid statistics.

4.4 Fluid and particle statistics for $St_{0}=100$

In this section we explore a set-up which is of greater relevance to practical gas–solid flows. Owing to the low mass loading for the $St_{0}=1$ case, there was small particle–turbulence interaction. In contrast, the present case exhibits significant coupling between the dispersed and carrier phases. Here the mass loading is increased by a factor of $100$ to $\unicode[STIX]{x1D6F7}_{m}=1.8$ . This is accomplished by keeping the particle diameter to initial Kolmogorov scale fixed at $d_{p}/\unicode[STIX]{x1D702}_{0}=1$ and increasing the density ratio to $\unicode[STIX]{x1D70C}^{(p)}/\unicode[STIX]{x1D70C}^{(f)}=1800$ .

Fluid kinetic energy and dissipation rate for the present case are shown in figure 2. In contrast with the previous case $(St_{0}=1)$ , it is clear in the present case $(St_{0}=100)$ that significant turbulence modulation is experienced in the fluid. Here we are showing a comparison of the particle-resolved simulation with four point-particle simulations; two of the point-particle simulations employ Stokes drag, with the difference being that one accounts for the undisturbed fluid velocity (St corrected) while the other does not (St uncorrected). The two other point-particle simulations employ the Schiller–Naumann Reynolds-number corrected steady drag formula; the difference in the Schiller–Naumann implementations is again that in one procedure the undisturbed fluid velocity is estimated (SN corrected), while in the other the disturbed fluid velocity is used (SN uncorrected).

Figure 2.

Comparison of (a) fluid kinetic energy and (b) fluid dissipation rate versus time from particle-resolved simulation with different point-particle models for the case $St_{0}=100$ . The fluid dissipation rate for the particle-resolved simulation, the right-hand side of (2.7), is compared against the sum of all terms on the right-hand side of (2.13). Both figures share the same line legends.

In the present case a clear difference is observed among the different point-particle simulations compared with the particle-resolved simulation. The point-particle simulations where the particles obey the Schiller–Naumann drag correlation whose slip velocity is calculated as the difference between the undisturbed fluid velocity and the particle velocity, show excellent agreement in both the fluid kinetic energy and dissipation rate compared with the particle-resolved simulation. The greatest discrepancy is seen between the particle-resolved simulation and the St uncorrected simulation. Interestingly, the St corrected implementation shows better agreement with the particle-resolved simulation than the SN uncorrected scheme. This is an interesting observation which indicates that, at least for the present set-up, it is more important to get the undisturbed fluid velocity right than to use the most applicable drag correlation. Notwithstanding, accounting for both the undisturbed fluid velocity in the drag correlation, and finite Reynolds number corrections to the drag formulation are necessary for agreement between the particle-resolved and point-particle simulation. Despite the particle size being borderline for the importance of unsteady effects (Ling et al. Reference Ling, Parmar and Balachandar2013; Olivieri et al. Reference Olivieri, Picano, Sardina, Iudicone and Brandt2014; Daitche Reference Daitche2015), good agreement between the point-particle and particle-resolved simulation is found without explicitly accounting for these effects in the point-particle model. We have also conducted tests of particles moving in simple flows using the analytical solutions provided by Coimbra & Rangel (Reference Coimbra and Rangel1998) and those results provide justification that the history term can be neglected for the present case.

Figure 3.

Comparison of particle-resolved simulation against different point-particle models for the case $St_{0}=100$ , (a) particle kinetic energy, and (b) solid-phase interphase kinetic energy transfer versus time.

The behaviour of the point-particle curves observed in figure 2 is also observed in the particle kinetic energy and dissipation rate shown in figure 3. The SN corrected point-particle scheme shows excellent agreement for particle kinetic energy compared with the PR-DNS. As before, St corrected is more accurate than SN uncorrected, and St uncorrected is the least accurate. To ensure the computation of $k^{(p)}$ is not subject to substantial statistical variability due to the finite number of particles ( $N_{p}=1689$ ) in the computational domain, we have calculated $95\,\%$ confidence intervals on particle kinetic energy at different times for PR-DNS and PP-DNS, as shown in table 1. The low statistical variability for both PR-DNS and PP-DNS results provide confidence that the observed differences in $k^{(p)}$ , especially at later times when particles have lost memory of their initial conditions, are primarily owing to the model form of drag and its implementation rather than sampling error.

Table 1.

The $95\,\%$ confidence intervals on particle kinetic energy at different times for PR-DNS and PP-DNS cases for $St_{0}=100$ . These quantities are normalized by the initial fluid kinetic energy $k_{0}^{(f)}$ .

Given that the particle kinetic energy in the particle-resolved simulation and SN corrected simulation are in good agreement, the source term that creates particle kinetic energy, namely the covariance of particle velocity and acceleration (see equations (2.5) and (2.12)), must also be in good agreement. Indeed, figure 3(b) reveals excellent agreement between the SN corrected scheme and the particle-resolved simulation. The SN corrected scheme predicts a slightly higher peak in particle velocity–acceleration covariance. This is consistent with the particle kinetic energy being slightly lower at late times for the SN corrected scheme compared with the PR-DNS. The well ordering of the point-particle schemes is also found in the particle kinetic energy source term, namely that St corrected is better than SN uncorrected. The least accurate, as before, is the St uncorrected scheme which under-predicts the peak particle velocity–acceleration covariance as seen in the PR-DNS by almost 50 %.

Thus far, we have demonstrated validation of the point-particle method via the SN corrected implementation for Eulerian integral quantities (fluid kinetic energy and dissipation rate) as well as Lagrangian integral quantities (particle kinetic energy and dissipation rate). Finer appraisal of the accuracy of the point-particle schemes can be studied by examining the distribution of particle acceleration events. Probability density functions (PDFs) of particle acceleration for each of the schemes at three times are shown in figure 4. The corresponding sample mean and standard deviation of particle acceleration are presented in table 2. In examination of the PDF shapes and their moments, the SN corrected particle acceleration PDFs are in good agreement with the acceleration PDFs derived from the particle-resolved simulation at the same non-dimensional times. It is clear that the particle acceleration PDF is skewed towards higher accelerations in the particle-resolved simulation than any of the point-particle methods except for SN corrected. The SN corrected scheme shows good agreement in predicting the PDF mode and maximum particle acceleration, although the SN corrected scheme slightly over-predicts the mean and standard deviation of acceleration compared with the PR-DNS. Though the other point-particle schemes have the correct qualitative shape, they considerably underestimate the mean and standard deviation of acceleration as well as the maximum acceleration compared with those predicted by the PR-DNS simulation. Again the same ordering of the point-particle schemes is apparent, with St uncorrected showing the lowest level of agreement with the PR-DNS, followed by SN uncorrected. The St corrected scheme does a good job at predicting the location of the peak in the PR-DNS particle acceleration PDF, however, the St corrected PDF falls off too quickly, underpredicting the lower probability high acceleration events at early ( $t\unicode[STIX]{x1D700}_{0}^{(f)}/k_{0}^{(f)}=0.54$ ) and intermediate ( $t\unicode[STIX]{x1D700}_{0}^{(f)}/k_{0}^{(f)}=2.70$ ) times.

Figure 4.

Particle acceleration PDFs at different times from particle-resolved simulation compared with different point-particle models.

Table 2.

Mean ( $|\boldsymbol{a}_{p}|$ ) and standard variation ( $\unicode[STIX]{x1D70E}_{|\boldsymbol{a}_{p}|}$ ) of the particle acceleration magnitude at different times for the PR-DNS compared with the PP-DNS cases for $St_{0}=100$ . These quantities are normalized by the initial Kolmogorov velocity scale $u_{\unicode[STIX]{x1D702}}$ and time scale $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ .