Direct numerical simulation of particle-laden flow in an open channel at $Re_{\tau }=5186$

Abstract

We perform two-way coupled direct numerical simulation of particle-laden flow in an open channel at a friction Reynolds number ($Re_{\tau }$) of 5186, which exhibits many characteristics of high-Reynolds-number wall-bounded turbulence, such as the distinct separation of scales in the inner and outer layers. Three representative cases, an unladen case and low- and high-Stokes-number particle-laden cases, are performed to investigate the turbulent modification by particles. To this end, we compare several statistical quantities to understand the particle effect on momentum exchange and interphasial energy transfer. The modulation of large-scale motions (LSMs) and very-large-scale motions (VLSMs) are analysed using spectral information, and we find that the LSMs and VLSMs are generally weakened in the inner and outer layers, which is qualitatively different from similar simulations at lower Reynolds numbers ($Re_{\tau } \approx 500$). The spatial structures are investigated with correlation analysis, and inclined VLSMs are observed in the near-wall region, with decreased inclination angles by particles. The particles tend to widen and shorten the spanwise and streamwise extent of coherent structures, respectively. Furthermore, we find that the vorticity vector displays a preferential alignment with the eigenvector corresponding to the intermediate eigenvalue of the strain-rate tensor, independent of the particle Stokes number.

1. Introduction

Wall-bounded particle-laden turbulent flows exist widely in nature and industry, including systems such as wind-blown sand storms in the atmosphere, fluidisation beds in chemical reactors and sediment transport in river and coastal systems. As such, the past decades have seen numerous experimental and numerical studies devoted to particle-laden turbulent flows (Crowe, Troutt & Chung 1996; Balachandar & Eaton 2010; Brandt & Coletti 2022). The complex multiscale interactions between fluid and particles can produce rich and complex physics, such as turbulence augmentation (Squires & Eaton 1990; Elghobashi & Truesdell 1993; Tanaka & Eaton 2008), modifications of near-wall coherent structures (Dritselis & Vlachos 2008; Richter & Sullivan 2014; Lee & Lee 2015), clustering and dispersion of inertial particles (Rouson & Eaton 2001; Bernardini, Pirozzoli & Orlandi 2013; Ahmadi, Sanders & Ghaemi 2021) and multiscale energy redistribution (Poelma & Ooms 2006; Richter 2015; Zheng, Wang & Zhu 2021b). From the perspective of direct numerical simulation (DNS) coupled with point-particle (Balachandar 2009) or particle-resolved methods (Tenneti & Subramaniam 2014) for representing the dispersed phase; however, most, if not all, of the investigations of two-way coupling are restricted to low and moderate Reynolds numbers.

It is well-known that the turbulence modification by inertial particles depends on many factors including the ratio of particle diameter to turbulent length scale, particle Reynolds number, Stokes number and mass loading (Gore & Crowe 1989; Tanaka & Eaton 2008). Lee & Lee (2015) examined the effect of Stokes number on turbulence modification in channel flow ($Re_{\tau }=180$) and found that the turbulence statistics (turbulence intensities, Reynolds stress, viscous dissipation, vorticity fluctuations) are augmented or suppressed by low and high inertial particles in a non-monotonic way, respectively. Richter (2015) investigated the spectral extent of the particle effect on the turbulent energy cascade in plane Couette flow ($Re_{\tau }=125,325,900$), and found that the aggregate behaviour of the turbulent kinetic energy (TKE) is a complex function of wall-normal height, Stokes number and Reynolds number.

In recent years, high-Reynolds-number wall turbulence (without particles) has been investigated extensively (Smits, McKeon & Marusic 2011), with $Re_{\tau }$ ($= u_{\tau } \delta / \nu$, where $u_{\tau }$ is the friction velocity, $\nu$ is the kinematic viscosity and $\delta$ is the boundary layer thickness) up to $O(10^6)$. At sufficiently high $Re_{\tau }$, distinct features emerge in the turbulence statistics which differ significantly from low $Re_{\tau }$. For example, the inner-scaled streamwise turbulence intensity has two distinct peaks, one in the inner layer corresponding to so-called large-scale motions (LSMs) (routinely observed at low Reynolds numbers), whereas a second peak in the outer layer which is closely related to so-called very-large-scale motions (VLSMs) that are unique to high-Reynolds-number, wall-bounded flows (Hutchins & Marusic 2007; Marusic et al. 2010). These distinct, multiscale turbulent motions can produce a wide range of particle-to-fluid timescale ratios varying in the wall-normal direction, and complicate simple descriptions of particle transport in wall turbulence (Fessler, Kulick & Eaton 1994; Rouson & Eaton 2001; Marchioli & Soldati 2002; Bernardini et al. 2013).

Therefore, although many studies have focused on turbulence–particle interactions in low-Reynolds-number wall turbulence, it remains unknown the degree to and mechanisms by which particles modulate turbulence in high-Reynolds-number wall turbulence. Wang & Richter (2020) analysed the corresponding interactions of inertial particles with LSMs and VLSMs in open channel flow ($Re_{\tau }=550$) by spatial filtering, and found that LSMs and VLSMs contribute separately or jointly to the particle dynamics, depending strongly on the Stokes number and wall-normal location, i.e. only in certain regions and for certain particle Stokes numbers, particle behaviour is dictated only by small-scale motions; generally speaking, particles interact with the full spectrum of turbulent motion, having strong implications on large-eddy simulation (LES) strategies.

In the very near-wall region ($z^+ < 15$), Lee & Moser (2019) found in particle-free channel flow at $Re_{\tau }=5200$ that there is an inverse scale transfer from the dominant LSMs to VLSMs, driven by interactions of the streaks and large-scale outer-layer structures; it is one of the driving hypotheses of the present study that this can be disrupted by the presence of a dispersed phase. To this end, Wang & Richter (2019) studied the effect of inertial particles on VLSMs at moderate Reynolds number ($Re_{\tau }=550,950$) in open channel flow, and found that the inertial particles can enhance the VLSMs by both direct TKE transportation and indirect upscale energy transfer from LSMs. This indirect mechanism is consistent with the inner–outer coupling described by Lee & Moser (2019).

Jie et al. (2022) revealed that large-scale particle streaks are formed for relatively high-inertia particles at high $Re_{\tau }$ ($=1000, 2000$), which are distinct from the typical small-scale particle streaks observed in low $Re_{\tau }$ channel flow. To the best of the authors’ knowledge, this is the highest-Reynolds-number particle-laden channel flow case simulated with DNS, presumably due to the rapidly increasing computational cost with higher Reynolds numbers. Liu, Feng & Zheng (2022) experimentally investigated the effects of collision bounces and impact splashing on turbulence statistics in an open-loop wind tunnel at $Re_{\tau }=3950$ with multiple modes of particle introduction (both released from above as well as eroded from below). They found that top-released particles attenuate the Reynolds stress and weaken the VLSMs, whereas bouncing and splashing generate ascending particles that transport the kinetic energy upwards, resulting in more energetic VLSMs. Zhu et al. (2021) experimentally investigated the preferential concentration of inertial particles in a particle-laden turbulent boundary layer at $Re_{\tau } = 5500$ under net sedimentation conditions, and proposed a conceptual model in which a large proportion of particles prefer to accumulate on the back ridges of low-momentum wall-attached structures.

