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 3–5 are discussed. Afterwards, the main simulation results are discussed in the light of experimental work found in the literature (table 6).

Table 4.

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

Table 5.

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

Table 6.

Comparison between experiments of particle-laden stationary homogeneous isotropic turbulence in the literature (Hwang & Eaton 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. (b–d) Radial profiles of the turbulence dissipation rate of point-particle simulations PP0 (b), PP2 (c) and PP4 (d), for which ${\it\epsilon}$ (black solid line and open symbols) has been decomposed into a resolved part ${\it\epsilon}^{(1)}$ (red solid) and a Schiller–Naumann part ${\it\epsilon}^{(2)}$ (red dashed). The ${\it\epsilon}$ value from the point-particle simulations should be compared to ${\it\epsilon}$ from the particle-resolved simulation L128 (blue dashed).

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

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

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

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

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

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

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

For a turbulent flow around a single particle with diameter $d\approx 2{\it\eta}$ , a case in which the global turbulence kinetic energy was not modified, Burton & Eaton (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.