Direct Numerical Simulation of Particle-Laden Swirling Flows on Turbulence Modulation

The modulation of turbulence by the laden particles in swirling flows is studied via direct numerical simulation. The statistical characteristics of turbulence modulation are investigated in detail under the effects of different mass loadings as well as Stokes numbers. It is found that the characteristics of turbulence modulation for different Stokes numbers are very similar to each other when the mass loading is light. As the mass loading increases, small particles seem to modulate turbulence more rapidly than large particles. The number concentration or the number flow rate of particles plays an important role in modulation of turbulence. It induces the preferential attenuation of turbulence for small particles in the near field region. Moreover, the trends of modulation of the axial/azimuthal fluctuations, the turbulent kinetic energy, and the Reynolds stress tenor as well as its invariants are similar in the near field region. However, when the turbulence is decayed sufficiently in the downstream region, the inverse turbulence modulation may occur especially for the regions with local intensive accumulation of small particles.


Introduction
The swirling flow occurs in a variety of natural phenomena and industrial applications.In the past decades, enormous researches were devoted to it, for example, the studies of the swirling flames [1], the large scale coherent structures [2,3], the velocity characteristics [4,5], and so forth.In general, swirling flow is always characterized by the remarkable vortical motions, for example, the vortex breakdown (VB) [6], the processing vortex core (PVC) [7], the helical instability [8,9], and the hysteresis [10].All of them are very interesting topics of turbulence research.
With regard to the modulation of turbulence in particle/droplet-laden flows, many researches were carried out too.But most of them are based on the isotropic turbulent flows [11][12][13][14][15]. Thus, it is still not well understood for the turbulence modulation characteristics of the anisotropic particle-laden swirling flows.For example, Anagnostopoulos and Bergeles [16] carried out a numerical study of the effects of droplet size, loading, and swirl intensity on the internal recirculation region of the carrier phase, indicating that the recirculating mass increases with the droplet size and decreases with loading, and the internal recirculation region breaks down into two subregions with high loadings, swirl, and so forth.
In addition, it still lacks good well-accepted criteria for turbulence modulation.An early work on this issue is carried out by Gore and Crowe [17] on the effect of particle size on turbulence modulation.However, it is still in controversy against other criteria.For example, Lucci et al. [18] carried out a numerical study on the indicator of turbulent modulation.With the same particle response time but different mass loadings or volume fractions, the modifications of the turbulent kinetic energy are different.Thus, they concluded that the Stokes number is not a good indicator for turbulence modulation.
In our opinion, the modulation of turbulence may not be determined by a sole factor.At least, it should be influenced by a group of parameters instead.As a result, it is necessary to carry out more investigations on the modulation of turbulence in gas-solid flows.For this purpose, the main objective of this study is to explore the mechanism of turbulence modulation in particleladen swirling flows.It focuses on the modification of statistical characteristics of the carrier phase by the laden particles.A direct numerical simulation of a gas-solid swirling jet is carried out.The turbulence statistics as well as the various influencing factors on turbulence modulation are analyzed, especially for the conditions with different Stokes numbers but the same mass loading.

Governing Equations.
The particle-laden swirling jet is issued through a nozzle into a rectangular tank on the basis of a previous experiment [19].The governing equations of the carrier phase are the time-dependent three-dimensional Navier-Stokes equations for incompressible viscous Newtonian fluids [20][21][22]: where   is the fluid velocity and  is the pressure.The Reynolds number is defined as Re =  0 ⋅ /], where  = 0.4 mm is the diameter of the jet inlet. 0 is the inlet mean velocity, and ] is the kinematic viscosity. , is the feedback force from particle to fluid.In simulation, it is the whole feedback forces from all the particles within a mesh volume to the fluid; that is,  , = ∑  =1  () , , where  is the number of particles in the mesh volume and  ()  , is the feedback force from the th particle to the fluid in the th direction.
To solve the above equations, the flow domain of 30 × 10 × 10 is discretized by 384 × 128 × 128 mesh grids (Figure 1).The finite volume method and the fractional-step projection technique are applied [23].The spatial discretization is of second order precision.An explicit low-storage, third-order Runge-Kutta scheme is used for time integration Density of particle,   (kg/m 3 ) 2 5 0 0 Time step, Δ (-) 0.005 Simulation time,   (-) 100 [24].A direct fast elliptic solver is used to solve the Poisson equation.
Generally, the smallest scale of vortices or the Kolmogorov length scale  is needed to be resolved in direct numerical simulation.However, it is not necessary to meet this requirement exactly since it is too stringent [25].In this study, the mesh spacing is of the same order of the Kolmogorov scale under the central difference scheme; that is, Δ ∼ () [3].Thus, the requirement of spatial resolution for DNS is guaranteed.Moreover, a fine time step of Δ = 0.005 is used to stabilize the simulation, and a total time of  = 100 is simulated (Table 1).In addition, the validation for the numerical methods has been done in a previous study [3].