In summary, there are many outstanding questions regarding particle–turbulence interactions at high Reynolds numbers which have not been explored systematically with DNS. In particular, will the turbulence modification of VLSMs at sufficiently high Reynolds number be same as moderate Reynolds number? Will the mechanisms of coherent structure modification remain the same? Furthermore, most LES work on particle-laden flow focus on low and moderate Reynolds numbers, due to the absence of reliable high-Reynolds-number particle-laden wall-bounded turbulence data (Kuerten 2016). To address these issues, we perform two-way coupled DNS of particle-laden flow in the pressure-driven open channel at $Re_{\tau } =5186$, with low and high Stokes numbers, which achieves the same high $Re_{\tau }$ as performed by Lee & Moser (2015) in DNS of single-phase channel flow. This hitherto is the highest $Re_{\tau }$ case in DNS of particle-laden channel flow, and would be informative for developing LES two-way coupled particle-laden turbulence parameterisations.

This paper is organised as follows. First, the numerical methods and simulation parameters are described in § 2. In § 3, we present the DNS results, including comparisons with Lee & Moser (2015) for the unladen case, statistically averaged quantities for the fluid and particle phases, turbulence modification and modulation analysis, particle clustering behaviour and correlation analysis. Finally, conclusions are drawn in § 4.

2. Numerical details

DNSs of the carrier phase are performed for an incompressible Newtonian fluid. For the carrier phase, periodic boundary conditions are employed in both the streamwise ($x$) and spanwise ($y$) directions, whereas a no-slip and shear-free boundary condition are imposed at the wall and free surface, respectively. Such boundary conditions have been implemented by Calmet & Magnaudet (2003), Adrian & Marusic (2012) and Wang & Richter (2019) in simulations of open channel flows. A pseudospectral method is employed in the periodic directions ($x$ and $y$) and second-order finite differences are used for spatial discretisation in wall-normal direction ($z$), where a single-sided stretched grid is used. The solution is advanced in time by a low-storage, third-order Runge–Kutta (RK) scheme. Incompressibility is enforced via the solution of a pressure Poisson equation, which updates the velocity field to guarantee a divergence free field at each RK stage.

The trajectories and dynamics of particles are based on the point-force approximation. Here the particle-to-fluid density ratio $\rho _p/\rho _f \gg 1$, and the particle diameter $d$ is of the order of the local Kolmogorov length scale, corresponding to inner-scaled value of $d^+ = 2.59$, which is similar to the choice by Bernardini (2014) and Motoori, Wong & Goto (2022) in point-particle DNS of particle-laden channel flow at $Re_{\tau }=1000$, i.e. $d^+ = 2.58$. The average particle Reynolds number $Re_p$ for the high- and low-$St$ cases are 1.95 and 0.17, respectively. This is far smaller than the suggested maximum $Re_p \approx 800$ for the Stokes drag correction (Schiller 1933). As a consequence of this and the low volume concentrations (maximum bulk volume fraction $\varPhi _v$ of $O(10^{-3})$) in our investigations, only the Stokes drag force is incorporated in the two-way coupling approach (Balachandar & Eaton 2010), and the particle feedback force is geometrically projected onto the eight surrounding computational nodes. There is evidence that terms such as the lift force in the point-particle equations need adjustments near the wall (Costa, Brandt & Picano 2020), but it remains unclear what is the best path forward because the lift force model depends in unknown ways on many parameters, such as Reynolds number, Stokes number, distance from wall and particle volume fraction (Arcen, Tanière & Oesterlé 2006). Presently, the lift force is neglected and we do not expect this to fundamentally change the inner–outer turbulence interactions that are our focus. The simplistic treatment of two-way coupling is validated by Wang et al. (2019), who made a side-by-side comparison between the current numerical model and the more sophisticated feedback treatment (kernel-based model) by Capecelatro & Desjardins (2013), where minimal differences are found in the two-way coupling regime. In order to highlight the effect of particle response time, the effect of gravitational settling is neglected. The collisions between particles are not taken into account, but purely elastic collisions between particles and the wall/free surface are applied. The elastic rebound at the free surface is mathematically analogous to a no-flux condition on the particle concentration (Wang & Richter 2019).

This in-house DNS code has been utilised to simulate particle-laden flows in open channels by Wang & Richter (2019, 2020) at $Re_{\tau } =550,950$ and in planar Couette flow at $Re_{\tau } = 125, 325$ and $900$ by Richter & Sullivan (2014) and Richter (2015). Grid convergence has been verified and the numerical scheme has been validated in Wang & Richter (2019). Presently, we perform two-way coupled DNS of particle-laden flows in open channel at $Re_{\tau }=5186$ for three different cases: unladen (C0), high $St$ (C1) and low $St$ (C2); these are summarised in table 1. The mass and volume fractions are the same as Wang & Richter (2019), who investigated the modulation mechanism of VLSMs with particles of low and high $St$ in the open channel flow.

Table 1. Summary of simulation cases. The particle relaxation time is $\tau _p = \rho _p d^2/ (18 \rho _f \nu )$, where $d$ is the particle diameter, and the ratio of $d$ to the wall-normal averaged Kolmogorov length scale $\eta _K$ is maintained at $d/ \eta _K \approx 0.42$. The particle Reynolds number remains $O(1)$ or lower. Here $\varPhi _m$ and $\varPhi _v$ are the particle mass and volume concentration, and $N_p$ is the total particle number. The dimensionless Stokes number $St^+$ is based on the viscous time scale, i.e. $\nu / u^2_{\tau }$. $St_K$ represents the dimensionless particle time scale based on the averaged Kolmogorov time scale in the inner layer $(z^+ < 100)$ and outer layer $(z^+ > 100)$. $St_{out}$ represents the dimensionless particle time scale based on the bulk time scale $h/U_b$, where $U_b$ is the bulk mean velocity of the channel. We use $T$ to denote the averaging time for flow and particle statistics.

A careful choice of the computational domain size is important to resolve physically reliable turbulent scales and the associated turbulent statistics. Lozano-Durán & Jiménez (2014) examined the effect of domain size on fluid one-point statistics in single-phase turbulent channel flows, and found that DNS with streamwise and spanwise sizes of $L_x =2 {\rm \pi}h$ and $L_y ={\rm \pi} h$ can reproduce identical one-point statistics as a larger domain size of $L_x =8 {\rm \pi}h$ and $L_y = 3 {\rm \pi}h$ at moderate to high Reynolds numbers ($547 \le Re_{\tau } \le 4179$). It is argued that the interactions between very large structures and well-resolved scales is correctly represented despite the periodic boundary conditions. Furthermore, Del Álamo et al. (2004) and Lozano-Durán & Jiménez (2014) found that the resolved part of the velocity spectrum is not strongly affected by the size of the domain. Wang, Park & Richter (2020b) investigated the effects of domain size on inertial particle one-point statistics in DNS of turbulent open channel flows at moderate Reynolds numbers ($Re_{\tau } =550$ and $950$), which are seeded with two-way coupled particles at low volume concentration ($\varPhi _v \le 1.5 \times 10^{-3}$) for a wide range of particle Stokes numbers ($2.42 \le St^+ \le 908$). They confirmed the same conclusion as Lozano-Durán & Jiménez (2014) for the unladen case by comparing the fluid one-point statistics in the medium domain ($L_x =6 h$, $L_y =2 h$) and larger domain ($L_x =6 {\rm \pi}h$, $L_y =2 {\rm \pi}h$). For the particle-laden case with low ($St^+ =2.41$) and very high Stokes numbers ($St^+ =908$), they found that all of the particle one-point statistics could be well captured by the medium domain size. In the present work, we perform DNS of the unladen flow (C0), to compare particle-laden flow at low ($St^+ =6.0$, C2) and very high Stokes numbers ($St^+ =448$, C1), as listed in table 1. Hence, the computational domain size adopted herein is similar to Wang, Park & Richter (2020b), i.e. $(L_x, L_y, L_z) = (2 {\rm \pi}h, {\rm \pi}h, h)$. Although the domain size is smaller than VLSMs in the streamwise direction, which typically scales as ${\lambda _x} = O(10h)$ (Kim & Adrian 1999), we believe it is large enough for the present purpose. We briefly remark that this domain size is the same as Jie et al. (2022) in DNS of particle-laden channel flow at $Re_{\tau } = 2000$, but smaller than the DNS of single-phase turbulent channel flow at $Re_{\tau } = 5186$ by Lee & Moser (2015), using $(L_x, L_y, L_z) = (8 {\rm \pi}h, 3{\rm \pi} h, h)$. However, the mesh resolution is same as Lee & Moser (2015), i.e. $(\Delta x^+,\Delta y^+, \Delta z^+)=(12.7, 6.4, 1)$ with $N_x \times N_y \times N_z =2560 \times 2560 \times 768$.

