Particle-resolved direct numerical simulation of homogeneous isotropic turbulence modified by small fixed spheres

Abstract

A statistically stationary homogeneous isotropic turbulent flow modified by 64 small fixed non-Stokesian spherical particles is considered. The particle diameter is approximately twice the Kolmogorov length scale, while the particle volume fraction is 0.001. The Taylor Reynolds number of the corresponding unladen flow is 32. The particle-laden flow has been obtained by a direct numerical simulation based on a discretization of the incompressible Navier–Stokes equations on 64 spherical grids overset on a Cartesian grid. The global (space- and time-averaged) turbulence kinetic energy is attenuated by approximately 9 %, which is less than expected. The turbulence dissipation rate on the surfaces of the particles is enhanced by two orders of magnitude. More than 5 % of the total dissipation occurs in only 0.1 % of the flow domain. The budget of the turbulence kinetic energy has been computed, as a function of the distance to the nearest particle centre. The budget illustrates how energy relatively far away from particles is transported towards the surfaces of the particles, where it is dissipated by the (locally enhanced) turbulence dissipation rate. The energy flux towards the particles is dominated by turbulent transport relatively far away from particles, by viscous diffusion very close to the particles, and by pressure diffusion in a significant region in between. The skewness and flatness factors of the pressure, velocity and velocity gradient have also been computed. The global flatness factor of the longitudinal velocity gradient, which characterizes the intermittency of small scales, is enhanced by a factor of six. In addition, several point-particle simulations based on the Schiller–Naumann drag correlation have been performed. A posteriori tests of the point-particle simulations, comparisons in which the particle-resolved results are regarded as the standard, show that, in this case, the point-particle model captures both the turbulence attenuation and the fraction of the turbulence dissipation rate due to particles reasonably well, provided the (arbitrary) size of the fluid volume over which each particle force is distributed is suitably chosen.

1 Introduction

It is intriguing that a low volume fraction of small particles with a high Stokes number is able to dampen turbulence significantly (Tsuji, Morikawa & Shiomi 1984; Gore & Crowe 1989; Hetseroni 1989). The phenomenon has been observed experimentally in pipe flows (Tsuji et al. 1984), channel flows (Kulick, Fessler & Eaton 1994; Kussin & Sommerfeld 2002) and stationary homogeneous turbulence (Hwang & Eaton 2006; Tanaka & Eaton 2010). These works clearly show that the turbulence kinetic energy of the carrier phase of a turbulent flow laden with small particles is lower than in the corresponding unladen flow. Compared to the amount of experimental data on macroscopic turbulence modification, there are only a few experimental results on the modification of the small-scale turbulence. High-resolution experiments of turbulent air flow laden with small solid particles were performed by Tanaka & Eaton (2010), who refined the experiments performed by Hwang & Eaton (2006), and observed that the turbulence dissipation rate was augmented in a roughly spherical region around the particles. The experiments reported in these two papers were performed for stationary homogeneous isotropic turbulence, forced by synthetic jet actuators (woofers) driven with sine waves with random frequencies.

A popular method to simulate particle-laden flows is the point-particle technique. In this approach, the scales in the direct vicinity of the particles are not resolved, but a simplified expression, such as Stokes’ law or the Schiller–Naumann correlation, models the forces exerted on the particles. The point-particle method or a related method has also been applied in simulations of particle-laden isotropic turbulence (see e.g. Squires & Eaton 1990; Elghobashi & Truesdell 1993; Ferrante & Elghobashi 2003). These simulations confirmed that small particles with large Stokes number tend to attenuate the turbulence. The attenuation was attributed to the particle-induced dissipation term, a sink term that appears explicitly in the equation of the turbulence kinetic energy derived from the point-particle or related formulation. In particle-laden pipe and channel flows, the turbulence attenuation can be further enhanced by a second effect, which is non-uniformity of the mean particle force (Vreman 2007, 2015).

In order to gain a more fundamental understanding of multiphase flows, the recently emerged particle-resolved direct numerical simulation technique (PR-DNS) looks very promising (Balachandar & Eaton 2010; Tryggvason, Scardovelli & Zaleski 2011; Prosperetti 2015). Bagchi & Balachandar (2003) and Burton & Eaton (2005) used PR-DNS to compute the force on a single fixed particle in decaying homogeneous isotropic turbulence. In both works, a body-fitted spherical grid was used, while in the latter work the body-fitted grid was overset on a Cartesian grid. Burton & Eaton (2005) included results for the local turbulence dissipation rate around a particle of $d_{p}\approx 2{\it\eta}$ , where $d_{p}$ is the particle diameter and ${\it\eta}$ the Kolmogorov length scale. However, since the particle volume fraction ${\it\alpha}$ was negligible (approximately $10^{-7}$ ), no significant modifications of the global turbulence kinetic energy ( $\overline{K}$ ) and global turbulence dissipation rate ( $\overline{{\it\epsilon}}$ ) were found.

More recently, PR-DNS studies of turbulent flows with multiple solid fixed and moving spherical particles have occurred in the literature: studies of turbulent channel flow (Uhlmann 2008; Picano, Breugem & Brandt 2015) and homogeneous isotropic turbulence (Lucci, Ferrante & Elghobashi 2010; Mehrabadi et al. 2015). In all these cases, the flow was not dilute ( ${\it\alpha}\geqslant 0.01$ ) and the particle diameter was relatively large (approximately $10$ wall units in the turbulent channel flows and at least $5{\it\eta}$ in the homogeneous isotropic cases). Furthermore, in all these cases the immersed boundary method was used, in which the no-slip boundary conditions on the surfaces of the particles were approximated on a uniform Cartesian grid, which did not conform to the particle bodies. This method is powerful; Picano et al. (2015) used it to simulate turbulence attenuation in plane channel flow with 10 000 moving neutrally buoyant spheres. Another method in which the grid is Cartesian and does not conform to the particle bodies is the lattice Boltzmann method, which has also been applied to particle-laden homogeneous isotropic turbulence (ten Cate et al. 2004), for ${\it\alpha}\geqslant 0.02$ and $d_{p}\approx 8{\it\eta}$ . Both the lattice Boltzmann method and the immersed boundary method have also been used to validate or improve drag and heat transfer correlations (Hill, Koch & Ladd 2001; Tenneti, Garg & Subramaniam 2011; Deen et al. 2012; Tang et al. 2014). Finally, the PHYSALIS method is mentioned, in which local analytical representations are used for the flow in the regions between particles and adjacent nodes of the Cartesian grid (Takagi, Oguz & Prosperetti 2003). This method has been used to simulate the flow around nine fixed particles with $d_{p}\approx 8{\it\eta}$ in decaying homogeneous isotropic turbulence (Botto & Prosperetti 2012).

According to the literature, there is a strong indication that point-particle models underpredict turbulence attenuation (Hwang & Eaton 2006; Eaton 2009; Balachandar & Eaton 2010). Hwang & Eaton (2006) concluded that the extra dissipation caused by particles with $d_{p}\approx {\it\eta}$ was greatly underestimated by conventional models. Furthermore, Burton & Eaton (2005) found that the instantaneous error of the modelled particle force varied between 15 and 30 % for $d_{p}\approx 2{\it\eta}$ . Indeed, point-particle models are formally not justified if the condition that the particle diameter $d_{p}\ll {\it\eta}$ is not satisfied (Lucci et al. 2010). Nonetheless, point-particle methods are frequently used to simulate particle-laden flows for $d_{p}\geqslant {\it\eta}$ . Because of the documented discrepancies between simulated and experimentally observed turbulence attenuation, Eaton (2009) remarked that ‘it is incumbent on the developers of these codes to prove that the models are valid, either through fully resolved simulations, or direct comparison to experiments’. The present paper shows such a validation, for one specific case.

The objective of the present work is to perform a study of stationary homogeneous isotropic turbulence modified by small solid particles ( $d_{p}\approx 2{\it\eta}$ ) at low volume fraction ( ${\it\alpha}\approx 0.001$ ), based on simulations in which all relevant scales of the flow are well resolved. This objective prompts several research questions. (1) How are basic quantities such as turbulence kinetic energy and the turbulence dissipation rate modified as a function of the distance to the nearest particle centre? (2) How much attenuation of the global (space- and time-averaged) turbulence kinetic energy is predicted by a particle-resolved simulation, and is the degree of attenuation comparable to what is known from experiments in the literature? (3) What does the budget of turbulence kinetic energy as a function of the distance to the nearest particle look like, and what does it teach us about the mechanics of turbulence attenuation? (4) How do small particles modify the statistical properties of the small-scale turbulence, for example, the skewness and flatness (intermittency) factors of individual velocity derivatives? (5) Is the point-particle model based on the Schiller–Naumann correlation able to predict the turbulence attenuation and the extra dissipation produced by particles?

In view of the objective and related research questions, it is necessary to perform a direct numerical simulation (DNS) that resolves the sub-Kolmogorov structures near small particles in a dilute flow. This means that the smallest grid size needs to be much smaller than the Kolmogorov length scale (approximately 13 times smaller in the present case). The immersed boundary method, as described above, is an efficient tool for particle-resolved simulation of turbulent flows with relatively large particles at relatively high particle volume fraction. The uniform Cartesian grid (in combination with at least two periodic directions) brings the advantage that the pressure Poisson equation and implicit discretizations of the viscous terms can be solved directly using fast Fourier transforms, which are much more efficient than iterative procedures. This advantage cannot be exploited for the present investigation since sub-Kolmogorov scales near particles need to be resolved in a flow with a low particle volume fraction. For such a case, the use of a uniform grid would become a heavy burden, since it would dramatically increase the number of grid points required. Therefore, the immersed boundary method is not chosen for the present investigation. Instead, the overset grid approach used by Burton & Eaton (2005) is adopted and extended to multiple particles: around each particle, a body-fitted (spherical) staggered non-uniform mesh is used, overset on a staggered Cartesian mesh. A finite difference technique that uses multiple overset curvilinear meshes has also been called a composite mesh difference technique (Starius 1977), a chimera grid-embedding technique (Benek, Buning & Steger 1985) and a composite overlapping grid technique (Chesshire & Henshaw 1990; Henshaw 1994). With the overset technique, not only can the number of grid points be limited, but also its body-fitted property facilitates a straightforward computation of the turbulent dissipation rate on the particle surfaces. A drawback is that the pressure Poisson equation can only be solved by iterative techniques, which are computationally much more expensive than fast direct techniques. It is remarked that, owing to the requirement of a locally high spatial resolution, the temporal resolution needs to be high as well or one has to resort to fully implicit methods.

To alleviate the burdens of technical complexity and computational cost to some extent, the Taylor Reynolds number is limited to 32, the number of particles is limited to 64, and the particles are fixed and ordered in a cubic pattern. Thus we consider stationary homogeneous isotropic turbulence modified by an array of stagnant spherical particles, solved by means of body-fitted PR-DNS. This flow case lies in the gap between, on the one hand, body-fitted PR-DNS of turbulent flow around a single fixed small particle (Bagchi & Balachandar 2003; Burton & Eaton 2005) and, on the other, immersed boundary PR-DNS of turbulent flow around many freely moving coarser particles (Uhlmann 2008; Lucci et al. 2010; Mehrabadi et al. 2015; Picano et al. 2015). Burton & Eaton (2005) remarked that the use of fixed particles in homogeneous turbulence could be justified as relevant for turbulence modification by heavy moving particles or particles in microgravity. Tanaka & Eaton (2010) concluded that the small heavy particles in their experiments of stationary homogeneous turbulence behaved as localized dampers of the fluid motion, like grids act in turbulence. From this perspective, it is a logical choice to perform a DNS of the idealized flow of homogeneous spatially periodic turbulence modified by an array of fixed spherical particles, in order to investigate in further detail how turbulence is modified by these spherical dampers. Thus, the present purpose is not to reproduce a practical flow as accurately as possible, but to obtain a numerically accurate solution of the Navier–Stokes equations that describes turbulent flow around small particles, a solution that provides fundamental understanding of the mechanics and small-scale behaviour of turbulent flow around small heavy spherical particles at low volume fractions.

The structure of the present paper is as follows. The governing equations and numerical method are described in § 2 and appendix A. The simulation cases are defined in § 3. The simulation results are shown in § 4, which has five subsections, one for each research question formulated above. Finally, the conclusions are summarized in § 5.

2 Particle-resolved simulation method

The numerical method used to simulate the flow around the spherical particles is a second-order staggered overset grid method, inspired by and similar to the method used by Burton & Eaton (2002, 2005) to simulate turbulent flow around a single spherical particle. However, the details are not the same. An important difference is that in the present method the pressure Poisson equations for the different domains are assembled into one large matrix that is iteratively solved, while a block-matrix technique was used in the references just mentioned. Another difference is that the communication among the Poisson equations is arranged via interpolation of the pressure instead of the pressure gradient. These two choices make the method more convenient for multiple particles. Furthermore, we use an alternative mathematically equivalent formulation of the viscous terms in spherical coordinates, which leads to a different spatial discretization. The treatment of the apparent singularities in the spherical frames of reference is also different. Finally, we use a new procedure to accurately compute integrals over the entire flow domain for variables that are defined on overlapping grids. Thus, in order to make the results of this paper reproducible, a complete specification of the numerical method is required.

The coordinate systems and required relations are introduced in § 2.1. The form of the Navier–Stokes equations used is defined in § 2.2. The overlapping grids are defined in § 2.3. The details of the staggered second-order discretization and third-order interpolation method are presented in appendix A. Definitions and implementations of the statistical averaging operators are presented in § 2.4. The results of four test cases are presented in § 2.5.

The summation convention over repeated indices is assumed, except when mentioned otherwise and except when the index is $t$ , $r$ , ${\it\theta}$ or ${\it\phi}$ . Partial derivatives are denoted using the comma notation, for example $u_{j,i}=\partial u_{j}/\partial x_{i}$ , $u_{j,t}=\partial u_{j}/\partial t$ and $u_{{\it\theta},r}=\partial u_{{\it\theta}}/\partial r$ .

2.1 Coordinate systems

The Cartesian base vectors are $\boldsymbol{e}_{1}$ , $\boldsymbol{e}_{2}$ and $\boldsymbol{e}_{3}$ . The Cartesian position vector is $\boldsymbol{x}=[x_{1},x_{2},x_{3}]^{\text{T}}$ , and the Cartesian velocity vector is $\boldsymbol{u}=[u_{1},u_{2},u_{3}]^{\text{T}}$ . We consider a rigid spherical particle, centred at position $\boldsymbol{x}^{p}$ . The particle velocity, denoted by $\boldsymbol{u}^{p}$ , is zero in the simulations, but, to make the formulation more general, the particle velocity is carried into the equations. The spherical coordinates around the particle are given by

(2.1a-c ) $$\begin{eqnarray}r=|\boldsymbol{x}-\boldsymbol{x}^{p}|,\quad {\it\theta}=\text{arccos}\left(\frac{x_{3}-x_{3}^{p}}{r}\right),\quad {\it\phi}=\text{atan2}(x_{2}-x_{2}^{p},x_{1}-x_{1}^{p}),\end{eqnarray}$$

where $\text{atan2}(y,x)$ is the function that provides the argument of the complex number with real part $x$ and complex part $y$ , such that $-{\rm\pi}<{\it\phi}\leqslant {\rm\pi}$ . Here $r$ denotes the radial, ${\it\theta}$ the polar and ${\it\phi}$ the azimuthal coordinate. The inverse transformation is

(2.2) $$\begin{eqnarray}\boldsymbol{x}=\boldsymbol{x}^{p}+\left[\begin{array}{@{}c@{}}r\sin {\it\theta}\cos {\it\phi}\\ r\sin {\it\theta}\sin {\it\phi}\\ r\cos {\it\theta}\end{array}\right].\end{eqnarray}$$

The standard base vectors of the spherical coordinate system are $\boldsymbol{e}_{r}$ , $\boldsymbol{e}_{{\it\theta}}$ and $\boldsymbol{e}_{{\it\phi}}$ , while the spherical velocity components are denoted by $u_{r}$ , $u_{{\it\theta}}$ and $u_{{\it\phi}}$ .

The base vector relations for an orthogonal coordinate transformation from the basis $\{\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}\}$ to the basis $\{\tilde{\boldsymbol{e}}_{1},\tilde{\boldsymbol{e}}_{2},\tilde{\boldsymbol{e}}_{3}\}$ are given by

(2.3a,b ) $$\begin{eqnarray}\tilde{\boldsymbol{e}}_{i}=\unicode[STIX]{x1D608}_{ij}\boldsymbol{e}_{j},\quad \boldsymbol{e}_{i}=\unicode[STIX]{x1D608}_{ji}\tilde{\boldsymbol{e}}_{j},\end{eqnarray}$$