Motion Equations.
For the dispersed phase, some assumptions are used: (1) the particles are rigid spheres with a uniform diameter   and a uniform density   ; (2) the particle density is far larger than the fluid density (  /  ≫ 1); (3) the particle phase is dilute, and the effect of particleparticle collision can be neglected.With regard to the forces of particles, for example, the drag force, gravity force, Basset force, Magnus force, Saffmann force, virtual mass force, and so forth, the drag force is of the leading order based on the above assumptions [26,27].Moreover, the effect of gravity force can also be neglected when the jets are injected in the horizontal direction.Thus, only the drag force is considered in this simulation.Following Newton's law of motion, the governing equations of particles are as follows: where  , , V , ,   , and  , are the positions, velocities, masses, and drag forces of particles, respectively.The drag force is as follows: is the drag factor [28].The dimensionless form of the motion equation is as follows: where St =   /  = (   2  /(18  ))/(/ 0 ) is the Stokes number (the  and  0 are the characteristic length scale and velocity, resp.).

Boundary and Operating Conditions.
For the inflow condition, the profiles of the azimuthal and axial inflow velocities at the inlet are the same with Figures 1b and 1c of [3], respectively.The velocities on the side walls of the flow domain are set to zero to satisfy the nonslip condition.At the outlet, the nonreflecting boundary condition is applied [29].Additionally, the sketch of simulation setup is shown in Figure 1, and parameters are listed in Table 1 for clarity.
The swirl number  is defined as the rate of the maximum azimuthal inflow velocity to the axial inflow velocity [3].The particles are issued through the nozzle into the flow domain at different flow rates.Three different mass loadings and three different Stokes numbers are simulated (Table 2).Thus, the number flow rates, that is, the number of particles issued in any time step, are calculated according to the mass loadings and the Stokes numbers.

Theoretical Parameters of Analysis.
To illustrate the effects of laden particles on modulation of fluid statistical characteristics, it is necessary to study the axial and azimuthal fluctuations, as well as the correlation between them (the Reynolds stress tensor       ), under the conditions with either the same mass loading or the same Stokes number.The fluctuations    (so-called the r.m.s velocity) and the Reynolds stress tensor are defined, respectively, as follows: where where  is the simulation time.As it is a strongly swirling jet which decays rapidly, the near field range is very short.Hence, in this study, the region of / ≤ 3 is regarded as the near field region, whereas the others are regarded as the downstream regions.In the following sections, three axial locations to the nozzle outlet, that is, / = 1, 2, 3 in the near field region, and an axial location of / = 5 in the downstream are analyzed in detail for case study.Moreover, the effects of the laden particles on modulation of the turbulence fluctuations, including the normal and shear components, are all illustrated above.Suppose that all the turbulent fluctuations       constitute a stress tensor which is related to turbulence fluctuations.It is necessary to study the turbulence modulation on this constituted stress tensor.Thus, this section focuses on the effect of modulation on the invariants of the tensor; that is, where  1 ,  2 , and  3 are the eigenvalues of the tensor       .Since the tensor invariants might indicate important geometrical and physical characteristics of turbulence fluctuations, to study the modulation effect on the invariants may be more general than that of any individual component.Based on Figure 2, it is found that the fluid vortices are rather complicated and fully turbulent.Thus, the following sections are going to study the turbulence modulation via statistical characteristics of fluid motion.By comparison, it is seen that the axial fluctuations of St = 0.5 (the solid lines) are attenuated with increasing mass loadings in the near field region (/ = 1, 2, 3, resp.).It is due to that the small particle can "suck" the turbulent kinetic energy immediately after the nozzle outlet.However, in the downstream region of / = 5, the axial r.m.s velocity of St = 0.5 is inversely increased in the center of the radial direction (|/| ≤ 1) when the mass loading is large (  = 0.252).It is probably interpreted as follows: the attenuation of the fluid kinetic energy is much faster than that of the particles in the central region of |/| ≤ 1.Thus, in the downstream, the kinetic energy is transferred backwardly from the particles to the fluid.When the mass loading is large, the backwardly transferred kinetic energy is of a great amount.Thus, the axial turbulent fluctuation is augmented.It is especially evident for the regions with large local concentration of particles, for example, in the central region of |/| ≤ 1.On the contrary, the r.m.s velocity is still attenuated in other regions.

