Hostname: page-component-7c8c6479df-nwzlb Total loading time: 0 Render date: 2024-03-28T01:43:43.097Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  28 April 2016

A. W. Vreman*
Affiliation:
AkzoNobel, Research Development and Innovation, Process Technology, PO Box 10, 7400 AA Deventer, The Netherlands
*
Email address for correspondence: bert@vremanresearch.nl

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.

Type
Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© 2016 Cambridge University Press

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 Reference Tsuji, Morikawa and Shiomi1984; Gore & Crowe Reference Gore and Crowe1989; Hetseroni Reference Hetseroni1989). The phenomenon has been observed experimentally in pipe flows (Tsuji et al. Reference Tsuji, Morikawa and Shiomi1984), channel flows (Kulick, Fessler & Eaton Reference Kulick, Fessler and Eaton1994; Kussin & Sommerfeld Reference Kussin and Sommerfeld2002) and stationary homogeneous turbulence (Hwang & Eaton Reference Hwang and Eaton2006; Tanaka & Eaton Reference Tanaka and Eaton2010). 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 (Reference Tanaka and Eaton2010), who refined the experiments performed by Hwang & Eaton (Reference Hwang and Eaton2006), 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 Reference Squires and Eaton1990; Elghobashi & Truesdell Reference Elghobashi and Truesdell1993; Ferrante & Elghobashi Reference Ferrante and Elghobashi2003). 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 Reference Vreman2007, Reference Vreman2015).

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 Reference Balachandar and Eaton2010; Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011; Prosperetti Reference Prosperetti2015). Bagchi & Balachandar (Reference Bagchi and Balachandar2003) and Burton & Eaton (Reference Burton and Eaton2005) 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 (Reference Burton and Eaton2005) 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 Reference Uhlmann2008; Picano, Breugem & Brandt Reference Picano, Breugem and Brandt2015) and homogeneous isotropic turbulence (Lucci, Ferrante & Elghobashi Reference Lucci, Ferrante and Elghobashi2010; Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015). 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. (Reference Picano, Breugem and Brandt2015) 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. Reference ten Cate, Derksen, Portela and van den Akker2004), 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 Reference Hill, Koch and Ladd2001; Tenneti, Garg & Subramaniam Reference Tenneti, Garg and Subramaniam2011; Deen et al. Reference Deen, Kriebitzsch, van der Hoef and Kuipers2012; Tang et al. Reference Tang, Kriebitzsch, Peters, van der Hoef and Kuipers2014). 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 Reference Takagi, Oguz and Prosperetti2003). 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 Reference Botto and Prosperetti2012).

According to the literature, there is a strong indication that point-particle models underpredict turbulence attenuation (Hwang & Eaton Reference Hwang and Eaton2006; Eaton Reference Eaton2009; Balachandar & Eaton Reference Balachandar and Eaton2010). Hwang & Eaton (Reference Hwang and Eaton2006) concluded that the extra dissipation caused by particles with $d_{p}\approx {\it\eta}$ was greatly underestimated by conventional models. Furthermore, Burton & Eaton (Reference Burton and Eaton2005) 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. Reference Lucci, Ferrante and Elghobashi2010). 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 (Reference Eaton2009) 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 (Reference Burton and Eaton2005) 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 Reference Starius1977), a chimera grid-embedding technique (Benek, Buning & Steger Reference Benek, Buning and Steger1985) and a composite overlapping grid technique (Chesshire & Henshaw Reference Chesshire and Henshaw1990; Henshaw Reference Henshaw1994). 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 Reference Bagchi and Balachandar2003; Burton & Eaton Reference Burton and Eaton2005) and, on the other, immersed boundary PR-DNS of turbulent flow around many freely moving coarser particles (Uhlmann Reference Uhlmann2008; Lucci et al. Reference Lucci, Ferrante and Elghobashi2010; Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015; Picano et al. Reference Picano, Breugem and Brandt2015). Burton & Eaton (Reference Burton and Eaton2005) 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 (Reference Tanaka and Eaton2010) 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 (Reference Burton and Eaton2002, Reference Burton and Eaton2005) 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 Reference Phan-Thien2013). 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 (Reference Bagchi and Balachandar2002), 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 (Reference Bagchi and Balachandar2002). In immersed boundary simulations of this flow in a smaller computational domain, $C_{D}=1.178$ (Tang et al. Reference Tang, Kriebitzsch, Peters, van der Hoef and Kuipers2014) and $C_{D}=1.175$ (Baltussen Reference Baltussen2015) 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 (Reference Zeng, Balachandar and Fischer2005), 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. (Reference Zeng, Balachandar and Fischer2005), $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. (Reference Tenneti, Garg and Subramaniam2011), 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. (Reference Tenneti, Garg and Subramaniam2011). 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. (Reference Mehrabadi, Tenneti, Garg and Subramaniam2015)). 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 (Reference Uhlmann2005), Mark & van Wachem (Reference Mark and van Wachem2008) and Breugem (Reference Breugem2012). In the test presented by Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011), 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 (Reference Hwang and Eaton2006) and Tanaka & Eaton (Reference Tanaka and Eaton2010) 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 Reference Djenidi and Antonia2015).

