1 Introduction

Interactions between particles and turbulent flows exist extensively in industrial processes and environmental applications, such as pneumatic transport, spray forming, pollutant dispersion, evaporating milk droplets in spray dryers and pulverized coal combustion^[1]. Due to fundamental importance and practical interest of these problems in nature and engineering, particle flows have received considerable attention in experiments and numerical studies over the last few decades^[2-7].

A crucial concept in the analysis of particle-laden turbulent flows is coupling. Elghobashi^[8] proposed regimes to classify the level of interactions between particles and turbulence. There are one-way, two-way, and four-way coupling forces in terms of particle volume fraction (Φ_p). In the past decades, particle-laden turbulent flows have been performed both experimentally and by means of numerical simulation. It is found that turbulence modulation by the presence of particles is relevant to two primary parameters: the particle Stokes number (St) and mass loading (φ_m). Kulick et al.^[2] conducted an experiment to investigate interactions between turbulent flows and small dense particles in a fully developed downward channel flow in air. They found that the fluid turbulence was attenuated by the presence of particles with the increasing Stokes number, the particle mass loading, and the distance from the wall. Paris^[9] also performed an experiment of particle-laden flows to study the momentum transfer in the two phases. He observed that the moderate particle mass loading decreased the level of the gas-phase turbulent kinetic energy, and regraded the modification of the gas-phase turbulence structure by the particles as a significant mechanism for turbulence attenuation in the particle-laden flow. Direct numerical simulation (DNS) and large eddy simulation (LES) have also been used to examine turbulence modulation in particle-laden channel flows. Pan and Banerjee^[10] numerically investigated the two-way coupled point-particle DNS in an open channel flow. They found that particles with the diameter smaller than the Kolmogorov scale in the near wall region tended to suppress turbulence, whereas larger particles enhance turbulence. Li et al.^[11] also used the DNS to study the effects of particle feedback on the gas phase and particle-particle collisions, and they observed that collisons between particles and fluids played a significant role in suppressing the gas-phase turbulence. The results of the LES combined with the Lagrangian particle tracking in a fully developed gas-solid vertical channel were presented by Dritselis and Vlachos^[12], and they were in good agreement with the DNS results. However, their LES prediction of particle-laden turbulent flows with momentum exchange between two phases was more sensitive to the subgrid-scale turbulence model relative to particle-free flows. Zhao et al.^[13] found that spherical particles can reduce the drag because of a larger bulk fluid velocity at the constant pressure gradient compared with particle-free flows. Then, they also found that particle dissipation will always give rise to a loss of mechanical energy from the fluid-particle system as long as a slip velocity exists. However, for particles of the size comparable to or larger than the Kolmogorov length scale, conclusions are always inconsistent. Shao et al.^[14] and Yu et al.^[15] performed fully resolved simulations of particle-laden turbulent flows in a horizontal channel, and their results showed that the flow friction drag increased first and then reduced, as particle inertia increased.

Although particle-laden channel flows have been investigated widely in the past few decades, the mechanisms of turbulence modulation and their parametric dependence are poorly understood and are wide open for fundamental investigation^[3]. The present study focuses on the ability of particles to modulate the dynamic characteristics of the fluid phase in wall-bounded turbulent flows. We investigate the turbulence modification in channel flows quantified by fluid velocities statistics and characterized by the existence of persistent structures, particularly in the near wall regions, such as low-and high-velocity bands and quasi-streamwise vortices. Different from the above mentioned computational method, this study is implemented with the lattice Boltzmann (LB) method to acquire the numerical solution of Navier-Stokes (NS) equations. The LB method was previously thought to be a very attractive model to simulate complex flows^[16], such as turbulent flows^[17], particle-laden flows^[18], and porous media flows^[19]. Compared with the conventional computational fluid dynamics (CFD), the LB method is an efficient alternative candidate for solving incompressible NS equations^[20]. In fact, the LB method is robust, and it is regarded as a convenient tool for handling boundary conditions in complex flow geometries. In addition, using the LB code can easily realize parallel computing.

