3.1 Validation

The code used in the present work has been already validated against several different cases in previous studies (Breugem Reference Breugem2012; Picano et al. Reference Picano, Breugem and Brandt2015; Fornari et al. Reference Fornari, Picano and Brandt2016c ; Kazerooni et al. Reference Kazerooni, Fornari, Hussong and Brandt2017). To further investigate the accuracy of the code, we calculate the friction factor $f=8(\bar{U}_{\ast }/U_{b})^{2}$ for the reference unladen case with $Re_{b}=5600$ , and compare it with the value obtained from the empirical correlation by Jones (Reference Jones1976)

(3.1) $$\begin{eqnarray}\displaystyle 1/f^{2}=2\log _{10}(1.125Re_{b}f^{1/2})-0.8. & & \displaystyle\end{eqnarray}$$

The same value of $f=0.035$ is obtained from the simulation and the empirical formula. This corresponds to a mean $Re_{\unicode[STIX]{x1D70F}}=185$ .

We also performed a simulation at lower $Re_{b}=4410$ and $\unicode[STIX]{x1D719}=0$ , see figure 2, where we report the profile of the streamwise velocity fluctuation at the wall bisector, normalized by the local friction velocity $U_{\ast }$ . This is compared to the results by Gavrilakis (Reference Gavrilakis1992) and Joung, Choi & Choi (Reference Joung, Choi and Choi2007) at $Re_{b}=4410$ and 4440. We see a good agreement with both works. For all components of the root-mean-square velocity, the relative difference between our results (at the wall bisector) and those by Gavrilakis (Reference Gavrilakis1992) is typically less than 1.5 % (and at most 4 % locally).

Figure 2.

Streamwise velocity fluctuation at the wall bisector from the present simulation at $Re_{b}=4410$ , and from the data by Gavrilakis (Reference Gavrilakis1992) at $Re_{b}=4410$ and Joung et al. (Reference Joung, Choi and Choi2007) at $Re_{b}=4440$ .

3.2 Mean velocities, drag and particle concentration

In this section we report and discuss the results obtained for the different solid volume fractions $\unicode[STIX]{x1D719}$ considered. The phase ensemble averages for the fluid (solid) phase have been calculated by considering only the points located outside (inside) of the volume occupied by the particles. The statistics reported are obtained by further averaging over the eight symmetric triangles that form the duct cross-section.

The streamwise mean fluid and particle velocities in outer units (i.e. normalized by the bulk velocity $U_{b}$ ), $U_{f/p}(y,z)$ , are illustrated in figure 3(a,b) for all $\unicode[STIX]{x1D719}$ . The contour plots are divided in four quadrants showing results for $\unicode[STIX]{x1D719}=0.0$ (top left), 0.05 (top right), 0.1 (bottom right) and $0.2$ (bottom left). The streamwise mean particle velocity contours closely resemble those of the fluid phase. In particular we observe that the maximum velocity at the core of the duct grows with $\unicode[STIX]{x1D719}$ . The increase with $\unicode[STIX]{x1D719}$ is similar to that reported for turbulent channel flows (Picano et al. Reference Picano, Breugem and Brandt2015), except for $\unicode[STIX]{x1D719}=0.2$ where the increase of $U_{f/p}(y,z)$ in the duct core is substantially larger. We observe that the convexity of the mean velocity contours also increases with the volume fraction up to $\unicode[STIX]{x1D719}=0.1$ . This is due to the increased intensity of secondary flows that convect mean velocity from regions of large shear along the walls towards regions of low shear along the corner bisectors (Prandtl Reference Prandtl1963; Gessner Reference Gessner1973; Vinuesa et al. Reference Vinuesa, Noorani, Lozano-Durán, Khoury, Schlatter, Fischer and Nagib2014). For $\unicode[STIX]{x1D719}=0.2$ the secondary flow intensity is substantially reduced and accordingly also the convexity of the contours of $U_{f/p}$ reduces.

Figure 3.

Contours of streamwise mean fluid (a) and particle velocity (b) in outer units. In each figure, the top left, top right, bottom right and bottom left quadrants show the data for $\unicode[STIX]{x1D719}=0.0,0.05,0.1$ and 0.2.

Next, we show in figure 4(a,b) the streamwise mean fluid and particle velocity profiles along the wall bisector in outer and inner units ( $U_{f/p}(y)$ and $U_{f/p}^{+}$ ). The local value of the friction velocity at the wall bisector, $U_{\ast }=\sqrt{\unicode[STIX]{x1D70F}_{w,bis}/\unicode[STIX]{x1D70C}_{f}}$ , is used to normalize  $U_{f/p}(y)$ . Solid lines are used for $U_{f}(y)$ while symbols are used for $U_{p}(y)$ . We observe that the mean velocity profiles of the two phases are almost perfectly overlapping at equal $\unicode[STIX]{x1D719}$ , except very close to the walls ( $y^{+}\leqslant 30$ ) where particles have a relative tangential motion (slip velocity). Note also that the mean particle velocity decreases with $\unicode[STIX]{x1D719}$ very close to the walls. We also observe that by increasing $\unicode[STIX]{x1D719}$ , the profiles of $U_{f/p}(y)$ tend towards the laminar parabolic profile with lower velocity near the wall and larger velocity at the centreline, $y/h=1$ . Concerning $U_{f/p}^{+}(y^{+})$ , we observe a progressive downward shift of the profiles with the volume fraction $\unicode[STIX]{x1D719}$ denoting a drag increase, at least up to $\unicode[STIX]{x1D719}=0.1$ .

The mean velocity profiles still follow the log law (Pope Reference Pope2000)

(3.2) $$\begin{eqnarray}\displaystyle U_{f/p}^{+}(y^{+})=\frac{1}{k}\log (y^{+})+B, & & \displaystyle\end{eqnarray}$$