The procedure used for the stochastic forcing is similar to the one used by Yeung & Pope (Reference Yeung and Pope1989) 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 (Reference Yeung and Pope1989), 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 (Reference Eswaran and Pope1988), 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 (Reference Yeung and Pope1989) 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. (Reference Jimenez, Wray, Saffman and Rogallo1993) 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 Reference Yeung and Pope1989; Jimenez et al. Reference Jimenez, Wray, Saffman and Rogallo1993; Ishihara et al. Reference Ishihara, Kaneda, Yokokawa, Itakura and Uno2007). 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 Reference Yeung and Pope1989; Jimenez et al. Reference Jimenez, Wray, Saffman and Rogallo1993; Ishihara et al. Reference Ishihara, Kaneda, Yokokawa, Itakura and Uno2007). 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 Reference Vreman and Kuerten2014b ). 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 Reference Vreman and Kuerten2014a ). 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.14.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 Reference Burton and Eaton2005; Tanaka & Eaton Reference Tanaka and Eaton2010). 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 (Reference Mansour, Kim and Moin1988) and Vreman & Kuerten (Reference Vreman and Kuerten2014b ) 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 (Reference Burton and Eaton2005), since only the part of the curve where the enhancement was less than a factor of two was shown. Tanaka & Eaton (Reference Tanaka and Eaton2010) 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. Reference Lucci, Ferrante and Elghobashi2010). 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. (Reference Lucci, Ferrante and Elghobashi2010) 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 (Reference Hwang and Eaton2006) and Tanaka & Eaton (Reference Tanaka and Eaton2010) 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 35 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 Reference Hwang and Eaton2006; Tanaka & Eaton Reference Tanaka and Eaton2010) 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 Reference Hwang and Eaton2006), glass beads of diameter 165  ${\rm\mu}$ m were used in air, while in the third (fourth) experiment (Tanaka & Eaton Reference Tanaka and Eaton2010), 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. Reference Kulick, Fessler and Eaton1994; Tanaka & Eaton Reference Tanaka and Eaton2010). 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 (Reference Tanaka and Eaton2010) ( $Re_{{\it\lambda}}=127$ ) was slightly less than in the experiments of Hwang & Eaton (Reference Hwang and Eaton2006) ( $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 (Reference Tanaka and Eaton2010) measured a much smaller reduction of $\overline{{\it\epsilon}}$ than Hwang & Eaton (Reference Hwang and Eaton2006) (in the fourth experiment listed in table 6, the turbulence dissipation rate did not reduce at all). Tanaka & Eaton (Reference Tanaka and Eaton2010) 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 (Reference Tanaka and Eaton2010) 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 (Reference Tanaka and Eaton2010) 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. Reference Mansour, Kim and Moin1988). 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. Reference Mansour, Kim and Moin1988).

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. Reference Ishihara, Kaneda, Yokokawa, Itakura and Uno2007, 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. Reference Jimenez, Wray, Saffman and Rogallo1993); 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 Reference Vedula and Yeung1999). 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. Reference Ishihara, Kaneda, Yokokawa, Itakura and Uno2007). 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 Reference Vreman and Kuerten2014b ).

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. (Reference Link, Cuypers, Deen and Kuipers2005) 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. (bd) 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 35 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 (Reference Burton and Eaton2005) 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. Reference Hill, Koch and Ladd2001). 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 (Reference Hwang and Eaton2006) 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 (Reference Tanaka and Eaton2010) measured significantly higher turbulence dissipation rates than Hwang & Eaton (Reference Hwang and Eaton2006) did. If the turbulence dissipation rate measured by Hwang & Eaton (Reference Hwang and Eaton2006) 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 (Reference Bagchi and Balachandar2003) 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 Reference Bagchi and Balachandar2003). 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 (Reference Bagchi and Balachandar2003), 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 (Reference Burton and Eaton2005), 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 (Reference Bagchi and Balachandar2003) also showed that the errors in the Schiller–Naumann correlation increase with increasing $d_{p}/{\it\eta}$ . Botto & Prosperetti (Reference Botto and Prosperetti2012) 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. (Reference Kulick, Fessler and Eaton1994) 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 Reference Vreman2015). 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.

Acknowledgement

The author is grateful to J. G. M. Kuerten and to the reviewers for useful comments on the manuscript.

Appendix A. Further specification of the numerical method