3. Results and discussion

In the subsequent sections, the streamwise ($x$), spanwise ($y$) and wall-normal ($z$) velocity components are represented by $u$, $v$ and $w$, respectively. The superscript ($+$) indicates that the velocity components are normalised by the friction velocity $u_{\tau }$, and the length scale is normalised by the viscous length scale $\delta _{\nu }$ ($=\nu / u_{\tau }$, $\nu$ is the kinematic viscosity). The subscripts $p$ and $f$ denote particle and fluid, respectively. The velocity fluctuations are indicated with primes ($\prime$), whereas the statistically (temporally and horizontally) averaged quantities are indicated with bars or overlines ($-$).

3.1. Statistically averaged quantities

Figure 1 shows the mean streamwise velocity profiles ($\bar {u}^+$) and profiles of root-mean-square (r.m.s.) of velocity fluctuations ($u'^+_{rms}$, $v'^{+}_{rms}$ and $w'^{+}_{rms}$), turbulent Reynolds stress ($-\overline {u'w'}^{+}$) and TKE of the carrier phase. For the unladen case, the present DNS results are compared with the DNS results of turbulent channel flow at $Re_{\tau } = 5186$ from Lee & Moser (2015). For $z^+ \le 1000$, the mean streamwise velocity profile in the open channel compares well with the closed channel, and a distinct wake region can be observed in both unladen and laden cases for $z^+ > 1000$. Here $u'^+_{rms}$ compares well in the inner layer (i.e. $z^+ \le 100$), whereas $v'^{+}_{rms}$ and $w'^{+}_{rms}$ compare well in a slightly larger region (i.e. $z^+ \le 500$). We find that $-\overline {u'w'}^{+}$ compares well across the channel height, but the TKE shows similar behaviour as the r.m.s. of velocity fluctuations. The peak values of r.m.s. velocity fluctuations and the corresponding locations matches well with the closed channel. Differences in the outer layer are attributed to the differences in channel versus open channel flow near the centreline/boundary, i.e. blocking effect of the shear-free condition on the free surface (Calmet & Magnaudet 2003).

Figure 1. (a) Mean streamwise velocity profile and (b–f) profiles of r.m.s. of velocity fluctuations: (b) $u'^+_{rms}$, (c) $v'^{+}_{rms}$, (d) $w'^{+}_{rms}$, (e) $\overline {u'w'}^{+}$ and ( f) TKE; ($\circ$, cyan), DNS results from Lee & Moser (2015); ——, present DNS results for the unladen case (C0); – – – –, high-$St$ case (C1); —$\cdot$—, low-$St$ case (C2). Values in all panels are normalised by viscous scales.

For the particle-laden cases, the mean streamwise velocity profiles tend to follow the log-law with a smaller Kármán constant for $30 \le z^+ \le 1000$, and a distinct wake region can also be observed for $z^+ > 1000$, which is noticeably flatter than the unladen case. The r.m.s. velocity fluctuations, turbulent Reynolds stress and TKE are all reduced compared with the unladen case, and this suppression is strongest for the high-$St$ particles. These modifications to first- and second-order statistics are generally consistent with the observations in both experimental (Righetti & Romano 2004) and numerical investigations (Dritselis & Vlachos 2008; Zhao, Andersson & Gillissen 2010; Lee & Lee 2015; Vreman 2015), with some exceptions (particularly in the changes to $u'^+_{rms}$). Meanwhile, the peak locations are the same as for the unladen case.

Figure 2 shows the distributions of the mean particle volume concentration ($\bar {\varPhi }_v$), mean streamwise particle velocity ($\bar {u}_p^+$), particle Reynolds stresses ($-\overline {u'_p w'_p}^+$) and r.m.s. particle streamwise velocity fluctuations ($u'^+_{p,rms}$). In addition to the near-wall preferential accumulation, the particles also exhibit a tendency to accumulate preferentially near the free surface, which is observed by Jie, Andersson & Zhao (2021) in particle-laden channel flow at $Re_{\tau } =600$ but missing in low-$Re_{\tau }$ flows, due to the existence of a quiescent core region characterised by relatively uniform velocity magnitude and weak turbulence levels. This is ultimately a result of turbophoresis due to the gradient of TKE. The near-surface particle accumulation is also observed by Fornari et al. (2016) and Yu et al. (2017) in channel flow laden with finite-size particles ($Re_{\tau } = 180$), but they attributed the particle migration towards the channel centre to particle collision and turbulence–particle interactions, respectively. For $z^+ < 30$, the streamwise particle velocity profile for the low-$St$ case transitions from a logarithmic to a linear behaviour, whereas for the high-$St$ case no linear law is observed; this is also reported by Righetti & Romano (2004). The particle Reynolds stress is important for particle deposition in wall-bounded turbulent flows (Young & Leeming 1997). It is found that $-\overline {u'_p w'_p}^+$ increases with decreasing $\bar {\varPhi }_v$ (see also Zhu, Yu & Shao 2018), and then decreases to zero at the free surface; however, the peak location of $-\overline {u'_p w'_p}^+$ is smaller than the location of lowest $\bar {\varPhi }_v$. The low-$St$ case exhibits larger $-\overline {u'_p w'_p}^+$, which is similar to r.m.s. velocity fluctuations of the carrier phase (see figure 1). R.m.s. particle streamwise velocity fluctuations $u'^+_{p,rms}$ for the low-$St$ case display a similar behaviour as the carrier phase, but the peak location of the high-$St$ case is closer to the wall, due to their more ballistic trajectories and collisions with the wall, which produce larger differences between the fluid and particle phases (i.e. slip velocity).

Figure 2. (a) Distributions of the mean particle volume concentration $\bar {\varPhi }_v$ along the wall-normal direction, scaled by the bulk mean concentration $\varPhi _v$. (b) Mean streamwise particle velocity, $\bar {u}_p^+$. (c) Particle Reynolds stress, $-\overline {u'_p w'_p}^+$. (d) R.m.s. of particle streamwise velocity fluctuations, $u'^+_{p,rms}$. Curves:– – – –, high-$St$ case (C1); —$\cdot$—, low-$St$ case (C2). Values in (b,c,d) are normalised by $u_{\tau }$.

3.2. Momentum exchange and interphasial energy transfer

The slip velocity describes the exchange of momentum between the fluid and particle phases. An accurate characterisation of the slip velocities is essential to the modelling of particle trajectories in particle-laden LES (Fede & Simonin 2006), Reynolds-averaged Navier–Stokes (RANS)-based particle dispersion models (Arcen & Tanière 2009), two-fluid model development (Simonin, Deutsch & Minier 1993), etc. The reaction force, directly related with the slip velocity, governs the particle trajectories and segregation (Marchioli & Soldati 2002), and subsequently modulates the turbulent flow (Tanaka & Eaton 2008; Zhao, Andersson & Gillissen 2013).

The streamwise and wall-normal slip velocity, i.e. $\Delta \bar {u}^+$ and $\Delta \bar {w}^+$, are shown in figure 3. In the streamwise direction, Zhao, Marchioli & Andersson (2012) found that particles tend to lead the fluid ($\Delta \bar {u}^+<0$) for $z^+< 20$, whereas lag behind the fluid ($\Delta \bar {u}^+>0$) for $z^+> 20$. This behaviour is confirmed in our high-$St$ case, however, the low-inertia particles always lag behind the fluid, except very close to the wall ($z^+<1.5$). In the wall-normal direction, Zhao et al. (2012) found that particles tend to lead the fluid ($\Delta \bar {w}^+<0$) for $z^+<50$, whereas lag behind the fluid ($\Delta \bar {w}^+>0$) for $z^+> 50$. Two distinct regions can also be observed in the high-$St$ case, and the transition is located at approximately $z^+=1000$. The low-inertia particles (C2) always tend to lag behind the fluid, similar to $\Delta \bar {u}^+$. The magnitude of the slip velocity increases with particle inertia, which indicates that the adaption of particles to the fluid motion is reduced with increasing inertia (Zhao et al. 2012). Muste & Patel (1997) experimentally found that $\Delta \bar {u}^+$ increases with the particle concentration (while referring to constant particle diameter). In the present DNS, the streamwise slip velocity seems to be associated with $\bar {\varPhi }_v$, i.e. for $z^+>200$, $\Delta \bar {u}^+$ tends to be zero when $\bar {\varPhi }_v$ approaches the bulk mean volume concentration.

Figure 3. (a) Streamwise slip velocity $\Delta \bar {u}^+ = (\bar {u}-\bar {u}_p)/{u_{\tau }}$ and (b) wall-normal slip velocity $\Delta \bar {w}^+ = (\bar {w}-\bar {w}_p)/{u_{\tau }}$. Curves: – – – –, high-$St$ case (C1); —$\cdot$—, low-$St$ case (C2).

To gain further insight, we examine the streamwise momentum balance. The total stress $\tau (z)$ can be given by the integral of the streamwise mean momentum equation from $h$ to $z$,

(3.1)\begin{equation} \tau (z) = \rho \nu \frac{\partial \bar{u}}{\partial z} - \rho \overline{u'w'} + \rho \int_h^z \bar{F}_x(z')\,{\rm d} z' = \rho u^2_{\tau} \left(1-\frac{z}{h} \right), \end{equation}

where $\bar {F}_x$ is the mean horizontal feedback force. The first and second terms on the right-hand side of (3.1) denote the fluid viscous and turbulent Reynolds stresses, and the third term is the particle stress, which represents horizontal momentum transferred in the wall-normal direction by the average motions of particles. Profiles of these stress components are shown in figure 4. The total stress shows a linear profile with a slope of $-1/h$ (Costa, Brandt & Picano 2021; Zheng, Feng & Wang 2021a). The viscous sub-layer is dominant by the viscous stress, above which the turbulent Reynolds stress achieves the maximum value and decrease monotonically up to the free surface, almost linear for the unladen case. The particles reduce the turbulent flux by momentum extraction, and this reduction is more pronounced for the high-$St$ case due to larger slip velocity (figure 3), which has been reported in numerous experimental and numerical investigations (Li et al. 2001; Righetti & Romano 2004; Richter & Sullivan 2014). In the wake region ($z/h>0.2$), the particle stress is dominant in the momentum budget, which makes the mean velocity profile more uniform (see figure 1a).

Figure 4. Profiles of the total (solid), viscous (dot-dot-dashed), turbulent (dashed) and particle (dot-dashed) stresses over the channel height. Different colours refer to different cases: black, unladen case (C0); red, high-$St$ case (C1); blue, low-$St$ case (C2).

The momentum exchange between the particle and fluid phases acts as a direct source/sink in the fluid Reynolds stress and TKE budgets (Richter 2015). Particle sources to the $\overline {u'u'}$, $\overline {w'w'}$, $\overline {u'w'}$ and TKE budget are given by

(3.2a–d)\begin{align} \bar{\varPsi}_{11}=\overline{F_{x}^{\prime} u^{\prime}},\quad \bar{\varPsi}_{22}=\overline{F_{z}^{\prime} w^{\prime}},\quad \bar{\varPsi}_{12}=\overline{F_{x}^{\prime} w^{\prime}+F_{z}^{\prime} u^{\prime}},\quad \bar{\varPsi}=\overline{F_{x}^{\prime} u^{\prime} +F_{y}^{\prime} v^{\prime} + F_{z}^{\prime} w^{\prime}}, \end{align}

where $F_{x}^{\prime }$, $F_{y}^{\prime }$ and $F_{z}^{\prime }$ are fluctuations of the particle feedback force on the carrier phase in the streamwise, spanwise and wall-normal directions, respectively. The particle sources are dependent on the characteristics of particle clusters (Capecelatro, Desjardins & Fox 2018) and also strongly related to the particle inertia (Richter 2015).

Figure 5 shows the particle source terms $\bar {\varPsi }_{11}^+$ ($\overline {u'u'}$ budget), $\bar {\varPsi }_{22}^+$ ($\overline {w'w'}$ budget), $\bar {\varPsi }_{12}^+$ ($\overline {u'w'}$ budget) and $\bar {\varPsi }^+$ (TKE budget). The feedback terms are maximum near the wall where the shear rate is the strongest, which generates the major evident signature of the turbulence modulation by particles. $\bar {\varPsi }_{11}^+$ is always positive (source) for the low-$St$ case, meaning that the streamwise velocity fluctuations are forced by the Stokes drag (Gualtieri et al. 2013). For the high-$St$ case, $\bar {\varPsi }_{11}^+$ is almost always negative (sink) except for $z^+<5$. Here $\bar {\varPsi }^+$ shows similar trend as $\bar {\varPsi }_{11}^+$. This is consistent with the idea that heavy, inertial particles inhibit strong turbulent fluctuations by simply acting as obstacles in the viscous and buffer layers (Richter 2015). We find that $\bar {\varPsi }_{22}^+$ is always negative (sink), meaning that Stokes drag intercepts momentum from the fluid turbulence and depletes the velocity fluctuations in the wall-normal direction. The $St$ dependence of $\bar {\varPsi }_{11}^+$ and $\bar {\varPsi }_{22}^+$ is similar to the DNS results of Wang & Richter (2020). We find that $\bar {\varPsi }_{12}^+$ is always negative (sink) for low $St$, but positive (source) for high $St$, which is actually opposite to $\bar {\varPsi }_{11}^+$. It can be noted that the particle contributions to both the Reynolds stress and TKE budgets tend to vanish for $z^+ \ge 300$. As argued in Richter (2015) and Wang & Richter (2019), however, it is important to recognise that the direct contribution $\varPsi$ is not the only mechanism by which particles can modify turbulence.

Figure 5. Profiles of the particle feedback terms to Reynolds stress budget: (a) particle sources to the $\overline {u'u'}$ budget, $\bar {\varPsi }_{11}^+$; (b) particle sources to the $\overline {w'w'}$ budget, $\bar {\varPsi }_{22}^+$; (c) particle sources to the $\overline {u'w'}$ budget, $\bar {\varPsi }_{12}^+$; (d) particle contributions to the TKE budget, $\bar {\varPsi }^+$. Curves: – – – –, high-$St$ case (C1); —$\cdot$—, low-$St$ case (C2). All terms are normalised by $u_{\tau }^3 / \delta _{\nu }$.

Owing to the slip velocity induced by particle inertia, the particles working on the fluid not only act as a direct source/sink in the Reynolds stress budgets, but also transfer energy between fluid and particles. The imbalance between the energy transferred from the fluid to the particles and the energy received from the particles to the fluid reflects energy dissipation, which may help describe the mechanism of drag reduction in particle-laden flow. The time rate of the work done by the local fluid to a particle $W_{p}$, the work done by a particle on the local fluid $W_{f}$, and the particle dissipation $\epsilon$ are expressed as (Zhao et al. 2013)

(3.3) \begin{equation} \left.\begin{gathered} W_{p} =6 {\rm \pi}\mu a(u_{f,i}-u_{p,i}) u_{p,i}, \\ W_{f} ={-}6 {\rm \pi}\mu a(u_{f,i}-u_{p,i}) u_{f,i}, \\ \epsilon=W_{p}+W_{f} ={-}6 {\rm \pi}\mu a(u_{f,i}-u_{p,i})(u_{f,i}-u_{p,i}), \end{gathered}\right\} \end{equation}

where $a$ is the particle radius, and $u_{p,i}$ and $u_{f,i}$ are the particle velocity and the fluid velocity seen by the particle, respectively.

Figure 6 shows the mean power transferred between fluid and particles. Due to larger magnitudes of slip velocities at high $St$, i.e. $\Delta \bar {u}^+$ and $\Delta \bar {w}^+$ (see figure 3), the intensity of energy transfer and dissipation at high $St$ is much stronger than that at low $St$. For the low-$St$ case, the particles exert work on the local fluid ($\bar {W}_{p}<0$, $\bar {W}_{f}>0$) for $z^+ < 20$, above which the energy transfer mechanism is negligible. For the high-$St$ case, the particles exert work on the local fluid for $z^+ < 500$, whereas receive energy from the fluid ($\bar {W}_{p}>0$, $\bar {W}_{f}<0$) for $z^+ > 20$. The peak locations of $\bar {W}_{p}$ and $\bar {W}_{f}$ are consistent with the streamwise slip velocity (figure 3). The non-zero $\bar {\epsilon }$ denotes extra energy dissipation caused by particles. The particle dissipation is most pronounced in the near-wall region but retain an appreciable level below the wake region for the high-$St$ case. The overall trend of $W_{p}$, $W_{f}$ and $\epsilon$ is qualitatively similar to Zhao et al. (2013), although their analysis is restricted to heavy particles at $Re_{\tau } =180$. Zhao et al. (2013) reported a transient location, i.e. $z^+ =36$, above which that particles receive energy from the fluid ($\bar {W}_{p}>0$, $\bar {W}_{f}<0$). The present DNS indicates that the transient location is not only a function of the Reynolds number, but also dependent on the Stokes number.

Figure 6. Profiles of the mean power transferred between fluid and particles: (a) from the fluid to the particle, $\bar {W}_{p}^+$; (b) from the particle to the fluid, $\bar {W}_{f}^+$; (c) particle dissipation, $\bar {\epsilon }^+$. Curves: – – – –, high-$St$ case (C1); —$\cdot$—, low-$St$ case (C2). The mean power is normalised by $6{\rm \pi} \mu a u_{\tau }^2$.

3.3. Spectral analysis

High-Reynolds-number wall-bounded turbulence is characterised by scale separations between the near-wall and outer-layer turbulence, which is at the heart of all scaling theories (McKeon & Sreenivasan 2007). This separation of scales can be seen in the one-dimensional velocity spectra. Figure 7 shows the wavenumber-premultiplied one-dimensional $u$-spectra $k_y \varPhi _{u'u'} (k_y)/u^2_{\tau }$, where $\varPhi _{u'u'} (k_y) = \langle \hat {u}'(k_y) \hat {u}^{* \prime }(k_y) \rangle$ ($\hat {u}'$ is the Fourier coefficient of $u'$), as functions of spanwise wavelength $\lambda ^+_y$. For the unladen case, the energy spectral density of $u$ in the spanwise direction contains two distinct peaks, i.e. one at $\lambda ^+_y =117$, $z^+=13$, and the other at $\lambda ^+_y =7762$, $z^+=1122$. Such distinct peaks have previously been observed in high-Reynolds-number experimental (Hutchins & Marusic 2007; Marusic et al. 2010; Smits et al. 2011) and DNS investigations (Lee & Moser 2015; Yamamoto & Tsuji 2018). In DNS of particle-laden open channel flows ($Re_{\tau }=550, 950, 2.42 \le St^+ \le 908$), Wang & Richter (2019) found that particles weaken LSMs in the inner layer whereas they enhance VLSMs in the outer layer. At the higher Reynolds number of the present simulations, however, both LSMs and VLSMs are generally weakened in the inner and outer layers, as shown in figure 7. The peaks in the spectra located at $\lambda _{y}^{+} \approx 150$ and $\lambda _{y}^{+} \approx 10{\,}000$ are significantly reduced in magnitude (particularly the VLSM peak), whereas there is a small region of energy enhancement around $\lambda _{y}^{+} \approx 3000$ for the high-$St$ particles. The vertical range containing VLSMs is also reduced by the particles, especially for the high-$St$ case. The weakened strength increases with increased $St$, and the contribution from VLSMs forms a bimodal spanwise spectrum in the outer layer for the high-$St$ case: this feature is muted in the low-$St$ case.

Figure 7. Wavenumber-premultiplied one-dimensional $u$-spectra $k_y \varPhi _{u'u'} (k_y)/u^2_{\tau }$ as functions of spanwise wavelength and wall-normal direction: (a) unladen case (C0), (b) high-$St$ case (C1) and (c) low-$St$ case (C2). The values are scaled by the maximum value in the inner layer of the unladen flow.

In Wang & Richter (2019), it was demonstrated that the particles tended to enhance VLSM energy at $Re_{\tau } = 550$ and $Re_{\tau } = 950$, and that this was a non-monotonic function of Stokes number. The DNS showed that VLSM energy was enhanced at two different $St$, but through two different mechanisms: direct interaction between particles and VLSM structures at high $St$, and indirect, nonlinear interactions between particles, near-wall LSM and VLSM structures at low $St$. In the present simulations, however, figure 7 indicates that, at least for the two $St$ considered presently, the peak LSM and VLSM energy is damped, and that the behaviour observed by Wang & Richter (2019) is somehow not present. Although this discrepancy may be simply due to the fact that the previous Reynolds numbers were only just entering the asymptotic regime for wall-bounded turbulence and therefore the VLSMs are only emergent (Smits et al. 2011), we turn to the spectral signature of energy production to understand the mechanism behind the changes in LSM and VLSM energy observed in figure 7.

In spectral space, the modulation of LSMs and VLSMs by particles is at least partially related to the direct influence on velocity fluctuations, which, in turn, can modify the production of TKE and Reynolds shear stress (see figure 5 and the discussion in Wang & Richter 2019). This is demonstrated in figure 8, where we present the $\overline {u'u'}$ spectral production term $\hat {P}_{11}=-\langle \hat {u}^{\prime } (k_{y}, z) \hat {w}^{\prime *} (k_{y}, z) \rangle \mathrm {d} \bar {u} / \mathrm {d} z$, as well as the $\overline {u'w'}$ spectral production term $\hat {P}_{12}=-\langle \hat {w}^{\prime } (k_{y}, z) \hat {w}^{\prime *} (k_{y}, z) \rangle \mathrm {d} \bar {u} / \mathrm {d} z$, as functions of spanwise wavelength $\lambda ^+_y$. Throughout the wall-normal direction, $\hat {P}_{11}$ is positive whereas $\hat {P}_{12}$ is negative. Comparing the premultiplied $u$-spectrum (figure 7) and $\hat {P}_{11}$-spectrum, we find that they have a similar overall shape. In the inner layer associated with LSMs, the $\overline {u'u'}$ production term $\hat {P}_{11}$ is weakened by particles, but enhanced in the outer layer associated with VLSMs. This means that in the inner layer, the presence of particles reduces the production of $\overline {u'u'}$ compared with the unladen case, whereas in the outer layer they act to promote VLSM energy. Indeed, the peak in $\hat {P}_{11}$ near $\lambda _{y}^{+} \approx 3000$ for the high-$St$ case seems to correspond directly to the local maximum of TKE observed in figure 7(b) – at these wavelengths at the large Stokes numbers, energy is being directly used to promote LSMs at a very narrow wavenumber band. Thus, the particles would seem to act as a direct source of streamwise kinetic energy at the VLSM scale based only on $\hat {P}_{11}$.

Figure 8. Spectral production contribution to the streamwise TKE budget and Reynolds shear stress budget: (a) $\hat {P}_{11}$, (b) $\hat {P}_{12}$, normalised by $u^3_{\tau }/{\delta _v}$. Panels from left to right correspond to the unladen case (C0), high-$St$ case (C1) and low-$St$ case (C2), respectively.

Figure 7, however, showed a considerable decrease in VLSM energy centred around $\lambda _{y}^{+} \approx 10{\,}000$. In figure 5(b), it is observed that the particles simultaneously act as a sink in the Reynolds shear stress budget at the VLSM scale. This is evident in figure 1(e), and indicates a severe reduction of TKE shear production. This reduction of energy production ultimately wins the competition for decreasing the energy contained in the VLSM range, over the direct contribution $\hat {P}_{11}$ from the particles. Figure 5 also shows that the primary effect of $St$ is to alter the magnitude of the impacts on $\hat {P}_{11}$ and $\hat {P}_{12}$ (similar to the energy in figure 7). This is in contrast to Wang & Richter (2019), who observed qualitatively different patterns in $\hat {P}_{11}$ and $\hat {P}_{12}$ for the low- and high-Stokes number cases (see their figure 12). The present results at the higher Reynolds number would seem to suggest that the results of Wang & Richter (2019) may unique to Reynolds numbers near the cusp of being considered asymptotic in the outer layer, but additional simulations at a wider range of $St$ at $Re_{\tau } = 5186$ would be needed to confirm this. It is worth noting that at much larger Reynolds numbers in the atmospheric surface layer ($Re_{\tau } = O(10^{6})$), Wang, Gu & Zheng (2020a) and Liu, He & Zheng (2021) have observed that the effect of airborne sand/dust is to increase the energy contained in the VLSM range as compared with dust-free conditions.

3.4. Instantaneous flow structure and particle behaviour

In order to visualise the turbulent structures associated with local enhancements of particle population, figure 9 shows the instantaneous streamwise velocity fluctuation and the corresponding snapshots of particle locations at four $x$–$y$ slabs, i.e. $z^+ = 13, 270, 1122, 3630$, where $z^+ = 13$ and $z^+ = 1122$ correspond to the peak locations of the premultiplied $u$-spectrum for the unladen case; $z^+ = 270$ corresponds to the outer peak locations of the premultiplied $u$-spectrum for the laden case; $z^+ = 3630$ corresponds to the wake region. Particles accumulating in the inner-flow low-speed streaks can be observed for the high-$St$ case (figure 9b,d), whereas low-inertia particles tend to distribute homogeneously in wall-normal planes (figure 9e). This is similar to many other studies (Pan & Banerjee 1996; Marchioli & Soldati 2002; Wang & Richter 2019).

Figure 9. Instantaneous streamwise velocity fluctuation $u'/u_{\tau }$ for the (a) unladen case (C0), (b) high-$St$ case (C1), (c) low-$St$ case (C2) and the corresponding instantaneous snapshots of particle locations for (d) high-$St$ case (C1), $(e)$ low-$St$ case (C2) in four $x$–$y$ planes (panels from left to right: $z^+ = 13$, $z^+ = 270$, $z^+ = 1122$, $z^+ = 3630$). The particles are coloured by the instantaneous streamwise slip velocity $\Delta u^+ = (u-u_p)/u_{\tau }$, which shares the colourmap of $u'/u_{\tau }$ within $[-0.03,0.03]$. For better visualisations, all the figures are stretched in the $y$-direction.

To characterise the anisotropic clustering of the particles in the ensemble-averaged sense, we employ the angular distribution functions (ADFs) (Gualtieri, Picano & Casciola 2009). Here we consider the two-dimensional angular distribution function, $\operatorname {ADF}(r, \theta )$, defined as

(3.4) \begin{equation} \operatorname{ADF}(r, \theta)=\frac{\displaystyle\sum_{i=1}^{n_{p}} \delta N_{i}(r, \theta)/ (\delta r \delta \theta r n_{p})}{N /(L_{x} L_{y})},\quad 0 \leqslant \theta \leqslant {\rm \pi}/ 2, \end{equation}

where particles are collected from a slab with thickness of $2d$; $\delta N_{i}(r, \theta )$ is the particle number in a sector between $r-\delta r/2$ and $r + \delta r/2$ in the radial direction and $r-\delta \theta /2$ and $r + \delta \theta /2$ in the angular direction from the centre of particle $i$; $\theta =0$ and $\theta = {\rm \pi}/2$ correspond to the spanwise and streamwise directions, respectively. Following Wang & Richter (2020), we set $\delta r = 0.008h$ ($\delta r^+ \approx 40$) and $\delta \theta = 0.025 {\rm \pi}$ to compute $\operatorname {ADF}(r, \theta )$. The mean value is from the average of $n_p$ particles from multiple snapshots in time. Finally, the distribution functions are normalised by the surface average particle number in the $x$–$y$ plane ($N/L_xL_y$ representing a randomly distributed particle number density, where $N$ denotes the mean particle number of the two-dimensional $x$–$y$ slab taken in the wall-normal direction). Periodic boundary conditions are used for particles in the streamwise and spanwise directions.

Figure 10 shows $\operatorname {ADF}(r, \theta )$ in the spanwise ($\theta =0$) and streamwise ($\theta = {\rm \pi}/2$) directions at two typical wall-normal heights: $z^+ = 13$ (inner layer) and $z^+ = 1122$ (outer layer). For the high-$St$ case, the particle density from a reference particle in the streamwise direction is higher than in the spanwise direction. For the low-$St$ case, the particles tend to distribute homogeneously without obvious clustering behaviour, consistent with figure 9.

Figure 10. The (a) spanwise and (b) streamwise ADF of particles in the $x$–$y$ slabs with thickness of $2d$ at two wall-normal heights: $z^+ = 13$ (left), $z^+ = 1122$ (right). Curves: – – – –, high-$St$ case (C1); —$\cdot$—, low-$St$ case (C2).

Figures 11 and 12 show the instantaneous wall-normal velocity fluctuation and the corresponding snapshots of particle locations in the $x$–$z$ and $y$–$z$ planes, respectively. In the near-wall region ($z/h \le 0.2$), we can observe inclined structures that high-speed regions alternate with low-speed regions, which is similar to the observation by Zhang, Hu & Zheng (2018) in the atmospheric boundary layer with suspended dust. In the wake region ($z/h>0.2$), vertical structures can be observed in the laden cases, i.e. regions of positive vertical velocity fluctuations are flanked on either side by regions of negative vertical velocity fluctuations in both the streamwise and spanwise directions. The size of these vertical structures are larger due to the particle–turbulence interaction, compared with the unladen case. The low-inertia particles act as scalar tracers near the free surface, in accordance with the observations in figure 9.

Figure 11. Instantaneous wall-normal velocity fluctuation $w'/u_{\tau }$ for the (a) unladen case (C0), (b) high-$St$ case (C1), (c) low-$St$ case (C2) and the corresponding instantaneous snapshots of particle locations for (d) high-$St$ case (C1), (e) low-$St$ case (C2) in the $x$–$z$ plane. The particles are coloured by the instantaneous wall-normal slip velocity $\Delta w^+ = (w-w_p)/u_{\tau }$, which shares the colourmap of $w'/u_{\tau }$ within $[-0.01,0.01]$.

Figure 12. Instantaneous wall-normal velocity fluctuation $w'/u_{\tau }$ for the (a) unladen case (C0), (b) high-$St$ case (C1), (c) low-$St$ case (C2) and the corresponding instantaneous snapshots of particle locations (right) in the $y$–$z$ plane. The particles are coloured by the instantaneous wall-normal slip velocity $\Delta w^+ = (w-w_p)/u_{\tau }$, which shares the colourmap of $w'/u_{\tau }$ within $[-0.01,0.01]$.

3.5. Correlation analysis

The spatial two-point correlations of streamwise velocity fluctuations are utilised to explore the effect of particles on the average spatial structure of the flow. The streamwise two-point autocorrelations of the streamwise velocity fluctuations at different reference planes are computed as

(3.5)\begin{equation} R_{uu} (\Delta x, z_r)= \frac{\overline{ \langle u'(x+\Delta x,y,z) u'(x,y,z_r) \rangle} }{\langle u_{rms}'(x+\Delta x,y,z) \rangle \langle u_{rms}'(x,y,z_r) \rangle}, \end{equation}

where the brackets denotes an averaging operation over the homogenous horizontal plane.

Figure 13 shows the contour plots of $R_{uu} (\Delta x, z_r)$ for both laden and unladen cases, which are computed at four reference planes, as illustrated in figure 9. The high-correlation region near the wall highlights the inclined VLSMs, as noted in figures 11 and 12, which is similar to the atmospheric surface layer at $Re_{\tau } \sim O(10^6)$ (Wang et al. 2020a). However, the particles tend to decrease the inclination angles of VLSMs, which is opposite to the observations by Wang et al. (2020a). The vertical structures are also observed in the wake region ($z/h>0.2$) for the laden cases. In the inner layer associated with LSMs ($z_r^+ =13$), the LSMs in the unladen case have stronger correlations with the VLSMs of the outer region compared with the laden case, which indicates stronger interactions between LSMs and VLSMs. In the wake region ($z_r^+ =1122$), the particles weaken the streamwise correlation, and the weakened strength is even stronger for the high-$St$ case, consistent with the $u$-spectrum shown in figure 7. However, this behaviour is opposite at $z_r^+ =3630$ for the laden cases, i.e. the streamwise correlation of high-$St$ case is slightly stronger than the low-$St$ case.

Figure 13. Contours of $R_{uu}(\Delta x,z_r)$ for the (a) unladen case (C0), (b) high-$St$ case (C1), (c) low-$St$ case (C2) at four $z_r$ locations $z_r^+ = 13$, $z_r^+ = 270$, $z_r^+ = 1122$ and $z_r^+ = 3630$.

A better understanding of the growth of coherent structures with wall distance can be gained from $R_{uu} (\Delta x, \Delta y)$, computed as

(3.6)\begin{equation} R_{uu} (\Delta x, \Delta y)= \frac{\overline{ \langle u'(x,y) u'(x+\Delta x,y+\Delta y) \rangle}}{\langle u_{rms}'(x,y) \rangle \langle u_{rms}'(x+\Delta x,y+\Delta y) \rangle}. \end{equation}

Figure 14 shows the contour plots of $R_{uu} (\Delta x, \Delta y)$ at four reference planes (same as figure 13). Here only positive correlated regions ($R_{uu} (\Delta x, \Delta y) \ge 0.05$) are shown, where the threshold value of $0.05$ indicates the coherent structure edge according to Hutchins & Marusic (2007). Positive correlation regions are flanked by anti-correlated behaviour (not shown) in the spanwise direction, reflecting the spanwise striping of high- and low-momentum regions as noted in figure 9 (Hutchins & Marusic 2007). For all cases, the spanwise width of coherent structures increases with wall distance, and the growth rate is larger in the log region whereas smaller in the wake region, which is in consistent with the observation by Monty et al. (2007). The spanwise scale in the log region is similar, but the particles widen these scales in the wake region, particularly for the low-$St$ case (see also figure 9). The streamwise size of the VLSMs is largest for the unladen case, and the presence of particles shortens their extent, especially for high $St$, which indicates the weakening effects of particles as noted in figure 7.

Figure 14. Contours of $R_{uu}(\Delta x, \Delta y)$ for the (a) unladen case (C0), (b) high-$St$ case (C1), (c) low-$St$ case (C2) in four $x$–$y$ planes $z^+ = 13$, $z^+ = 270$, $z^+ = 1122$ and $z^+ = 3630$. Here $R_{uu} < 0.05$ are cut off.

3.6. Alignment of the vorticity and the strain-rate tensor

In order to gain a better insight into the dynamics at small scales, we analyse the alignment of the vorticity vector, ${\boldsymbol {\omega }}$ with the eigenvectors of the strain-rate tensor $s_{ij}$, whose associated eigenvalues are $\lambda _1 < \lambda _2 < \lambda _3$, also referred to as the principal strain rates. Due to incompressibility, the trace of $s_{ij}$ must be zero, which results in $\lambda _1 < 0$ and $\lambda _3 > 0$. The eigenvectors associated with $\lambda _1$, $\lambda _2$ and $\lambda _3$ represent the principal compressing, intermediate and stretching stress directions, respectively. The orientation of the vorticity vector ${\boldsymbol {\omega }}$ with the principal eigenvectors is commonly expressed by the cosine of the angle, i.e. $\cos (\boldsymbol {\omega }, \lambda _i) (i=1,2,3)$, which is determined from the inner products of vectors.

Figure 15 shows the probability density functions (p.d.f.s) of $\cos (\boldsymbol {\omega }, \lambda _i)$. The vorticity displays a preferential alignment with the intermediate eigenvector, showing as an increased probability of $| \cos (\omega, \lambda _2) | =1$, for both laden and unladen cases. This is in agreement with previous observations in isotropic turbulence (Ashurst et al. 1987), free and contracted shear flow (Mullin & Dahm 2006; Mugundhan et al. 2020), turbulent boundary layer flow (Elsinga & Marusic 2010) and turbulent channel flow (Lozano-Durán, Holzner & Jiménez 2016, $Re_{\tau }=932$). It should be noted that these references only include the spatially averaged p.d.f.s, which lack the detailed structures along the wall-normal direction. For the unladen case, the p.d.f.s of $\cos (\boldsymbol {\omega }, \lambda _i)$ are almost uniform in the $z$-direction, but for the particle-laden cases, the alignment of ${\boldsymbol {\omega }}$ with $\lambda _2$ intensifies in the near-wall and free-surface regions, showing increased probability of $| \cos (\omega, \lambda _1) | =0$ and $| \cos (\omega, \lambda _3) | =0$. The wall-normal averaged p.d.f.s of $| \cos (\omega, \lambda _2) |$ (see figure 15d) also show this difference, i.e. slightly higher p.d.f.s close to $| \cos (\omega, \lambda _2) | =1$. Furthermore, the vorticity is mostly oriented normal to the compressing eigenvector, and particles slightly intensify this. The particle effects become apparent for the stretching eigenvector, more of which are perpendicular to the vorticity. Overall, the alignment preference seems to be independent of the particle Stokes number.

Figure 15. Probability density functions of the cosine of the angle between the vorticity vector and the three eigenvectors of the strain-rate tensor (panels from left to right: compressing vector, intermediate vector and stretching vector): (a) unladen case (C0), (b) high-$St$ case (C1), (c) low-$St$ case (C2) and (d) wall-normal averaged p.d.f.s. Curves: ——, unladen case; – – – –, high-$St$ case; —$\cdot$—, low-$St$ case.

These findings are essential particularly in the context of LES and the development of subgrid-scale (SGS) model, especially the structure-based SGS models (Meneveau & Katz 2000). One such physical model is the stretched-spiral vortex SGS model originally developed by Misra & Pullin (1997). They introduced a stretched vortex to model the SGS effect for LES in the following way: embedded within each computational cell, there exists a superposition of stretched vortices, each having orientation taken from a p.d.f. It is essentially a local alignment model, which assumes a fraction of the local dominant subgrid structures tend to align with the stretching eigenvector and the remainder align with the intermediate eigenvector. Chung & Pullin (2009) orients the SGS vortices along the stretching eigenvector, and this idea has been applied successfully in LES of wall-bounded turbulent flows: turbulent channel flow (Chung & Pullin 2009; Gao, Cheng & Samtaney 2020), turbulent boundary layer flow (Inoue & Pullin 2011) and airfoil flow (Gao et al. 2019). The alignment of the vortical structures in particle-laden flows is, in a sense, a validation of the assumption of orientation along stretching eigenvectors in the stretched-spiral vortex SGS model. Our DNS results in the particle-laden open channel flow may serve as a guide to other SGS models that rely on the vortex alignment in LES of particle-laden flows.

4. Conclusions

We have presented results from two-way coupled DNS of particle-laden flow in an open channel at $Re_{\tau } = 5186$. Three cases, including an unladen case and particles of low and high inertia, were performed with the Eulerian–Lagrangian point-particle approach. The results are compared with the DNS of channel flow at the same $Re_{\tau }$ by Lee & Moser (2015), and we find that the statistically averaged quantities match well below the wake region, and the differences can be attributed to the applied shear-free boundary condition at the free surface. The r.m.s. velocity fluctuations, turbulent Reynolds stress and TKE are suppressed by the particles, and even stronger for the high-inertia particles.

The effects of particles on the momentum exchange are evaluated with the slip velocity. The low-inertia particles tend to lead the fluid in both streamwise and wall-normal directions. The high-inertia particles tend to lead the fluid in the streamwise direction for $z^+ > 20$, whereas lag behind the fluid in the wall-normal direction for $z^+ < 20$, respectively. In the wall-normal direction, the high-inertia particles show opposite behaviour. From the perspective of streamwise momentum balance, the particles reduce the turbulent flux by momentum extraction and this reduction is more pronounced for the high-$St$ case. In the wake region ($z^+>1000$), the particle stress is dominant in the momentum budget. Furthermore, we examined the particle feedback terms to Reynolds stress and TKE budgets. The particle effects are strongest in the near-wall region, and almost vanish for $z^+ > 300$. The low-inertia particles act as source to the $\overline {u'u'}$ and TKE budget, and a sink to the $\overline {u'w'}$ budget, which is opposite for the high-$St$ case. For the $\overline {w'w'}$ budget, the particles always act as sink. The momentum exchange between the particle and fluid phases also induces energy transfer and extra energy dissipation. In the near-wall region, the particles exert work on the local fluid, whereas above some location, the particles receive energy from the fluid. This location is affected by both Reynolds number and Stokes number. The imbalance of energy transferred between fluid and particles causes extra energy dissipation, which is more significant for the high-$St$ case. We furthermore examine the preferential accumulation patterns of particles in the inner and outer layers, and find that the high-inertia particles are clustering in the low-speed streaks, but this behaviour is missing for the low-inertia particles.

Based on the spectral analysis of $u$-spectra, we find the LSMs and VLSMs are generally weakened in the inner and outer layers. The magnitude and vertical range of VLSM are reduced, especially for the high-$St$ case. This is different from the moderate $Re_{\tau }$ cases ($Re_{\tau } =550, 950$) which show VLSM enhancement with particles, as investigated by Wang & Richter (2019). Correspondingly, we find the qualitatively different behaviour of particle effect on the spectral production of TKE and Reynolds shear stress, as compared with the moderate $Re_{\tau }$ cases by Wang & Richter (2019). The $\overline {u'u'}$ spectral production term, $\hat {P}_{11}$ is positive, which is weakened by particles in the inner layer, but enhanced in the outer layer with increasing $St$; the $\overline {u'w'}$ spectral production term, $\hat {P}_{12}$ is negative. The reason may lie in the fact that non-sufficient high $Re_{\tau }$ just enters the asymptotic regime for wall-bounded turbulence and the VLSMs are only emergent.

The effect of particles on the average spatial structure of the flow are analysed with spatial two-point correlations. In the near-wall region, inclined VLSMs are observed, and the particles tend to decrease the inclination angles. In the inner layer, the unladen case shows stronger interactions of LSMs and VLSMs, whereas the particles weaken the streamwise correlation in the wake region. The spanwise width of coherent structures is growing with wall distance, and the particles tend to widen these scales, especially for the low-$St$ case. However, the streamwise extent of the VLSMs are decreased with increasing $St$.

The particle effects on the angle between the vorticity and the eigenvectors of the rate-of-strain tensor is also studied. The results show that the vorticity, ${\boldsymbol {\omega }}$, tends to align with the second eigenvector of the rate-of-strain tensor, $\lambda _2$, which is slightly intensified by particles in the near-wall and free-surface regions, independent of the particle Stokes number. These findings can help further development of SGS models for LES of particle-laden flows.

Finally, we recall the scientific questions posed at the beginning of this paper, namely whether the turbulence modification of VLSMs will persist at high Reynolds number, and whether the mechanisms which modify them are similar to those found in Wang & Richter (2019). Substantial modification of VLSMs (and turbulence statistics more generally) was observed in the present study, and evidence of similar behaviour is being observed in the atmospheric surface layer (Liu et al. 2021). The complexity of the nonlinear coupling of inertial particles with high-$Re$ wall turbulence and the high cost of the simulations, however, precludes a full description at this stage. At such high Reynolds numbers, a wide range of turbulence time and length scales exist from which multiple particle Stokes numbers can be defined. This confounds a clear description of whether the qualitative differences seen presently compared with the lower Reynolds numbers of Wang & Richter (2019) are due simply to a strong Reynolds number dependence, or if the relevant Stokes numbers have not been properly matched between the two (e.g. Jie et al. 2022). Future work would benefit from a more complete picture of how these modification statistics behave across a denser range of $St$ and $Re$, in addition to the dependence on mass fraction.