Axial Fluctuation
For St = 1, the trends are similar to that of St = 0.5.However, for St = 5, the situation becomes a bit changed.The degree of attenuation of r.m.s velocity with increased mass loading is not as evident as before.Recall that with the same mass loading the number flow rates for St = 0.5, St = 1, and St = 5 are different.It seems to be that the degree of modulation of turbulent fluctuations is related to not only the mass loading but also the number flow rates.

With Different Stokes Numbers.
On the other hand, it is necessary to study the effect of particle size (  ∼ St 1/2 ) on modulation of fluid fluctuation, as shown in Figure 3 for analysis.
It is observed that the discrepancies in the r.m.s velocity profiles of different Stokes numbers but the same mass loading are not very evident when the loading is light (Figure 3(a)), although minor differences exist in the central regions at / ≤ 3.However, when the loading is increased (Figures 3(b) and 3(c)), the effects of particle sizes on modulation of turbulent fluctuations become evident.With the same mass loading, the axial r.m.s velocity is augmented with increased Stokes numbers/particle sizes (Figure 3(b)).Moreover, with a large mass loading (Figure 3(c)), the r.m.s velocity is augmented only in the near field region (/ ≤ 3) and in the downstream region (/ = 5) for |/| > 1,  whereas it is relatively attenuated at / = 5 in the central region around the jet axis (|/| < 1).
As aforementioned, the modulations for the same large particle of St = 5 but different mass loadings (or different number flow rates) are similar, whereas they are evidently different for different particle sizes (or different number flow rates) but the same mass loading.Thus, the situation seems to be rather complicated.
In our opinion, it is considered that the modulation should be determined by a group of parameters, for example, at least by particle size (  , and   ∼ St 1/2 ), loading (  ), and number flow rate (  ).These parameters are closely related to each other (remember St ∼  2  ,   ∼  3    ∼     St).Thus, it seems to be true that the particle number flow rate is a very important influencing factor.It can affect the effect of particle size as well as the effect of particle mass loading on turbulent modulation.In other words, the effect of mass loading on modulation is appropriately due to the common contribution of particle number flow rate and particle size.For this way, the situation becomes complicated.Thus, when the number flow rate is low, the effect of mass loading on modulation of turbulent fluctuation is very weak (Figure 3(a)).With the number flow rates increased, the effect of particle size on turbulence modulation becomes evident (Figure 3(b)).Moreover, to reach the same mass loading, the increase in the number flow rate of the small particle is larger than that of the large particle.Thus, the modifications on r.m.s velocity are different.The r.m.s velocity for St = 0.5 is attenuated evidently due to the largest increase in the number flow rate, whereas it remains almost the same for St = 5 due to the least increase in the number flow rate (from Figure 3(a) and Figure 3(b)).In addition, the profile of r.m.s velocity is augmented for St = 0.5 at / = 5 in Figure 3(c) compared to Figures 3(a) and 3(b).Thus, it seems that with increasing number flow rate the small particles attenuate the turbulence fluctuation in the near field (/ ≤ 3) and augment the turbulence fluctuation in the central of the downstream region (/ > 3).In contrast, with the same mass loading, large particle produces the weaker modulation due to the less increase in number flow rate.