where $\unicode[STIX]{x1D63C}$ is an orthogonal matrix. The transformation to the spherical coordinate system is obtained by defining $\tilde{\boldsymbol{e}}_{1}=\boldsymbol{e}_{r}$ , $\tilde{\boldsymbol{e}}_{2}=\boldsymbol{e}_{{\it\theta}}$ and $\tilde{\boldsymbol{e}}_{3}=\boldsymbol{e}_{{\it\phi}}$ , specified by

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D63C}=[\boldsymbol{e}_{r},\boldsymbol{e}_{{\it\theta}},\boldsymbol{e}_{{\it\phi}}]^{\text{T}}=\left[\begin{array}{@{}ccc@{}}\sin {\it\theta}\cos {\it\phi} & \sin {\it\theta}\sin {\it\phi} & \cos {\it\theta}\\ \cos {\it\theta}\cos {\it\phi} & \cos {\it\theta}\sin {\it\phi} & -\!\sin {\it\theta}\\ -\!\sin {\it\phi} & \cos {\it\phi} & 0\end{array}\right].\end{eqnarray}$$

The inverse of the coordinate transformation (2.1) becomes $x_{j}=x_{j}^{p}+r\unicode[STIX]{x1D608}_{1j}$ .

For the velocity components in the spherical coordinate system, we define the vector $\tilde{\boldsymbol{u}}=[\tilde{u} _{1},\tilde{u} _{2},\tilde{u} _{3}]^{\text{T}}=[u_{r},u_{{\it\theta}},u_{{\it\phi}}]^{\text{T}}$ . Since $\boldsymbol{u}=u_{i}\boldsymbol{e}_{i}=\boldsymbol{u}^{p}+\tilde{u} _{i}\tilde{\boldsymbol{e}}_{i}$ , we obtain the velocity transformations

(2.5a,b ) $$\begin{eqnarray}\tilde{u} _{j}=\unicode[STIX]{x1D608}_{ji}(u_{i}-u_{i}^{p}),\quad u_{j}=u_{j}^{p}+\unicode[STIX]{x1D608}_{ij}\tilde{u} _{i}.\end{eqnarray}$$

The gradient of the velocity is the tensor

(2.6) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{u}=u_{j,i}\boldsymbol{e}_{i}\boldsymbol{e}_{j}=\unicode[STIX]{x1D60E}_{kl}\tilde{\boldsymbol{e}}_{k}\tilde{\boldsymbol{e}}_{l},\end{eqnarray}$$

where $u_{j,i}=\partial u_{j}/\partial x_{i}$ and

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D642}=\left[\begin{array}{@{}ccc@{}}u_{r,r} & u_{{\it\theta},r} & u_{{\it\phi},r}\\ {\displaystyle \frac{u_{r,{\it\theta}}}{r}}-{\displaystyle \frac{u_{{\it\theta}}}{r}} & {\displaystyle \frac{u_{{\it\theta},{\it\theta}}}{r}}+{\displaystyle \frac{u_{r}}{r}} & {\displaystyle \frac{u_{{\it\phi},{\it\theta}}}{r}}\\ {\displaystyle \frac{u_{r,{\it\phi}}}{r\sin {\it\theta}}}-{\displaystyle \frac{u_{{\it\phi}}}{r}} & {\displaystyle \frac{u_{{\it\theta},{\it\phi}}}{r\sin {\it\theta}}}-{\displaystyle \frac{u_{{\it\phi}}\cot {\it\theta}}{r}} & {\displaystyle \frac{u_{{\it\phi},{\it\phi}}}{r\sin {\it\theta}}}+{\displaystyle \frac{u_{r}}{r}}+{\displaystyle \frac{u_{{\it\theta}}\cot {\it\theta}}{r}}\end{array}\right]\end{eqnarray}$$

(Phan-Thien 2013). Equation (2.6) implies a useful expression for the spatial velocity derivatives

(2.8) $$\begin{eqnarray}u_{j,i}=\unicode[STIX]{x1D60E}_{kl}\unicode[STIX]{x1D608}_{ki}\unicode[STIX]{x1D608}_{lj}.\end{eqnarray}$$

2.2 Governing equations

The Navier–Stokes equations read

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{,t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{u})=-\boldsymbol{{\rm\nabla}}p+{\it\nu}{\rm\nabla}^{2}\boldsymbol{u}+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$

where $t$ denotes time, $p$ the pressure divided by the constant density ${\it\rho}$ , ${\it\nu}$ the kinematic viscosity and $\boldsymbol{f}$ the forcing vector, which is a function of time and space. This is the form of the Navier–Stokes equations used in the Cartesian domain.

The form of the Navier–Stokes equations used in the spherical frame of reference of a particle is given by

(2.11) $$\begin{eqnarray}\tilde{\boldsymbol{{\rm\nabla}}}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}~=\frac{(r^{2}u_{r})_{,r}}{r^{2}}+\frac{(u_{{\it\theta}}\sin {\it\theta})_{,{\it\theta}}}{r\sin {\it\theta}}+\frac{u_{{\it\phi},{\it\phi}}}{r\sin {\it\theta}}=0,\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & & \displaystyle u_{r,t}+\frac{(r^{2}u_{r}u_{r})_{,r}}{r^{2}}+\frac{(u_{r}u_{{\it\theta}}\sin {\it\theta})_{,{\it\theta}}}{r\sin {\it\theta}}+\frac{(u_{r}u_{{\it\phi}})_{,{\it\phi}}}{r\sin {\it\theta}}-\frac{u_{{\it\theta}}^{2}+u_{{\it\phi}}^{2}}{r}\nonumber\\ \displaystyle & & \displaystyle \quad =-p_{,r}+{\it\nu}\left(\tilde{{\rm\nabla}}^{2}u_{r}+2\frac{u_{r}}{r^{2}}+2\frac{u_{r,r}}{r}\right)-\unicode[STIX]{x1D608}_{1j}u_{j,t}^{p}+\unicode[STIX]{x1D608}_{1j}\,f_{j},\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & & \displaystyle u_{{\it\theta},t}+\frac{(r^{2}u_{{\it\theta}}u_{r})_{,r}}{r^{2}}+\frac{(u_{{\it\theta}}u_{{\it\theta}}\sin {\it\theta})_{,{\it\theta}}}{r\sin {\it\theta}}+\frac{(u_{{\it\theta}}u_{{\it\phi}})_{,{\it\phi}}}{r\sin {\it\theta}}+\frac{u_{{\it\theta}}u_{r}-u_{{\it\phi}}^{2}\cot {\it\theta}}{r}\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{p_{,{\it\theta}}}{r}+{\it\nu}\left(\tilde{{\rm\nabla}}^{2}u_{{\it\theta}}+2\frac{u_{r,{\it\theta}}}{r^{2}}-\frac{u_{{\it\theta}}+2u_{{\it\phi},{\it\phi}}\cos {\it\theta}}{r^{2}\sin ^{2}{\it\theta}}\right)-\unicode[STIX]{x1D608}_{2j}u_{j,t}^{p}+\unicode[STIX]{x1D608}_{2j}\,f_{j},\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & & \displaystyle u_{{\it\phi},t}+\frac{(r^{2}u_{{\it\phi}}u_{r})_{,r}}{r^{2}}+\frac{(u_{{\it\phi}}u_{{\it\theta}}\sin {\it\theta})_{,{\it\theta}}}{r\sin {\it\theta}}+\frac{(u_{{\it\phi}}u_{{\it\phi}})_{,{\it\phi}}}{r\sin {\it\theta}}+\frac{u_{{\it\phi}}u_{r}+u_{{\it\phi}}u_{{\it\theta}}\cot {\it\theta}}{r}\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad =-\frac{p_{,{\it\phi}}}{r\sin {\it\theta}}+\,{\it\nu}\left(\tilde{{\rm\nabla}}^{2}u_{{\it\phi}}+\frac{2u_{r,{\it\phi}}\sin {\it\theta}+2u_{{\it\theta},{\it\phi}}\cos {\it\theta}-u_{{\it\phi}}}{r^{2}\sin ^{2}{\it\theta}}\right)-\unicode[STIX]{x1D608}_{3j}u_{j,t}^{p}+\unicode[STIX]{x1D608}_{3j}\,f_{j}.\qquad\end{eqnarray}$$

In these equations, the Laplace operator in spherical coordinates appears. Applied to a scalar $v$ , it is defined by

(2.15) $$\begin{eqnarray}\tilde{{\rm\nabla}}^{2}v=\frac{(r^{2}v_{,r})_{,r}}{r^{2}}+\frac{(v_{,{\it\theta}}\sin {\it\theta})_{,{\it\theta}}}{r^{2}\sin {\it\theta}}+\frac{v_{,{\it\phi}{\it\phi}}}{r^{2}\sin ^{2}{\it\theta}}.\end{eqnarray}$$

The divergence operator in spherical coordinates, denoted by $\tilde{\boldsymbol{{\rm\nabla}}}\boldsymbol{\cdot }$ , is specified in (2.11). We denote the particle radius by $r_{0}=d_{p}/2$ , the particle surface by $S_{p}$ and the outward normal vector by $\boldsymbol{n}$ . The particle force divided by the fluid density is defined by

(2.16) $$\begin{eqnarray}F_{j}^{p}=\int _{S_{p}}\boldsymbol{n}\boldsymbol{\cdot }(-p\boldsymbol{e}_{j}+{\it\nu}u_{j,i}\boldsymbol{e}_{i})\,\text{d}S.\end{eqnarray}$$

After substituting $\boldsymbol{n}=\boldsymbol{e}_{r}=\tilde{\boldsymbol{e}}_{1}=\unicode[STIX]{x1D608}_{1m}\boldsymbol{e}_{m}$ , equation (2.8), $\boldsymbol{e}_{m}\boldsymbol{\cdot }\boldsymbol{e}_{i}={\it\delta}_{mi}$ and $\unicode[STIX]{x1D608}_{1i}\unicode[STIX]{x1D608}_{ki}={\it\delta}_{1k}$ into (2.16), where ${\it\delta}$ is the Kronecker delta, we find

(2.17) $$\begin{eqnarray}\displaystyle F_{j}^{p} & = & \displaystyle \frac{d_{p}^{2}}{4}\int _{0}^{2{\rm\pi}}\!\int _{0}^{{\rm\pi}}(-p\unicode[STIX]{x1D608}_{1j}+{\it\nu}\unicode[STIX]{x1D60E}_{1i}\unicode[STIX]{x1D608}_{ij})\sin {\it\theta}\,\text{d}{\it\theta}\,\text{d}{\it\phi}\nonumber\\ \displaystyle & = & \displaystyle \frac{\text{d}_{p}^{2}}{4}\int _{0}^{2{\rm\pi}}\!\int _{0}^{{\rm\pi}}(-p\unicode[STIX]{x1D608}_{1j}+{\it\nu}(u_{r,r}\unicode[STIX]{x1D608}_{1j}+u_{{\it\theta},r}\unicode[STIX]{x1D608}_{2j}+u_{{\it\phi},r}\unicode[STIX]{x1D608}_{3j}))\sin {\it\theta}\,\text{d}{\it\theta}\,\text{d}{\it\phi}.\end{eqnarray}$$

The boundary conditions on the surfaces of the particles follow from impermeability and no slip. If each spherical coordinate system moves and rotates with the same translative and angular velocities as the corresponding particle, the boundary conditions are that $u_{r}$ , $u_{{\it\theta}}$ and $u_{{\it\phi}}$ are zero on particle surfaces. Then (2.11) implies $u_{r,r}=0$ on particle surfaces. For moving particles, the system is closed by an ordinary differential equation for the particle velocity vector (based on a force balance which contains $\boldsymbol{F}^{p}$ ) and an ordinary differential equation for the particle angular velocity vector (based on a torque balance). In the present paper, we will consider stagnant non-rotating particles only. Thus we define $\boldsymbol{u}^{p}=\mathbf{0}$ , and we define $u_{r}=u_{{\it\theta}}=u_{{\it\phi}}=0$ on particle surfaces, while $u_{r,r}$ in (2.17) is replaced by zero.

2.3 Overlapping grids

The entire flow domain is denoted by ${\it\Omega}$ . It excludes the regions inside the particles. The flow domain is a cube with length $L$ and spherical holes of diameter $d_{p}=2r_{0}$ . Periodic boundary conditions are applied in each direction.

Staggered grids are used to discretize the equations on the different meshes. Thus the pressure is defined at cell centres, and the velocity components are defined at cell faces (faces of $p$ -cells). The Cartesian grid is uniform on the cubic domain of size $L_{1}^{3}$ and contains $N_{1}^{3}$ cubic cells of size $h_{1}$ . However, not each cell appears in the equations; the Cartesian grid contains spherical holes centred at particle locations. Indicator functions are used to cut out each cell whose centre lies within a distance $r_{a}$ from a particle. The cells that remain are called the interior Cartesian cells. The faces shared by these cells and the cells cut out of the grid are called the Cartesian boundary faces. The non-interior cells with one or more Cartesian boundary faces are called virtual Cartesian cells. For virtual Cartesian cells, the pressure at the centre and the velocity components at the faces (including the Cartesian boundary faces) are obtained by interpolation from the corresponding spherical grid.

Around each particle, a spherical domain is defined, meshed by a grid of $N_{r}\times N_{{\it\theta}}\times N_{{\it\phi}}$ interior spherical cells. The boundaries in the radial direction are located at the radii $r_{0}=d_{p}/2$ and $r_{b}$ . The faces at radius $r_{b}$ are called spherical boundary faces. For each particle, there is a layer of virtual spherical cells, exterior cells with one face located at radius $r_{b}$ . The pressure at the centre and the velocity components at the faces of each spherical virtual cell (including the spherical boundary face) are obtained by interpolation from the Cartesian grid. To limit the complexity of the discretization method to some extent, spherical grids do not overlap each other. Each spherical grid has its own distinct overlap with the Cartesian grid, the overlap region being the region between $r_{a}$ and $r_{b}$ ( $r_{a}=3r_{b}/5$ in the PR-DNS and the Stokes flow test case, while $r_{a}=(r_{0}+r_{b})/2$ in the other test cases in § 2.5). The grid sizes are ${\rm\Delta}r$ and ${\rm\Delta}{\it\theta}={\rm\Delta}{\it\phi}={\rm\pi}/N_{{\it\theta}}$ ( $N_{{\it\phi}}=2N_{{\it\theta}}$ ). The grid is stretched in the radial direction. We define $r_{j}^{s}=r_{0}\,\text{exp}(j{\rm\pi}/N_{{\it\theta}})$ for the locations of $u_{r}$ ( $j=0,1,\ldots ,N_{r}$ ) and $r_{j}^{c}=r_{0}\,\text{exp}((j-1/2){\rm\Delta}{\it\theta})$ for the locations of $p$ , $u_{{\it\theta}}$ and $u_{{\it\phi}}$ ( $j=0,1,\ldots ,N_{r}+1$ ). The superscripts $s$ and $c$ refer to staggered and cell-centre locations, respectively. The stretching is such that at each radial location the cell size in the wall-normal direction and the maximum cell sizes in the tangential direction are approximately the same: ${\rm\Delta}r\approx r{\rm\Delta}{\it\theta}=r{\rm\Delta}{\it\phi}$ . See figure 1 for an illustration of the overlapping grids. The discretization schemes on these grids and the interpolations from one grid to another are specified in appendix A.

Figure 1. Illustration of the concentric spheres, $r=r_{a}$ , $r=r_{e}$ (defined in § 2.4) and $r=r_{b}$ , drawn on the overlapping grids used to mesh the region around one of the 64 particles in simulation L128 introduced in § 3. Part of the plane $x_{3}=4$ is shown.

2.4 Statistical operators

For a given quantity $q$ we introduce two means: (1) the global mean $\overline{q}$ , obtained by averaging over time and the entire flow domain ${\it\Omega}$ , and (2) the radial mean $\langle q\rangle =\langle q\rangle (r)$ , which is a function of $r$ and obtained by averaging over time and all points with distance $r$ to the nearest particle centre. The corresponding standard deviations or root-mean-square values are defined by $\overline{\text{RMS}}(q)=\overline{q^{2}}-\overline{q}^{2}$ and RMS $(q)=\langle q^{2}\rangle -\langle q\rangle ^{2}$ . Owing to the nature of isotropic turbulence, different variables can be statistically equivalent, for example $\langle u_{1}^{2}\rangle =\langle u_{2}^{2}\rangle =\langle u_{3}^{2}\rangle$ . This feature has been exploited in the computation of the statistics by averaging over components with equivalent statistical behaviour. Thus, by definition, statistically equivalent quantities are the same, also numerically.