where $k$ is the von Kármán constant and $B$ is an additive coefficient. For the unladen case with $Re_{b}=4410$ , Gavrilakis (Reference Gavrilakis1992) fitted the data between $y^{+}=30$ and 100 to find $k=0.31$ and $B=3.9$ . In the present simulations, the extent of the log region is larger due to the higher bulk Reynolds number and we hence fit our data between $y^{+}=30$ and 140. The results are reported in table 2. For the unladen case we obtain results in agreement with those of Gavrilakis (Reference Gavrilakis1992) although our constant $B$ is slightly smaller. Increasing $\unicode[STIX]{x1D719}$ , $k$ decreases and the additive constant $B$ decreases becoming negative at $\unicode[STIX]{x1D719}=0.1$ . Values for $k$ and $B$ are also reported for $\unicode[STIX]{x1D719}=0.2$ in table 2, although a log layer cannot be clearly identified.

Figure 4.

Streamwise mean fluid and particle velocity along the wall bisector at $z/h=1$ in outer (a) and inner units (b). Lines are used for fluid velocity profiles while symbols are used for particle velocity profiles.

Figure 5.

(a) Mean friction Reynolds number $Re_{\unicode[STIX]{x1D70F},mean}$ estimated from the pressure gradient, friction Reynolds number at the wall bisector $Re_{\unicode[STIX]{x1D70F},bis}$ and mean friction Reynolds number for the channel flow, $Re_{\unicode[STIX]{x1D70F},2D}$ (Picano et al. Reference Picano, Breugem and Brandt2015), for all volume fractions $\unicode[STIX]{x1D719}$ . (b) Profile of $Re_{\unicode[STIX]{x1D70F}}$ along the duct wall.

Table 2.

The von Kármán constant and additive constant $B$ of the log law at the wall bisector estimated from the present simulations for different volume fractions $\unicode[STIX]{x1D719}$ . Here $k$ and $B$ have been fitted in the range $y^{+}\in [30,140]$ . The friction Reynolds number calculated at the wall bisector $Re_{\unicode[STIX]{x1D70F},bis}$ , the mean friction Reynolds number $Re_{\unicode[STIX]{x1D70F},mean}$ estimated via the pressure gradient, and the corresponding friction Reynolds number found by Picano et al. (Reference Picano, Breugem and Brandt2015) for turbulent channel flow, $Re_{\unicode[STIX]{x1D70F},2D}$ , are also reported.

At the wall bisector, we also calculated the local friction Reynolds $Re_{\unicode[STIX]{x1D70F},bis}$ . This is reported in figure 5(a) and in table 2 as function of the volume fraction $\unicode[STIX]{x1D719}$ . The results of Picano et al. (Reference Picano, Breugem and Brandt2015) are also reported for comparison in table 2. We see that although the initial value of $Re_{\unicode[STIX]{x1D70F},bis}$ for $\unicode[STIX]{x1D719}=0.0$ is substantially larger than that of the corresponding channel flow at $Re_{b}=5600$ , the increase with the volume fraction is smaller than for $Re_{\unicode[STIX]{x1D70F},2D}$ . Indeed for $\unicode[STIX]{x1D719}=0.2$ we find that $Re_{\unicode[STIX]{x1D70F},bis}\simeq Re_{\unicode[STIX]{x1D70F},2D}$ . In the absence of particles, the increase of $Re_{\unicode[STIX]{x1D70F},bis}$ with respect to $Re_{\unicode[STIX]{x1D70F},2D}$ is due to the fact that here the secondary flows are directed towards the core, transporting streamwise momentum in the wall-normal direction, therefore inducing a steeper gradient, $\text{d}U_{f}/\text{d}y$ .

It is also interesting to observe the behaviour of the mean friction Reynolds number $Re_{\unicode[STIX]{x1D70F},mean}$ as a function of $\unicode[STIX]{x1D719}$ , as this directly relates to the overall pressure drop along the duct. For each case, we estimate the wall shear stress directly via the pressure gradient needed to drive the flow ( $\langle \unicode[STIX]{x1D70F}_{w}\rangle =-(\text{d}P/\text{d}x)L_{z}/4$ ). From the shear stress we then compute the mean friction Reynolds number. From figure 5(a) and table 2 we see that $Re_{\unicode[STIX]{x1D70F},mean}$ strongly increases with the volume fraction up to $\unicode[STIX]{x1D719}=0.1$ . The increase in $Re_{\unicode[STIX]{x1D70F},mean}$ with $\unicode[STIX]{x1D719}$ is similar to that observed for $Re_{\unicode[STIX]{x1D70F},2D}$ in channel flow. However, by further increasing the volume fraction to $\unicode[STIX]{x1D719}=0.2$ , $Re_{\unicode[STIX]{x1D70F},mean}$ remains approximately constant in duct flow (it actually increases by 0.5 % as $\unicode[STIX]{x1D719}$ increases from 0.1 to 0.2), while it increases by ${\sim}6\,\%$ in channel flow. It is also interesting to note that $Re_{\unicode[STIX]{x1D70F},mean}$ would be underestimated if only computed via the mean fluid streamwise velocity gradient at the walls. Indeed, the wall friction is altered by the shear exerted by the particles on the walls (through lubrication and collisions). This is particularly important for $\unicode[STIX]{x1D719}=0.2$ . At this $\unicode[STIX]{x1D719}$ , $Re_{\unicode[STIX]{x1D70F},mean}$ would be underestimated by approximately 3 %. In figure 5(b) we report the profiles of $Re_{\unicode[STIX]{x1D70F}}$ , displaying the signature of the near-wall structures (estimated by the velocity gradient) along one wall. We note that the friction Reynolds number increases with $\unicode[STIX]{x1D719}$ , especially towards the corners. For $\unicode[STIX]{x1D719}=0.2$ , instead, the profile exhibits a sharp change at approximately $z/h=0.1\sim 2a$ , and the maxima move towards the wall bisector ( $z/h\sim 0.65$ ). This is probably due to the clear change in local particle concentration in the cross-section as we will explain in the following.

Figure 6.

Mean particle concentration $\unicode[STIX]{x1D6F7}(y,z)$ in the duct cross-section for $\unicode[STIX]{x1D719}=0.05$ (a), $\unicode[STIX]{x1D719}=0.1$ (b) and $\unicode[STIX]{x1D719}=0.2$ (c).