In this paper, we use the DNS-LB method to investigate the turbulence and skin friction modulation in channel flows for different Stokes numbers and different mass loadings. In the next section, we describe the physical model of the Eulerian-Lagrangian approach and the numerical methodology, and we present the details on how the LB method can solve the three-dimensional particle-laden turbulent channel flows. Then, in Section 3, we describe the simulation overview of flows and make the numerical simulation for validation. The results are presented in Section 4. Finally, in Section 5, some conclusions are stated.

2 Physical model of Eulerian-Lagrangian approach and numerical methodology 2.1 LB method of particle-laden channel flows

In this study, we consider that fully developed particle-laden flows are bounded by two infinite parallel walls at y=-H and y=H, where y is in the wall-normal direction. We assume that the fluid is incompressible and Newtonian with constant physical properties.

The continuity and momentum equations of the flow driven by a constant pressure gradient can be written as

(1)

(2)

where u_i represents the velocity components in the streamwise (x), wall-normal (y), and spanwise directions (z), respectively. The Kronecker function δ_{1, i} (δ_{i, j}=1 for i=j, and δ_{i, j}=0 for i≠j) is the constant pressure gradient that drives the flow in the streamwise direction. The effect of particles on fluids is represented with f_{p, i}, a two-way momentum coupling force. The equations are non-dimensionalized by the channel half-width H and the friction velocity u_τ. The friction Reynolds number is defined as Re_τ=u_τH/ν, where ν is the kinematic viscosity of the fluid.

The standard lattice Bhatnagar-Gross-Krook (LBGK)^[21] model with only a single relaxation parameter is thought to be the most widely used LB equation model. It is well known that it has second-order accuracy in the time and space if the flow is weakly compressible. In the present study, the LBGK model is used to simulate the turbulent flow.

The three-dimensional particle-laden wall-bounded turbulent flow is numerically solved with the LB method for computing the time evolution of a probability distribution function f(x, u, t). According to the LBGK model including the forcing term^[22], the Boltzmann equation is written as

(3)

and the relaxation time τ and the equilibrium distribution function f_i^eq can be adjusted to recover the incompressible NS equation. For the nineteen-velocity (D3Q19) model, the discrete velocity is defined by

(4)

and the equilibrium distribution function can be expressed as

(5)

where ω_i is the weight coefficient, and c_s is the sound speed of the model. Equation (3) can be divided into two processes: the collision step, i.e., the distribution functions are updated according to the collision rules, and the streaming step, i.e., the evolved distribution functions are moved to new sites according to the D3Q19 discrete velocity model.

The particles feedback forcing term is

(6)

where F is the sum of the particles feedback force and the driving force. The macroscopic density and the velocity can be obtained from the probability distribution function,

(7)

From the LB equation model as shown in Eq.(3), we can use Chapman-Enskog expansion multi-scale technology and transform to the dimensionless form to get the corresponding NS equations, i.e., Eqs.(1) and (2).

2.2 Calculation of particle trajectories

The carrier phase can be simulated with the LB method, while the Lagrangian tracking method is used to calculate the particle dynamics in turbulent flows^[23]. In the current study, the density ratio of the dispersed phase to the carrier phase is large enough, and the particle diameter is smaller than the Kolmogorov length scale. Thus, the most important force acting on the particle is the drag force exerted by the fluid. The governing equations for the location and velocity of each particle can be written as

(8)

(9)

where x_i represents the particle position, and v_i is the particle velocity. u_i is the local fluid velocity at the particle position. Re_p is defined as the particle Reynolds number, τ_p=ρ_pd_p²/(18μ) is the particle response time which stands for the particle inertia. ρ_p and d_p denote the density and diameter of a particle, respectively, and the dynamic viscosity is denoted by μ.