The governing equations and the overlapping grid structure were defined in §§ 2.2 and 2.3. When the simulation is started, the metric terms are determined and stored for each staggered location in each mesh. In addition, the interpolation data structure that arranges the communication among the different meshes is set up. The grid points for which interpolation is required are identified, and all interpolation coefficients are computed and stored. Then the initial velocity field is set, and the boundary conditions and interpolations are applied to the initial velocity field. The time level $n$ is set to zero, and the time step is set to ${\rm\Delta}t$ . The projection method is then used to solve the velocity and pressure at the next time level $(n+1)$ . It consists of two steps: the intermediate velocity step and the projection step. In the following, the intermediate velocity step is described first, then the projection step is described, and finally the spatial discretization and interpolation procedures are described.

A.1 Intermediate velocity step

The intermediate velocity step of the projection method consists of four substeps. First the auxiliary velocities are computed according to

(A 1a,b ) $$\begin{eqnarray}\boldsymbol{u}^{n+1/2}=3\boldsymbol{u}^{n}/2-\boldsymbol{u}^{n-1}/2,\quad \tilde{\boldsymbol{u}}^{n+1/2}=3\tilde{\boldsymbol{u}}^{n}/2-\tilde{\boldsymbol{u}}^{n-1}/2,\end{eqnarray}$$

except in the first time step ( $n=0$ ), when these auxiliary velocities are set equal to the initial condition.

The second substep is the computation of the intermediate velocity $\boldsymbol{u}^{\ast }$ on the Cartesian grid:

(A 2) $$\begin{eqnarray}\boldsymbol{u}^{\ast }=\boldsymbol{u}^{n}+{\rm\Delta}t(-\boldsymbol{H}^{n+1/2}+\boldsymbol{J}^{n}),\end{eqnarray}$$

where $\boldsymbol{H}^{n+1/2}$ represents the nonlinear terms applied to the auxiliary velocity $\boldsymbol{u}^{n+1/2}$ and $\boldsymbol{J}^{n}$ represents the viscous terms applied to the velocity $\boldsymbol{u}^{n}$ plus the forcing terms at level $n$ .

The third substep is the computation of the intermediate velocities on the spherical grids. Two terms in the equation for $\tilde{u} _{j}$ are treated implicitly: the first-order convective derivative with respect to ${\it\phi}$ and the second-order viscous derivative with respect to ${\it\phi}$ . All other terms are treated explicitly. On each spherical grid the following system is solved for the intermediate velocity:

(A 3) $$\begin{eqnarray}(\unicode[STIX]{x1D644}-{\textstyle \frac{1}{2}}{\rm\Delta}t(-\tilde{\unicode[STIX]{x1D643}}^{{\it\phi}}+\tilde{\unicode[STIX]{x1D645}}^{{\it\phi}}))\tilde{\boldsymbol{u}}^{\ast }=\tilde{\boldsymbol{u}}^{n}+{\textstyle \frac{1}{2}}{\rm\Delta}t(-\tilde{H}^{{\it\phi}}+\tilde{J}^{{\it\phi}})\tilde{\boldsymbol{u}}^{n}+{\rm\Delta}t(\tilde{\boldsymbol{H}}^{n+1/2}+\tilde{\boldsymbol{J}}^{n}),\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the identity matrix, and $\tilde{\unicode[STIX]{x1D643}}^{{\it\phi}}$ and $\tilde{\unicode[STIX]{x1D645}}^{{\it\phi}}$ are the coefficient matrices for the implicit convective term and implicit viscous term, respectively. The matrix $\tilde{\unicode[STIX]{x1D645}}^{{\it\phi}}$ is the same in each time step, but the matrix $\tilde{\unicode[STIX]{x1D643}}^{{\it\phi}}$ needs to be computed in each time step since it depends on the auxiliary velocity component $\tilde{u} _{3}^{n+1/2}=u_{{\it\phi}}^{n+1/2}$ . The vector $\tilde{\boldsymbol{H}}^{n+1/2}$ represents the explicit nonlinear terms (applied to the auxiliary velocity $\tilde{\boldsymbol{u}}^{n+1/2}$ ), and $\tilde{\boldsymbol{J}}^{n}$ represents all explicit viscous terms (applied to the velocity $\tilde{\boldsymbol{u}}^{n}$ ) plus the forcing terms at level $n$ . The forcing terms are given by $\unicode[STIX]{x1D63C}\boldsymbol{f}$ after interpolation of the Cartesian $\boldsymbol{f}$ to the spherical grid. For moving particles, the particle velocity term that appears in the Navier–Stokes equation in the spherical frame of reference should also be included into $\tilde{\boldsymbol{J}}^{n}$ . The implicit systems for the spherical intermediate velocities are solved by a direct method, Gauss elimination, which is fast since only tri-diagonal matrices need to be inverted.

The last substep is that the same boundary condition and interpolation routines used for the velocity are executed for the intermediate velocity fields.

A.2 Projection step