Radial mean quantities were computed only for $r<r_{b}\approx 7r_{0}$ , except when indicated otherwise. The evaluation of the global mean causes complications, due to the integral over ${\it\Omega}$ involved. For the approximation of the latter integral, radii $r_{e}^{-}$ , $r_{e}$ and $r_{e}^{+}$ are defined, such that $r_{a}<r_{e}^{-}<r_{e}<r_{e}^{+}<r_{b}$ . More specifically $r_{e}=r_{N_{r}-2}^{s}$ ( $r_{i}^{s}$ is defined in § 2.3), $r_{e}^{-}=r_{e}-h_{1}\sqrt{3}$ and $r_{e}^{+}=r_{e}+h_{1}\sqrt{3}$ ( $h_{1}$ is the Cartesian grid size). The domain ${\it\Omega}$ is split into three disjoint parts: ${\it\Omega}_{1}$ , ${\it\Omega}_{2}$ and ${\it\Omega}_{3}$ . The set ${\it\Omega}_{1}$ contains all points that have a distance larger than $r_{0}$ but smaller than $r_{e}^{-}$ from the centre of a particle. The set ${\it\Omega}_{3}$ is composed of all Cartesian cells whose centre has a distance larger than $r_{e}^{+}$ from each particle. The set ${\it\Omega}_{2}$ is defined by ${\it\Omega}\setminus ({\it\Omega}_{1}\cup {\it\Omega}_{3})$ . The integral of a quantity $q$ over ${\it\Omega}_{1}$ is obtained by numerical integration in spherical coordinates, using the spherical grid cells. The integral of a quantity $q$ over ${\it\Omega}_{3}$ is obtained by numerical integration in Cartesian coordinates, using the Cartesian grid cells. To integrate $q$ over ${\it\Omega}_{2}$ , ${\it\Omega}_{2}$ is meshed by small cubic cells of size $h_{1}/5$ . All centres of the small cubic cells are inside ${\it\Omega}_{2}$ . The small cubic cells do not overlap ${\it\Omega}_{3}$ . The discrete values of $q$ in ${\it\Omega}_{2}$ are based on the solution in the spherical domains. Owing to the approximation of ${\it\Omega}_{2}$ by cubic cells, the integral over ${\it\Omega}_{2}$ is only first-order accurate, but, since the size of these cells is five times smaller than $h_{1}$ , the error in the approximation is sufficiently small.

2.5 Results of test cases

In this subsection, we present the results of four validation studies for steady flows past one or multiple fixed spheres. In each case $d_{p}=1$ and ${\it\nu}=1$ . The forcing function was zero and the particle Reynolds number $Re_{p}$ was based on the particle diameter and the velocity at infinity $(\boldsymbol{U}^{\infty })$ , except in the fourth validation. The Cartesian grid was uniform and cubic, except in the third validation. The polar axis of each spherical grid was parallel to the $x_{3}$ -direction, while the radial stretching function was the one described in § 2.3. The ratio $d_{p}/h$ represents the number of grid points by which the diameter was resolved. In the overset simulations, $h={\rm\pi}d_{p}/N_{{\it\phi}}$ (which is the minimum grid size in the plane ${\it\theta}={\rm\pi}/2$ ), such that $d_{p}/h=N_{{\it\phi}}/{\rm\pi}$ . In each grid refinement, the resolution was doubled in each direction. Each test was run until the steady state was achieved.

Table 1. Simulation results for steady Stokes flow past a sphere in an infinite domain.

First, the solver was validated against the analytical solution for steady Stokes flow past a single spherical particle in an infinite domain ( $|\boldsymbol{U}^{\infty }|=1$ ). Numerically, the particle was placed at the centre of a cubic domain with length 16. Thus one stretched spherical grid was overset on a uniform Cartesian grid. The analytical velocity vector was prescribed on the faces of the cube, and the convective terms were multiplied by zero. Tests were performed for $\boldsymbol{U}^{\infty }$ pointing in the $x_{1}$ -direction, the $x_{2}$ -direction and the $x_{3}$ -direction, respectively. Since the differences between the results of the first two cases were negligible, only the first and third case are shown in table 1. Three resolutions were used. The coarsest resolution was given by $N_{r}=15$ , $N_{{\it\theta}}=24$ , $N_{{\it\phi}}=48$ and $N_{1}=32$ ( $d_{p}/h=15$ ). Table 1 shows the maximum norm of the difference between the numerical and analytical velocity vectors. It also shows the drag force $F^{p}$ on the particle, decomposed into a pressure part $F^{p,pres}$ and a viscous part $F^{p,visc}$ . The deviations from the analytical values are small (less than 2 % for the coarse grid) and converge to zero upon grid refinement.

As a second validation, the drag coefficient $C_{D}$ of a non-Stokesian particle in an infinite domain was computed for three resolutions and compared to the value computed by Bagchi & Balachandar (2002), who used a spectral method. In this case $Re_{p}=100$ , which is significantly higher than the maximum particle Reynolds number occurring in the turbulent flow considered in this paper. The velocity at infinity was pointing in the $x_{3}$ -direction and prescribed as an inflow condition at one face of the cube. On the other five faces, outflow conditions were prescribed (the pressure and the first-order velocity derivatives were set to zero). The sphere was placed in the centre of a cubic domain with length 30. The coarsest resolution was given by $N_{r}=21$ , $N_{{\it\theta}}=24$ , $N_{{\it\phi}}=48$ and $N_{1}=30$ ( $d_{p}/h=15$ ). For this resolution $C_{D}=1.16$ was obtained. After one grid refinement ( $d_{p}/h=31$ ) $C_{D}=1.096$ was found, while after two grid refinements ( $d_{p}/h=61$ ) $C_{D}=1.091$ was obtained. The latter two values of the drag coefficient are within 1 % of $1.09$ , the value reported by Bagchi & Balachandar (2002). In immersed boundary simulations of this flow in a smaller computational domain, $C_{D}=1.178$ (Tang et al. 2014) and $C_{D}=1.175$ (Baltussen 2015) were found, both for $d_{p}/h=20$ .

As a third validation, the drag $C_{D}$ and lift $C_{L}$ coefficients of a particle near a moving flat wall were computed and compared to the values reported by Zeng, Balachandar & Fischer (2005), who used a spectral element method. In this case $Re_{p}=10$ , and the wall, located at $x_{3}=0$ , was moving with constant velocity in the $x_{1}$ -direction. The computational domain size was $24\times 14\times 8$ and the centre of the particle was put at position $x_{1}=8$ , $x_{2}=7$ and $x_{3}=1$ . The flow was simulated for two resolutions; the coarsest one was given by $N_{r}=4$ , $N_{{\it\theta}}=24$ , $N_{{\it\phi}}=48$ , $N_{1}=88$ , $N_{2}=68$ and $N_{3}=45$ $(d_{p}/h=15)$ . In this case the Cartesian grid was also stretched. For the coarse resolution, $C_{D}=4.64$ and $C_{L}=0.342$ ; while after one grid refinement ( $d_{p}/h=31$ ) $C_{D}=4.70$ and $C_{L}=0.348$ were obtained. The latter two values are within 1 % of the values reported by Zeng et al. (2005), $C_{D}=4.72$ and $C_{L}=0.351$ .

Figure 2. Comparison between overset method (circles) and an immersed boundary method (triangles). Convergence characteristics of the pressure part (open symbols) and viscous part (filled symbols) of the drag force for a face-centred cubic arrangement of particles in a periodic domain ( $Re_{p}=40$ and ${\it\alpha}=0.2$ ). The results of the immersed boundary method were taken from figure 4(a) in Tenneti et al. (2011), where they were also denoted by triangles.

As a fourth validation, the pressure and viscous parts of the drag force on a sphere in the flow past a structured array of spheres were computed and compared to results obtained by Tenneti et al. (2011). These authors used an immersed boundary method (second-order accurate direct forcing embedded into a pseudo-spectral computation; the same method was recently used by Mehrabadi et al. (2015)). It is remarked that many variants of immersed boundary methods exist and have been validated for flows past multiple particles – see also, for example, Uhlmann (2005), Mark & van Wachem (2008) and Breugem (2012). In the test presented by Tenneti et al. (2011), a face-centred cubic (FCC) array of particles was used with a particle volume fraction ${\it\alpha}=0.2$ . To simulate this case with the present method, four particles were placed in the FCC pattern in a cubic periodic domain of length $(10{\rm\pi}/3)^{1/3}$ . The flow was driven by a time-dependent uniform forcing in the $x_{1}$ -direction, chosen such that the flow converged to a steady state in which $U_{1}$ , the volume average of the $u_{1}$ velocity, was $40$ , such that $Re_{p}$ based on $U_{1}$ was also 40. Simulations were performed for three resolutions; the coarsest one was given by $N_{r}=3$ , $N_{{\it\theta}}=24$ , $N_{{\it\phi}}=48$ and $N_{1}=24$ ( $d_{p}/h=15$ ). The pressure and viscous parts of the drag force, normalized by $3{\rm\pi}{\it\nu}d_{p}U_{1}(1-{\it\alpha})$ , are shown in figure 2. The pressure part of the normalized force on the three consecutive grids was 3.363, 3.541 and 3.573, respectively. The viscous part of the normalized force on the three consecutive grids was 5.715, 5.817 and 5.846, respectively. Apparently the errors in the two components of the drag force are approximately 1 % if the present numerical method is used for $d_{p}/h=31$ and $Re_{p}=40$ . According to figure 2, the slopes of the curves for the overset method are nearly horizontal for $d_{p}/h\geqslant 30$ , while this is clearly not the case for the immersed boundary method. According to this and the second validation, the overset method is apparently at least as accurate as some immersed boundary methods.

3 Definition of the direct numerical simulations

We consider two flow cases: an unladen statistically stationary homogeneous isotropic turbulent flow, and the same flow laden with fine solid fixed particles. In the following three subsections, the stochastic forcing of the stationary turbulence is defined, the simulation cases are specified, and the accuracy of the simulations is discussed.

3.1 Stochastic forcing

The forcing in the experiments by Hwang & Eaton (2006) and Tanaka & Eaton (2010) was caused by jet actuators driven by sine waves at random frequencies. In view of these experiments, it is a logical step to simulate (particle-laden) homogeneous isotropic turbulence forced by a number of Fourier modes at random frequencies. The stochastic forcing is applied to the large spatial (and temporal) scales only, such that the energy is not injected into scales that correspond to strong velocity gradients (scales where most of the dissipation due to turbulence and particles occurs). Thus the small spatial and temporal scales in the flow are generated by the turbulence cascade process and the particles, and not by the stochastic forcing, at least not directly. Nonetheless, the statistical results of forced homogeneous isotropic turbulence for given $Re_{{\it\lambda}}$ are expected to depend on the specific forcing, at least to some extent. However, the statistical state of decaying homogeneous isotropic turbulence is also not unique; the power-law exponent depends on the initial condition in a simulation or the upstream condition in an experiment (Djenidi & Antonia 2015).

The procedure used for the stochastic forcing is similar to the one used by Yeung & Pope (1989) in their simulations of homogeneous isotropic turbulence. More specifically, the forcing function $\boldsymbol{f}$ is defined by

(3.1) $$\begin{eqnarray}\boldsymbol{f}(\boldsymbol{x},t)=\int _{0}^{t}\left(\frac{1}{t_{fil}}\text{e}^{-(t-t^{\prime })/t_{fil}}\mathop{\sum }_{j=1}^{18}\left[\boldsymbol{g}_{j}(t^{\prime })-\frac{(\boldsymbol{g}_{j}(t^{\prime })\boldsymbol{\cdot }\boldsymbol{k}_{j})\boldsymbol{k}_{j}}{|\boldsymbol{k}_{j}|^{2}}\right]\text{e}^{\sqrt{-1}\,\boldsymbol{k}_{j}\boldsymbol{\cdot }\boldsymbol{x}}\right)\text{d}t^{\prime }\end{eqnarray}$$

(the summation convention is not used). For each $j$ , the symbol $\boldsymbol{k}_{j}$ refers to one of the 18 three-dimensional wavevectors that satisfy $0<|\boldsymbol{k}|\leqslant \sqrt{2}(2{\rm\pi}/L_{1})$ and fit into the periodic domain. Forcing based on a larger number of modes requires a larger computational box to achieve the same $Re_{{\it\lambda}}$ for given ${\it\eta}$ . It was not feasible to simulate the particle-resolved case in a larger domain without making compromises on the numerical and statistical accuracy of the results. In some of the simulations reported by Yeung & Pope (1989), forcing was applied to the same $18$ modes.

The expression within square brackets in (3.1) represents the projection on the space of divergence-free periodic functions, such that $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{f}=0$ . For each $j$ , each one of the three components of $\boldsymbol{g}_{j}(t)$ , say $X(t)$ , is an independent stochastic process with the same properties. More specifically, $X(t)$ is an Ornstein–Uhlenbeck process, with parameters chosen such that the mean of $X(t)$ is zero, the standard deviation of $X(t)$ is ${\it\sigma}$ , and the correlation between $X(t)$ and $X(t+t^{\prime })$ is $\text{exp}(-t^{\prime }/t_{OU})$ . The Ornstein–Uhlenbeck process was first applied to the forcing of homogeneous isotropic turbulence by Eswaran & Pope (1988), and since then by many others. It is the only stochastic process that is stationary, continuous, Gaussian and Markovian. However, the Ornstein–Uhlenbeck process is not differentiable and may therefore introduce unphysical effects at small time scales. Yeung & Pope (1989) mention that for this reason they replaced the Ornstein–Uhlenbeck by an integrated Ornstein–Uhlenbeck process to obtain a continuously differentiable $X(t)$ . The exponential filter applied in (3.1) has the same effect: it also makes $\boldsymbol{f}$ continuously differentiable in time. It can be proven that, despite the filtering, each component of $\boldsymbol{f}(\boldsymbol{x},t)$ is a stationary Gaussian process for each $\boldsymbol{x}$ . The Markov property means that, at time $t$ , future states (at times larger than $t$ ) depend only on the present state (at time $t$ ) and not on the past states (at times smaller than $t$ ). The consequence of filtering is that $\boldsymbol{f}$ is formally not Markovian; however, it can be considered as approximately Markovian for time scales larger than $t_{fil}$ . The simulations were performed for ${\it\sigma}=0.60$ , $t_{OU}=2$ and $t_{fil}=0.50$ .

3.2 Direct numerical simulation cases

Table 2. Unladen and laden DNS. Particle volume fraction ${\it\alpha}$ , grid parameters ( $N_{1}$ , $N_{r}$ , $N_{{\it\theta}}$ and $N_{{\it\phi}}$ ) and time step  ${\rm\Delta}t$ .

The set-up of four DNS cases is described: two simulations of the unladen flow (U128 and U64), and two particle-resolved simulations of the laden flow (L128 and L64). Here U and L denote unladen and laden. The number 128 or 64 refers to $N_{1}$ . The main cases are the fine grid simulations (U128 and L128). The two coarser cases (U64 and L64) have been included to show the effect of the numerics when the grid is coarsened.

The fluid viscosity and the length of the periodic cubic flow domain are ${\it\nu}=1$ and $L_{1}=32$ in each case, while the particle diameter and volume fraction in the laden simulations are $d_{p}=2r_{0}=1$ and ${\it\alpha}=0.001$ . The outer boundaries of the spherical domains are located at $r_{b}\approx 7.124r_{0}$ . The resolution and time step of each case are listed in table 2. Since the turbulence enters the system through the stochastic forcing term $\boldsymbol{f}$ , the initial condition of simulation U128 is equal to zero. All other simulations were started at $t=18$ from the velocity field of U128 at the same time. The time interval for statistical averaging is defined by $[t_{1},t_{2}]$ , where $t_{1}=20$ and $t_{2}=300$ , unless mentioned otherwise. The Taylor Reynolds number of the unladen flow is approximately 32 and the Kolmogorov length scale ${\it\eta}$ is approximately $0.46$ . The integral length scale based on the three-dimensional energy spectrum (see e.g. Jimenez et al. (1993) for a definition) is approximately 9. The eddy turnover time, defined by the ratio of the integral length scale and the square root of two-thirds of the turbulence kinetic energy, is approximately 1.4, such that the statistical averaging time, $t_{2}-t_{1}$ , is approximately 200 eddy turnover times.

The laden flow contains $N_{p}=64$ fixed particles, arranged in a simple cubic array. Denoting one of the corners of the cubic domain by the coordinates $(0,0,0)$ , the centres of the particles are located at $(4+8i,4+8j,4+8k)$ , where $i$ , $j$ and $k$ are elements of $\{0,1,2,3\}$ . The particles were injected into the velocity field of U128 at $t=18$ . Since the number of time steps was large in each case and roughly 100 iterations per time step were needed to solve the pressure Poisson equation within the prescribed tolerance, the simulations were computationally expensive. For example, simulation L128 (3 million time steps, 11 million grid points) took around 200 days (wall-clock time) on a single modern computer node with 20 threads (the code was parallelized with OpenMP). Thus approximately $10^{5}$  CPU hours in total were required for this run. These numbers do not imply that the overset method is less efficient than immersed boundary methods, which are typically used on uniform grids to allow the pressure Poisson equation to be solved by a fast direct method. The number of grid points of a uniform grid with the same resolution near particles $(d_{p}/h=31)$ would be approximately $10^{9}$ , approximately 90 times more than in case L128. Thus the local grid refinement facilitated by the overset method can clearly be an advantage.

3.3 Accuracy of the direct numerical simulations