The particle velocity and position are simultaneously calculated with the NS equations of carrier phase, by means of Lagrangian particle tracking coupled with the DNS by the LB method. The velocity and position of each particle are advanced in the time using the second-order Adams-Bashforth approach^[24]. In Eq.(9), the fluid velocities at the particle position are interpolated by the fourth-order Lagrangian interpolation polynomial. After the flow has fully developed to turbulence, the particles are released at random locations within the computational box. The initial particle velocity is set to be equal to that of the fluid at the particle initial position. The particles impinging the smooth wall are treated as perfectly-elastic collisions. To satisfy the periodic boundary conditions, particles are reintroduced into the computational domain when they move outside the domain from the streamwise or spanwise direction. Then, the subsequent new particle position and velocity can be computed from Eqs.(8) and (9), respectively.

The direct particle-particle interaction is disregarded under the assumption that the particle volume fraction is smaller than 10^-3. However, the effects of the coupling between the two phases in momentum equations are considered. The dimensionless formulation of feedback terms in Eq.(2) is as follows^[5]:

(10)

with the sum taken in the computational cell of the volume V over the number of the particles n. In the LB model, since all variables are defined in the lattice nodes, we project the feedback term onto the eight nearest lattice nodes using a volume-weighting method.

3 Numerical simulation 3.1 Simulation setup

A schematic of the physical domain and thousands of particles is shown in Fig. 1. For efficiency, the length of domain in the streamwise direction is selected a little smaller than the classical size of a single-phase turbulent channel flow used by the LES or DNS from the other researchers, e.g., Kim et al.^[25]. In this simulation, the Reynolds number based on the friction velocity and the half channel height equals 180. We use 282×90×141 lattice nodes to discretize the computational domain in the stremwise, wall-normal, and spanwise directions, respectively.

As for the fluid phase, no-slip and no-penetration velocity conditions are imposed on the two parallel walls, and periodic boundary conditions for the physical quantities are used in both the streamwise and spanwise directions. When the velocity field has developed to a statistically stationary situation, we randomly release the particles whose velocities are equal to the fluid velocity at the particle locations. For channel flows, the superscript "+" indicates the non-dimensional quantity by the viscous scales. The non-dimensional velocity is given by u⁺=u/u_τ, and the distance from the wall measured in the viscous length (δ_ν=ν/u_τ) is denoted by y⁺=y/δ_ν.

The program for the simulation is parallelized with 64 threads with the OpenMP command in FORTRAN 90.

3.2 Validation of particle-laden flow

To validate quantitatively the present numerical method and in-house LB code, a particle-laden turbulent channel flow is calculated, and the results are compared with the previous results given by Dritselis and Vlachos^[12]. Table 1 shows the detailed parameters of the two cases. In the verification simulation, the computational box is 2πH, 2H, and πH, and it is discretized using N_x× N_y× N_z=282×90×141 gird points in the streamwise, wall-normal, and spanwise directions, respectively. The grid nodes are distributed uniformly in all directions. While in the study of Dritselis and Vlachos^[12], the DNS of turbulent channel flows based on the finite difference method is performed on the 128×129×128 grid points which are uniform in the homogeneous directions and uneven in the wall-normal direction.

In the two cases, the Reynolds number based on the friction velocity and the half channel height is equal to 180. The density ratio of the particle to the fluid (ρ_p/ρ_f) is 7333. The non-dimensional particle diameter d_p⁺ is 0.7, and the particle response time τ_p⁺(=τ_pu_τ²/ν) is 200, where τ_p=ρ_pd_p²/(18μ). The number of particles (N_p) is 35000, and the particle volume fraction is 2.7×10^-5. Therefore, the effects of particle collisions can be neglected. Figure 2 shows the fluid mean velocity (u⁺) and root-mean-square (RMS) velocity fluctuation (u_{i, RMS}) of the present LB method, in comparison with the DNS results of Ref.[12], which are in good agreement as expected. With the analysis above, we confirm that our strategy and codes are reliable for two-phase turbulent flow predictions.

4 Results and discussion