The projection step of the projection method consists of three substeps. In the first substep, the discrete pressure field (at level $n+1$ ) is computed from the following system of equations (no summation convention over index $i$ ):

(A 4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}cl@{}}\displaystyle a_{i}({\rm\nabla}^{2}p)_{i}+b=a_{i}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\ast })_{i}/{\rm\Delta}t\quad & \text{at all internal points }i\\ & \text{of the Cartesian domain},\\ \displaystyle a_{i}(\tilde{{\rm\nabla}}^{2}p)_{i}+b=a_{i}(\tilde{\boldsymbol{{\rm\nabla}}}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}^{\ast })_{i}/{\rm\Delta}t\quad & \text{at all internal points }i\\ & \text{of all spherical domains,}\\ \qquad \qquad \qquad \displaystyle \mathop{\sum }_{i=1}^{m}p_{i}=0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The system contains $m+1$ equations and $m+1$ unknowns, where $m$ is the total number of internal points in all domains. The first two equations represent the discretized Poisson equations, where $(\cdot )_{i}$ denotes the discretization of the operator inside the brackets at point $i$ . The equations are normalized such that the diagonal elements of the coefficient matrix are equal to 1; $a_{i}$ is the reciprocal of the coefficient of $p_{i}$ in the stencil of the discrete Poisson operator. This normalization is effectively Jacobi preconditioning. The augmented matrix method proposed by Henshaw (Reference Henshaw1994) leads to the appearance of the last equation and the additional unknown $b$ . This technique is used because, due to the interpolations, the global conservation property of the continuous Laplace and divergence operators is not exactly satisfied by the discrete system, although the property is satisfied in the limit of zero grid size ( $b$ converges to zero in the limit of zero grid size). Without this technique or another treatment, the system would be singular and the iterative procedure would not converge. The complication is inherent to pressure Poisson problems without a Dirichlet boundary condition.

The total linear system to be solved is written as $\unicode[STIX]{x1D648}\boldsymbol{y}=\boldsymbol{z}$ , where $\unicode[STIX]{x1D648}$ is the $(m+1)\times (m+1)$ coefficient matrix, $\boldsymbol{y}$ the vector of unknowns and $\boldsymbol{z}$ the vector of the $m+1$ right-hand sides of the system. The linear system $\unicode[STIX]{x1D648}\boldsymbol{y}=\boldsymbol{z}$ is iteratively solved by the BiCGStab(1) method proposed by van der Vorst (Reference van der Vorst1992). The implementation of the method is based on a Fortran subroutine written by M. A. Botchev, who used the code written by Fokkema (Reference Fokkema1996) and incorporated the enhancements proposed by Sleijpen & van der Vorst (Reference Sleijpen and van der Vorst1995, Reference Sleijpen and van der Vorst1996). The subroutine was parallelized by the author of the present work. When the iterative procedure is started, the vector $\boldsymbol{y}$ contains the pressure at the previous time level and $b$ is set to zero. The iterations continue until the maximum of the absolute elements of the residual vector $\boldsymbol{z}-\unicode[STIX]{x1D648}\boldsymbol{y}$ is less than the prescribed tolerance, $10^{-4}$ . For each interior cell that is not adjacent to a virtual cell or a particle surface, the matrix row contains eight non-zero elements (the standard seven elements plus one that represents the coefficient of $b$ ). For interior cells adjacent to virtual cells, the matrix row contains many more non-zero elements, since the virtual points that occur in the discrete definitions of ${\rm\nabla}^{2}p$ and $\tilde{{\rm\nabla}}^{2}p$ are replaced by the interpolation relations that relate the pressure at the virtual points to the pressure at interior points. For interior cells adjacent to the particle surface, the matrix row contains one element less since in these cells $\tilde{{\rm\nabla}}^{2}p$ does not represent the physical Laplace operator, but a discrete pseudo-Laplace operator, after $\partial p/\partial r$ on the wall has been replaced by zero (see Vreman Reference Vreman2014, and references therein), which makes the staggered projection method equivalent to the approach originally proposed by Harlow & Welch (Reference Harlow and Welch1965). Thus the method does not need the pressure gradient on the wall. Indeed, it can be shown that the physical boundary condition respected by the staggered discretization is that the Laplacian of the velocity divergence is zero in the near-wall limit, which is equivalent to imposing the continuity equation at the solid boundary (Vreman Reference Vreman2014).

After the pressure field has been obtained, the second substep is applied, which is the actual projection step. The pressure at the virtual points is set using the same interpolation relations that have been used in each matrix evaluation (iteration) of the BiCGStab method. The velocities at the internal points are obtained by projecting the intermediate velocities on the space of functions that is (approximately) divergence-free,

(A 5) $$\begin{eqnarray}\boldsymbol{u}^{n+1}=\boldsymbol{u}^{\ast }-{\rm\Delta}t\boldsymbol{{\rm\nabla}}p\end{eqnarray}$$