Azimuthal Fluctuation. In fact, above analyses on the axial fluctuations are also valid for the azimuthal fluctuations and the turbulent kinetic energy 𝐸
Taking / = 3, for example, Figures 4(a) and 4(b) show the azimuthal r.m.s velocities for different mass loadings (St = 0.5) and different Stokes numbers (  = 0.125), respectively.Figure 4(c) shows the turbulent kinetic energy  for different Stokes numbers (  = 0.252).Similar to Figure 3, it is found that the azimuthal fluctuation is attenuated with increased mass loading in the near field (/ ≤ 3) (Figure 4(a)).Moreover, compared to large particles with the same mass loading, the azimuthal r.m.s velocity is relatively attenuated with the greater increment of the number flow rate for small particles (Figure 4(b)).In addition, the profiles of turbulent kinetic energy in Figure 4(c) for   = 0.252 are also similar to that of Figure 3(c) at / = 3.Small particles attenuate the turbulent kinetic energy in the near filed region, and large particles relatively augment the turbulent kinetic energy under the same mass loading.Thus, the conclusions of the effect of laden particles on modulation of the axial fluctuation can be extended to the azimuthal fluctuation as well as the turbulent kinetic energy.

Reynolds Stress
Tensor       .The components of the Reynolds stress tensor can be categorized into the normal component       and the shear component       ( ̸ = ).The normal component has already been studied above.This section focuses on the shear component, taking   1   3 or     for case study (Figure 5).
For example, Figures 5(a), 5(b), and 5(c) show the shear component     of the Reynolds stress tensor at / = 3 for   = 0.015, 0.125, and 0.252, respectively.It seems that the modifications of the shear component     are very complicated too.For example, when the mass loading is low (Figure 5(a)), the profiles of different Stokes numbers are of similar trends generally.The large peak values occur around |/| = 2.However, as aforementioned, the number flow rate of the small particles (St = 0.5) increases more rapidly than that of the large particles (St = 5) when the mass loading is increased.As a result, the     of St = 0.5 is attenuated more evidently than that of St = 5 (Figure 5(b)).Moreover, for the even large mass loading (St = 5, Figure 5(c)), the     of St = 1 is suppressed greatly too.Hence, it seems that the increase of mass loading suppresses the shear component of the Reynolds stress tensor.To say this more precisely, the effect of mass loading might work at least via the approach of changing the particle number flow rate, which plays an important role in modulation of the fluctuations as well as the correlated fluctuations of turbulence.Thus, the small particles increase more rapidly in the number flow rate than the large particles when the mass loading is increased.It induces the preferential attenuation of the turbulence fluctuations as well as the correlated fluctuations for small particles in the near field region.

Stress Tensor Invariant
Figures 6(a), 6(b), and 6(c) illustrate the invariants for   = 0.015 at / = 2, for   = 0.125 at / = 3, and for   = 0.252 at / = 4, respectively.In Figure 6(a), the invariants  1 for different Stokes numbers are very similar to each other when the mass loading is light.It is also true for the other axial locations and the other two invariants.In Figure 6(b), when the mass loading is increased, it is seen that the second invariant is decreased for the small particles compared to that of the large particles.Once again, the trends are similar for other invariants and axial locations except / = 5.In addition, the third invariant is smaller for the small particles than that for the large particles (Figure 6(c)).
These observations extend the conclusion derived from the above sections; that is, the turbulent fluctuations, no matter what kind of combinations of them, for example,   or  1 ,  2 ,  3 , and so forth, are lower for small particles than for large particles under the same mass loadings in the near fields.