In this section, we study the turbulence modulation by particles. There are mainly two important parameters, i.e., the particle Stokes number (St=τ_pu_τ²/ν) and the particle mass loading (φ_m). We simulate six different cases listed in Table 2. Cases 1, 2, and 3 are the results by varying the Stokes number (at the constant mass loading), and Cases 2, 4 and 5 stand for different mass loadings at the constant Stokes number.

4.1 Velocity statistics

The mean streamwise velocity profiles of the fluid and particle are shown in Fig. 3. In the presence of particles, the mean fluid velocity has different changes. In Fig. 3(a), the lower inertia particle (St=36) closely follows the carrier fluid in the channel center region, but lags the fluid in the near wall region, and the mean streamwise fluid velocity has a significant increase in the channel center region. The intermediate inertia particle (St=74) follows the fluid in the channel center and the near wall region, but lags the fluid in the logarithmic region, and the mean flow velocity has a slightly increase relative to the particle-free case. However, both the higher inertia particle (St=270) and the relevant carrier fluid are smaller than the particle-free case. As expected, the higher inertia particle has a larger difference from the fluid velocity, and it is slightly lower than the fluid at the center region of channel but ahead of the fluid close to the wall. As shown in Fig. 3(b), with the mass loading increasing from 0.1 to 0.2 and then to 0.4 (at the constant Stokes number), particles and the fluid velocity have a consistent increase. The mean fluid velocity is effectively modulated by the higher mass loading cases (φ_m=0.2 and 0.4) but has barely been affected by the lower mass loading case (φ_m=0.1).

Figure 4 shows the simulated RMS of turbulent velocity fluctuation components in the case of two-way coupling at different particle inertias and different mass loadings relative to the particle-free flow. The turbulence modulation by spherical particles can be easily observed.

In most cases (except Cases 1 and 2), the addition of particles augments the streamwise turbulence intensity in almost the whole channel. However, the wall-normal and spanwise turbulence intensities v_RMS and w_RMS, respectively, are attenuated in all cases. This tendency is consistent with earlier experimental data and numerical results. With the increasing Stokes number, the streamwise turbulence intensity (u_RMS) does not show an evident change relative to the base line of the particle-free flow, but the wall-normal and spanwise turbulence intensities become more and more weak.

And for the increasing particle mass loading, the streamwise turbulence intensity is getting larger, and the peak value of each curve is shifted somewhat further away from the wall than that in the particle-free flow. The increase of the streamwise fluctuation is most obvious in Case 5. The dampenings of the spanwise and wall-normal directions are monotonic with the mass loading and Stokes number.

Figure 5 exhibits the effect of feedback on the Reynolds shear stress. The results clearly show that, compared with the particle-free flow, the Reynolds shear stress is reduced by the presence of particles. Obviously, the particle mass loading has a larger impact on the Reynolds shear stress than the particle inertia. This is in accord with the depression of turbulent spanwise and wall-normal intensities.

To obtain contributions to the Reynolds shear stress from different events in particle-laden flows, and to investigate more information of the carrier phase modulated by the particles on each contribution, we perform a quadrant analysis of the Reynolds shear stress. The results are presented in Fig. 6. We can see that the particle addition has great effects on the quadrants two and four (Q2 and Q4) but has few effects on the quadrants one and three (Q1 and Q3). It means that the particle addition can significantly suppress the sweep and ejection motions occurring near the wall, which is responsible for the reduction of the Reynolds shear stress.

4.2 Turbulent structures

In order to confirm that the particle addition can suppress the ejection and sweep motions, we plot the instantaneous visualizations of the vortices in Fig. 7. We identify the effect of the particles on the vortex structures by the Q-criterion which is named after the second invariant of velocity gradient tensor ▽ u by Hunt et al.^[26]. The criterion is defined as Q=, where W_ij and S_ij are the symmetric and antisymmetric parts of velocity gradient tensors, respectively, given by

(11)

We have chosen a value of Q=0.3 to exhibit the vortex structures of particle-free flow and Case 1 (St=26, φ_m=0.1). It is found that the addition of particles reduces the near wall vortex structures compared with the particle-free flow. This is the evidence of the intensities of sweep, and ejection motions near the wall are reduced by the particles.