The mean particle concentration over the duct cross-section $\unicode[STIX]{x1D6F7}(y,z)$ is displayed in figure 6 for all $\unicode[STIX]{x1D719}$ , whereas the particle concentration along the wall bisector ( $z/h=1$ ) and along a segment at $z/h=0.2$ is shown in figure 7. Finally, we report in figure 8 the secondary (cross-stream) velocities of both phases, defined as $\sqrt{V_{f/p}^{2}+W_{f/p}^{2}}$ . We shall now discuss these three figures together.

The particle concentration distribution is defined as

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}(y,z)=\frac{1}{N_{t}N_{x}}\mathop{\sum }_{m=1}^{N_{t}}\mathop{\sum }_{i=1}^{N_{x}}\unicode[STIX]{x1D713}(x_{ijk},t^{m}), & & \displaystyle\end{eqnarray}$$

where $N_{t}$ is the number of time steps considered for the average, $t^{m}$ is the sampling time and $\unicode[STIX]{x1D713}(x_{ijk},t^{m})$ is the particle indicator function at the location $x_{ijk}$ and time $t^{m}$ . The particle indicator function is equal to 1 for points $x_{ijk}$ contained within the volume of a sphere, and 0 otherwise. Two interesting observations are deduced from figure 6: (i) particle layers form close to the walls, and (ii) for $\unicode[STIX]{x1D719}=0.05$ and 0.1, the local particle concentration $\unicode[STIX]{x1D6F7}(y,z)$ is higher close to the duct corners. We have recently reported a similar result for laminar duct flow at $Re_{b}=550$ and the same duct-to-particle size ratio, $h/a=18$ , and volume fractions, $\unicode[STIX]{x1D719}$ (Kazerooni et al. Reference Kazerooni, Fornari, Hussong and Brandt2017). At those $Re_{b}$ and $h/a$ , particles undergo an inertial migration towards the walls and especially towards the corners, while the duct core is fully depleted of particles. Clearly, turbulence enhances mixing and the depletion of particles at the duct core disappears. It is also interesting to observe that the presence of particles further enhances the fluid secondary flow around the corners for $\unicode[STIX]{x1D719}\leqslant 0.1$ . This can be easily seen from figure 9, where both the maximum and the mean value of the secondary fluid velocity are shown as a function of the volume fraction $\unicode[STIX]{x1D719}$ . The maximum value of the secondary cross-stream velocity increases from approximately 2 to 2.5 % of $U_{b}$ . The relative increase of the mean $\sqrt{V_{f}^{2}+W_{f}^{2}}$ in comparison to the unladen case, is even larger than the increase in the maximum value at equal $\unicode[STIX]{x1D719}$ . Similarly, in laminar ducts, as particles migrate towards the corners, secondary flows are generated.

Figure 7.

Mean particle concentration along a line at $z/h=0.2$ and at the wall bisector, $z/h=1$ .

Figure 8.

Contours and vector fields of the secondary flow velocity $\sqrt{V_{f/p}^{2}+W_{f/p}^{2}}$ of the fluid (a) and solid phases (b).

Figure 9.

Mean and maximum value of the secondary flow velocity of the fluid phase as a function of the solid volume fraction $\unicode[STIX]{x1D719}$ . Results are normalized by the values of the single-phase case.

As shown in figure 6, the particle concentration close to the corners increases with the volume fraction. However, the mean particle distribution in the cross-section changes at the highest volume fraction considered, $\unicode[STIX]{x1D719}=0.2$ : the highest values of $\unicode[STIX]{x1D6F7}(y,z)$ are now found at the centre of the duct (see figure 7 b). This is not the case in turbulent channel flows and hence it can be related to the additional confinement of the suspension given by the lateral walls. As previously mentioned, the streamwise mean fluid and particle velocities are also substantially higher at the duct centre for $\unicode[STIX]{x1D719}=0.2$ .

To better quantify this effect, we analyse the numerical data by Picano et al. (Reference Picano, Breugem and Brandt2015) and calculate the number of particles crossing the spanwise periodic boundaries per unit time $h/U_{b}$ . For $\unicode[STIX]{x1D719}=20\,\%$ , we find that in 1 unit of $h/U_{b}$ approximately 1 % of the total number of particles cross the lateral boundaries. Inhibiting this lateral migration with lateral walls has therefore important consequences on the flow structure.

Concerning the secondary fluid velocity (figure 8 a), we note that for $\unicode[STIX]{x1D719}=0.2$ both the maximum and the mean values decrease below those of the unladen case. The presence of the solid phase leads to an increase of the secondary fluid velocity, which however saturates for volume fractions between 0.1 and 0.2. As shown by these mean velocity profiles, the turbulence activity is substantially reduced at $\unicode[STIX]{x1D719}=0.2$ .

Vector and contour plots of the secondary motions of the solid phase are reported in figure 8(b). These closely resemble those of the fluid phase. However, the velocities at the corner bisectors are almost half of those pertaining the fluid phase (in agreement with the fact that the particle concentration is high at these locations). Indeed, as discussed above there is a high particle concentration along the diagonals; particles reside for long times on these bisectors and, as a consequence they have lower cross-stream velocities than the fluid. Conversely, the cross-stream particle velocity is similar to that of the fluid phase at the walls, away from the corners. For $\unicode[STIX]{x1D719}=0.2$ we have previously noticed that along the duct walls, the highest particle concentration $\unicode[STIX]{x1D6F7}(y,z)$ is exactly at the corners. In agreement, we also observe from figure 8(b) that the secondary particle velocity is negligible at these locations (i.e. particles tend to stay at these locations for long times).

3.3 Velocity fluctuations

Next, we report the contours of the root-mean-square (r.m.s.) of the fluid and particle velocity fluctuations in outer units, see figure 10. Due to the symmetry around the corner bisectors, we show only the contours of $v_{f/p,rms}^{\prime }(y,z)$ in panels (c) and (d). Corresponding r.m.s. velocities along two lines at $z/h=0.2$ and 1 are depicted in figure 11. Results in inner units are shown at the wall bisector, $z/h=1$ in figure 12.

Figure 10.