in the Cartesian domain, and

(A 6) $$\begin{eqnarray}\tilde{\boldsymbol{u}}^{n+1}=\tilde{\boldsymbol{u}}^{\ast }-{\rm\Delta}t\tilde{\boldsymbol{{\rm\nabla}}}p\end{eqnarray}$$

in all spherical domains, where $\tilde{\boldsymbol{{\rm\nabla}}}$ represents the gradient operator applied to a scalar field in spherical coordinates.

The last substep is the execution of the boundary condition and interpolation routines for the velocity fields at the new time level $n+1$ . Whereas the discretization of the convective terms and the implicitly treated viscous terms is second order in time, the discretization of the pressure and the explicitly treated viscous terms is only first order in time. However, the simulations did not suffer from excessive temporal discretization errors, as has been confirmed by grid refinement, which included a reduction of the time step (see §§ 3 and 4). Since the intermediate velocity step is explicit (apart from the ${\it\phi}$ direction), a refinement of the mesh by a factor  $2$ requires the time step to be smaller by a factor  $4$ . Thus in combination with second-order spatial accuracy, the relative reduction of temporal and spatial errors are the same. To alleviate the time step restriction, one could consider using the Crank–Nicolson scheme for all viscous terms. However, a fast implicit treatment of all viscous terms in the spherical coordinate system is not straightforward. In addition, the computational effort per time step would increase.

A.3 Spatial discretization and interpolations

All spatial partial derivatives in the Navier–Stokes equations in the forms specified in § 2.2 are discretized by second-order central differences. In the computation of the nonlinear fluxes, the velocities are averaged to the faces (using weights $1/2$ and $1/2$ ) before they are multiplied. Each metric quantity in the spherical form is computed by direct substitution of the $r$ and ${\it\theta}$ coordinate of the location where the metric quantity is defined (centre or face of a pressure or velocity cell). This means that no averaging of any metric quantity is used (averaging of metric quantities may lead to inconsistent discretization near the poles).

The grid of $u_{{\it\theta}}$ has points at the poles ( ${\it\theta}=0$ or ${\it\theta}={\rm\pi}$ ), where the $u_{{\it\theta}}$ momentum equation is not defined. To approximate $u_{{\it\theta}}$ at the poles, fourth-order interpolation across the poles is used,

(A 7) $$\begin{eqnarray}\displaystyle u_{{\it\theta}}(r,0,{\it\phi},t) & = & \displaystyle {\textstyle \frac{2}{3}}[u_{{\it\theta}}(r,{\rm\Delta}{\it\theta},{\it\phi},t)-u_{{\it\theta}}(r,{\rm\Delta}{\it\theta},{\it\phi}+{\rm\pi},t)]\nonumber\\ \displaystyle & & \displaystyle -\,{\textstyle \frac{1}{6}}[u_{{\it\theta}}(r,2{\rm\Delta}{\it\theta},{\it\phi},t)-u_{{\it\theta}}(r,2{\rm\Delta}{\it\theta},{\it\phi}+{\rm\pi},t)],\end{eqnarray}$$

and a similar equation is used for ${\it\theta}={\rm\pi}$ . The minus sign of the terms inside the rectangular brackets occurs because at the poles $u_{{\it\theta}}(r,{\it\theta},{\it\phi}+{\rm\pi},t)=-u_{{\it\theta}}(r,{\it\theta},{\it\phi},t)$ should hold.

An essential feature of the discretization is the interpolation routines applied to the variables in the virtual cells. These interpolations are third-order interpolations. For example, the pressure in virtual cell $i$ of a spherical grid is defined as an average over $3\times 3\times 3$ neighbouring points of the Cartesian grid,

(A 8) $$\begin{eqnarray}p_{i}=\mathop{\sum }_{j\in W_{i}}w_{ij}p_{j},\end{eqnarray}$$

where $W_{i}$ is the set of Cartesian grid points $j$ on the Cartesian grid and the $27$ weight coefficients $w_{ij}$ are based on quadratic Lagrange polynomials for each Cartesian direction. The interpolation of a velocity component requires more effort. For example, to obtain $u_{r}$ in virtual cell $i$ of a spherical grid, all three Cartesian velocity components are needed. Each of the three Cartesian velocity components is interpolated to virtual cell $i$ , using expressions similar to (A 8), then the transformation to $u_{r}$ is performed, according to the first expression in (2.5). This means that $3\times 27$ Cartesian grid points and as many weight coefficients are used for each velocity component in a virtual cell. Burton & Eaton (Reference Burton and Eaton2002) also used third-order interpolations, but they interpolated the pressure gradient instead of the pressure. Interpolation of the pressure was used by Chesshire & Henshaw (Reference Chesshire and Henshaw1990) and Henshaw (Reference Henshaw1994), who developed discretization methods for collocated overlapping grids.