Since the near wall vortices are affected by the particles, we also expect the structure of the near wall streaks to be modified. Figure 8 shows the instantaneous streamwise velocity fluctuation contours in the wall region (y⁺=16). As seen, the near wall streaks spacing is increased in the presence of particles. The alternating high-and low-velocity bands in Fig. 8 appear longer and more regular than those in the particle-free flow. It implies that the particles can lead to drag reduction.

4.3 Skin friction

To investigate the influence of particles on the dynamic performance of the two-phase turbulent flows, the skin friction coefficient C_f of turbulent channel flows is analyzed. Fukagata et al.^[27] derived the first theoretical formula for the constituent to the frictional drag in the turbulent channel, pipe, and plane boundary layer flows. They showed the direct relationship between the Reynolds stress and the skin friction coefficient for fully developed turbulent channel flows as follows:

(12)

where Re_b=2U_m^*H/ν is the bulk Reynolds number. In this paper, we extend the theoretical analysis of Fukagata et al.^[27] to a fully developed particle-laden turbulent channel flow and derive the skin friction coefficient of a two-phase flow as shown below.

An analytical expression for the skin friction coefficient C_f is derived, which includes three contributions in a particle-laden turbulent channel flow,

(13)

where is the non-dimensional bulk mean velocity. The equation indicates that the skin friction coefficient in particle-laden channel flows includes three parts: the laminar contribution (C_L), the turbulent contribution (C_T), which is contributed by the average Reynolds stress, and the particles feedback contribution (C_P).

Table 3 shows the contributions of three terms to the skin friction for particle-free and particle-laden flows. The values of these terms are averaged over a sufficiently long period of time. C_f0 is the skin friction coefficient of particle-free flow, which is close to the results of Kim et al.^[25](C_f0=0.00818).

It can be observed that the skin-friction coefficient (C_f/C_f0) is reduced by the presence of particles. And the main reason is the considerable reduction of the turbulent contribution term. For Cases 1, 2, and 3 with the increase in the particle inertia (St =36, 74, and 270), the reduction of the skin friction decreases. As shown in Fig. 9(a), this is owing to the fact that the particles with the larger Stokes number make more contributions to skin friction, but not significant change through the terms (C_L) and (C_T). On the contrary, for the flow with the highest particle mass loading, its skin friction coefficient is the smallest. As shown in Fig. 9(b), in Case 5 (φ_m=0.4, St=74), the turbulent contribution term (C_T) decreases obviously than those in Case 4 (φ_m=0.1, St=74) and Case 2 (φ_m=0.2, St=74), and this is consistent with turbulent intensities which are dramatically suppressed by the presence of the highest mass loading particles as discussed above.

5 Conclusions

We have presented a detailed investigation of turbulence and skin friction modification in particle-laden channel flows. The model has been implemented in the DNS-LB method for the carrier phase, where the dispersed phase is modeled in a Lagrangian way. We simulate two groups of cases by designing different particle Stokes numbers and mass loadings to study the dynamic performance of particle-laden channel flows in terms of mean velocity, turbulent intensities, Reynolds stress, and turbulent structures.

The two-way coupling model in turbulent particle-laden channel flows is validated by comparing the velocity statistics with existing numerical simulation results. This validation confirms that our code is accurate and reliable. Compared with results of particle-free flows, it is found that all cases with presence of particle inertial depress the intensities of the spanwise and wall-normal components of velocity fluctuations and the Reynolds shear stress. Except Case 1, all the mean fluid streamwise velocities have a slight increase by the presence of particles. For a given constant pressure gradient, an increase of bulk mean velocity means the drag reduction. The near wall vortex structures are reduced by the presence of particles, and the streamwise low-and high-speed streaks appear longer and more regular than those of the unladen flow. The skin-friction coefficient is reduced by particles, especially with the particle St=270 and φ_m=0.4, it is reduced by 31.3%.