(a) Root-mean-square of the streamwise fluid velocity fluctuations in outer units, $u_{f,rms}^{\prime }$ , for all $\unicode[STIX]{x1D719}$ . (b) Root-mean-square of the streamwise particle velocity fluctuations in outer units, $u_{p,rms}^{\prime }$ , for all $\unicode[STIX]{x1D719}$ . (c) Root-mean-square of the vertical fluid velocity fluctuations in outer units, $v_{f,rms}^{\prime }$ , for all $\unicode[STIX]{x1D719}$ . (d) Root-mean-square of the vertical particle velocity fluctuations in outer units $v_{p,rms}^{\prime }$ , for all $\unicode[STIX]{x1D719}$ .

The contours of the streamwise fluid velocity fluctuations reveal that r.m.s. values are stronger near the walls (close to the wall bisector), while minima are found along the corner bisectors. In this region, $u_{f,rms}^{\prime }(y,z)$ is substantially reduced for $\unicode[STIX]{x1D719}=0.2$ , as also visible from the profiles in figure 11(a). From figure 11(a) we also see that the local maxima of $u_{f,rms}^{\prime }(y)$ close to the corners is of similar magnitude for $\unicode[STIX]{x1D719}=0.0,0.05$ and 0.1. On the other hand, at the wall bisectors the maxima of $u_{f,rms}^{\prime }(y,z)$ decrease with $\unicode[STIX]{x1D719}$ , well below the value of the unladen case (see figure 11 b). For $\unicode[STIX]{x1D719}=0.2$ the profile of $u_{f,rms}^{\prime }(y)$ deeply changes. In particular, the streamwise r.m.s. velocity $u_{f,rms}^{\prime }(y)$ is substantially smaller than in the unladen case in the core region, where the particle concentration and the mean velocity are high while the secondary flows are negligible. After the maximum value, $u_{f,rms}^{\prime }(y)$ initially decreases smoothly and then sharply for $y/h>0.6$ .

From figures 10(b) and 11(a,b) we see that the streamwise r.m.s. particle velocity, $u_{p,rms}^{\prime }(y,z)$ , resembles that of the fluid phase. However, $u_{p,rms}^{\prime }(y,z)$ is typically smaller than $u_{f,rms}^{\prime }(y,z)$ in the cross-section. Exceptions are the regions close to the walls where the particle velocity does not vanish (unlike the fluid velocity).

Figure 11.

Root-mean-square of fluid and particle velocity fluctuations: (a,b) streamwise, (c,d) wall-normal and (e,f) spanwise ( $z$ -direction) components in outer units at $z/h=0.2$ (a,c,e) and $z/h=1$ (b,d,f), for all $\unicode[STIX]{x1D719}$ . Lines and symbols are used for the fluid and solid phase statistics, respectively.

From the contour of $v_{f,rms}^{\prime }(y,z)$ (see figure 10 c), we see that r.m.s. velocities are larger in the directions parallel to the walls, rather than in the wall-normal direction. Close to the wall bisectors, the peak values of the wall-normal and parallel fluid r.m.s. velocities increase with $\unicode[STIX]{x1D719}$ , again except for $\unicode[STIX]{x1D719}=0.2$ when turbulence is damped. From figure 11(c,e) we see instead that close to the corners, the local maxima of both  $v_{f,rms}^{\prime }(y)$ (wall normal) and $w_{f,rms}^{\prime }(y)$ (wall parallel) increase with respect to the unladen case, for all $\unicode[STIX]{x1D719}$ . At the wall bisector (see figure 11 d), wall-normal velocity fluctuations are slightly larger than the single-phase case for $\unicode[STIX]{x1D719}\lesssim 0.1$ .

Finally, figure 11(f) shows profiles at the wall bisector of the parallel component of the fluid velocity r.m.s., $w_{f,rms}^{\prime }(y)$ . Note that the peak value of $w_{f,rms}^{\prime }(y)$ increases with the volume fraction up to $\unicode[STIX]{x1D719}=0.1$ and moves closer to the wall. There is hence a clear redistribution of energy due to the particle presence towards a slightly more isotropic state.

Figure 12.

Root-mean-square of fluid and particle velocity fluctuations: (a) streamwise, (b) wall-normal and (c) spanwise ( $z$ -)direction components in inner units at $z/h=1$ for all $\unicode[STIX]{x1D719}$ . Lines and symbols are used for the fluid- and solid-phase statistics as in figure 11.

Concerning the solid phase, the wall-normal particle r.m.s. velocity, $v_{p,rms}^{\prime }$ , exhibits a peak close to the walls where particle layers form, see figures 10(d) and 11(d) (and also Costa et al. Reference Costa, Picano, Brandt and Breugem2016). Differently from the channel flow analysed by Picano et al. (Reference Picano, Breugem and Brandt2015), in the square duct flow $v_{p,rms}^{\prime }$ decreases slowly with increasing volume fraction $\unicode[STIX]{x1D719}$ in the wall particle layer ( $y^{+}\leqslant 20$ ) and the maxima are similar for all $\unicode[STIX]{x1D719}$ . This may be related to the existence of the secondary flows that pull the particles away from the walls towards the duct core. As for the fluid phase, the particle r.m.s. velocity parallel to the wall is greater than the perpendicular component close to the walls. From the profiles at the wall bisector, see figure 11(f), we also see that the maximum of $w_{p,rms}^{\prime }(y)$ is similar for $\unicode[STIX]{x1D719}=0.05$ and 0.1, while it clearly decreases for the densest case simulated. More generally, we observe that particle r.m.s. velocities are similar to those of the fluid phase towards the core of the duct.

Both fluid- and solid-phase r.m.s. velocities are shown in inner units at the wall bisector in figure 12(a–c). The local friction velocity has been used to normalize the velocity fluctuations. The fluid-phase r.m.s. velocity increases with $\unicode[STIX]{x1D719}$ in the viscous sublayer. Clearly, the presence of solid particles introduces additional disturbances in the fluid increasing the level of fluctuations in regions where these are typically low (in the unladen case). On the other hand, particle velocity fluctuations are typically smaller than the corresponding fluid r.m.s. velocity, except in the inner-wall region, $y^{+}<20$ , where they are one order of magnitude larger.

Figure 13.

Mean turbulence intensity $I=\langle u^{\prime }\rangle /U_{b}$ for all $\unicode[STIX]{x1D719}$ , normalized by the value for the unladen case.

