4.1. Enhanced drag coefficients by turbulence

In this section, the effects of the turbulence environments on the particle drag are examined in detail. Nine cases are simulated, and their turbulence intensity $I = \langle u'\rangle /U$ ranges from 0.132 to 0.411, where $\langle u'\rangle$ is the turbulent r.m.s. velocity and $U$ is the mean flow velocity computed from the given $Re_p$ in each case. This range of turbulence intensity is chosen to cover the conditions in most experimental studies where the turbulence effect on the particle drag was observable. As mentioned, once the turbulence freely streams into region B, its intensity decays due to the viscous dissipation, and its isotropy is no longer maintained due to the redistribution of TKE among three components. The latter change can be seen from the time-averaged component-wise TKE budget equations

(4.1a)$$\begin{gather} \frac{{\rm d}}{{\rm d}\kern0.06em x}\left[\left\langle\frac{1}{2}\left(U + u_x'\right)u_x'^2\right\rangle + \frac{1}{\rho}\Big\langle p'u_x' \Big\rangle - 2\nu\Big\langle s_{xx}'u_x' \Big\rangle \right] = \Big\langle p's_{xx}'\Big\rangle - 2\nu\Big\langle s_{xj}'s_{xj}'\Big\rangle, \end{gather}$$

(4.1b)$$\begin{gather}\frac{{\rm d}}{{\rm d}\kern0.06em x}\left[\left\langle\frac{1}{2}\left(U + u_x'\right)u_y'^2\right\rangle - 2\nu \left\langle s_{xy}'u_y' \right\rangle \right] = \left\langle p's_{yy}'\right\rangle - 2\nu \left\langle s_{yj}'s_{yj}'\right\rangle, \end{gather}$$

(4.1c)$$\begin{gather}\frac{{\rm d}}{{\rm d}\kern0.06em x}\left[\left\langle\frac{1}{2}\left(U + u'\right)u_z'^2\right\rangle - 2\nu \left\langle s_{xz}'u_z' \right\rangle \right] = \left\langle p's_{zz}'\right\rangle - 2\nu\left\langle s_{zj}'s_{zj}'\right\rangle, \end{gather}$$

where $u_x'$, $u_y'$, $u_z'$ are the fluctuation velocities in the $x$, $y$ and $z$ directions, and $p'$ is the fluctuation pressure.

To better assess the flow information at the particle location, unladen cases corresponding to the nine examined cases are simulated, and the statistics of these unladen simulations are tabulated in table 4. These statistics (spatial and time averaged) are measured in a $y$–$z$ plane sized $L^2$ centred at $(x_c - 0.5d_p,y_c,z_c)$, corresponding to the flow stagnation point in the particle-laden cases. A similar approach was used by Botto & Prosperetti (Reference Botto and Prosperetti2012), except that they measured the statistics from the $y$–$z$ plane corresponding to the particle centre. As shown in table 4, the turbulent incoming flows seen by the particle are indeed anisotropic. To assist the comparisons between the drag enhancements from the present study and those reported by Homann et al. (Reference Homann, Bec and Grauer2013), an overall r.m.s. velocity $\langle u'\rangle = \sqrt {(\langle u_x'\rangle ^2 + \langle u_y'\rangle ^2 + \langle u_z'\rangle ^2)/3}$ is defined for the present cases, and the turbulence intensity $I = \langle u'\rangle /U$ can be computed accordingly, as listed in table 4. It should be made clear that the purpose of conducting these unladen cases is not to obtain the instantaneous undisturbed fluid velocity to assess the reliability of the particle equations of motion, as intended in some previous studies (e.g. Bagchi & Balachandar Reference Bagchi and Balachandar2003; Burton & Eaton Reference Burton and Eaton2005; Zeng et al. Reference Zeng, Balachandar, Fischer and Najjar2008). In the present study, we are only interested in how the turbulence modifies the statistical average of the particle drag and how these modifications can be predicted from the statistics of the turbulent flows. The probability distribution functions (p.d.f.s) of $u_x'$, $u_y'$ and $u_z'$ in three selected cases out of the nine cases are shown in figure 4. It can be clearly seen that the distributions of the incoming flow velocities in all three directions are still close to Gaussian.

Table 4.

Information of the incoming flow.

Figure 4.

The probability distribution functions (p.d.f.s) of the fluctuation velocity components of the incoming turbulent flows: (a) case C, $I = 0.410$, (b) case E, $I = 0.265$, (c) case G, $I = 0.132$. The black solid line in each plot represents the Gaussian distribution.

The drag enhancement coefficient, defined as $\langle F_D^T/F_D^L\rangle - 1$, where $\langle F_D^T\rangle$ and $F_D^L$ are the time-averaged drag forces in turbulent and laminar flows, respectively, are plotted in figure 5, together with the results reported in the Homann et al. (Reference Homann, Bec and Grauer2013) study. To increase the accessibility of our data, they are also tabulated in table 5. It is evident that turbulence significantly enhances the drag forces on the particle in all cases. The empirical correlation suggested by Uhlherr & Sinclair (Reference Uhlherr and Sinclair1971) for drag coefficients of spherical particles in turbulent flows,

(4.2)\begin{equation} C_D^{US} = 0.133\left(1 + \frac{150}{Re_p}\right)^{1.565} + 4I,\quad 50< Re_p<700,\ 0.05< I<0.5, \end{equation}

is also shown in figure 5 by the dashed lines. Since this correlation does not apply to laminar flow, its corresponding results of $\langle F_D^T/F_D^L\rangle - 1$ are computed as $C_D^{US}/C_{D}^{SN} - 1$, where $C_D^{SN}$ is the drag coefficient predicted by the Schiller & Naumann correlation. Uhlherr & Sinclair (Reference Uhlherr and Sinclair1971) predicted linear increases of drag enhancements with the turbulence intensity, which roughly agrees with the numerical data reported by Homann et al. (Reference Homann, Bec and Grauer2013) and this study. However, the quantitative values of the drag enhancement coefficients are significantly overestimated by the Uhlherr–Sinclair empirical correlation.

Figure 5.

The drag enhancement coefficients under different turbulence intensities: $\circ$, data from the present study; $*$, data from the Homann et al. (Reference Homann, Bec and Grauer2013) study; solid lines, nonlinear drag enhancements; dashed lines, the empirical correlation suggested by Uhlherr & Sinclair (Reference Uhlherr and Sinclair1971). Symbols and curves with the same colour have the same $Re_p$.

Table 5.

Comparison of nonlinear drag enhancement coefficients based on different models, and the actual drag enhancement coefficients from the simulations (last column).

The drag enhancement by turbulence, as analysed in the previous studies (e.g. Homann et al. Reference Homann, Bec and Grauer2013; Fornari et al. Reference Fornari, Picano, Sardina and Brandt2016b), is partially due to the nonlinear dependency of the particle drag on the flow velocity. Given that the incoming flow velocity is Gaussian distributed in each direction (confirmed in figure 4), Homann et al. (Reference Homann, Bec and Grauer2013) computed the averaged drag force as

(4.3a,b)\begin{equation} F_{D}^{H} \equiv \int f^{iso}\left({\boldsymbol{u}}\right) F_{D}^{SN}\left({\boldsymbol{u}}\right) {\rm d}{\boldsymbol{u}},\quad {\boldsymbol{u}}=(U+u_x',u_y',u_z'), \end{equation}

where

(4.4)\begin{equation} F_{D}^{SN}\left({\boldsymbol{u}}\right) = 3{\rm \pi} d_p \mu \left[1 + 0.15\left(\frac{\vert{\boldsymbol{u}}\vert d_p}{\nu}\right)^{0.687}\right]\left({ \boldsymbol{u}}\boldsymbol{\cdot} {\boldsymbol{e}}_x\right) \end{equation}

is the drag force predicted by the Schiller & Naumann correlation and

(4.5a,b)\begin{equation} f^{iso}({\boldsymbol{u}}) = \frac{1}{(\sqrt{2{\rm \pi}}\langle u'\rangle)^3 } \exp\left(-\frac{u_x'^2+u_y'^2+u_z'^2}{2 \langle u'\rangle^2}\right),\quad\langle u'\rangle = IU\end{equation}

is the p.d.f. of the instantaneous fluctuation velocity of the incoming flow assuming an isotropic distribution of TKE. The integral in (4.3a,b) can be numerically obtained, and the relative changes of $\langle F_D^{H}/F_D^L\rangle - 1$ with the turbulence intensity are shown by the solid lines in figure 5. It is clear that the enhancements due to the nonlinear drag effect are below the actual drag enhancements and are $Re_p$ insensitive after $Re_p > 200$, which is different from the actual drag enhancement where $Re_p$ plays a much more important role. For the present cases where the incoming turbulent flows are anisotropic, the nonlinear drag $F_D^{P}$ can be computed by replacing $f^{iso}$ with a more general p.d.f.

(4.6)\begin{equation} f^{aniso}\left({\boldsymbol{u}}\right) = \frac{1}{(\sqrt{2{\rm \pi}})^3 \langle u_x'\rangle \langle u_y'\rangle \langle u_z'\rangle} \exp\left(-\frac{u_x'^2}{2 \langle u_x'\rangle^2}-\frac{u_y'^2}{2\langle u_y'\rangle^2}-\frac{u_z'^2}{2\langle u_z'\rangle^2}\right), \end{equation}

and with the information of $\langle u_x'\rangle$, $\langle u_y'\rangle$ and $\langle u_z'\rangle$ in table 4. We also note that, in the earlier work of Bagchi & Balachandar (Reference Bagchi and Balachandar2003), two definitions of the particle drag force accounting for the nonlinear effect were proposed, and they are

(4.7a)$$\begin{gather} F_{D}^{BB1} = \frac{1}{8}\rho_f {\rm \pi}d_p^2 \langle C_D^{SN}(Re_p)\vert {\boldsymbol{u}}\vert^2 \rangle,\quad Re_p = Re_p\left({\boldsymbol{u}}\right) = \frac{\vert {\boldsymbol{u}}\vert d_p}{\nu}, \end{gather}$$

(4.7b)$$\begin{gather}F_{D}^{BB2} = \tfrac{1}{8}\rho_f {\rm \pi}d_p^2 C_D^{SN}\left(\langle Re_p\rangle\right) \vert\langle{\boldsymbol{u}}\rangle\vert^2 , \end{gather}$$

where

(4.8a,b)\begin{align} \langle C_D^{SN}(Re_p)\vert {\boldsymbol{u}}\vert^2 \rangle = \int f^{aniso}\left({\boldsymbol{u}}\right)C_D^{SN}\left(Re_p\right)\vert {\boldsymbol{u}}\vert^2 \,{\rm d}{\boldsymbol{u}}, \quad \langle Re_{p}\rangle = \int f^{aniso}\left({\boldsymbol{u}}\right) Re_p\left({\boldsymbol{u}}\right) {\rm d}{\boldsymbol{u}}. \end{align}

Given the nonlinear dependency of the instantaneous drag force on the flow velocity, it would make more sense that the drag depends on the average of the velocity squared, i.e. $\langle \vert {\boldsymbol {u}}\vert ^2\rangle$, rather than the square of the averaged velocity, i.e. $\vert \langle {\boldsymbol {u}}\rangle \vert ^2$ in (4.7b). With the Gaussian distribution of the fluctuation velocity, we shall also have $\langle \boldsymbol {u}\rangle = U\boldsymbol {e}_x$, which is constant. Based on this argument, a more reasonable definition for the nonlinear drag is formulated as

(4.9a,b) \begin{equation} F_{D}^{BB3} = \tfrac{1}{8}\rho_f {\rm \pi}d_p^2 C_D^{SN}\left(\langle Re_p\rangle\right)\langle\vert{\boldsymbol{u}}\vert^2\rangle,\quad\langle\vert{\boldsymbol{u}}\vert^2\rangle = \int f^{aniso}\left({\boldsymbol{u}}\right) \vert{\boldsymbol{u}}\vert^2 \,{\rm d}{\boldsymbol{u}}. \end{equation}

In table 5, these nonlinear drag enhancements are compared with the actual drag enhancements obtained from the numerical simulations. With the Homann et al. (Reference Homann, Bec and Grauer2013) definition of the nonlinear drag force, no matter whether the anisotropy in the incoming turbulent flows is considered, the nonlinear drag forces are significantly lower than their actual counterparts for all examined cases. With the Bagchi & Balachandar (Reference Bagchi and Balachandar2003) definition, the nonlinear drag forces are also considerably smaller than the actual drag forces, with only two exceptions in case B and case C. In particular, the original definition of the nonlinear drag in (4.7b) always shows drag reduction rather than drag enhancement, which could indicate this is an inappropriate estimation for the nonlinear drag. Compared with the Homann et al. (Reference Homann, Bec and Grauer2013) definition of nonlinear drag, the nonlinear drag of Bagchi & Balachandar (Reference Bagchi and Balachandar2003) is significantly higher. The difference is due to the different definitions of drag in the two works. In the Homann et al. (Reference Homann, Bec and Grauer2013) work, the drag was defined as the force component aligned with the mean flow direction, whereas in Bagchi & Balachandar (Reference Bagchi and Balachandar2003), it was the force aligned with the instantaneous relative velocity between the flow and the particle. Hence, the magnitudes of the latter are certainly higher. This detail partially explains the conflicting conclusions in the two works, i.e. Homann et al. (Reference Homann, Bec and Grauer2013) reported the inadequacy of using the nonlinear drag to explain the drag enhancement due to turbulence, while Bagchi & Balachandar (Reference Bagchi and Balachandar2003) claimed that the standard drag formula was already sufficient to predict the drag coefficient.

Based on the numerical results, this study adopts the Homann et al. (Reference Homann, Bec and Grauer2013) point and seeks mechanisms resulting in the differences between the nonlinear drag enhancements and the actual drag enhancements. It should be noted that, although the r.m.s. velocities of the turbulent flow are used in the calculation of $\langle F_D^{H}\rangle$, the nonlinear drag enhancement is essentially a laminar effect, which also exists when the incoming flow is laminar but time dependent.

4.2. Additional mechanisms of drag enhancements

Homann et al. (Reference Homann, Bec and Grauer2013) attributed the gaps between the actual drag enhancements and the nonlinear drag enhancements to the interactions between the small turbulent eddies and the particle boundary layer. This deduction is reasonable, but the associated mechanisms have not been clarified. In order to better understand the mechanisms of how turbulence affects the particle drag, the flow fields around the particle is visualized. In figure 6, the pressure fields around the particle are visualized for case E and its corresponding laminar case with the same $Re_p$. To ensure the two pressure fields can be compared at the same scale, we plot the deviation of the pressure from its instantaneous field mean, then average over 1000 time frames covering approximately 45 eddy turnover times, i.e. $\langle p(\boldsymbol {x}) \rangle - \langle P \rangle$. Although the wake structure is planar symmetric at this particle Reynolds number ($Re_p = 240$ in case E), the averaged pressure field still shows symmetry in the lower and top halves, so only half of the region is shown. To facilitate the following discussion, we locate $(x_c - 0.5d_p,y_c,z_c)$, $(x_c + 0.5d_p,y_c,z_c)$ and $(x_c,y_c+0.5d_p,z_c)$ as the front edge, trailing edge and the zenith, respectively. Compared with the laminar case, when the incoming flow is turbulent, the low-pressure region around the zenith is significantly extended, and the magnitudes of the pressure near the trailing edge are also reduced. These two changes certainly result in increased pressure drag. The pressure gradients between the high- and low-pressure regions are also enlarged by turbulence, as shown in the contours around the top-left corner. At the fixed particle Reynolds number $Re_p = 100$, the pressure contours from the laminar and turbulent flow cases with different intensities are compared in figure 7. It is clear that when the turbulence intensity increases, pressure gradients between the front edge and the zenith are further enhanced. As a result, the area of the low-pressure regions expands with the increase of the turbulence intensity, and the pressure in the low-pressure regions further decreases.

Figure 6.

The pressure fields around the fixed particle in (a) case E and (b) the corresponding laminar flow.

Figure 7.

The pressure contours around the particle in (a) the laminar flow, and turbulent flows with different intensities; (b) case A, (c) case B and (d) case C.

The strengthened low-pressure region in the turbulent cases also modifies the wake structures. Contours of the velocity magnitude and $Q$-criterion (Jeong & Hussain Reference Jeong and Hussain1995) are compared in figure 8 for case E and its corresponding laminar case. More high-speed fluids outside the boundary layer are transported towards the symmetry axis in response to the strengthened low-pressure region around the trailing edge, which not only suppresses the circulation region but also results in faster restorations of the velocity deficiency, as reported in the previous studies (e.g. Bagchi & Balachandar Reference Bagchi and Balachandar2004; Homann et al. Reference Homann, Bec and Grauer2013).

Figure 8.

Contours of the velocity magnitude and $Q$-criterion in case E (a,c) and the corresponding laminar case (b,d). (a,b) Velocity magnitude, (c,d) $Q$-criterion.

As turbulence reduces the pressure behind the zenith, the region with adverse pressure gradients shrinks. This change postpones the separation of the boundary layer, as can be partially seen from figure 8. In the classic boundary layer theory, postponed separation usually results in drag reductions (Kundu, Cohen & Dowling Reference Kundu, Cohen and Dowling2015, Chap. 9). However, this is not the case in the present study. This exception is because, unlike the uniform incoming flow in the classic boundary layer theory, the incoming flows in the present cases are already turbulent. In the classic boundary layer theory, turbulent fluctuations mainly affect the pressure after the zenith by resisting the boundary layer separation, which leads to better pressure recovery to reduce the pressure drag (Kundu et al. Reference Kundu, Cohen and Dowling2015, Chap. 9). In the present study, however, a significant portion of the turbulence effect functions before the zenith, as shown by the amplified pressure gradients around the top-left corner in the turbulent cases, i.e. figures 6(a) and 7(b)–(d). In order to show this phenomenon more clearly, streamwise distributions of the planar fluid-phase-averaged pressure around the particle are computed in some selected laminar and turbulent cases. The averaging is processed over $y$–$z$ planes with a size of $2.5d_p\times 2.5d_p$ that centre on the streamwise axis of the particle. As confirmed in figure 9, compared with the corresponding laminar case, the streamwise pressure gradient in the turbulent case is significantly larger, and most of the difference occurs before the zenith.

Figure 9.

Additional pressure drop in the streamwise direction around the particle in a turbulent flow relative to that in a laminar flow. Here, $P_0$ is the planar fluid-phase-averaged pressure at the front edge.

The amplified pressure gradients between the front edge and the zenith can be qualitatively understood by referring to the solution of the pressure field of the inviscid flow (Kundu et al. Reference Kundu, Cohen and Dowling2015, Chap. 4), where

(4.10)\begin{equation} p = p_{\infty} + \frac{1}{2}\rho U^2\left[\left(\frac{3a^3}{r^3} - \frac{3a^6}{4r^6}\right)\cos^2\theta - \left(\frac{a^3}{r^3} + \frac{a^6}{4r^6}\right) \right]. \end{equation}

The pressure gradient $\partial p/\partial \theta$ yields

(4.11)\begin{equation} \frac{\partial p}{\partial \theta} ={-}\frac{1}{2}\rho U^2\left(\frac{3a^3}{r^3} - \frac{3a^6}{4r^6}\right)\sin\left(2\theta\right) . \end{equation}

The region between the front edge and the zenith is where $\partial p/\partial \theta < 0$. The above solution does not apply rigorously for the present cases with viscous turbulent incoming flows. However, experimental measurements (Maxworthy Reference Maxworthy1969; Achenbach Reference Achenbach1972) found that the lowest pressure in a turbulent boundary layer around a sphere was still $\theta = 90$ degrees from the stagnation point and the pressure distribution was similar to the inviscid flow theory before the lowest pressure point. Despite the fact that these experimental measurements were conducted at much higher $Re_p$ than the present study, considering the similarity in their wake structures, it is reasonable to expect that the pressure gradient follows

(4.12)\begin{equation} \frac{\partial p}{\partial \theta} \propto \langle\vert {\boldsymbol{u}}\vert^2\rangle\sin\left(2\theta\right), \end{equation}

before the boundary layer separation. Equation (4.12) still predicts the lowest pressure at the zenith, which has been confirmed by the contour plots in figures 6 and 7. With the above assumption, it is clear that the more kinetic energy the flow has, the larger the pressure gradient expected. Compared with its laminar counterpart, a turbulent incoming flow has extra kinetic energy in the turbulent motions, i.e. TKE. Thus a greater pressure change is created, and also a greater resulting pressure drag. For the same reason, the pressure drop would further increase with the turbulence intensity. Since the large-scale motions mainly carry TKE in turbulent flows, this mechanism of pressure drag amplification is more associated with the large-scale eddies.

Different from large-scale eddies, the main effect of small turbulent eddies on the drag enhancement is to increase the velocity gradient in the particle boundary layer. This mechanism is the same one that increases the surface shear stress in the flat turbulent boundary layer. Since the particle has a finite size, the curvature of the particle might become unimportant for small scales in the turbulent flow. Unlike a laminar boundary layer that mainly transports momentum in the wall-normal direction via viscous diffusion, when turbulence is present, the small-scale velocity fluctuations in the wall-normal direction can directly bring high-speed fluids outside the boundary layer to the near-wall region and increase the velocity gradients there. As a result, the viscous drag is enhanced. With uniform incoming flows, this mechanism only occurs at very high Reynolds numbers after the transition from a laminar to a turbulent boundary layer occurs. In the present cases, the incoming flows are already turbulent, so the enhancement of the viscous drag by small turbulent eddies could happen at much lower Reynolds numbers. To confirm this mechanism, the streamwise velocity profiles along two lines, ${\boldsymbol {e}}_1 = (0,1,0)$ and ${\boldsymbol {e}}_2 = (-1/\sqrt {2},1/\sqrt {2},0)$ across the particle centre are plotted in figure 10. The turbulent cases generally have larger streamwise velocity gradients than corresponding laminar cases, especially with smaller $Re_p$. When $Re_p$ increases, the relative increase of the streamwise velocity gradient due to turbulence reduces. Case E with $Re_p = 240$ has a similar turbulence intensity to that in case B with $Re_p = 100$, but its enhancements of streamwise velocity gradient are much less evident.

Figure 10.

The distribution of the streamwise velocity along two selected lines across the particle centre: (a) a vertical line, (b) an inclined line.

4.3. Scaling of the drag enhancement

In the previous study of Homann et al. (Reference Homann, Bec and Grauer2013), it was reported that the drag enhancement coefficient $F_D^T/F_D^L - 1$ scales well with $0.18t_U/t_\kappa$, where $t_U = d_p/U$ is the convective time, and $t_\kappa = (\nu /\varepsilon )^{0.5}$ is the Kolmogorov time. The data points from the present study are also plotted using this scaling in figure 11 together with the data reported in the Homann et al. (Reference Homann, Bec and Grauer2013) study. It is worth mentioning that the dissipation rates used to evaluate $t_\kappa = (\nu /\varepsilon )^{0.5}$ in the present study are estimated indirectly from $\langle \varepsilon \rangle = \langle u'\rangle ^{3}\langle \varepsilon _{HIT}\rangle /\langle u'_{HIT}\rangle ^{3}$, where $\langle u'_{HIT} \rangle$ and $\langle \varepsilon _{HIT} \rangle$ are measured from the sustained HIT in region A (table 1). This essentially assumes that the integral length scales $L$ are unchanged while the turbulent flows travel downstream, i.e. $L = \langle u'\rangle ^3/\langle \varepsilon \rangle$ is unchanged.

Figure 11.

The drag enhancement coefficients under the Homann et al. (Reference Homann, Bec and Grauer2013) scaling.

As shown in figure 11, although the drag enhancement coefficients in the present study still qualitatively increase with $t_U/t_\kappa$, their values are considerably higher than the predictions with the Homann et al. (Reference Homann, Bec and Grauer2013) scaling. These quantitative deviations are likely due to the difference in flow configurations. Specifically, as the turbulent flow loses isotropy while travelling downstream, the anisotropy of the turbulence could play an important role in the drag enhancement.

In order to see this point more clearly, a volume-averaged streamwise momentum balance analysis is conducted here. In recent years, such analyses have been widely employed to enhance the understanding of how particles modify the mean flow velocity and TKE in particle-laden turbulent channel and pipe flows (e.g. Picano, Breugem & Brandt Reference Picano, Breugem and Brandt2015; Yu et al. Reference Yu, Lin, Shao and Wang2017; Vreman & Kuerten Reference Vreman and Kuerten2018; Peng et al. Reference Peng, Ayala and Wang2019c). By volume averaging the streamwise momentum equation, the difference between the particle drag in turbulent and laminar flows is related to flow properties at the front and trailing edges of the particle as

(4.13)\begin{equation} \langle F_D^T - F_D^L\rangle \approx{-}\frac{1}{\Delta x}\left.\left( \frac{1}{\rho}P^T - \frac{1}{\rho}P^L + k_x^{T} - k_x^{L}\right)\right\rvert_{x_c-({1}/{2})d_p}^{x_c+({1}/{2})d_p}. \end{equation}

Detailed derivation of this equation is provided in Appendix A. Equation (4.13) shows that, in an ideal case, the particle drag enhancement due to turbulence is only related to two changes, the amplified pressure drops in the turbulent cases (as confirmed in figure 9) and the modification of the streamwise kinetic energy distribution. Compared with laminar boundary layers, turbulent boundary layers could provide more viscous dissipation to the kinetic energy, further increasing its drop in the streamwise direction. Since the modification of streamwise kinetic energy dominates the drag enhancement, in anisotropic turbulence, the influence of flow anisotropy should be considered when modelling the drag enhancement based on the fluid and particle parameters.

We start to seek an appropriate scaling by conducting a dimensional analysis. When the turbulent flow is homogeneous and isotropic, it is natural to think that the drag enhancement coefficient is a function of the turbulent r.m.s. velocity $\langle u'\rangle$, dissipation rate $\langle \varepsilon \rangle$, viscosity $\nu$, the mean flow velocity $U$ and the particle diameter $d_p$, i.e.

(4.14)\begin{equation} \left\langle\frac{F_D^T}{F_D^L}\right\rangle - 1 = f\left(\langle u'\rangle, \langle\varepsilon\rangle, \nu, U, d_p\right). \end{equation}

The first three parameters define a HIT, while the last three define the particle Reynolds number. There are five dimensional parameters in this function and two basic units. Thus we can define three non-dimensional parameters, and they are

(4.15)\begin{equation} \left\langle\frac{F_D^T}{F_D^L}\right\rangle - 1 = f\left(\frac{\langle u'\rangle}{U}, \frac{d_p}{\eta}, \frac{U d_p}{\nu}\right) = f(I,d_r,Re_p), \end{equation}

where $d_r = d_p/\eta$ is short-hand notation for the relative particle size. Note that, since the function is undetermined, we use the same $f$ in (4.14) and (4.15) for convenience. With some mathematical manipulation, the Homann et al. (Reference Homann, Bec and Grauer2013) scale $t_U/t_\kappa$ can be recast as $Re_p^{-1}d_r^2$, i.e.

(4.16)\begin{equation} \frac{t_U}{t_\eta} = \frac{d_p}{U}\left(\frac{\varepsilon}{\nu} \right)^{0.5}= \frac{\nu}{Ud_p}\frac{d_p^2}{\left(\nu^3/\varepsilon\right)^{0.5}} = Re_p^{{-}1}\left(\frac{d_p}{\eta}\right)^2 = Re_p^{{-}1}d_r^2, \end{equation}

which is a special form of (4.15). For a more general anisotropic turbulent flow, a single r.m.s. velocity $\langle u'\rangle$ is no longer sufficient to describe the turbulence intensity, we therefore formulate

(4.17)\begin{equation} \left\langle\frac{F_D^T}{F_D^L}\right\rangle - 1 = f_{new}\left(\langle u_1'\rangle, \langle u_2'\rangle, \langle u_3'\rangle, \langle\varepsilon\rangle, \nu, U, d_p\right), \end{equation}

where $\langle u_1'\rangle$, $\langle u_2'\rangle$ and $\langle u_3'\rangle$ are the component-wise r.m.s. velocities in three directions. Since the small scales in a turbulent flow are nearly isotropic, we keep $\langle \varepsilon \rangle$ rather than separating it into component-wise quantities. An alternative argument for this choice is that even the viscous dissipation is affected by anisotropic large-scale turbulent fluctuations as long as the time scales $\langle u_i\rangle ^2/\langle \varepsilon _i\rangle$ for different velocity components are comparable, thus there is still no need to divide $\langle \varepsilon \rangle$ into component-wise quantities. Compared with (4.15), the new function $f_{new}$ contains two additional non-dimensional parameters, $\beta _1 = \langle u_1'\rangle /\langle u_2'\rangle$ and $\beta _2 = \langle u_1'\rangle /\langle u_3'\rangle$ that measure the anisotropy of the turbulent flow. For convenience, we can always make $\langle u_1'\rangle$ the r.m.s. velocity in the streamwise direction that usually contains the largest TKE, i.e. $\beta _1 > 1$ and $\beta _2 > 1$. In the present study, the two cross-streamwise directions are symmetric, so only one parameter $\beta = (\beta _1\beta _2)^{0.5}$ is needed to measure the anisotropy. Equation (4.17) then becomes

(4.18)\begin{equation} \left\langle\frac{F_D^T}{F_D^L} \right\rangle- 1 = f_{new}(\beta, I,d_r,Re_p). \end{equation}

Taking $\langle F_{D}^{T}/F_{D}^{L}\rangle - 1 \propto Re_p^{-1}d_r^2$ in Homann et al. (Reference Homann, Bec and Grauer2013) as the reference, it is reasonable to assume $f_{new}$ has power function dependencies on the four non-dimensional parameters, i.e.

(4.19)\begin{equation} \left\langle\frac{F_D^T}{F_D^L}\right\rangle - 1 = K\beta^{a}I^{b}d_r^{c}Re_p^{d}. \end{equation}

We can obtain the coefficients $a$, $b$, $c$, $d$ and $K$ in two ways. The first is to assume the most general form in (4.19) and obtain the optimal values of the coefficients to fit the DNS data, which are tabulated for the nine simulated cases in table 6, via multi-dimensional linear regression. This method gives $a = 1.4$, $b = 0.5$, $c = 0.9$, $d = -0.1$ and $K = 0.06$. The comparison between the predictions of this scaling ($A = \beta ^{1.4}I^{0.5}d_r^{0.9}Re_p^{-0.1}$) and the numerical simulation results is shown in figure 12(a). The second method is to introduce the anisotropic correction on the basis of the Homann et al. (Reference Homann, Bec and Grauer2013) scale, i.e. $b = 0$, $c = 2$, $d = -1$ and $K = 0.18$. The optimal value for $a$ is obtained from linear regression as 1.8. The comparison between the model prediction of the second scaling ($B = \beta ^{1.8}d_r^{2}Re_p^{-1}$) and the numerical data is presented in figure 12(b). Good agreements between the prediction and the numerical results are observed in both plots. Scaling $A$ has a slightly better match to the numerical data than scaling $B$, but its formula is also slightly more complicated. Under both scalings, the anisotropy and the ratio between the particle size and the small turbulent eddies positively affect the relative drag enhancement, while the particle Reynolds number delivers negative impacts. For comparative purposes, the drag enhancement predicted by the nonlinear drag effect, Homann, Bec & Grauer Reference Homann, Bec and Grauer2013's scaling, and the present modified scalings are shown in table 6.

Figure 12.

The drag enhancement coefficients under the revised scalings accounting for the anisotropy of the incoming turbulence flow.

Table 6.

Non-dimensional parameters in the modelling of drag enhancement. The fourth, sixth and eighth columns of the table represent the drag enhancements predicted by the nonlinear drag, the Homann et al. (Reference Homann, Bec and Grauer2013) scaling and the modified scaling. The DNS drag enhancements are shown in the last column.

4.4. The turbulent effect on the drag coefficient of particle arrays

Besides the single-particle cases discussed in the last section, cases with arrays of fixed particles in laminar and turbulent incoming flows are also investigated. The purpose is to understand how the lateral interactions among particles change the particle drag force. In order to reduce the computational cost, only three set-ups, cases B, D and E in table 4 are chosen. Case B and case E have similar turbulence intensities, while case D and case E have the same particle Reynolds number. Under each set-up, the number of particles varies to create particle arrays (see the sketch in figure 13) with different particle–particle relative gap distances, which are quantified by a non-dimensional parameter $r_\rho = d_p/d$, where $d$ is the gap distance between the two nearest particles. The range of $r_\rho$ is chosen between $1/9$ and 1.5, corresponding to a layer of 1 to 36 uniformly distributed particles in region B. Particles in each simulation are placed at the same streamwise location as the corresponding single-particle case. Given a fixed number of 256 grid points in each cross-streamwise direction, further increasing $r_\rho$ may lead to the flow in the throat between two particles being unresolved. Therefore, to avoid sabotaging the reliability of the numerical results, very dense particle cases with $r_\rho > 1.5$ are not covered.

Figure 13.

The drag coefficients of laminar and turbulent flows passing over arrays of particles.

The averaged drag coefficient in both the laminar and selected turbulent cases are shown in figure 13 and tabulated in table 7. With laminar incoming flows, the averaged drag force increases monotonically with $r_\rho$, whereas with turbulent incoming flows, the drag coefficient decreases slightly at the beginning, then increases with $r_\rho$. The overall increased particle drag with $r_\rho$ also partially explains why the settling speed of multiple particles is lower compared with the sedimentation speed of individual particles, as reported in Fornari et al. (Reference Fornari, Picano and Brandt2016a). Besides the hindering effect and the nonlinear drag effect (Yin & Koch Reference Yin and Koch2007; Fornari et al. Reference Fornari, Picano and Brandt2016a), the lateral interactions between particles could also play an important role in the reduced settling speed for multiple particles.

Table 7.

The averaged drag coefficients for laminar and turbulent flow passing over particle arrays.

The increasing drag coefficient with $r_\rho$ is not difficult to understand. As shown by the pressure contours in figure 14, when their distribution becomes denser, particles form contraction passages that force the flow to accelerate and lead to significantly amplified pressure drops between the particle front edge to the zenith location. In addition to this change, as particles approach each other, the low-pressure regions on the back of the particles merge and make the pressure recovery more difficult, which can be seen in figure 14. These two effects result in a significant increase in the pressure drag. Moreover, the flow acceleration in the contracted passages would create larger velocity gradients around particles and enhance the viscous drag. This latter mechanism is confirmed by the comparison of the streamwise velocity profiles along the two lines ${\boldsymbol {e}}_1 = (0,1,0)$ and ${\boldsymbol {e}}_{2} = (-1/\sqrt {2},1/\sqrt {2},0)$ across the particle centre, as shown in figure 15. When turbulence is present, the fluctuation motions of the turbulent eddies may partially counter these mechanisms and result in slight drag reduction. When $r_\rho$ further increases, these inverse effects eventually get overwhelmed, and enhancements of the drag coefficient are again observed. Nevertheless, more in-depth and decisive analyses are certainly needed to further confirm and reveal the subtle physics behind this inverse impact of turbulence on particle drag, which could be pursued in our future studies.

Figure 14.

The comparison between the pressure contours of laminar and turbulent cases passing over arrays of particles with different particle–particle relative gap distances; (a,b) $r_\rho = 1/9$, (c,d) $r_\rho = 1/4$, (e, f) $r_\rho = 1$; (a,c,e) laminar cases, (b,d, f) turbulent cases.

Figure 15.

The distribution of the streamwise velocity along (a) a vertical line and (b) an inclined line across the particle centre.

The results of relative drag enhancement due to turbulence are shown in figure 16(a). In all examined cases, the enhancements of the drag coefficient due to turbulence decrease monotonically with $r_\rho$ and eventually disappear around $r_\rho = 1.5$. This change essentially results from the fact that the drag enhancement due to the contracted passages dominates the contribution due to the turbulence, which is observed from the comparison of the pressure contours in figures 14(e)–14( f) and the comparison of the streamwise velocity distribution around the particle in figure 15. It is also observed that the changes in the drag enhancement coefficient with $r_\rho$ in different cases follow the same pattern.

Figure 16.

The enhancement of the mean drag force as a function of $r_\rho$: (a) relative drag enhancements compared with the laminar cases, (b) an empirical correlation for the relative drag enhancements.

Based on this observation, an empirical correlation to consider the effect of the particle–particle relative gap distance on the drag enhancement due to turbulence is proposed. Defining $r_{\rho 0}$ as the $r_\rho$ with which the particle interaction can be fully ignored (cases with a single particle are roughly this type), the relative drag enhancement compared with its baseline value $\varTheta (r_{\rho 0}) = \langle F_D^T(r_{\rho 0})/F_D^L(r_{\rho 0})\rangle - 1$ approximately satisfies,

(4.20)\begin{equation} \varTheta(r_\rho) = \varTheta(r_{\rho0})\frac{-0.147(r_\rho - r_{\rho0})}{r_\rho - r_{\rho0} + 0.228}. \end{equation}

We emphasize that this correlation is generalized from our simulation results under relatively narrow ranges of parameter settings, i.e. the particle–particle relative gap distances $1/9\leqslant r_\rho \leqslant 1.5$, particle Reynolds numbers $Re_p = 100$ and $240$ and turbulence intensities $0.18< I<0.3$. Applying the correlation for parameters beyond these ranges should proceed with caution.