In case U128, the grid size $h_{1}=L_{1}/N_{1}=0.25$ and the Kolmogorov length scale ${\it\eta}=0.46$ are such that $k_{max}={\rm\pi}/h_{1}=5.8/{\it\eta}$ , where $k_{max}={\rm\pi}/h_{1}$ is the maximum wavenumber (the corresponding wave is $\text{exp}(\text{i}k_{max}x_{1})$ ). For spectral methods, $k_{max}=2/{\it\eta}$ is commonly regarded to provide sufficient resolution to capture the velocity gradient (and thus the turbulence dissipation rate) in isotropic turbulence (Yeung & Pope 1989; Jimenez et al. 1993; Ishihara et al. 2007). Although spectral simulations using $k_{max}=2/{\it\eta}<{\rm\pi}/{\it\eta}$ do not resolve ${\it\eta}$ , they do resolve the small scales that significantly contribute to the dissipation of the turbulence (Yeung & Pope 1989; Jimenez et al. 1993; Ishihara et al. 2007). This indicates that the Kolmogorov scale ${\it\eta}$ tends to be a low estimate of the small-scale dynamics. It is remarked that the accuracy of any simulation depends on which variable is looked at. If the fall of a spectrum is exponential, then it is theoretically possible to find a high-order quantity (derivative) which is not accurately represented in any simulation that includes a finite number of wavenumbers (Vreman & Kuerten 2014b ). Of course, a second-order finite difference method will require a higher $k_{max}$ than a spectral method, but this is not necessarily an order of magnitude higher. For a staggered method second-order accurate in wall-normal and fourth-order accurate in the other directions, the small-scale turbulence of a turbulent channel flow was simulated with similar accuracy as in a fully spectral simulation, while $k_{max}$ was only $1.33$ times higher than in the spectral case (Vreman & Kuerten 2014a ). Thus it is reasonable to assume that if $k_{max}$ is 2.9 times higher (which is the case in U128 and L128), the standard second-order staggered method is sufficiently accurate to resolve the turbulence. The three-dimensional energy spectra shown in figure 3(a) confirm this. The decay of the spectra is large and the two spectra nearly coincide. This demonstrates that virtually all eddies are captured and that the flow is well resolved, especially in simulation U128, which shows a spectral fall-off of approximately 15 orders of magnitude.

Figure 3. Demonstration of the effect of resolution on the energy spectra. (a) The three-dimensional energy spectrum as a function of $k$ (the length of the three-dimensional wavevector) is shown for U128 (solid line) and U64 (symbols). (b) The energy spectrum as a function of the azimuthal wavenumber ( $k_{{\it\phi}}$ ) at $r=r_{b}$ and ${\it\theta}={\rm\pi}/2$ is shown for L128 (solid line) and L64 (symbols).

It is not straightforward to compute energy spectra for the laden flow, except in the periodic ${\it\phi}$ -direction. For both L128 and L64, ${\it\phi}$ spectra are shown in figure 3 b). These spectra represent the energy in the spherical components of the velocity field at $r=r_{b}$ and ${\it\theta}={\rm\pi}/2$ . This is a critical region in the flow since there the grid size of the spherical grid is maximum and the spherical velocity is directly computed by interpolation from the Cartesian grid. The interpolation errors show up as a horizontal tail in the spectrum. Fortunately, the level of the tail is low, in particular in case L128, in which case the horizontal part of the tail is nine orders of magnitude lower than the peak of the spectrum.

The accuracy of L128 was further confirmed by verification of the energy balance and a grid refinement study. With respect to the energy balance, the global balance error in L128 is less than 0.5 %, while the radial balance error is less than 1 % at almost all radial locations. These results will be shown in detail in the next section. With respect to the grid refinement study, L128 is a refinement of case L64, by a factor of two for the resolution in each spatial direction and by a factor of four for the temporal resolution. If a particular deviation between a result of simulation L128 and L64 is dominated by discretization errors, the effect of discretization errors on the quantity corresponding to the deviation is expected to be three times smaller than the deviation (because the method is second-order accurate). However, if the deviation is dominated by statistical errors, it provides just an estimate of the error in simulation L128 (and L64). Most deviations between L128 and L64 were found to be quite small (as demonstrated by the tables and figures in the next section). No deviation affected the conclusions.

As indicated above, the discrepancy between simulations at different resolutions is caused not only by discretization, but also by statistical errors. It is not easy to determine which type of error is dominant. All results in the next section correspond to $t_{2}=300$ ( $t_{2}-t_{1}=280$ ), but the dependence on the end time of the averaging $t_{2}$ has been investigated for fixed $t_{1}=20$ . Thus, the statistical error of the $L_{2}$ -norm of the stochastic forcing term was estimated to be approximately 2 % for $t_{2}=300$ . Since we are interested in the effect of the particles on the turbulence, the ratios of laden and corresponding unladen values are particularly relevant. A comparison of the modification ratios L128/U128 for $t_{2}=300$ and $t_{2}=150$ showed that the differences between the ratios for the two averaging end times were a maximum 0.01 in table 3 and a maximum 2 % in table 7. The differences in the quantities themselves were somewhat larger, but less than 2 % for all quantities, except for the skewness and flatness factors of the pressure, which showed a difference of approximately 5 %. It is mentioned that, in order to reduce the effect of the uncertainty due to the statistical errors in the forcing term on the comparison between laden and unladen simulations, all simulations were run for the same realization of the stochastic forcing $\boldsymbol{f}(\boldsymbol{x},t)$ . Thus, for a given $\boldsymbol{f}$ and $t_{2}$ , any statistical error between two simulations can only be due to the nonlinear chaotic behaviour of the turbulence (and not to the stochastic forcing). The effect of $t_{2}$ upon the turbulence kinetic energy is shown in figure 4. It is observed that for $t_{2}>150$ the dependence on $t_{2}$ is indeed sufficiently small to conclude that the ratio between the turbulence kinetic energies of the laden and unladen flow reduces to approximately 0.91 with an error estimate of 0.01.

Figure 4. The effect of the end time ( $t_{2}$ ) of the time interval for statistical averaging on the global turbulence kinetic energy, for simulations U128 (squares), L128 (circles), U64 (red dashed) and L64 (black dash-dotted).

Table 3. Global turbulence kinetic energy and dissipation rate and other basic quantities from simulations U128, L128 and L64. The ratios in the last five columns express the modifications of the global quantities by the particles in ratios of laden to unladen quantities for simulations L128, L64 and the point-particle simulations discussed in § 4.5. Ratios obtained for half the averaging time differed by maximum 0.01 from those obtained for the full averaging time (see § 3.3).

4 Results

The following five subsections are ordered according to the five research questions posed in § 1. The first four subsections are entirely devoted to the results of the DNS introduced in the previous section. From § 3.3, we conclude that U128 and L128 are sufficiently accurate for the present purposes. The unladen reference simulation is simulation U128. The coarse unladen simulation (U64) will not be used in the present section. Case L128 should be regarded as the primary laden simulation. The resolution used in that simulation is the resolution recommended for the flow with particles. However, since the reader may be interested to see the comparison between the fine and coarse laden cases in detail, the results of both L128 and L64 have been included in the figures and tables in §§ 4.1–4.4. In the fifth subsection, three point-particle simulations (PP0, PP2 and PP4) will be introduced and the results will be compared to the particle-resolved simulation L128. For the sake of conciseness, results of the point-particle simulations have been included in the tables in the first four subsections, but the discussion of these results is postponed to the fifth subsection.

Figure 5. A snapshot of the turbulence around one of the particles at $x_{2}=4$ and $t=150$ , taken from simulation L128. (a) Contours of the magnitude of the velocity vector, $(u_{r}^{2}+u_{{\it\theta}}^{2}+u_{{\it\phi}}^{2})^{1/2}$ ; the contour increment is 1.5. (b) Contours of the square root of the local dissipation rate, $({\it\nu}\unicode[STIX]{x1D60E}_{ij}\unicode[STIX]{x1D60E}_{ij})^{1/2}$ ; the contour increment is 8. (c) The projection of the velocity vector on the plane, $(u_{1},u_{2})$ , and contours of the radial velocity $u_{r}$ ; the contour increment is 2. (d) Contours of the pressure, $p$ ; the contour increment is 5. The vertical lines observed at $x_{1}=12$ denote the locations of the apparent singularities of the Navier–Stokes equations in spherical coordinates at ${\it\theta}=0$ and ${\it\theta}={\rm\pi}$ .

4.1 Contours and radial profiles of basic quantities

Before we consider the statistical radial profiles of basic variables, we show a snapshot from simulation L128 as an illustration of the modification of an intense structure of kinetic energy by a particle. Figure 5 zooms into the region around a particle in the $x_{1}$ – $x_{3}$ plane at $t=150$ . The velocity is roughly in the plane, i.e. the $u_{2}$ velocity is relatively small in this particular region at this particular time. The region of low kinetic energy around the stagnant particle is elongated in the main direction of the velocity field (figure 5 a), which points upwards and somewhat to the right (figure 5 c). Thus the velocity field in figure 5 is dominated by a positive $u_{3}$ velocity. For the discussion, it is convenient to introduce a front side and a rear side in terms of the direction of the flow surrounding the particle. The front side is the region on the particle surface where the outward-directed surface normal vector $\boldsymbol{n}$ is approximately opposite to the direction of the surrounding flow, while the rear side is the region on the particle surface where $\boldsymbol{n}$ is aligned with the direction of the surrounding flow. Since the surrounding flow velocity is an instantaneous vector, the front and rear sides are time-dependent regions. Figure 5(d) shows that the maximum pressure occurs at the front and the minimum pressure at the rear side, as expected. As has been reported by others, the local dissipation rate is enhanced in the vicinity of the particles (Burton & Eaton 2005; Tanaka & Eaton 2010). It is remarked that the square root of the dissipation rate is shown in figure 5, and thus the local enhancement can be very strong. The dissipation rate on the surface is due to shear caused by the tangential velocity components. The surface locations where this shear is large are relatively far away from the rear and front sides, further away from the rear than from the front side. We will refer back to figure 5 a number of times, since several features will be recognized in the statistics of the simulation.

We denote the global turbulence kinetic energy and dissipation rate by $\overline{K}$ and $\overline{{\it\epsilon}}$ , and the radial profiles by $K$ and ${\it\epsilon}$ . We recall that the global statistics are based on the average over time and the entire fluid volume, while the radial statistics are based on $\langle \cdot \rangle$ , the average over time and all locations with distance $r$ to the nearest particle centre (see § 2.4). Since the means of all velocity components are zero, we define

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{K}=\overline{u_{j}u_{j}}/2, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{{\it\epsilon}}={\it\nu}\overline{u_{j,i}u_{j,i}}, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle K=K(r)=\langle u_{j}u_{j}\rangle /2=\langle u_{r}^{2}+u_{{\it\theta}}^{2}+u_{{\it\phi}}^{2}\rangle /2, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\epsilon}={\it\epsilon}(r)={\it\nu}\langle u_{j,i}u_{j,i}\rangle ={\it\nu}\langle \unicode[STIX]{x1D60E}_{ij}\unicode[STIX]{x1D60E}_{ij}\rangle . & \displaystyle\end{eqnarray}$$

Results for the radial profiles $K(r)$ and ${\it\epsilon}(r)$ and other basic statistical quantities are shown in figure 6. All profiles have been normalized by appropriate unladen references values from U128, which are denoted by $K_{0}$ (unladen $\overline{K}$ ), ${\it\epsilon}_{0}$ (unladen $\overline{{\it\epsilon}}$ ) and $p_{0}$ (unladen $\overline{\text{RMS}}(p)$ ). The numbers can be found in table 3, discussed in the next subsection.

Figure 6. Radial profiles of basic quantities for cases L128 (blue, thick solid) and L64 (red, thin dashed), normalized by reference values of the unladen simulation U128. The thin dotted horizontal lines represent the normalized unladen values. (a) Turbulence kinetic energy $K(r)$ . (b) Variances of the radial velocity component (circles) and polar and azimuthal velocity components (squares), $\langle u_{r}^{2}\rangle$ , $\langle u_{{\it\theta}}^{2}\rangle$ and $\langle u_{{\it\psi}}^{2}\rangle$ , scaled by $2K_{0}/3$ . (c) RMS pressure (circles) and mean pressure $\langle p\rangle$ , scaled by $p_{0}$ (squares; $\overline{p}$ has been subtracted). (d) Turbulence dissipation rate ${\it\epsilon}(r)$ . The black dash-dotted line represents $2{\it\nu}\langle s_{ij}s_{ij}\rangle$ for case L128.

Figure 6(a) shows the attenuation of the turbulence kinetic energy as a function of $r$ , the distance to the centre of the nearest particle. For each $r$ , $K(r)$ is lower than the unladen value (indicated by the dotted horizontal line). Not surprisingly, the strongest attenuation occurs on particle surfaces, $K(r_{0})=0$ , and the weakest attenuation relatively far away from the particles. The largest distance to the nearest particle centre, $r_{max}=8\sqrt{3}r_{0}$ , is not shown on the figure since the radial profiles were computed on the spherical domains, i.e. for $r_{0}<r<r_{b}\approx 7r_{0}$ . However, $K(r_{max})$ was separately computed: $K(r_{max})\approx 0.93K_{0}$ . It is concluded that the turbulence kinetic energy is significantly attenuated in the entire flow domain.

The profiles of the individual contributions of the three Cartesian velocity components to $2K$ , $\langle u_{1}^{2}\rangle$ , $\langle u_{2}^{2}\rangle$ and $\langle u_{3}^{2}\rangle$ , are the same and equal to $2K/3$ . However, the individual contributions of the three spherical velocity components are not the same (figure 6 b). It appears that the particles suppress the variance of the radial component $\langle u_{r}^{2}\rangle$ more than the variances of the other two components, $\langle u_{{\it\theta}}^{2}\rangle =\langle u_{{\it\phi}}^{2}\rangle$ . Figure 6(c) shows that the RMS of the pressure is reduced, except near particles ( $r<1.8r_{0}$ ; at $r_{0}$ it is enhanced by approximately 40 %). Unlike the mean profiles of the velocity components, the mean profile of the pressure is not zero. Figure 6(c) shows that the mean pressure is relatively low on the particle surfaces.

Another basic quantity is the turbulence dissipation rate profile ${\it\epsilon}(r)$ , which is shown in figure 6(d). Near the surfaces of the particles, ${\it\epsilon}(r)$ has increased by more than a factor 100 (note the logarithmic scaling of the vertical axis, ${\it\epsilon}(r_{0})=2610=118{\it\epsilon}_{0}$ in case L128 and $2650$ in L64). However, further away from the particles ( $r>3.6r_{0}$ ), the turbulence dissipation rate is attenuated, like the turbulence kinetic energy and the pressure fluctuation. Figure 6(d) shows that the profile ${\it\epsilon}(r)={\it\nu}\langle u_{j,i}u_{j,i}\rangle$ is marginally different from the profile $2{\it\nu}\langle s_{ij}s_{ij}\rangle$ , where $s_{ij}=(u_{j,i}+u_{i,j})/2$ is the rate of strain. The second expression is precisely the viscous source term in the internal energy equation and is therefore regarded as the thermodynamically correct definition of the energy dissipation rate. However, since the difference between the two quantities is only a transport term (the volume integral of the difference vanishes), it is not uncommon to base the definition of the profile of the turbulence dissipation rate on ${\it\nu}u_{j,i}u_{j,i}$ (see, for example, Mansour, Kim & Moin (1988) and Vreman & Kuerten (2014b ) and references therein). It simplifies the computation and leads to a more compact form of the turbulence kinetic energy transport equation than the definition based on $2{\it\nu}s_{ij}s_{ij}$ . It is remarked that the two quantities are formally equal at $r=r_{0}$ , since $2s_{ij}s_{ij}=(\unicode[STIX]{x1D60E}_{ij}+\unicode[STIX]{x1D60E}_{ji})(\unicode[STIX]{x1D60E}_{ij}+\unicode[STIX]{x1D60E}_{ji})/2$ and $u_{j,i}u_{j,i}=\unicode[STIX]{x1D60E}_{ij}\unicode[STIX]{x1D60E}_{ij}$ , while $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ are the only non-zero components of $\unicode[STIX]{x1D642}$ at $r=r_{0}$ .

That the local dissipation rate is enhanced near the particles is known, but that the factor of the increase on the surfaces is this large (of the order of 100) is a surprise. The high peak of the dissipation rate profile could not be inferred from Burton & Eaton (2005), since only the part of the curve where the enhancement was less than a factor of two was shown. Tanaka & Eaton (2010) reported that the dissipation rate close to the particle surface is enhanced by a factor of three. In simulations performed with another numerical method (the immersed boundary method) and for much larger moving particles, the dissipation rate on particle surfaces was found to be approximately three times larger than the dissipation rate far away from surfaces of the particles (Lucci et al. 2010). A theoretical estimate of the turbulence dissipation rate on particles, presented in § 4.5, shows that it strongly increases with decreasing particle diameter. Thus the main reason for the large discrepancy with the results of Lucci et al. (2010) is probably that in that work $d_{p}/{\it\eta}$ was roughly eight times larger.

Figure 7. Fraction of the global turbulence dissipation rate that occurs within distance $r$ to a particle centre: (a)  $c_{{\it\epsilon}}$ versus $r/r_{0}$ and (b)  $c_{{\it\epsilon}}$ versus $c_{v}$ . Cases L128 (blue, thick solid) and L64 (red, thin dashed).