Finally, we also computed the mean turbulence intensity defined as $I=\langle u^{\prime }\rangle /U_{b}$ , with $u^{\prime }=\sqrt{1/3(u_{f,rms}^{\prime 2}+v_{f,rms}^{\prime 2}+w_{f,rms}^{\prime 2})}$ . This is shown is figure 13, where results are normalized by the value obtained for the unladen case. As for secondary flows, we see that the turbulence intensity $I$ increases up to $\unicode[STIX]{x1D719}=0.1$ , for which $I$ is approximately 5 % larger than the unladen value. For $\unicode[STIX]{x1D719}=0.2$ , instead, the mean turbulence intensity decreases well below the single-phase value, and $I(\unicode[STIX]{x1D719}=0.2)\sim 0.9\,I(\unicode[STIX]{x1D719}=0)$ . It can be expected that when the volume fraction is larger than a threshold value ( $0.1<\unicode[STIX]{x1D719}\leqslant 0.2$ ), the particles are almost closely packed in the near-wall region. As a consequence, the quasi-coherent turbulent structures are quickly disrupted (as can be seen in figure 14) and the turbulence regenerating cycle is modified. As we will see in the next section, there is a strong reduction of the primary Reynolds stresses at $\unicode[STIX]{x1D719}=0.2$ , especially in the regions close to the wall bisector. Since the mean fluid velocity gradients are also reduced in the near-wall region with respect to $\unicode[STIX]{x1D719}=0.1$ , it can be expected that the associated turbulence production ( $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle \text{d}U_{f}/\text{d}y+\langle u_{f}^{\prime }w_{f}^{\prime }\rangle \text{d}U_{f}/\text{d}z$ ) will be strongly reduced, leading to the observed reduction in turbulence intensity.

In addition, we also calculated the fraction of the global kinetic energy that is contained in each phase, and the overall energy dissipation (per unit density) of the fluid phase, $\unicode[STIX]{x1D716}_{T}=2\unicode[STIX]{x1D707}\langle \unicode[STIX]{x1D60C}_{ij}\unicode[STIX]{x1D60C}_{ij}\rangle$ , averaged over time and volume, and normalized by $\unicode[STIX]{x1D708}U_{b}^{2}/(2h)^{2}$ ( $\unicode[STIX]{x1D60C}_{ij}=0.5(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}+\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i})$ is the strain-rate tensor). The global kinetic energy is defined as $K=K_{f}+K_{p}=0.5(U_{f,i}^{2}+u_{f,rms,i}^{\prime 2})+0.5(U_{p,i}^{2}+u_{p,rms,i}^{\prime 2})$ , where $K_{f}$ and $K_{p}$ are the total energies of the fluid and solid phase. We find that 6 %, 12 % and 24 % of $K$ is contained in the particle phase for $\unicode[STIX]{x1D719}=0.05,0.1$ and 0.2. Clearly, these values are proportional to the solid volume fraction.

Concerning the mean values of the overall energy dissipation, we see that they increase with the volume fraction, up to $\unicode[STIX]{x1D719}=0.1$ . In particular, $\unicode[STIX]{x1D716}_{T}=23,28$ and 30 for $\unicode[STIX]{x1D719}=0.0,0.05$ and 0.1. Instead, for $\unicode[STIX]{x1D719}=0.2$ , $\unicode[STIX]{x1D716}_{T}=28$ . Generally, a reduction in $\unicode[STIX]{x1D716}_{T}$ is related to turbulence attenuation as discussed, for example, by Spandan, Lohse & Verzicco (Reference Spandan, Lohse and Verzicco2016). As shown in figure 14, it is difficult to identify quasi-coherent structures at $\unicode[STIX]{x1D719}=0.2$ (due to the large number of particles near the walls), and as we will later explain, ejection events that are typically (very) intense close to the wall bisectors, substantially reduce in this case. The turbulence activity is hence reduced as shown by the reduction in energy dissipation, turbulence intensity and, as we will discuss in the next sections, of the primary Reynolds stresses and mean streamwise vorticity.

At $\unicode[STIX]{x1D719}=0.2$ , the stronger reduction in turbulence activity in comparison to channel flow may be due to the geometrical constraint given by the lateral walls. From the data of Picano et al. (Reference Picano, Breugem and Brandt2015), it can be found that on average 7.5 % of the particles reside near the walls (between $y/h=0$ and $y/h\sim 0.12$ ). In duct flow instead, we find that 13.1 % of the particles are located close to one of the 4 walls. There are hence almost twice as many particles interacting with the near-wall turbulence, leading to its stronger attenuation. It is also interesting to note that for all $\unicode[STIX]{x1D719}$ , the amount of particles located close to the walls is always approximately 13 % of the total. This observation can be relevant for modelling purposes as done by Costa et al. (Reference Costa, Picano, Brandt and Breugem2016) who divide a particle-laden turbulent flow into a near-wall particle layer and a homogeneous suspension bulk flow.

Figure 14.

The magnitude of the fluid velocity (normalized by $U_{b}$ ) at $y^{+}\simeq 20$ for $\unicode[STIX]{x1D719}=0.0,0.1$ and 0.2.

3.4 Reynolds stress and mean fluid streamwise vorticity

We now turn to the discussion of the primary Reynolds stress in the duct cross-section, as this plays an important role in the advection of mean streamwise momentum. We show in figure 15(a) the $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ component of the fluid Reynolds stress in the duct cross-section, for all $\unicode[STIX]{x1D719}$ . The component $\langle u_{f}^{\prime }w_{f}^{\prime }\rangle$ is not shown as it is the $90^{\circ }$ rotation of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ . We see that the maximum of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ , located close to the wall bisector, increases with the volume fraction up to $\unicode[STIX]{x1D719}=0.1$ . The maximum $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ then progressively decreases with $\unicode[STIX]{x1D719}$ , denoting a reduction in turbulent activity. For $\unicode[STIX]{x1D719}=0.2$ we observe that the maximum of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ reaches values even lower than in the unladen case. The contour of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ also changes with increasing $\unicode[STIX]{x1D719}$ . We find that $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ increases towards the corners, and also for $\unicode[STIX]{x1D719}=0.2$ it is larger than in the unladen case. However, the mean value of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ in one quadrant slightly increases up to $\unicode[STIX]{x1D719}=0.1$ , while for $\unicode[STIX]{x1D719}=0.2$ the mean is 26 % smaller than for the unladen case.