The tangential boundary conditions on the surfaces of the particles are discretized by linear extrapolation of $u_{{\it\theta}}$ and $u_{{\it\phi}}$ at radius $r_{0}^{c}$ from the values at the radii $r_{0}$ and $r_{1}^{c}$ . The approximation of the radial $u_{{\it\theta},r}$ at the wall is approximated by $(u_{{\it\theta}}(r_{1}^{c},{\it\theta},{\it\phi},t)-u_{{\it\theta}}(r_{0}^{c},{\it\theta},{\it\phi},t))/(r_{1}^{c}-r_{0}^{c})$ , and similarly $u_{{\it\phi},r}$ is approximated. As indicated, the pressure at location $r_{0}^{c}$ is not required in the discretized partial differential equations. However, the pressure at $r_{0}^{s}=r_{0}$ is required, to compute the force $F_{i}$ , exerted on the particle surface. This pressure is approximated by linear extrapolation: $(3/2)p(r_{1}^{c},{\it\theta},{\it\phi},t)-(1/2)p(r_{2}^{c},{\it\theta},{\it\phi},t)$ .

The radial statistical profiles are defined at the cell centres (radial locations $r_{j}^{c}$ , defined in § 2.3). In order to obtain the radial statistics based on velocity components (not on velocity derivatives), the spherical velocity components are first transferred to pressure locations, with use of second-order interpolation for $u_{{\it\theta}}$ and $u_{{\it\phi}}$ and third-order interpolation for $u_{r}$ .

The statistics of the small-scale turbulence are expressed in first-order velocity derivatives. The Cartesian finite difference formula in the post-processing routine is, like the first-order derivatives inside the Laplacian, based on finite differences applied to a small stencil. For example, $u_{2,1}$ is defined at the $x_{1}$ faces of the staggered $u_{2}$ -cells, using a central difference based on two points with distance $h_{1}$ . For the spherical grid, it is more convenient to compute the nine components of the gradient of the velocity, $\unicode[STIX]{x1D642}$ , at the cell centres of the $p$ -cells. This means that, for example, the numerical expression for $u_{r,r}$ is based on two points whose coordinates differ by ${\rm\Delta}r$ , while, for example, the expression for $u_{r,{\it\theta}}$ is based on two points whose coordinates differ by $2{\rm\Delta}{\it\theta}$ .

References