The very large dissipation rate near particles supports the view expressed by Hwang & Eaton (2006) and Tanaka & Eaton (2010) for heavy particles, namely that such particles, like grids, act as localized dampers of the turbulence motion. To quantify how much of the total amount of dissipation in the flow occurs near the particles, we define $V_{r}$ as the set of all points in the flow domain within distance $r$ of a particle centre, and we define

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle c_{v}(r)=(4/3){\rm\pi}N_{p}(r^{3}-r_{0}^{3})/(1-{\it\alpha})L_{1}^{3}, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle c_{{\it\epsilon}}(r)=4{\rm\pi}N_{p}\int _{r_{0}}^{r}{\it\epsilon}(r^{\prime }){r^{\prime }}^{2}\,\text{d}r^{\prime }\bigg/\int _{{\it\Omega}}{\it\nu}u_{j,i}u_{j,i}\,\text{d}{\it\Omega}. & \displaystyle\end{eqnarray}$$

The function $c_{v}$ is the fraction between the volume of $V_{r}$ and the total flow volume, while $c_{{\it\epsilon}}$ is the fraction between the dissipation that occurs in $V_{r}$ and the total dissipation. The function $c_{{\it\epsilon}}$ is shown in figure 7, plotted against $r/r_{0}$ and $V_{r}$ . Figure 7 shows that 10 % of the total dissipation occurs in 0.5 % of the total volume, and 5 % occurs in only 0.1 % of the volume.

4.2 Global turbulence attenuation

In this subsection we investigate the modification of global (domain- and time-averaged) statistics of the flow, most of them related to $\overline{K}$ and $\overline{{\it\epsilon}}$ . First, the simulation results shown in tables 3–5 are discussed. Afterwards, the main simulation results are discussed in the light of experimental work found in the literature (table 6).

Table 4. Global energy balance for the simulations U128, L128, L64 and the point-particle simulations discussed in § 4.5.

Table 5. Particle force and particle-induced dissipation in simulations L128, L64 and the point-particle simulations discussed in § 4.5.

Table 6. Comparison between experiments of particle-laden stationary homogeneous isotropic turbulence in the literature (Hwang & Eaton 2006; Tanaka & Eaton 2010) and the present particle-resolved DNS L128.

The simulation results for the global turbulence kinetic energy and dissipation rate and various other basic quantities are summarized in table 3. The last five columns show the ratios of laden to unladen quantities and thus quantify the suppression or enhancement of the quantities due to particles. The bold ratios quantify the global turbulence modification in the most accurate particle-resolved case (L128). Table 3 shows that the global turbulence kinetic energy $\overline{K}$ in simulation L128 reduces to 91 % of the unladen value, which corresponds to an attenuation of 9 %. The modification of the global turbulence dissipation rate $\overline{{\it\epsilon}}$ due to the particles is much smaller, an increase of 3 %. The next quantity in table 3 is the $L_{2}$ -norm of the stochastic forcing term (the $L_{2}$ -norm of any vector field $\boldsymbol{q}$ is defined by $\overline{|\boldsymbol{q}|^{2}}^{1/2}$ in this paper). The fourth quantity shown in table 3 is the $L_{2}$ -norm of the velocity, which is equal to $(2\overline{K})^{1/2}$ . It is also the $L_{2}$ -norm of the particle Reynolds number $Re_{p}$ since the velocity of the particles is zero, $d_{p}=2r_{0}=1$ and ${\it\nu}=1$ . Thus the $L_{2}$ -norm of $Re_{p}$ is approximately $11$ . The results of $\overline{\text{RMS}}(p)$ show that not only $\overline{K}$ but also the pressure fluctuations are significantly attenuated, even somewhat more than $\overline{K}$ . Owing to the effect of particles upon global turbulence kinetic energy and to a smaller extent upon the global turbulence dissipation rate, quantities constructed from $\overline{K}$ and $\overline{{\it\epsilon}}$ are also modified. The Taylor Reynolds number, $Re_{{\it\lambda}}=\overline{K}(20/3{\it\nu}\overline{{\it\epsilon}})^{1/2}$ , is attenuated by approximately 10 %. The Kolmogorov length scale, ${\it\eta}=({\it\nu}^{3}/\overline{{\it\epsilon}})^{1/4}$ , is not significantly modified. In contrast, the length scale of the large-scale turbulence, $L_{{\it\epsilon}}=\overline{K}^{3/2}/\overline{{\it\epsilon}}$ , is relatively strongly reduced (16 %). The relatively strong decrease of $L_{{\it\epsilon}}$ indicates that the particles diminish the size or the strength of large eddies in particular. The corresponding time scale of large eddies, $t_{{\it\epsilon}}=\overline{K}/\overline{{\it\epsilon}}$ , is also significantly reduced (approximately 12 %).

In summary, table 3 confirms that fine particles significantly attenuate the turbulence kinetic energy in stationary homogeneous isotropic turbulence. The global attenuation of the turbulence kinetic energy is found to be approximately 9 %. That a small volume fraction of fine particles is able to attenuate turbulence is well known, but it is the first time that global turbulence attenuation has been confirmed by particle-resolved DNS for a case in which the particle diameter is only $2.2{\it\eta}$ and the particle volume fraction is only $0.1\,\%$ .

In table 4, the global energy balance is shown for all simulations. The global production of energy by the forcing term is defined by $\overline{P}=\overline{u_{j}\,f_{j}}$ . Theoretically, we should find $\overline{P}-\overline{{\it\epsilon}}=0$ . The numerical evaluation of the energy balance is shown in table 4. It appears that in each case the global balance error, $\overline{E}=\overline{P}-\overline{{\it\epsilon}}$ , is small. The last row shows the absolute value of the relative error, $|\overline{E}|/\overline{{\it\epsilon}}$ . The balance error is very small in case U128 (only 0.1 %). The balance error in the particle-resolved laden case L128 is somewhat larger, but still small ( $0.4\,\%$ ). The balance error in case L64 is 1.6 %, which is not large, but four times larger than in case L128. Since the numerical method is second-order, it indicates that, at least in case L64, the balance error is mainly caused by discretization errors. This is supported by the finding that the difference between the domain average of $u_{i}u_{i}/2$ at $t_{2}$ and $t_{1}$ divided by $t_{2}-t_{1}$ is only 0.15 % of $\overline{{\it\epsilon}}$ in case L64 and 0.07 % in case L128, considerably smaller than the balance error.

In the particle-resolved simulations, the particle force is a numerical evaluation of the exact expression in equation (2.17). To compute this $L_{2}$ -norm, the domain integration in the definition of the global mean has been replaced by the average over all particles in the domain. The $L_{2}$ -norms of the particle forces in the particle-resolved cases differ by only 0.5 % and are 188.9 and 189.6 in cases L128 and L64, respectively (table 5). The magnitudes of the distinct pressure and viscous contributions to the force are also shown.

Table 7. Global skewness and flatness factors and ratio $\overline{{\it\chi}}$ from simulations U128, L128 and L64. The last five columns express the modifications of the global quantities by the particles in ratios of laden to unladen quantities for simulations L128, L64 and the point-particle simulations discussed in § 4.5. Ratios obtained for half the averaging time differed by a maximum 2 % from those obtained for the full averaging time (see § 3.3).

The forces of the particles on the fluid lead to a strong increase of the turbulence dissipation rate ${\it\epsilon}(r)$ near the particles, as we saw in figure 6(d). It is interesting to quantify the extra dissipation due to particles by a global quantity, which we call the particle-induced dissipation $\overline{{\it\epsilon}}^{p}$ . The quantity multiplied by the total fluid volume should estimate the energy dissipated by all the particles in the flow domain, like $\overline{{\it\epsilon}}$ multiplied by the total fluid volume represents the total amount of dissipated energy in the flow domain. Intuitively, the energy dissipated by the particles is related to a volume integral of the increase of ${\it\epsilon}(r)$ in figure 6(d). Therefore, we define a radius $R$ sufficiently far away from the particle, where ${\it\epsilon}(r)$ has become approximately flat. We integrate ${\it\epsilon}(r)-{\it\epsilon}(R)$ (the local increase of the radial profile of the turbulence dissipation rate) over the volume between $r_{0}$ and $R$ . The result multiplied by the number of particles and divided by the total fluid volume is used as the definition for the particle-induced dissipation in the particle-resolved simulation:

(4.7) $$\begin{eqnarray}\overline{{\it\epsilon}}^{p}=\frac{N_{p}}{(1-{\it\alpha})L_{1}^{3}}\int _{r_{0}}^{R}4{\rm\pi}r^{2}({\it\epsilon}(r)-{\it\epsilon}(R))\,\text{d}r.\end{eqnarray}$$

In the present paper, we choose the radius $R=r_{b}\approx 7r_{0}$ . If the profile ${\it\epsilon}(r)$ is perfectly horizontal for $r\geqslant R$ then ${\it\epsilon}(r)-{\it\epsilon}(R)=0$ for $r\geqslant R$ and $\overline{{\it\epsilon}}^{p}=\overline{{\it\epsilon}}-{\it\epsilon}(R)$ . However, for the sake of clarity, we use (4.7) to compute $\overline{{\it\epsilon}}^{p}$ for the two particle-resolved simulations. The results are listed in table 5. According to the ratio $\overline{{\it\epsilon}}^{p}/\overline{{\it\epsilon}}$ , approximately 14 % of the turbulence dissipation rate can be attributed to the particles. Also the values of ${\it\epsilon}(r_{b})$ , the turbulence dissipation rate far away from the particles, are included. These are significantly smaller than the unladen dissipation. This is perhaps not unexpected because the turbulence kinetic energy far away from particles is also reduced. The turbulence dissipation rate far away from particles is probably set by the turbulence kinetic energy far away and the outer scale of the turbulence far away. The latter can be estimated by $(K(r_{b}))^{3/2}/{\it\epsilon}(r_{b})$ , which is equal to $20.8$ in case L128, in fact very close to the global outer scale of the turbulence ( $L_{{\it\epsilon}}$ ) in the unladen case. Indeed, the outer scale far away from particles is set by the stochastic forcing, which is the same in the laden and unladen cases.

At this point it is appropriate to discuss the results obtained for the modification of $\overline{K}$ and $\overline{{\it\epsilon}}$ in relation to the two experimental works mentioned. Table 6 shows an overview of four experiments selected from these works. In the first two experiments listed (Hwang & Eaton 2006), glass beads of diameter 165  ${\rm\mu}$ m were used in air, while in the third (fourth) experiment (Tanaka & Eaton 2010), glass (polystyrene) beads of 250  ${\rm\mu}$ m were used in air. The parameter ${\rm\Delta}r$ denotes the minimum spatial resolution in the experiments and the simulation. The material density of the particles does not play a role in the simulation, since the particles are fixed. Thus the mass loading ( ${\it\psi}$ ) and the Stokes number with respect to Kolmogorov time ( $St$ ) are not defined in the simulation, while the terminal velocity $v_{t}$ is zero since gravity is ignored.

According to the particle-resolved simulation, $\overline{K}$ is reduced by 9 %, while $\overline{{\it\epsilon}}$ is increased by 3 %. However, in the experiments, $\overline{K}$ reduced more (except in the fourth experiment) and $\overline{{\it\epsilon}}$ also reduced (except in the fourth experiment), although the particle volume fraction in the experiments was smaller than in the simulation. Thus, it is surprising that the simulated attenuation of $\overline{K}$ is not much larger than 9 % and that the dissipation is not reduced. According to experiments, the degree of turbulence attenuation tends to increase with Stokes number (Kulick et al. 1994; Tanaka & Eaton 2010). Larger Stokes number usually implies a weaker response of the motion of the particles to the fluid and a larger particle–fluid slip velocity. Therefore, we do not expect that the simulated attenuation of $\overline{K}$ would have been stronger if the particles had not been fixed. However, the simulated attenuation of $\overline{K}$ could have been stronger if the terminal velocity had not been zero, i.e. if the effect of gravity had been included. The particle Reynolds number $Re_{p}$ based on the total velocity in the experiments was on average probably larger than $Re_{p}$ based on the relative fluctuating velocity and $Re_{p}$ based on the terminal velocity. Provided $Re_{p}$ remains small enough to avoid vortex shedding in the wakes of the particle, larger $Re_{p}$ is expected to lead to stronger turbulence attenuation. However, it seems unlikely that the quantitative discrepancy between experiments and simulations is only caused by $Re_{p}$ . Perhaps, the relatively weak attenuation of $\overline{K}$ in the simulations is due to the relatively low $Re_{{\it\lambda}}$ . Indeed, the maximum attenuation of $\overline{K}$ in the experiments of Tanaka & Eaton (2010) ( $Re_{{\it\lambda}}=127$ ) was slightly less than in the experiments of Hwang & Eaton (2006) ( $Re_{{\it\lambda}}=230$ ). Although the simulated flow at $Re_{{\it\lambda}}=32$ can be regarded as turbulent (see the discussion of the higher-order statistics of the flow in § 4.4), the number 32 is low and much lower than in the experiments. In addition, some quantitative effect of the limited size of the computational box cannot be ruled out. However, if $Re_{{\it\lambda}}$ is not increased simultaneously, a larger box size is not expected to reduce the discrepancy between simulation and experiments significantly.

With respect to the discrepancies in the turbulence dissipation rate, Tanaka & Eaton (2010) measured a much smaller reduction of $\overline{{\it\epsilon}}$ than Hwang & Eaton (2006) (in the fourth experiment listed in table 6, the turbulence dissipation rate did not reduce at all). Tanaka & Eaton (2010) attributed the higher $\overline{{\it\epsilon}}$ to the increase of the spatial resolution. The values ${\rm\Delta}r/{\it\eta}$ listed in table 6 show that in the high-resolution experiments performed by Tanaka & Eaton (2010) the turbulence dissipation rate was much better resolved indeed ( ${\rm\Delta}r=0.55{\it\eta}$ ). However, figure 7 shows that, if in the simulations ${\rm\Delta}r=0.55{\it\eta}\approx 0.5r_{0}$ had been chosen, 8 % of $\overline{{\it\epsilon}}$ would not have been captured (the amount of the dissipation in between $r_{0}$ and $1.5r_{0}$ ). Whether this means that even the high-resolution measurements by Tanaka & Eaton (2010) did not entirely capture $\overline{{\it\epsilon}}$ is difficult to say, since we do not know how figure 7 scales with decreasing ${\it\alpha}$ and increasing $Re_{p}$ and $Re_{{\it\lambda}}$ .

4.3 The radial turbulence kinetic energy budget

For spherical particles in statistically stationary homogeneous isotropic turbulence, the radial budget of turbulence kinetic energy is defined by

(4.8) $$\begin{eqnarray}P+T+{\it\Pi}+D-{\it\epsilon}=0,\end{eqnarray}$$

for all $r$ , where

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle P(r)=\langle u_{r}f_{r}+u_{{\it\theta}}f_{{\it\theta}}+u_{{\it\phi}}f_{{\it\phi}}\rangle , & \displaystyle\end{eqnarray}$$

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle T(r)=-\frac{1}{r^{2}}\frac{\text{d}}{\text{d}r}(r^{2}\langle (u_{r}^{2}+u_{{\it\theta}}^{2}+u_{{\it\phi}}^{2})u_{r}\rangle ), & \displaystyle\end{eqnarray}$$

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\Pi}(r)=-\frac{1}{r^{2}}\frac{\text{d}}{\text{d}r}(r^{2}\langle pu_{r}\rangle ), & \displaystyle\end{eqnarray}$$

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle D(r)=\frac{1}{r^{2}}\frac{\text{d}}{\text{d}r}\left(r^{2}{\it\nu}\frac{\text{d}K}{\text{d}r}\right), & \displaystyle\end{eqnarray}$$

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\epsilon}(r)={\it\nu}\langle \unicode[STIX]{x1D60E}_{ij}\unicode[STIX]{x1D60E}_{ij}\rangle . & \displaystyle\end{eqnarray}$$

The first term ( $P$ ) represents the production of $K=K(r)$ due to $\boldsymbol{f}$ , the last term ( ${\it\epsilon}$ ) the sink of $K$ due to dissipation of $K$ , while the three terms in between, $T$ , ${\it\Pi}$ and $D$ , are transport terms, which are called the turbulent transport term, the pressure diffusion term and the viscous diffusion term, respectively. It is remarked that the conventional production term in turbulent shear flows, the product of the mean shear and minus the Reynolds shear stress, does not appear in the balance because the mean velocity and the Reynolds shear stress are both zero in this flow.

Figure 8. Radial turbulence kinetic energy budget for cases L128 (blue, thick solid) and L64 (red, thin dashed). (a,b) Production $P$ (circles), turbulent transport $T$ (squares), pressure diffusion term ${\it\Pi}$ (stars), viscous diffusion $D$ (upward-pointing triangles), and minus the turbulence dissipation rate $-{\it\epsilon}$ (downward-pointing triangles): (a) normalized by the unladen dissipation rate ( ${\it\epsilon}_{0}$ ) and (b) divided by ${\it\epsilon}$ . (c) The balance error $E$ (diamonds) divided by ${\it\epsilon}$ .