Figure 15.

Contours of the primary Reynolds stress $\langle u_{f/p}^{\prime }v_{f/p}^{\prime }\rangle$ of (a) the fluid and (b) particle phase for all volume fractions considered. The profiles along the wall bisector at $z/h=1$ are shown in outer and inner units in panels (c) and (d). Lines and symbols are used for the fluid- and solid-phase statistics as in figure 11.

The profiles of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ at the wall bisector ( $z/h=1$ ) are shown in figure 15(c). While the profiles for $\unicode[STIX]{x1D719}=0.05$ and 0.1 are similar and assume larger values than the reference case, we see that for $\unicode[STIX]{x1D719}=0.2$ the profile is substantially lower for all $z/h>0.2$ . This is also different to what is found in channel flow as close to the core, $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ is found to be similar for $\unicode[STIX]{x1D719}=0.0$ and 0.2. This is probably related to the high particle concentration at the duct core, see figure 6(c). Looking at the profiles in inner units, see figure 15(d), we also see that the peak values of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle ^{+}$ are similar for the unladen case and for the laden cases with $\unicode[STIX]{x1D719}=0.05$ and 0.1. This is in contrast to what is found in channel flows, where a reduction of the peak with $\unicode[STIX]{x1D719}$ is observed (Picano et al. Reference Picano, Breugem and Brandt2015). This shows that up to $\unicode[STIX]{x1D719}=0.1$ , the turbulence activity is not significantly reduced by the presence of particles. On the other hand, for $\unicode[STIX]{x1D719}=0.2$ we observe a large reduction in the maximum $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle ^{+}$ , denoting an important reduction in turbulent activity.

The Reynolds stress of the solid phase, $\langle u_{p}^{\prime }v_{p}^{\prime }\rangle$ , is also shown for comparison in figure 15(b–d). The profiles are similar to those of the fluid phase, although $\langle u_{p}^{\prime }v_{p}^{\prime }\rangle <\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ for all $\unicode[STIX]{x1D719}$ , except close to the walls for $\unicode[STIX]{x1D719}=0.2$ . These local maxima may be related to the higher local particle concentration in layers close to the walls.

In figure 15(c) we also show the profile of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ obtained from an additional simulation. Here, we have used the Eilers fit (Stickel & Powell Reference Stickel and Powell2005) to estimate the effective viscosity of the suspension at $\unicode[STIX]{x1D719}=0.2$ and found $\unicode[STIX]{x1D708}_{e}/\unicode[STIX]{x1D708}=[1+2.5\unicode[STIX]{x1D719}/(1-\unicode[STIX]{x1D719}/0.6)]^{2}=1.89$ . With this value we can define an effective bulk Reynolds number, $Re_{e}=Re_{b}\unicode[STIX]{x1D708}/\unicode[STIX]{x1D708}_{e}=2962$ , ideally the bulk Reynolds number of a suspension of uniformly distributed monodispersed spheres. The results of a new unladen simulation with $Re_{b}=Re_{e}=2962$ are compared in figure 15(c) with those of the case with $\unicode[STIX]{x1D719}=0.2$ . We see that for $\unicode[STIX]{x1D719}=0.2$ , the maximum of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ is lower than that in the simulation with $Re_{e}=2962$ . This indicates that the turbulence attenuation is larger than what can be expected by just modelling the suspension rheological properties (i.e. $\unicode[STIX]{x1D708}_{e}/\unicode[STIX]{x1D708}$ ). Therefore the specific size of the particles ( ${\sim}20^{+}$ units), and their high concentration at the walls (13.1 % of $N_{p}$ ) contribute strongly to the turbulence attenuation. This clearly shows, that the effects of a non-uniform particle distribution should be considered when trying to model such flows.

While for turbulent channel flow the cross-stream component of the Reynolds stress tensor, $\langle v_{f}^{\prime }w_{f}^{\prime }\rangle$ , is identically zero, in duct flow it is finite and contributes to the production or dissipation of mean streamwise vorticity. Hence it is directly related to the origin of mean secondary flows (Gessner Reference Gessner1973; Gavrilakis Reference Gavrilakis1992). This can be seen from the Reynolds-averaged streamwise vorticity equation for a fully developed single-phase duct flow,

(3.4) $$\begin{eqnarray}\displaystyle V_{f}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}}{\unicode[STIX]{x2202}y}+W_{f}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}^{2}(\langle w_{f}^{\prime 2}\rangle -\langle v_{f}^{\prime 2}\rangle )}{\unicode[STIX]{x2202}y\unicode[STIX]{x2202}z}+\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}-\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\langle v_{f}^{\prime }w_{f}^{\prime }\rangle -\unicode[STIX]{x1D708}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\unicode[STIX]{x1D6FA}_{f}=0, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $V_{f}$ and $W_{f}$ are the mean velocities in the two wall-normal directions, and the mean vorticity is defined as

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}_{f}=\frac{\unicode[STIX]{x2202}W_{f}}{\unicode[STIX]{x2202}y}-\frac{\unicode[STIX]{x2202}V_{f}}{\unicode[STIX]{x2202}z}. & & \displaystyle\end{eqnarray}$$

The first two terms of (3.4) represent the convection of mean vorticity by the secondary flow itself and have been shown to be almost negligible (Gavrilakis Reference Gavrilakis1992). The third term is a source of vorticity in the viscous sublayer due to the gradients in the anisotropy of the cross-stream normal stresses. The fourth term, involving the secondary Reynolds stress, acts as source or sink of vorticity. Finally, the last term represents the viscous diffusion of vorticity.

Figure 16.

Contours of the secondary Reynolds stress $\langle v_{f/p}^{\prime }w_{f/p}^{\prime }\rangle$ of (a) the fluid and (b) the particle phase for all volume fractions $\unicode[STIX]{x1D719}$ under investigation.