Inverse Modulation.
In Figure 3(c), it is noticed that the turbulent axial fluctuation is augmented for St = 0.5 and St = 1 with the heavy mass loading   = 0.252, which seems to be in contradiction with the results of / ≤ 3 in the near field region.In fact, not only the axial fluctuation, but also the azimuthal fluctuation (Figure 7(a)), the Reynolds stress tensor (e.g.,     in Figure 7(b)), and the invariants (e.g.,  2 in Figure 7(c)) are contrary to the results of / ≤ 3, respectively.It seems that the turbulence fluctuations are decaying and the kinetic energy of the fluid is dissipating in the downstream region.Then, the small particles of large local number concentrations tend to augment the turbulence anomalously.We call it the "inverse" modulation.
In our opinion, the "inverse" modulation is mainly due to the abnormal local accumulation of particles in the central of the downstream region, for example, around |/| < 2 at / = 5.For validation, Figure 8 illustrates the total absolute number of particles within a mesh volume at / = 5.It is found that the small or moderate particles are intensively accumulated around the jet axis (|/| < 2) at / = 5 (Figure 8(a) for St = 0.5 and Figure 8(b) for St = 1).But it is not true for large inertia particles (Figure 8(c) for St = 5).With formation of the local accumulation of particles, the modulation is greatly influenced and changed.Remember that the attenuation of the kinetic energy of the fluid phase is greater than that of the particle phase in the downstream region.Hence, the net kinetic energy is transferred inversely from the particle phase to the fluid phase.With the large number flow rates or concentrations, the quantity of the inversely transferred kinetic energy is obviously large.Then, the turbulence within the local regions of intensive accumulation of particles is augmented considerably.In addition, the mass loading might indicate the order of the overall interphase momentum exchange.However, the spatial distribution of the interphase exchange of the momentum or the kinetic energy might depend to a great extent on the distribution of the particle number concentrations.Thus, although the mass loading is the same, the small particles might cause the inverse augmentation of turbulence through the local accumulation mechanism.

Brief Conclusion
Finally, based on the above analysis, the main findings of the present study are briefly summarized in the following.
(1) In our opinion, the turbulence modulation is determined by a group of parameters, for example, particle size, mass loading, number concentration, local accumulation, interphase energy equilibrium and transfer, and so forth.Thus, there might be no sole universal indicator for turbulence modulation.

3. 1 .
Particles and Vortices.For example, the typical results of fluid vortices and particles are shown in Figure 2. It shows the snapshots of fluid vortices (Figure 2(a)) and particles (Figure 2(b)), respectively, for   = 0.015 and St = 0.5, and the particle-vortices system for   = 0.252 and St = 1.0 (Figure 2(c)).In swirling flows, the vortices are expanded suddenly in the flow container in the radial direction and disturbed disorderly by the particles (Figure 2(a)).The light

( 2 )
Within these influencing factors, the mass loading is of the leading order.The trends of modulation are similar for the light loading.(3)The number flow rate or the concentration plays an important role in modulation.It changes the number of "patterns" of interphase energy transfer.Small particles take modulation more rapidly than large particles as the loading increased, leading to a preferential modulation in the near field.(4) The modulation depends also on the local interphase energy balance.The inverse modulation may take place when the equilibrium of interphase energy transfer is reversed due to the local accumulation of particles.(5) The modulation trends are similar for all of the combinations of turbulence fluctuation.Nomenclature   : Coefficient of drag force of particles : Diameter of the jet inlet (m)   : Particlediameter(m)   : Turbulent kinetic energy (J) : Dragfactorofparticles  , : Components of feedback force from particles to fluids (Pa)  , : Components of drag force of particles (Pa)  1 ,  2 ,  3 : Invariants of the strain tensor   : Mass loading of particles (kg/kg) : N um bero fpa rtic les : e( s )  0 : Mean inflow velocity (m/s)   ,   : Component of fluid velocity (m/s)    ,    : Component of fluid fluctuating velocity (m/s) V , : Component of particle velocity (m/s)  , : Component of particle displacement (m) Δ: Mesh spacing (m) Δ: Time step (s) : Kolmogorov length scale (m)   : Kinetic viscosity of fluids (Pa⋅s) ]: Kinematic viscosity (m 2 /s)   : Fluid density (kg/m 3 )   : Particle density (kg/m 3 )   : Characteristic time of flow (s)   : Characteristic response time of particle (s).

Table 1 :
Parameters used in the numerical simulation.

Table 2 :
The number flow rate for different mass loadings and Stokes numbers.