The budget normalized with the unladen global dissipation rate ( ${\it\epsilon}_{0}$ ) is shown in figure 8(a). Both $D$ and ${\it\epsilon}$ are very large at $r_{0}$ ; the values in the grid cell next to the surface are: ${\it\epsilon}(1.033r_{0})=96{\it\epsilon}_{0}$ in case L128 and ${\it\epsilon}(1.068r_{0})=82{\it\epsilon}_{0}$ in case L64, which is far outside the range of figure 8(a). Therefore, the terms are also shown after division by ${\it\epsilon}(r)$ (figure 8 b). The production by the forcing term is zero at $r_{0}$ , since $\boldsymbol{u}=\mathbf{0}$ there. At the further locations $r>3.8r_{0}$ , we observe $P/{\it\epsilon}>1$ . Thus at most locations in the flow (the volume of a spherical shell is proportional to $r^{2}$ ), the production is not locally balanced by the dissipation. The sum of the transport terms, $T+{\it\Pi}+D={\it\epsilon}-P$ , is negative for $r>3.8r_{0}$ and positive for $r<3.8r_{0}$ . Therefore, the sum of the three transport terms transfers energy from the region $r>3.8r_{0}$ to the regions $r<3.8r_{0}$ , which surround the particles. The regions $r<3.8r_{0}$ correspond to 5.5 % of the flow domain. Thus, in the direct vicinity of the particles, turbulence attenuation is caused by the enhanced dissipation rate, but in the major part of the flow domain, the attenuation is caused by transport of kinetic energy towards the particles.

There are three transport terms. Sufficiently far away from the nearest particle surface, all transport terms are negative and take kinetic energy. However, moving towards the nearest particle surface, each transport term becomes positive at some point and starts to give kinetic energy. Turbulent transport ( $T$ ) is negative for $r>5.5r_{0}$ and positive for $r_{0}<r<5.5r_{0}$ ; pressure diffusion term ( ${\it\Pi}$ ) is negative for $r>3.1r_{0}$ and positive for $r<3.1r_{0}$ ; while viscous diffusion ( $D$ ) is negative for $r>1.5r_{0}$ and positive for $r<1.5r_{0}$ . The relative importance of the pressure diffusion term among the transport fluxes is surprising since the pressure diffusion transport term is considered to be the least relevant contribution in the turbulent kinetic energy budget of near-wall turbulence (Mansour et al. 1988). The maximum attained by the pressure diffusion term is 13 % of the maximum turbulence dissipation rate in the present flow. This is 6 % in turbulent channel flow at friction Reynolds number 180 (Mansour et al. 1988).

An interesting check on the statistical and numerical accuracy of the results is whether the balance error,

(4.14) $$\begin{eqnarray}E(r)=P+T+{\it\Pi}+D-{\it\epsilon},\end{eqnarray}$$

is sufficiently small, since it should theoretically vanish. Profiles of the relative balance error, $E(r)/{\it\epsilon}(r)$ , are shown in figure 8(c). The figure shows that L128 is much more accurate than L64, although the accuracy of L64 is not very poor since the relative balance error less than 5 % for all $r$ locations. However, the balance error of L128 is less than 1 %, except for a few $r$ locations. The maximum error (2 %) is attained at the first grid point off the particle surface (the error at the second grid point is only $-0.5\,\%$ ).

Figure 9. (a) Energy fluxes related to the three transport terms: turbulent transport flux (squares), the pressure diffusion flux (stars) and the viscous diffusion flux (triangles). The pressure diffusion flux is dominant over the other two in the region between the black thin dotted vertical demarcation lines. (b) Pressure velocity correlation coefficient ${\it\beta}$ . Both panels show results for cases L128 (blue, thick solid) and L64 (red, thin dashed).

To further investigate the energy transport due to particles, we integrate the three transport terms over the volume between $r_{0}$ and $r$ , multiply by $-1$ and divide by the surface area $4{\rm\pi}r^{2}$ to obtain the corresponding energy fluxes: the turbulent transport flux $\langle u_{i}u_{i}u_{r}\rangle$ , the pressure diffusion flux $\langle pu_{r}\rangle$ and the viscous diffusion flux $-{\it\nu}\,\text{d}K/\text{d}r$ . Figure 9(a) shows that the energy flux profiles are negative everywhere, which means that, for each transport term and for all $r$ , the radial flux of kinetic energy is on average directed towards the surface of the nearest particle. The demarcation lines separate three regions: (1)  $r_{0}<r<1.9r_{0}$ , dominated by the flux of viscous diffusion; (2)  $1.9r_{0}<r<4.7r_{0}$ , dominated by the flux of pressure diffusion; and (3)  $r>4.7r_{0}$ , dominated by the flux of turbulent transport. The volumes of the three regions correspond to 0.6 %, 9.9 % and 89.4 % of the total volume (0.1 % is occupied by particles).

Figure 9(b) shows the correlation coefficient

(4.15) $$\begin{eqnarray}{\it\beta}(r)=\frac{\langle pu_{r}\rangle }{\text{RMS}(p)~\text{RMS}(u_{r})}.\end{eqnarray}$$

The coefficient is interesting because it is the only non-zero coefficient of correlation between two different basic variables ( $\langle u_{r}u_{{\it\theta}}\rangle =\langle u_{r}u_{{\it\phi}}\rangle =\langle u_{{\it\theta}}u_{{\it\phi}}\rangle =0$ and $\langle pu_{{\it\theta}}\rangle =\langle pu_{{\it\phi}}\rangle =0$ ). The negative sign is consistent with the qualitative behaviour in, for example, turbulent channel flow and is intuitively not difficult to understand: owing to the impermeability of the surface, the pressure is increased at times when the instantaneous velocity vector is directed towards the surface. Figure 5 provides some further insight into why the correlation is likely to be negative: near the surface at the front side $u_{r}$ is negative and pressure fluctuation positive, while near the surface at the rear side $u_{r}$ is positive and pressure fluctuation negative. Thus the product of pressure fluctuation and $u_{r}$ tends to be negative at both sides, such that $\langle pu_{r}\rangle$ and ${\it\beta}$ become negative.

4.4 Higher-order statistics

Standardized higher-order moments, in particular the skewness and flatness factors of velocity, pressure and their derivatives, have been a subject of extensive study in the field of homogeneous isotropic turbulence (see e.g. Ishihara et al. 2007, and references therein). It is therefore interesting to investigate how these quantities are influenced by the particles. The global skewness $\overline{S}(q)$ is defined by $\overline{(q-\overline{q})^{3}}/(\overline{\text{RMS}}(q))^{3}$ and the radial skewness by $S(q)=\langle (q-\langle q\rangle )^{3}\rangle /(\text{RMS}(q))^{3}$ . The definitions of global flatness $\overline{F}(q)$ and radial flatness $F(q)$ are obtained if the third powers in the skewness expressions are replaced by fourth powers. The skewness and flatness factors characterize the shape of the probability density function (p.d.f.) of the variable $q$ . The skewness quantifies the asymmetry of the p.d.f. of $q$ , while the flatness quantifies the importance of the tails of the p.d.f. Higher flatness corresponds to stronger intermittency. A large flatness of a variable corresponds to a relatively large probability that extreme events occur, events in which the absolute value of the variable is much larger than the standard deviation of its probability distribution. The standardized third and fourth moments can also be used to quantify the non-Gaussianity of a variable. The skewness of a Gaussian variable is zero and the flatness equals three. The skewness and flatness of velocity derivatives are commonly used to describe the properties of the small-scale turbulence. Negative skewness of the longitudinal velocity derivative is related to the generation of small scales by vortex stretching and the flatness of the components of the velocity gradient to the intermittency of small-scale motions.

Global skewness and flatness factors of the forcing term, velocity, pressure and velocity derivatives are shown in table 7. In addition, the ratio $\overline{{\it\chi}}=\overline{u_{1,1}^{2}}~/~\overline{u_{1,2}^{2}}$ is shown, which according to the theory of single-phase incompressible homogeneous isotropic turbulence should be equal to $1/2$ . To limit the length of the table, the skewness factors of $f_{1}$ , $u_{1}$ and $u_{1,2}$ have not been included since these are theoretically zero and numerically small ( $|\overline{S}(f_{1})|<0.01$ , $|\overline{S}(u_{1})|<0.01$ and $|\overline{S}(u_{1,2})|<0.03$ ).

The results for $\overline{F}(f_{1})$ are very close to 3, which confirms that the temporal filtering of the Ornstein–Uhlenbeck process does not affect the Gaussianity of the stochastic process used to force the large scales. The derivative ratio $\overline{{\it\chi}}$ is equal to 0.474 in case U128, reasonably close to the theoretical isotropic value. In the literature on DNS of single-phase isotropic turbulence, skewness and flatness factors have been reported for single-phase isotropic turbulence at $Re_{{\it\lambda}}$ somewhat higher than 32: $\overline{F}(u_{1})=2.80$ , $-\overline{S}(u_{1,1})=0.49$ , $\overline{F}(u_{1,1})=4.2$ and $\overline{F}(u_{1,2})=5.7$ at $Re_{{\it\lambda}}=35$ (Jimenez et al. 1993); and $-\overline{S}(p)=0.88$ , $\overline{F}(p)=5.6$ and $\overline{F}(u_{1,1})=4.0$ at $Re_{{\it\lambda}}=38$ (Vedula & Yeung 1999). The numbers for U128 in table 7 are sufficiently close to these values to conclude that this case provides an acceptable description of homogeneous isotropic turbulence at $Re_{{\it\lambda}}=32$ . The quantitative discrepancies between table 7 and the values cited from the literature can perhaps be attributed to differences in $Re_{{\it\lambda}}$ and the forcing. Although the forcing has been applied to large scales only, some sensitivity of the small-scale turbulence to the forcing probably cannot be avoided if $Re_{{\it\lambda}}$ is low. For this reason, the forcing has fully been specified in § 3.

We proceed to discuss the effect of particles on the global statistics listed in table 7. The value of $\overline{F}(u_{1})$ is close to 3 in laden and unladen cases; the nearly Gaussian behaviour of the velocity is not significantly affected by particles. Likewise, the ratio $\overline{{\it\chi}}$ is not significantly affected by particles. However, the skewness and flatness factors of the pressure and the velocity derivatives are modified by particles. Both $-\overline{S}(p)$ and $\overline{F}(p)$ appear to be reduced. Thus not only is the size of the pressure field affected, but also the structure of the pressure field changes due to particles: it becomes less skewed and less intermittent. The negative skewness of the pressure in a turbulent flow means that local pressure minima, which usually occur at the cores of vortices, tend to be stronger (but also narrower) than the pressure maxima. The reduction of $-\overline{S}(p)$ is rather strong and seems to be larger than can be explained from the reduction of $Re_{{\it\lambda}}$ from 32 in the unladen to 29 in the laden cases. The reverse trend is visible in the skewness of the longitudinal velocity derivative: $-\overline{S}(u_{1,1})$ is slightly increased by particles, while this skewness normally decreases with decreasing $Re_{{\it\lambda}}$ (Ishihara et al. 2007). The largest effect of particles is observed in the global flatness factors of the velocity derivatives: $\overline{F}(u_{1,1})$ becomes 6 and $\overline{F}(u_{1,2})$ almost 8 times larger. These derivatives and thereby also the turbulence dissipation rate become much more intermittent in a particle-laden field. It means that if we sampled a Cartesian velocity gradient at a random location in the flow, the probability that we would find a very large velocity gradient (much larger than the global standard deviation) is much higher than in an unladen turbulent flow. The shape of the turbulence dissipation rate profile shows that these very strong gradients preferentially occur in the direct vicinity of the particles.

Figure 10. (a) Skewness of $u_{r}$ (circles) and $p$ (downward-pointing triangles) and (b) flatness of $u_{1}$ (stars), $u_{r}$ (circles), $u_{{\it\theta}}$ and $u_{{\it\phi}}$ (squares) and $p$ (downward-pointing triangles). (c) Coefficient ${\it\chi}$ (diamonds) and skewness of $u_{1,1}$ (plus signs). (d) Flatness of $u_{1,1}$ (plus signs) and $u_{1,2}$ (upward-pointing triangles). Results are from simulations L128 (blue, thick solid) and L64 (red, thin dashed). The corresponding unladen quantities from simulation U128 are denoted by black filled symbols. The black thin dotted horizontal lines in (b) and (d) represent the Gaussian flatness.

The question arises how skewness and flatness factors are modified if the p.d.f. is conditioned on the distance to the nearest particle. That information is provided by the radial skewness and flatness factors, shown in figure 10. It appears that the skewness of $u_{r}$ , approximately zero far away from particles, becomes more and more negative if a particle is approached. This means that an extremely large radial velocity tends to be negative (the meaning of extremely large is very large compared to the local standard deviation, which is a function of $r$ ). This tendency is probably due to the asymmetry between the flow at the front and rear sides of the non-Stokesian particles, with the result that an extremely large velocity directed towards the surface becomes more probable than an extremely large velocity directed away from the surface. The radial flatness factors (figure 10 b) show that the intermittency of all velocity components is enhanced near the surfaces of the particles, while the intermittency of the radial component is enhanced most. In contrast, the intermittency of the pressure reduces, at all $r$ locations, but near $r_{0}$ in particular. The radial profile ${\it\chi}(r)$ defined by $\langle u_{1,1}^{2}\rangle /\langle u_{1,2}^{2}\rangle$ is shown in figure 10(c). Although ${\it\chi}$ displays a significant variation around the isotropic value ( $1/2$ ), the local anisotropy in the Cartesian velocity gradient tensor seems to be surprisingly weak. Furthermore, $S(u_{1,1})$ is shown in figure 10(c) ( $S(u_{1,2})$ is approximately zero). The curve displays a remarkably strong oscillation: from approximately zero at $r_{0}$ it decreases sharply to approximately $-0.45$ , then it rises to almost zero, before it slowly decreases, until it has reached a level slightly below the unladen value.

The radial flatness factors of the Cartesian velocity derivatives are shown in figure 10(d). For all $r$ locations, the radial flatness appears to be much lower than the corresponding global flatness values listed in table 7. This is due to the fact that extremely large velocity gradients occur at preferential locations, namely near particles. Nonetheless, not only the global but also the radial intermittency of the velocity gradient tensor are stronger than in unladen turbulence, at least near particle surfaces.

Figure 11. (a) RMS ( $\unicode[STIX]{x1D60E}_{ij}$ ), (b) zoomed RMS ( $\unicode[STIX]{x1D60E}_{ij}$ ), (c) skewness $S(\unicode[STIX]{x1D60E}_{ij})$ and (d) flatness $F(\unicode[STIX]{x1D60E}_{ij})$ of the components of the gradient of the velocity in spherical coordinates: $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ (squares), $\unicode[STIX]{x1D60E}_{11}$ (circles), $\unicode[STIX]{x1D60E}_{22}$ and $\unicode[STIX]{x1D60E}_{33}$ (stars), $\unicode[STIX]{x1D60E}_{21}$ and $\unicode[STIX]{x1D60E}_{31}$ (upward-pointing triangles), and $\unicode[STIX]{x1D60E}_{23}$ and $\unicode[STIX]{x1D60E}_{32}$ (downward-pointing triangles). Results are from simulations L128 (blue, thick solid) and L64 (red, thin dashed). The corresponding unladen quantities from simulation U128 are denoted by black filled symbols; Cartesian unladen values are denoted by squares. The black thin dotted horizontal line in (d) represents the Gaussian flatness.

Additional insight into the structure of the small-scale turbulence around particles is obtained from the statistics of $\unicode[STIX]{x1D642}$ , the gradient of the velocity in spherical coordinates. These statistics are shown in figure 11. Whereas all components of the Cartesian gradient of the velocity are enhanced at surfaces of the particles ( ${\it\chi}(r_{0})\approx 1/2$ ), this is not the case for the components of $\unicode[STIX]{x1D642}$ . Figures 11(a) and 11(b) show that only $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ are enhanced; the other components are zero on the surfaces of the particles. According to (2.7), the components $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ are $u_{{\it\theta},r}$ and $u_{{\it\phi},r}$ , respectively. The boundary conditions imply that the other components of $\unicode[STIX]{x1D642}$ are zero at $r_{0}$ ( $\unicode[STIX]{x1D60E}_{11}=u_{r,r}$ is zero due to the incompressibility constraint). This is true for stagnant particles, but also for moving particles, provided each spherical coordinate system moves and rotates with the same translative and angular velocities as the corresponding particle. However, slightly off the particles surfaces ( $r\approx 1.3r_{0}$ ), all components are significantly enhanced compared to the unladen values, except $\unicode[STIX]{x1D60E}_{23}$ and $\unicode[STIX]{x1D60E}_{32}$ . At $r\approx 1.3r_{0}$ , the ordering from large to small is (1)  $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ , (2)  $\unicode[STIX]{x1D60E}_{11}$ , (3)  $\unicode[STIX]{x1D60E}_{22}$ and $\unicode[STIX]{x1D60E}_{33}$ , (4)  $\unicode[STIX]{x1D60E}_{21}$ and $\unicode[STIX]{x1D60E}_{31}$ , and (5)  $\unicode[STIX]{x1D60E}_{23}$ and $\unicode[STIX]{x1D60E}_{32}$ . Between $r\approx 1.55r_{0}$ and $r\approx 3.1r_{0}$ , the longitudinal component $\unicode[STIX]{x1D60E}_{11}$ is the largest one. It should be recalled that the sum of the squares of the components constitutes the local turbulence dissipation rate divided by the viscosity. The dominance of $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ near $r_{0}$ is caused by the viscous friction due to the no-slip boundary condition imposed on the tangential velocity components. This is also illustrated in figure 5(b): the largest dissipation rate occurs in thin shear layers near the particle, away from the front and rear sides (in the plane shown, $\unicode[STIX]{x1D60E}_{12}=u_{{\it\theta},r}$ is very large). Figure 5 also suggests that, at somewhat larger radius, the radial compression and expansion regions near the front and rear sides of a particle lead to large turbulence dissipation rate, which is consistent with the dominance of $\text{RMS}(\unicode[STIX]{x1D60E}_{11})=\text{RMS}(u_{r,r})$ for $1.55r_{0}<r<3.1r_{0}$ in figure 11(b).