The contours of $\langle v_{f}^{\prime }w_{f}^{\prime }\rangle$ are shown in figure 16 for all $\unicode[STIX]{x1D719}$ . These are non-negligible along the corner bisectors (being directly related to mean secondary flows), and approximately one order of magnitude smaller than the primary Reynolds stress. Interestingly, we see that the maxima of $\langle v_{f}^{\prime }w_{f}^{\prime }\rangle$ strongly increases with $\unicode[STIX]{x1D719}$ up to $\unicode[STIX]{x1D719}=0.1$ and also for $\unicode[STIX]{x1D719}=0.2$ , the maximum value is still larger than that for $\unicode[STIX]{x1D719}=0.05$ . We also notice that regions of finite $\langle v_{f}^{\prime }w_{f}^{\prime }\rangle$ become progressively broader with increasing $\unicode[STIX]{x1D719}$ .

The secondary Reynolds stress of the solid phase, $\langle v_{p}^{\prime }w_{p}^{\prime }\rangle$ , resembles qualitatively that of the fluid phase, except at the corners were a high value of opposite sign is encountered. As a consequence, close to the corners this term may contribute in opposite way to the production or dissipation of vorticity with respect to $\langle v_{f}^{\prime }w_{f}^{\prime }\rangle$ .

We conclude this section by showing in figure 17 the mean streamwise fluid vorticity $\unicode[STIX]{x1D6FA}_{f}$ for all $\unicode[STIX]{x1D719}$ . The region of maximum vorticity at the wall is found between $z/h=0.2$ and $z/h=0.5$ for the unladen case and it extends closer to the corner for increasing $\unicode[STIX]{x1D719}$ . We also find that the maximum $\unicode[STIX]{x1D6FA}_{f}$ initially increases with the solid volume fraction (for $\unicode[STIX]{x1D719}=0.1$ the maximum $\unicode[STIX]{x1D6FA}_{f}$ is just slightly smaller than for $\unicode[STIX]{x1D719}=0.05$ ). This is expected as also the intensity of the secondary motions increase. The contours of $\unicode[STIX]{x1D6FA}_{f}$ become more noisy for $\unicode[STIX]{x1D719}=0.2$ , with a maximum value below that of the single-phase case. Also, for $\unicode[STIX]{x1D719}=0.2$ the location of maximum $\unicode[STIX]{x1D6FA}_{f}$ is similar to that found for the unladen case.

Figure 17.

Contours of the mean fluid streamwise vorticity $\unicode[STIX]{x1D6FA}_{f}$ (scaled by $h/U_{b}$ ) for all $\unicode[STIX]{x1D719}$ under investigation.

The mean streamwise vorticity attains a positive or negative sign in the near-wall region, depending on the specific octant considered. As we move further away from the walls, the mean streamwise vorticity changes sign. For the unladen case, the magnitude of the maximum vorticity at the walls is 2.5 times larger than the magnitude of the maximum vorticity of opposite sign found closer to the corner bisector. As reference and for clarity, we will consider the vorticity to be positive at the wall, and negative further away, as found for the top and bottom walls (for the lateral walls the situation is reversed). We find that the vorticity minimum approaches progressively the walls as $\unicode[STIX]{x1D719}$ increases up to $\unicode[STIX]{x1D719}=0.1$ , and that the ratio between the magnitudes of the maximum and minimum of vorticity increases up to ${\sim}3$ . On the other hand, for $\unicode[STIX]{x1D719}=0.2$ we clearly notice several local vorticity minima. One minimum is further away from the walls, at a location similar to what is found for the unladen case. Another minimum is further away towards the corner bisector. The other minimum is instead close to the corners. As we have previously shown, at the largest $\unicode[STIX]{x1D719}$ there is a significant mean particle concentration exactly at the corners and correspondingly we see two intense spots of vorticity around this location, antisymmetric with respect to the corner bisector.

Although (3.4) is valid for single-phase duct flow, we have calculated the convective, source/sink and diffusive terms for the cases of $\unicode[STIX]{x1D719}=0,0.05$ and 0.1. For $\unicode[STIX]{x1D719}=0.2$ there is a strong coupling between the dynamics of both fluid and solid phases, and it is therefore difficult to draw conclusions only by estimating the terms in (3.4).

As discussed by Gavrilakis (Reference Gavrilakis1992), the production of vorticity within the viscous sublayer is the main factor responsible for the presence of vorticity in the bulk of the flow. The main contribution to the production of vorticity is given by the term involving the gradients of the cross-stream normal stresses. For the unladen case, the maxima of this term are located close to the corners (at $y/h=0.016$ , $z/h=0.16$ from the bottom-left corner, similarly to what was found by Gavrilakis (Reference Gavrilakis1992)). Another positive, almost negligible contribution to the production of mean streamwise vorticity is given by the convective term. On the other hand, the diffusive term and the term involving the gradients of the secondary Reynolds stress give a negative contribution to the generation of vorticity in this inner-wall region. Note that the largest negative contribution due to the latter term is found close to the maximum production ( $(y/h,z/h)=(0.016,0.15)$ for $\unicode[STIX]{x1D719}=0$ ). Generally, we observe that the maximum production and dissipation increase and approach the corner as the volume fraction increases up to $\unicode[STIX]{x1D719}=0.1$ (not shown). In order to understand how the presence of particles contributes to the generation of vorticity, we calculate the ratio between the overall production and dissipation at the location of the maximum of the normal stress term (i.e. the summation of the convective term and the normal stress term, divided by the absolute value of the summation between the diffusive and secondary Reynolds stress terms). For the unladen case we have a good balance and the ratio is approximately 1. For $\unicode[STIX]{x1D719}=0.05$ and 0.1, the ratio is 0.96 and 0.95. Hence, the presence of the solid phase gives an additional contribution of approximately 5 % to the generation of mean streamwise vorticity in the near-wall region at the location of maximum production. The increased production due to larger gradients of the normal stress difference and due to the additional contribution by particles, leads globally to the larger $\unicode[STIX]{x1D6FA}_{f}$ observed at $\unicode[STIX]{x1D719}=0.05$ and 0.1.

3.5 Quadrant analysis