Bagchi, P. & Balachandar, S. 2002 Steady planar straining flow past a rigid sphere at moderate Reynolds number. J. Fluid Mech. 466, 365407.CrossRefGoogle Scholar
Bagchi, P. & Balachandar, S. 2003 Effect of turbulence on the drag and lift of a particle. Phys. Fluids 15, 34963513.CrossRefGoogle Scholar
Balachandar, S. & Eaton, J. K. 2010 Dispersed turbulent multiphase flow. Annu. Rev. Fluid Mech. 42, 111133.CrossRefGoogle Scholar
Baltussen, M. W.2015 Bubbles on the cutting edge – direct numerical simulations of gas–liquid–solid three-phase flows. PhD thesis, Eindhoven University of Technology.Google Scholar
Benek, J. A., Buning, P. G. & Steger, J. L.1985 A 3-D chimera grid embedding technique, AIAA Paper 85-1523.CrossRefGoogle Scholar
Botto, L. & Prosperetti, A. 2012 A fully resolved numerical simulation of turbulent flow past one or several spherical particles. Phys. Fluids 24, 013303.CrossRefGoogle Scholar
Breugem, W. P. 2012 A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows. J. Comput. Phys. 231, 44694498.CrossRefGoogle Scholar
Burton, T. M. & Eaton, J. K. 2002 Analysis of a fractional-step method on overset grids. J. Comput. Phys. 177, 336364.CrossRefGoogle Scholar
Burton, T. M. & Eaton, J. K. 2005 Fully resolved simulations of particle–turbulence interaction. J. Fluid Mech. 545, 67111.CrossRefGoogle Scholar
ten Cate, A., Derksen, J. J., Portela, L. M. & van den Akker, H. E. A. 2004 Fully resolved simulations of colliding monodisperse spheres in forced isotropic turbulence. J. Fluid Mech. 539, 233271.CrossRefGoogle Scholar
Chesshire, G. & Henshaw, W. D. 1990 Composite overlapping meshes for the solution of partial differential equations. J. Comput. Phys. 90, 164.CrossRefGoogle Scholar
Deen, N. G., Kriebitzsch, S. H. L., van der Hoef, M. A. & Kuipers, J. A. M. 2012 Direct numerical simulation of flow and heat transfer in dense fluid–particle systems. Chem. Engng Sci. 81, 329344.CrossRefGoogle Scholar
Djenidi, L. & Antonia, R. A. 2015 A general self-preservation analysis for decaying homogeneous isotropic turbulence. J. Fluid Mech. 773, 345365.CrossRefGoogle Scholar
Eaton, J. K. 2009 Two-way coupled turbulence simulations of gas–particle flows using point-particle tracking. Intl J. Multiphase Flow 35, 792800.CrossRefGoogle Scholar
Elghobashi, S. & Truesdell, G. C. 1993 On the two-way interaction between homogeneous turbulence and dispersed solid particles. I: Turbulence modification. Phys. Fluids A 5, 17901801.CrossRefGoogle Scholar
Eswaran, V. & Pope, S. B. 1988 An examination of forcing in direct numerical simulations of turbulence. Comput. Fluids 16, 257278.CrossRefGoogle Scholar
Ferrante, A. & Elghobashi, S. 2003 On the physical mechanisms of two-way coupling in particle-laden isotropic turbulence. Phys. Fluids 15, 315329.CrossRefGoogle Scholar
Fokkema, D. R.1996 Subspace methods for linear, nonlinear, and eigen problems. PhD thesis, University of Utrecht, The Netherlands.Google Scholar
Gore, R. A. & Crowe, C. T. 1989 Effect of particle size on modulating turbulence intensity. Intl J. Multiphase flow 15, 279285.CrossRefGoogle Scholar
Harlow, F. E. & Welch, J. E. 1965 Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Phys. Fluids 8, 21822189.CrossRefGoogle Scholar
Henshaw, W. D. 1994 A fourth-order accurate method for the incompressible Navier–Stokes equations on overlapping grids. J. Comput. Phys. 113, 1325.CrossRefGoogle Scholar
Hetseroni, G. 1989 Particles–turbulence interaction. Intl J. Multiphase Flow 15, 735746.CrossRefGoogle Scholar
Hill, R., Koch, D. L. & Ladd, A. J. C. 2001 Moderate-Reynolds-number flows in ordered and random arrays of spheres. J. Fluid Mech. 448, 243278.CrossRefGoogle Scholar
Hwang, W. & Eaton, J. K. 2006 Homogeneous and isotropic turbulence modulation by small heavy (St ∼ 50) particles. J. Fluid Mech. 564, 361393.CrossRefGoogle Scholar
Ishihara, T., Kaneda, Y., Yokokawa, M., Itakura, K. & Uno, A. 2007 Small-scale statistics in high-resolution direct numerical simulation of turbulence: Reynolds number dependence of one-point velocity gradient statistics. J. Fluid Mech. 592, 335366.CrossRefGoogle Scholar
Jimenez, J., Wray, A. A., Saffman, P. G. & Rogallo, R. S. 1993 The structure of intense vorticity in isotropic turbulence. J Fluid Mech. 255, 6590.CrossRefGoogle Scholar
Kulick, J. D., Fessler, J. R. & Eaton, J. K. 1994 Particle response and turbulence modification in fully developed channel flow. J. Fluid Mech. 277, 109134.CrossRefGoogle Scholar
Kussin, J. & Sommerfeld, M. 2002 Experimental studies on particle behaviour and turbulence modification in horizontal channel flow with different wall roughness. Exp. Fluids 33, 143159.CrossRefGoogle Scholar
Link, J. M., Cuypers, L. A., Deen, N. G. & Kuipers, J. A. M. 2005 Flow regimes in a spout-fluid bed: a combined experimental and simulation study. Chem. Engng Sci. 60, 34253442.CrossRefGoogle Scholar
Lucci, F., Ferrante, A. & Elghobashi, S. 2010 Modulation of isotropic turbulence by particles of Taylor length-scale size. J. Fluid Mech. 650, 555.CrossRefGoogle Scholar
Mansour, N. N., Kim, J. & Moin, P. 1988 Reynolds-stress and dissipation-rate budgets in a turbulent channel flow. J. Fluid Mech. 194, 1544.CrossRefGoogle Scholar
Mark, A. & van Wachem, B. G. M. 2008 Derivation and validation of a novel implicit second-order accurate immersed boundary method. J. Comput. Phys. 227, 66606680.CrossRefGoogle Scholar
Mehrabadi, M., Tenneti, S., Garg, R. & Subramaniam, S. 2015 Pseudo-turbulent gas-phase velocity fluctuations in homogeneous gas–solid flow: fixed particle assemblies and freely evolving suspensions. J. Fluid Mech. 770, 210246.CrossRefGoogle Scholar
Phan-Thien, N. 2013 Understanding Viscoelasticity. Springer.CrossRefGoogle Scholar
Picano, F., Breugem, W.-P. & Brandt, L. 2015 Turbulent channel flow of dense suspensions of neutrally-buoyant spheres. J. Fluid Mech. 764, 463487.CrossRefGoogle Scholar
Prosperetti, A. 2015 Life and death by boundary conditions. J. Fluid Mech. 768, 14.CrossRefGoogle Scholar
Sleijpen, G. & van der Vorst, H. 1995 Maintaining convergence properties of BiCGstab methods in finite precision arithmetic. Numer. Algor. 10, 203223.CrossRefGoogle Scholar
Sleijpen, G. & van der Vorst, H. 1996 Reliable updated residuals hybrid BiCG methods. Computing 56, 141163.CrossRefGoogle Scholar
Squires, K. D. & Eaton, J. K. 1990 Particle response and turbulence modification in isotropic turbulence. Phys. Fluids A 2, 11911203.CrossRefGoogle Scholar
Starius, G. 1977 Composite mesh difference methods for elliptic boundary value problems. Numer. Math. 28, 243258.CrossRefGoogle Scholar
Takagi, S., Oguz, H. Z. Z. & Prosperetti, A. 2003 PHYSALIS: a new method for particle simulations. Part II: two-dimensional Navier–Stokes flow around cylinders. J. Comput. Phys. 187, 371390.CrossRefGoogle Scholar
Tanaka, T. & Eaton, J. K. 2010 Sub-Kolmogorov resolution particle image velocimetry measurements of particle-laden forced turbulence. J. Fluid Mech. 643, 177206.CrossRefGoogle Scholar
Tang, Y., Kriebitzsch, S. H. L., Peters, E. J. A. F., van der Hoef, M. A. & Kuipers, J. A. M. 2014 A methodology for highly accurate results of direct numerical simulations: drag force in dense gas–solid flows at intermediate Reynolds number. Intl J. Multiphase Flow 62, 7386.CrossRefGoogle Scholar
Tenneti, S., Garg, R. & Subramaniam, S. 2011 Drag law for monodisperse gas–solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres. Intl J. Multiphase Flow 37, 10721092.CrossRefGoogle Scholar
Tryggvason, G., Scardovelli, R. & Zaleski, S. 2011 Direct Numerical Simulations of Gas–Liquid Flows. Cambridge University Press.Google Scholar
Tsuji, Y., Morikawa, Y. & Shiomi, H. 1984 LDV measurements of an air–solid two-phase flow in a vertical pipe. J. Fluid Mech. 139, 417434.CrossRefGoogle Scholar
Uhlmann, M. 2005 An immersed boundary method with direct forcing for the simulation of particulate flows. J. Comput. Phys. 209, 448476.CrossRefGoogle Scholar
Uhlmann, M. 2008 Interface-resolved direct numerical simulation of vertical particulate channel flow in the turbulent regime. Phys. Fluids 20, 053305.CrossRefGoogle Scholar
Vedula, P. & Yeung, P. K. 1999 Similarity scaling of acceleration and pressure statistics in numerical simulations of isotropic turbulence. Phys. Fluids 11, 12081220.CrossRefGoogle Scholar
van der Vorst, H. A. 1992 Bi-CGSTAB – a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear-systems. SIAM J. Sci. 13 (2), 631644.Google Scholar
Vreman, A. W. 2007 Turbulence characteristics of particle-laden pipe flow. J. Fluid Mech. 584, 235279.CrossRefGoogle Scholar
Vreman, A. W. 2014 The projection method for the incompressible Navier–Stokes equations: the pressure near a no-slip wall. J. Comput. Phys. 263, 353374.CrossRefGoogle Scholar
Vreman, A. W. 2015 Turbulence attenuation in particle-laden flow in smooth and rough channels. J. Fluid Mech. 773, 103136.CrossRefGoogle Scholar
Vreman, A. W. & Kuerten, J. G. M. 2014a Comparison of direct numerical simulation databases of turbulent channel flow at Re 𝜏 = 180. Phys. Fluids 26, 015102.CrossRefGoogle Scholar
Vreman, A. W. & Kuerten, J. G. M. 2014b Statistics of spatial derivatives of velocity and pressure in turbulent channel flow. Phys. Fluids 26, 085103.CrossRefGoogle Scholar
Yeung, P. K. & Pope, S. B. 1989 Lagrangian statistics from direct numerical simulations of isotropic turbulence. J. Fluid Mech. 207, 531586.CrossRefGoogle Scholar
Zeng, L., Balachandar, S. & Fischer, P. 2005 Wall-induced forces on a rigid sphere at finite Reynolds number. J. Fluid Mech. 536, 125.CrossRefGoogle Scholar
Figure 0

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.

Figure 1

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

Figure 2

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.

Figure 3

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

Figure 4

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

Figure 5

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

Figure 6

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

Figure 7

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

Figure 8

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 9

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

Figure 10

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

Figure 11

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

Figure 12

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.

Figure 13

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

Figure 14

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

Figure 15

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

Figure 16

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.

Figure 17

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.

Figure 18

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. (bd) 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).