The non-zero skewness factors and all flatness vectors of the components of $\unicode[STIX]{x1D642}$ are shown in figure 11(c,d). Near $r_{0}$ , $S(\unicode[STIX]{x1D60E}_{11})$ is strongly negative, while $S(\unicode[STIX]{x1D60E}_{22})$ and $S(\unicode[STIX]{x1D60E}_{33})$ are strongly positive. This indicates that, near particle surfaces, a strongly negative $\unicode[STIX]{x1D60E}_{11}$ (relatively strong radial compression) is more likely than a strongly positive $\unicode[STIX]{x1D60E}_{11}$ (relatively strong radial expansion). Owing to the continuity equation, the reverse must be true for the sum of $\unicode[STIX]{x1D60E}_{22}$ and $\unicode[STIX]{x1D60E}_{33}$ , and, owing to symmetry, the reverse is true for both $\unicode[STIX]{x1D60E}_{22}$ and $\unicode[STIX]{x1D60E}_{33}$ . The skewness factors of the non-diagonal components of $\unicode[STIX]{x1D642}$ are zero and not shown. The radial flatness factors of the components of $\unicode[STIX]{x1D642}$ (figure 11 d) seem to be lower than those of the Cartesian components of the velocity gradient, at least near $r_{0}$ (figure 10 d). All curves show a minimum, and the minimum of $F(\unicode[STIX]{x1D60E}_{11})$ is even lower than the Gaussian value. Very close to the surfaces of the particles, the intermittency is relatively large for the diagonal components. Apart from the $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ , the components of $\unicode[STIX]{x1D642}$ are small near $r_{0}$ , but in an intermittent way, the diagonal components in particular. Since the particle surfaces are solid impermeable no-slip walls, it is interesting to mention that also in turbulent channel flow the flatness factors of several spatial velocity derivatives peak at the walls (Vreman & Kuerten 2014b ).

At first glance, it is surprising that at $r_{0}$ the intermittency of the transverse Cartesian derivative $u_{1,2}$ is much larger than the intermittency of $\unicode[STIX]{x1D60E}_{12}=u_{{\it\theta},r}$ and $\unicode[STIX]{x1D60E}_{13}=u_{{\it\phi},r}$ , in particular when we realize that, at $r_{0}$ , the only non-zero components of $\unicode[STIX]{x1D642}$ are $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ . This is probably due to preferential concentration of large values of $u_{1,2}$ at specific locations on the particle surface, namely in the neighbourhood where the normal of the surface is approximately parallel to the $x_{2}$ -direction. This leads to a larger radial flatness of $u_{1,2}$ at $r_{0}$ since the radial averaging operator includes averaging over spherical surfaces. In contrast, the local p.d.f. of $\unicode[STIX]{x1D60E}_{12}$ (and $\unicode[STIX]{x1D60E}_{13}$ ) is nearly independent of the position on the particle surface.

We finish this subsection with analytical expressions for the radial flatness of Cartesian derivatives and the value of ${\it\chi}$ at $r_{0}$ . We denote the time averaging operator at $r_{0}$ by $\langle \cdot \rangle (r_{0},{\it\theta},{\it\phi})$ , where ${\it\theta}$ and ${\it\phi}$ are the coordinates of the spherical reference frame of the nearest particle. The radial averaging operator at $r_{0}$ is denoted by $\langle \cdot \rangle (r_{0})$ . We assume that the local p.d.f. of $\unicode[STIX]{x1D60E}_{ij}$ does not depend on the position on the surface. Thus the $k$ th moment of $\unicode[STIX]{x1D60E}_{ij}$ satisfies

(4.16) $$\begin{eqnarray}\langle \unicode[STIX]{x1D60E}_{ij}^{k}\rangle (r_{0},{\it\theta},{\it\phi})=\langle \unicode[STIX]{x1D60E}_{ij}^{k}\rangle (r_{0}).\end{eqnarray}$$

In the present case the assumption is only approximately valid since the particles are ordered in a lattice. The components of $\unicode[STIX]{x1D642}$ are zero at $r_{0}$ , except $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ . Since $\unicode[STIX]{x1D60E}_{12}$ and $\unicode[STIX]{x1D60E}_{13}$ behave statistically in the same way, we define

(4.17) $$\begin{eqnarray}M_{k}=\langle \unicode[STIX]{x1D60E}_{12}^{k}\rangle (r_{0})=\langle \unicode[STIX]{x1D60E}_{13}^{k}\rangle (r_{0}).\end{eqnarray}$$

Note that $M_{1}=0$ . With the use of (2.8), we derive for the moments of the Cartesian derivatives

(4.18) $$\begin{eqnarray}\displaystyle \langle u_{j,i}^{k}\rangle (r_{0}) & = & \displaystyle \langle (\unicode[STIX]{x1D60E}_{12}\unicode[STIX]{x1D608}_{1i}\unicode[STIX]{x1D608}_{2j}+\unicode[STIX]{x1D60E}_{13}\unicode[STIX]{x1D608}_{1i}\unicode[STIX]{x1D608}_{3j})^{k}\rangle (r_{0})\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{4{\rm\pi}}\int _{0}^{2{\rm\pi}}\!\int _{0}^{{\rm\pi}}\langle (\unicode[STIX]{x1D60E}_{12}\unicode[STIX]{x1D608}_{1i}\unicode[STIX]{x1D608}_{2j}+\unicode[STIX]{x1D60E}_{13}\unicode[STIX]{x1D608}_{1i}\unicode[STIX]{x1D608}_{3j})^{k}\rangle (r_{0},{\it\theta},{\it\phi})\sin {\it\theta}\,\text{d}{\it\theta}\,\text{d}{\it\phi}.\end{eqnarray}$$

The integrals can be evaluated after substitution of (2.4) and with use of (4.16) and (4.17) and the notion that any component of $\unicode[STIX]{x1D63C}$ is not affected by the time average $\langle \cdot \rangle (r_{0},{\it\theta},{\it\phi})$ . The integral over ${\it\theta}$ reduces to an integral over a polynomial in $y$ after the substitution $y=\cos {\it\theta}$ . For the evaluation of the integral over ${\it\phi}$ , the recursive relation

(4.19) $$\begin{eqnarray}\displaystyle \int _{0}^{2{\rm\pi}}\sin ^{n}{\it\phi}\cos ^{m}{\it\phi}\,\text{d}{\it\phi} & = & \displaystyle \frac{n-1}{n+m}\int _{0}^{2{\rm\pi}}\sin ^{n-2}{\it\phi}\cos ^{m}{\it\phi}\,\text{d}{\it\phi}\nonumber\\ \displaystyle & = & \displaystyle \frac{m-1}{n+m}\int _{0}^{2{\rm\pi}}\sin ^{n}{\it\phi}\cos ^{m-2}{\it\phi}\,\text{d}{\it\phi}\end{eqnarray}$$

is useful. This integral is zero when both $m$ and $n$ are odd. We thus obtain

(4.20a,b ) $$\begin{eqnarray}\langle u_{1,1}^{2}\rangle (r_{0})={\textstyle \frac{2}{15}}M_{2},\quad \langle u_{1,2}^{2}\rangle (r_{0})={\textstyle \frac{4}{15}}M_{2},\end{eqnarray}$$

and therefore ${\it\chi}(r_{0})=1/2$ . Figure 7(c) is consistent with this result. Similarly the fourth moments can be rewritten as

(4.21) $$\begin{eqnarray}\displaystyle & \displaystyle \langle u_{1,1}^{4}\rangle (r_{0})=\langle u_{3,3}^{4}\rangle (r_{0})={\textstyle \frac{8}{315}}M_{4}, & \displaystyle\end{eqnarray}$$

(4.22) $$\begin{eqnarray}\displaystyle & \displaystyle \langle u_{1,2}^{4}\rangle (r_{0})=\langle u_{3,2}^{4}\rangle (r_{0})={\textstyle \frac{16}{105}}M_{4}. & \displaystyle\end{eqnarray}$$

These expressions and those in (4.20) imply that the corresponding radial flatness factors of the Cartesian derivatives satisfy

(4.23) $$\begin{eqnarray}\displaystyle & \displaystyle F(u_{1,1})(r_{0})={\textstyle \frac{10}{7}}F(\unicode[STIX]{x1D60E}_{12})(r_{0}), & \displaystyle\end{eqnarray}$$

(4.24) $$\begin{eqnarray}\displaystyle & \displaystyle F(u_{1,2})(r_{0})={\textstyle \frac{15}{7}}F(\unicode[STIX]{x1D60E}_{12})(r_{0}). & \displaystyle\end{eqnarray}$$

Figures 10(d) and 11(d) are consistent with these results since very close to $r_{0}$ the intermittency (flatness factor) of $u_{1,2}$ is indeed one and a half times higher than the intermittency of $u_{1,1}$ and more than twice as high as the intermittency of $\unicode[STIX]{x1D60E}_{12}$ .

4.5 A posteriori testing of the point-particle model

Point-particle simulations, also called Eulerian–Lagrangian simulations, are based on an empirical model for the particle force, usually the standard drag law, also called Schiller–Naumann correlation, which is defined by

(4.25) $$\begin{eqnarray}\boldsymbol{F}^{p,SN}=3{\rm\pi}d_{p}{\it\rho}{\it\nu}(1+0.15Re_{p}^{0.687})(\boldsymbol{u}-\boldsymbol{v}),\end{eqnarray}$$

where $Re_{p}=d_{p}|\boldsymbol{u}-\boldsymbol{v}|/{\it\nu}$ and $\boldsymbol{v}$ is the particle velocity. The point-particle model for the fluid motion reads

(4.26) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(4.27) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{,t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{u})=-\boldsymbol{{\rm\nabla}}p+{\it\nu}{\rm\nabla}^{2}\boldsymbol{u}+\boldsymbol{f}+\boldsymbol{f}^{p}, & \displaystyle\end{eqnarray}$$

(4.28) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}^{p}=-\mathop{\sum }_{p}{\it\delta}(\boldsymbol{x}-\boldsymbol{x}^{p})\boldsymbol{F}^{p,SN}, & \displaystyle\end{eqnarray}$$

where ${\it\delta}$ is the delta function. A straightforward implementation of the feedback force $\boldsymbol{f}^{p}$ on a staggered grid consists of four steps: (1) the fluid velocity $\boldsymbol{u}$ is linearly interpolated to each particle location, (2) each particle force is computed, (3) each particle force is distributed to the eight surrounding centres of the pressure cells, and (4) the Eulerian force field is then linearly interpolated from the pressure locations to the staggered velocity locations. The simulation following this procedure using the grid of U128 is called PP0.

Some modellers will object that the grid size of $d_{p}/2$ in simulation PP0 is too small because the point-particle method should only be applied if the grid size is at least a few times larger than the particle size. However, if the grid size is increased, the turbulence will be resolved less accurately, which is undesirable since we wish to use simulation U128 in our comparison. A logical approach is then to modify step 3 in the described point-particle procedure and to distribute each particle force over more cells than just the eight surrounding cells. This has been reported in the literature: Link et al. (2005) distributed each particle force over a cube of $5^{3}$ cubic cells around each particle. This approach can be formalized through the application of a spatial convolution filter $H(\boldsymbol{x}-\boldsymbol{y})$ with filter width ${\it\Delta}$ to equation (4.28). We choose the top-hat filter, which essentially means that the particle force is distributed over a cube with volume ${\it\Delta}^{3}$ around each particle. For example, the top-hat filter kernel is defined by $H(\boldsymbol{x}-\boldsymbol{y})=1/{\it\Delta}^{3}$ if $|x_{i}-y_{i}|<{\it\Delta}/2$ ( $i=1,2,3$ ) and $H(\boldsymbol{x}-\boldsymbol{y})=0$ otherwise. Application of the top-hat filter to (4.28) leads to

(4.29) $$\begin{eqnarray}\boldsymbol{f}^{p}=-\mathop{\sum }_{p}H(\boldsymbol{x}-\boldsymbol{x}^{p})\boldsymbol{F}^{p,SN}.\end{eqnarray}$$

Two point-particle simulations were performed with (4.28) replaced by (4.29): PP2 ( ${\it\Delta}=2d_{p}$ ) and PP4 ( ${\it\Delta}=4d_{p}$ ). Since PP0 represents a discretization of the limit ${\it\Delta}\rightarrow 0$ , we say that ${\it\Delta}=0$ in PP0. However, the effective ${\it\Delta}$ cannot be smaller than the grid size $h_{1}=d_{p}/4$ . Owing to interpolation to eight surrounding pressure locations in case PP0, PP0 in fact corresponds to ${\it\Delta}\approx d_{p}/2$ .

The three point-particle simulations were performed on a uniform grid with $128^{3}$ cells (like U128), but, to reduce the computation time, the iterative Poisson solver was replaced by a direct solver based on fast Fourier transforms. In addition, the time step was increased to 0.002, while the first-order explicit temporal discretization of the viscous terms was replaced by the second-order explicit Adams–Bashforth method. The forcing function and the statistical averaging time were the same as in simulations U128 and L128. It was verified that an unladen simulation using the same time step, time discretization and Poisson solver as the point-particle simulations produced results very similar to those of simulation U128 ( $\overline{K}$ became 59.33 instead of 59.03).

In the comparison presented below, the results of the particle-resolved simulation L128 are considered as the standard. Using the terminology common in the literature of large-eddy simulations, we could call a comparison of results of a point-particle simulation against results of a corresponding particle-resolved DNS an a posteriori test. An a priori test would then be a comparison in which the result of a model for the drag force (or another quantity of interest) is evaluated by inserting the fields of the particle-resolved DNS into the model expression and then compared to the corresponding particle-resolved result.

Figure 12. (a) Radial profiles of turbulence kinetic energy $K$ , normalized by reference values of the unladen simulation U128: L128 (blue dashed), PP0 (black solid and circles), PP2 (red solid) and PP3 (red dash-dotted). The thin dotted horizontal line represents the normalized unladen turbulence kinetic energy. (b–d) Radial profiles of the turbulence dissipation rate of point-particle simulations PP0 (b), PP2 (c) and PP4 (d), for which ${\it\epsilon}$ (black solid line and open symbols) has been decomposed into a resolved part ${\it\epsilon}^{(1)}$ (red solid) and a Schiller–Naumann part ${\it\epsilon}^{(2)}$ (red dashed). The ${\it\epsilon}$ value from the point-particle simulations should be compared to ${\it\epsilon}$ from the particle-resolved simulation L128 (blue dashed).

We define the global turbulence dissipation rate in a point-particle simulation by

(4.30) $$\begin{eqnarray}\overline{{\it\epsilon}}=\overline{{\it\epsilon}}^{(1)}+\overline{{\it\epsilon}}^{(2)},\quad \text{where }\overline{{\it\epsilon}}^{(1)}={\it\nu}\overline{u_{j,i}u_{j,i}}\text{ and }\overline{{\it\epsilon}}^{(2)}=-\overline{u_{j}\,f_{j}^{p}}.\end{eqnarray}$$

The turbulence dissipation rate thus consists of a resolved part ( $\overline{{\it\epsilon}}^{(1)}$ ) and a Schiller–Naumann part ( $\overline{{\it\epsilon}}^{(2)}$ ). The corresponding radial profiles are defined analogously:

(4.31) $$\begin{eqnarray}{\it\epsilon}(r)={\it\epsilon}^{(1)}(r)+{\it\epsilon}^{(2)}(r),\quad \text{where }{\it\epsilon}^{(1)}(r)={\it\nu}\langle u_{j,i}u_{j,i}\rangle \text{and }{\it\epsilon}^{(2)}(r)=-\langle u_{j}\,f_{j}^{p}\rangle .\end{eqnarray}$$

They are shown in figure 12, together with the profiles of the turbulence kinetic energy. It is observed that, with increasing ${\it\Delta}$ , the minimum of $K(r)$ increases, the maximum of ${\it\epsilon}(r)$ reduces, while ${\it\epsilon}^{(1)}(r)$ becomes nearly flat and the Schiller–Naumann contribution to the turbulence dissipation rate, ${\it\epsilon}^{(2)}(r)$ , increases. The negativity of ${\it\epsilon}^{(2)}(r)$ in case PP0 is remarkable since it means that in the limit of ${\it\Delta}\rightarrow 0$ the Schiller–Naumann contribution does not dissipate turbulence kinetic energy but generates it.