In this section we employ the quadrant analysis to identify the contribution from so-called ejection and sweep events to the production of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ and $\langle u_{f}^{\prime }w_{f}^{\prime }\rangle$ . In single-phase wall-bounded turbulent flows, it is known that these phenomena are associated with pairs of counter-rotating streamwise vortices that exist in the shear layer near the wall. These force low-momentum fluid at the wall towards the high-speed core of the flow. This is typically referred to as a Q2 or ejection event, with negative $u_{f}^{\prime }$ and positive $v_{f}^{\prime }$ . On the other hand, events with positive $u_{f}^{\prime }$ and negative $v_{f}^{\prime }$ are associated with the inrush of high-momentum fluid towards the walls and are known as sweeps or Q4 events. Events Q1 (positive $u_{f}^{\prime }$ and positive $v_{f}^{\prime }$ ) and Q3 (negative $u_{f}^{\prime }$ and negative $v_{f}^{\prime }$ ) are not associated with any particular turbulent structure when there is only one inhomogeneous direction. However, as discussed by Huser & Biringen (Reference Huser and Biringen1993) in turbulent duct flows these also contribute to the total turbulence production.

Figure 18.

Maps of probability of Q1, Q2, Q3, Q4 events in panels (a), (c), (e), (g) and (b), (d), (f), (h) for particle volume fractions $\unicode[STIX]{x1D719}=0.0$ and 0.05, respectively.

Figure 19.

Maps of probability of Q1, Q2, Q3, Q4 events in panels (a), (c), (e), (g) and (b), (d), (f), (h) for particle volume fractions $\unicode[STIX]{x1D719}=0.1$ and 0.2, respectively.

We first show contours of the probability of finding Q1, Q2, Q3, Q4 events for $\unicode[STIX]{x1D719}=0.0$ and 0.05 in figure 18, and for $\unicode[STIX]{x1D719}=0.1$ and 0.2 in figure 19. The maximum probability of an event is 1, so that the sum $\text{Q1}+\text{Q2}+\text{Q3}+\text{Q4}=1$ at each point. Since these events are typically important close to the walls, the $y$ -coordinate is reported in inner units (the viscous length at the wall bisector is chosen), while the wall-parallel $z$ -coordinate spans half of the duct. The contours of the probability of the different events directly enables us to easily compare and discuss all cases.

Results for the unladen case are substantially in agreement with those by Huser & Biringen (Reference Huser and Biringen1993) and Joung et al. (Reference Joung, Choi and Choi2007). Ejection (Q2) events are important at the wall bisector ( $y/h=1$ ) and around $y/h\sim 0.8$ . Joung et al. (Reference Joung, Choi and Choi2007) report strong Q2 events already $y/h=0.62$ , and this is probably due to the fact that their $Re_{b}$ was smaller than ours ( $Re_{b}=4410$ instead of 5600). Q4 events are instead dominant closer to the vertical walls, around $y/h=0.3$ (as also reported by Joung et al. (Reference Joung, Choi and Choi2007)). Q3 events are found when both $u_{f}^{\prime }$ and $v_{f}^{\prime }$ are negative. As also shown by Huser & Biringen (Reference Huser and Biringen1993), we see that these increase from the wall bisector towards the vertical walls, being dominant at the corners. Consequently, these events are created by ejections from the vertical wall (and particularly from the corners), resulting in the increase of $\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$ in this region.

The general picture is similar for $\unicode[STIX]{x1D719}=0.05$ and 0.1. The probability maps of Q1, Q2, Q3 and Q4 events are only slightly changed with respect to those of the unladen case. Concerning the probability of Q3 events we see that it is still higher at the corner and that it increases with $\unicode[STIX]{x1D719}$ along the vertical wall and along the corner bisector. It is clear from figure 19 that the turbulence activity is strongly reduced for $\unicode[STIX]{x1D719}=0.2$ . The probability of ejection and sweep events close to the walls is drastically reduced. The probability of sweep (Q4) events is around 35 % only in a small region close to the wall bisector ( $y/h\sim 0.85$ ). Instead, the probabilities of Q2 (ejections) and Q3 events are above 35 % in a region around the corner bisector, where we recall there is high mean particle concentration ( $y/h\in (0;0.4)$ ). At this volume fraction the probability of these events is therefore mostly correlated to fluid–particle interactions. The probability contours for $\unicode[STIX]{x1D719}=0.2$ are indeed very different from those at lower $\unicode[STIX]{x1D719}$ and from those of the unladen case. The absence of high probability regions of Q2 and Q4 events around $z/h\sim 0.8$ and $z/h\sim 0.3$ denotes the disruption of the coherent streamwise vortical structures.

The presence of two inhomogeneous directions allows us to extend the quadrant analysis to the secondary shear stress $\langle v_{f}^{\prime }w_{f}^{\prime }\rangle$ , which contributes to the production, dissipation and transport of mean streamwise vorticity, as discussed above. As in Huser & Biringen (Reference Huser and Biringen1993), to distinguish between the events related to the primary and secondary Reynolds stresses, we will refer to the latter as $\text{Q}1_{s},\text{Q}2_{s},\text{Q}3_{s}$ and $\text{Q}4_{s}$ events. These authors showed that as the corners are approached, contributions of the $\text{Q}1_{s}$ and $\text{Q}3_{s}$ events progressively decrease. On the other hand, $\text{Q}2_{s}$ and $\text{Q}4_{s}$ events are found to be stronger than the latter close to the corners. This is due to the interaction between ejections from both walls that tilt the ejection stem toward the perpendicular wall. We looked at probability maps of these events as previously done for the events related to the primary Reynolds stress. For the unladen case we see indeed that close to the corners the probability of having $\text{Q}2_{s}$ and $\text{Q}4_{s}$ events is substantially larger than that of $\text{Q}1_{s},\text{Q}3_{s}$ events (not shown). As we increase the volume fraction $\unicode[STIX]{x1D719}$ the probability of these events remains high until the highest $\unicode[STIX]{x1D719}=0.2$ is reached. In this case, we find that there is an approximately equal share of probability between all type of events indicating again that there is a substantial alteration of the turbulence and less structured flow in the presence of many particles.