The results of the global statistics for PP0, PP2 and PP4 are shown in tables 3–5 and 7. The attenuation of $\overline{K}$ is 4 %, 7 % and 8 % in cases PP0, PP2 and PP4, respectively, which should be compared to 9 %, the attenuation in the PR-DNS (L128). The results in the tables show that the effect of the particles on the turbulence is not properly captured by PP0, but most predictions improve with increasing ${\it\Delta}$ . In fact, the overall agreement between PP4 and L128 is reasonably good. An exception are the flatness factors of the velocity derivatives, which are strongly increased in both PP0 and L128 ( $\overline{F}(u_{1,2})$ in PP0 is even larger than the corresponding PR-DNS value). However, no increase is found in PP4 (and PP2), due to the smoothing effect of the top-hat filter on the peaks created by the delta functions in the feedback force.

The differences between PP0, PP2 and PP4 show that the point-particle method is sensitive to the chosen approximation of the delta function. Without any filter the point-particle method depends on the grid size because the effective volume over which each particle force is distributed is then set by the volume of the grid cell. This introduces arbitrariness into the method because the optimum value of the grid size, or the optimum value of ${\it\Delta}$ if the delta function force field is explicitly filtered, is not known a priori and may vary from case to case. Nonetheless, the fact that one of the point particle simulations (PP4) is able to deliver reasonable predictions of the turbulence is surprising since the point-particle model is formally valid for $d_{p}\ll {\it\eta}$ only, and this condition is not satisfied in the present case. It is therefore interesting to investigate the accuracy of the Schiller–Naumann correlation and its contribution to the dissipation rate in further detail.

For a turbulent flow around a single particle with diameter $d\approx 2{\it\eta}$ , a case in which the global turbulence kinetic energy was not modified, Burton & Eaton (2005) reported that the Schiller–Naumann correlation underpredicts the magnitude of the true particle force. In addition, they found that the relative error of the force, defined as the relative difference between the model signal and the signal of the true force, varied between 15 % and 30 %. Table 5 shows that the magnitude ( $L_{2}$ -norm) of the Schiller–Naumann force is also underpredicted in the present assessment of the point-particle method: the values are 97.3, 148.2 and 171.2 in simulations PP0, PP2 and PP4, respectively, while the PR-DNS result is 188.9. Thus the particle force is underpredicted by 48 %, 22 % and 9 % in cases PP0, PP2 and PP4. The degree of underprediction of the turbulence attenuation roughly corresponds to the degree of underprediction of the particle force in each point-particle simulation.

In the point-particle simulations, the fluid velocity at the particle location is used to compute the Schiller–Naumann force. The reduction of the velocity at the particle location (see figure 12 a) is the reason why the force reduces with decreasing ${\it\Delta}$ and one of the reasons why the magnitude of the Schiller–Naumann force is too low in cases PP0 and PP2. The magnitude of the Schiller–Naumann force based on the global magnitude of the velocity from the PR-DNS, $(2\overline{K})^{1/2}$ from simulation L128, is equal to 170.6, very close to the prediction by PP4. The magnitude of the Schiller–Naumann force based on the unladen velocity is higher, of course, but it seems questionable to base the Schiller–Naumann force upon the unladen velocity if the velocity of the laden flow is globally modified by the turbulence.

No evidence was found that the underprediction of the particle force by the Schiller–Naumann correlation is caused by hydrodynamic particle–particle interactions. For uniform flow over a regular particle array, the hydrodynamic actions can be very strong (Hill et al. 2001). Figure 12 in that paper is particularly interesting in the present context, as it shows the drag force obtained from particle-resolved simulations of uniform flows over a simple cubic particle array for a volume fraction of 0.001 and a particle Reynolds number around 10. That figure applies to uniform flows pointing in five different directions, and it shows that, when the uniform flow was directed along the unit vector $\boldsymbol{e}_{1}$ , the force was roughly $30\,\%$ lower, while when the uniform flow was directed along $\boldsymbol{e}_{1}+\boldsymbol{e}_{2}$ , the force was roughly 15 % lower than in the other three cases. To investigate whether a similar dependence on the angle of the velocity over each particle occurs in the present turbulent flow, the magnitude and the cosine of the angle of the instantaneous particle force vectors from case L128 were binned into 20 uniform bins, spanning the range of the cosine from $-1$ to 1. The angle was defined as the angle between the force vector and a Cartesian unit vector, first $\boldsymbol{e}_{1}$ , and then $\boldsymbol{e}_{2}$ and $\boldsymbol{e}_{3}$ (to improve the statistics). For each bin the $L_{2}$ -norm of the particle force deviated less than 1 % from the overall magnitude of the force ( $188.9$ ). Thus particle–particle hydrodynamic interactions in the present turbulent flow seem to be much weaker than in some uniform flows over simple cubic particle arrays.

The particle-induced dissipation in the particle-resolved case is defined by (4.7) and is approximately 14 % of the global turbulence dissipation rate. The particle-induced dissipation in the point-particle simulations is defined as the sum of two effects, the global Schiller–Naumann dissipation term, $\overline{{\it\epsilon}}^{(2)}=-\overline{u_{j}\,f_{j}^{p}}$ , and the integral of the extra dissipation near particles observed in the profile ${\it\epsilon}^{(1)}(r)$ . Thus in the point-particle simulations

(4.32) $$\begin{eqnarray}\overline{{\it\epsilon}}^{p}=\overline{{\it\epsilon}}^{(2)}+\frac{N_{p}}{(1-{\it\alpha})L_{1}^{3}}\int _{r_{0}}^{R}4{\rm\pi}r^{2}({\it\epsilon}^{(1)}(r)-{\it\epsilon}^{(1)}(R))\,\text{d}r,\end{eqnarray}$$

where $R=7r_{0}$ . The Schiller–Naumann term, $\overline{{\it\epsilon}}^{(2)}$ , is equal to $-0.76$ , 2.56 and 2.99 in cases PP0, PP2 and PP4, respectively. The second term on the right-hand side of (4.32) is equal to 1.50, $-0.06$ and $-0.09$ in cases PP0, PP2 and PP4, and thus significant only in case PP0. The particle-induced dissipation itself (the sum of the two terms) is approximately 3 % of the turbulence dissipation rate in case PP0, 11 % in case PP2, and 13 % in case PP4 (table 5). It underpredicts the PR-DNS value in each case, but the underprediction in case PP4 is small.

We conclude that, provided the size of fluid volume over which each point-particle is distributed is suitably chosen, the point-particle model based on the Schiller–Naumann correlation provides reasonable predictions for the effect of two-way coupling on the turbulence, at least for the flow case investigated in this paper. The differences between point-particle and particle-resolved global quantities are noticeable, but most of them are much smaller than expected, at least in case PP4. In particular, the underpredictions of the turbulence attenuation and the particle-induced dissipation are no more than 10 % in case PP4. This is much better than expected since, as we mentioned in the introduction, Hwang & Eaton (2006) concluded from measurements that not only the attenuation of turbulence kinetic energy but also the extra dissipation due to particles were greatly underpredicted by conventional models. However, in the present simulations, a strong underprediction of the extra dissipation is only observed in case PP0, in which each particle force is fed back to a relatively small fluid volume, much smaller than in common Eulerian–Lagrangian simulations. The extra dissipation in the experiments was not directly measured but computed from the overall energy balance, by subtracting the measured turbulence dissipation rate from the net input of potential energy and energy from the synthetic jet actuators. However, using the same set-up with a high-resolution particle image velocimetry system, Tanaka & Eaton (2010) measured significantly higher turbulence dissipation rates than Hwang & Eaton (2006) did. If the turbulence dissipation rate measured by Hwang & Eaton (2006) was too low, the estimate of the extra dissipation due to particles was perhaps too high.

The fact that replacing the true particle forces by point-particle forces applied to a sufficiently large volume of fluid around the particle does not alter the most relevant turbulence modification results by more than 10 % suggests that each particle experiences the ambient turbulent flow basically as a (time-varying) uniform flow in the case that $d_{p}/{\it\eta}\approx 2$ and $Re_{{\it\lambda}}$ is low. To further investigate this hypothesis, we derive an expression for the turbulence dissipation rate on particle surfaces, assuming uniform flow around each particle. For uniform Stokes flow past a sphere, the analytical solution is known. It is straightforward to derive that the corresponding dissipation rate averaged over the particle surface is equal to $6{\it\nu}(U_{\infty }/d_{p})^{2}$ , where $U_{\infty }$ is the ambient (relative) velocity far away from the sphere. For fixed $U_{\infty }$ , the dissipation rate increases if $d_{p}$ reduces. Compared to the Stokes drag force, the drag force for a uniform flow at non-zero Reynolds number increases by a factor $1+0.15Re_{p}^{0.687}$ , according to the Schiller–Naumann correlation. Since the shear is expected to increase by approximately the same factor, the dissipation rate on the surface is expected to increase by the square of that factor. After identifying $U_{\infty }$ with $(2\overline{K})^{1/2}$ , we arrive at the following (rough) estimate of the turbulence dissipation rate on the surface of a fixed particle:

(4.33) $$\begin{eqnarray}{\it\epsilon}(r_{0})\approx 12{\it\nu}\frac{\overline{K}}{d_{p}^{2}}(1+0.15Re_{p}^{0.687})^{2}.\end{eqnarray}$$

For the second test case in § 2.5, uniform flow over a particle at $Re_{p}=100$ , the expression was found to estimate the average dissipation rate on the surface with an error of 1 %. For simulation L128, in which $\overline{K}\approx 59$ and $Re_{p}\approx 10.4$ , the expression leads to ${\it\epsilon}(r_{0})\approx 1970$ , approximately 25 % smaller than $2610$ , the simulated value. Although the underprediction is significant, the estimate based on the assumption of uniform flow is able to capture the correct order of magnitude of the strongly enhanced turbulence dissipation rate on the particle surfaces.

Although it is outside the scope of this paper to test many different options of the point-particle formulation (the scope was set by the research questions in the introduction), a few other options are briefly discussed. The top-hat filter was chosen because of its compact support and its separability into three one-dimensional filters (properties relevant for computational efficiency). Although the effect of the filter size is expected to be more important than the effect of the shape of the filter kernel, the selection of another shape (for example, spherically Gaussian) could lead to further improvement. Second, one could include forces other than the drag force into the model. In fact, a variant of simulation PP4 was performed with inclusion of the lift force. Using the common lift coefficient of 0.5, the results were hardly different from those of the original PP4. Inertial, mass and history forces were not included, since Bagchi & Balachandar (2003) reported that the best approximation to their DNS result was obtained without these forces. Third, instead of the fluid velocity at the particle location, one could insert into the Schiller–Naumann expression the fluid velocity averaged over a volume around the point particle (Bagchi & Balachandar 2003). To test the latter option for the present flow, another variant of simulation PP4 was performed, in which the averaging operator (top-hat filter) was applied not only to the force but also to the fluid velocity inserted into the force. The results were slightly but not significantly different from those of the original PP4.

5 Conclusions

A statistically stationary homogeneous isotropic turbulent flow modified by 64 small fixed non-Stokesian spherical particles was investigated by means of particle-resolved DNS. The particle diameter $d_{p}$ was approximately twice the Kolmogorov length scale ${\it\eta}$ , while the particle volume fraction was approximately 0.001. The Taylor Reynolds number of the corresponding unladen flow was approximately 32. The DNS of the particle-laden flow was based on a discretization of the incompressible Navier–Stokes equations on 64 spherical grids overset on a Cartesian grid. The numerical method was described in detail because non-standard procedures were required to obtain an overset grid method that could be applied to a turbulent flow with multiple small spherical particles. The sensitivity of the numerical results to resolution and statistical averaging time was addressed. In the following five paragraphs, following the order of the five research questions formulated in the introduction, we summarize the conclusions of this work.

Radial statistics, statistics that are functions of the distance to the nearest particle, showed that the kinetic energy is attenuated in the entire domain. The radial velocity variance is attenuated more than the azimuthal and polar velocity variances. The turbulence dissipation rate on particle surfaces is enhanced by two orders of magnitude, much larger than any local enhancement of the turbulence dissipation rate previously reported. More than 5 % of the total dissipation occurs in only 0.1 % of the flow domain. An equation that provides a rough estimate of the turbulence dissipation rate on particle surfaces was derived.

The global (space- and time-averaged) turbulence kinetic energy is attenuated by approximately 9 %, which is less than expected. The global turbulence dissipation rate is slightly increased by the particles. The particle-induced dissipation, the extra dissipation attributed to particles, was estimated to be 14 % of the global turbulence dissipation rate. The quantitative discrepancies between the results of the present particle-resolved simulations and the results of experiments on statistically steady homogeneous isotropic turbulence modified by particles were discussed.

In addition, the budget of the turbulence kinetic energy was computed as a function of the distance to the nearest particle centre. The budget illustrates how energy relatively far away from particles is transported towards the surfaces of the particles, where it is dissipated by the (locally enhanced) turbulence dissipation rate. The energy flux towards the particles is dominated by turbulent transport relatively far away from particles, by viscous diffusion very close to the particles and by pressure diffusion in a significant region in between. The volume where the pressure diffusion flux dominates over the other two transport fluxes is larger than expected and covers approximately 10 % of the entire flow volume.

To investigate the modification of the small-scale turbulence, the variances and higher-order statistics of the components of the gradient of the velocity were computed. The global flatness factor of the longitudinal velocity gradient, which characterizes the intermittency of small scales, is enhanced by a factor of six. The radial profiles of the higher-order turbulence statistics were also investigated. Significant modification of the statistical properties of the unladen small-scale turbulence was found, in particular near particles. Analytical relations were derived to explain the unexpected difference between moments of Cartesian and moments of spherical velocity derivatives on particle surfaces.

The same flow was also simulated by three point-particle simulations based on the Schiller–Naumann drag correlation. Three different sizes for the fluid volume over which each particle volume was distributed were used. The conclusion of the a posteriori tests, comparisons in which the particle-resolved results are regarded as the standard, is that, in this case, the point-particle model captures both the turbulence attenuation and the fraction of the turbulence dissipation rate due to particles reasonably well, provided the (arbitrary) size of the fluid volume over which each particle force is distributed is sufficiently large (for example $64d_{p}^{3}$ ). This is a surprising conclusion since $d_{p}\approx 2.2{\it\eta}>{\it\eta}$ . It was expected from the literature that for this particle size all three point-particle simulations would severely underpredict the turbulence attenuation and the part of the turbulence dissipation rate due to the particles. In the point-particle simulation in which each particle volume was distributed over a fluid volume of $64d_{p}^{3}$ , the magnitude of the true particle force was underpredicted by 9 %. This is in line with the conclusion of Bagchi & Balachandar (2003), who showed, for a single fixed particle subjected to a stationary frozen turbulent flow, that the Schiller–Naumann correlation provides a reasonably accurate prediction of the particle force if $d_{p}/{\it\eta}=1.5$ . It is also in line with the conclusion of Burton & Eaton (2005), who showed, for a single fixed particle in isotropic turbulence, that the Schiller–Naumann correlation tends to underpredict the magnitude of the particle force if $d_{p}/{\it\eta}=2$ . Bagchi & Balachandar (2003) also showed that the errors in the Schiller–Naumann correlation increase with increasing $d_{p}/{\it\eta}$ . Botto & Prosperetti (2012) performed particle-resolved simulations of turbulent flow past one or nine fixed particles for $6<d_{p}/{\it\eta}<10$ . They found that, for this range of particle sizes, the dissipation induced by a particle can instantaneously be very different from the work done by the particle force (in the reference frame of the fluid).

We finish the paper with a few remarks on possible extensions of the research, extensions that fall outside the scope of the present paper since they increase the computational complexity and demand. In this paper we discussed the quantitative discrepancies between the results of the present particle-resolved simulations and the results of high-resolution experiments on statistically steady homogeneous isotropic turbulence modified by particles. In view of this discussion, future particle-resolved simulations of turbulent flows laden with small particles could focus on the effects of (1) gravity, (2) higher $Re_{{\it\lambda}}$ and (3) particles that do respond to the flow.

Assuming that computer resources continue to increase, the long-term goal could be to use particle-resolved DNS to study the effect of a small concentration of small heavy particles on a more realistic turbulent flow, for example a turbulent channel flow. Kulick et al. (1994) measured strong turbulence attenuation in a gas–solid turbulent channel flow with a particle volume fraction of 0.0001, particle size $d_{p}^{+}=2.3$ and $Re_{{\it\tau}}=640$ . Recently, a point-particle DNS of this case was performed in a configuration that contained roughly $10^{5}$ particles (Vreman 2015). However, significant differences from the experimental results were still observed. To simulate such a configuration with the present particle-resolved method, roughly $3\times 10^{10}$ points for the uniform Cartesian and $1.5\times 10^{10}$ points for the overset grids would be required to achieve a resolution with $d_{p}/h=30$ near the particles. If the same case with the same $d_{p}/h$ were simulated on a uniform grid (for example by an immersed boundary method equipped with a pressure solver based on fast Fourier transforms), many more grid points would be required (roughly $10^{13}$ ). Despite the disadvantage that the pressure Poisson equation needs to be solved iteratively, the overset method seems to be a promising method for future particle-resolved simulation of (dilute) particle-laden flows.