Investigations of Gas–Particle Two-Phase Flow in Swirling Combustor by the Particle Stokes Numbers

: Gas turbulence modulations and particle dispersions of swirling gas–particle two-phase ﬂow in the combustor is investigated under the large spans of the particle Stokes numbers. To fully consider the preferential concentrations and anisotropic dispersions of a particle, a kinetic frictional stress model coupled with a second-order moment two-phase turbulent model and granular temperature equation is improved. The proposed modeling and simulations are in good agreement with the experimental validations. Results show turbulent modulations and particle dispersions exhibit strongly anisotropic characteristics, keeping a close relationship with ﬂow structure. The axial gas velocity and RMS ﬂuctuation velocity of 45.0- µ m EGP was approximately 5.0 times and 3.0 times greater than 1000.0 µ m Copper particles, and their axial particle velocity was 0.25 times and twice greater than those of 45.0 µ m EGP. The degree of modulation in the axial–radial direction is larger than those of radial–tangential and axial–tangential direction. Particle dispersions are sensitive to particle diameter parameters and intensiﬁed by higher Stokes


Introduction
Gas-particle two-phase turbulent flow in swirling combustor have been employed widely in the fields of industrial processes, e.g., chemical engineering, aerospace engineering, and biotechnology engineering. Gas turbulent flow and particle dispersions play very important roles for the mass and momentum transfer and energy exchange between gas and particle phases. In order to take an effective control for the two-phase hydrodynamics and explore the effects of particle Stokes numbers on gas turbulence modulations and particle dispersions, it is indispensable to obtain better a better understanding of the complicated multiphase mechanism by means of modeling, simulation strategy, and experimental techniques [1][2][3][4][5].
In recent years, computational fluid dynamics (CFD) methods are being utilized substantially compared to hydrodynamic modelling for the single/multicomponent gas-particle two-phase turbulent flow, e.g., the direct numerical simulation (DNS) algorithm [6][7][8][9][10][11], large eddy simulation (LES) SGS model [12][13][14][15][16][17][18], and Reynolds time-averaged Navier-Stokes simulation (RANS) [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. Regarding gas turbulence modulations, Gore et al. (1989) [35] found that turbulence modifications generally have two kind of effects-that is, the augmentation and attenuation activities-using a dimensionless parameter based on the ratio of particle diameter and length of maximum eddy. Elghobashi (1994) [36] established threshold values of the volume fraction of the particle phase to investigate these kind of effects, which are set to α < 10 −6 for one-way couplings, 10 −6 < α < 10 −3 for two-way couplings, and α > 10 −3 for four-way couplings. With respect to the DNS algorithm, Eaton at al. (1994) [7] revealed that fine particles reduce gas turbulence in the boundary layer in developing regions, as well as at the developed flow regions when the mass loading ratio is larger than 0.1. Meanwhile, gas turbulent modulations are highly sensitive to mass loadings rather than particle response, which is the primary order from a group of parameters: particle size, density, swirling number, etc. Moreover, turbulent kinetic energy of gas is determined by intermediate, smaller-scale vortices [8][9][10]. As for LES SGS modeling, Apte et al. (2003) [13] and Oefelein et al. (2007) [14] simulated the swirling particle-laden two-phase turbulent flow in a coaxial-jet combustor based on the Sommerfield's experiment [37]. The hydrodynamic parameters such as particle dispersions, residence time under different diameters, and evolution of coherent structures were predicted. However, the hydrodynamics and turbulent modulation affected by particle Stokes number or particle-density-diameter need to be further investigated. Compared with the abovementioned direct numerical simulation and the large eddy simulation, the Reynolds time-averaged simulation is a somewhat applicable compromised strategy due to its greater economic computation and robust turbulent models, though it fails to capture instantaneous turbulent flow information and two-phase coherent structures. Anagnostopoulos et al. (1992) [20] numerically simulated the effects of particle diameter and the swirling intensities located at corner recirculation regions by a k-ε two-equation turbulent model. Sommerfield et al. (1992) [38] also solved the gas turbulence and tracked particle trajectory using a discrete particle model for swirling gas-particle two-phase turbulence flow. An improved k-ε particle collision model demonstrated that the interparticle length scale had a great effect on gas turbulence by heavy particles [22]. In addition, the turbulence attenuation phenomenon was found by means of the dimensionless Reynolds-particle momentum numbers proposed by Tanaka et al. (2008) [23]. Pakhomov et al. (2015) [24] proposed a Reynolds stresses model for gas-particle turbulence flow and revealed that gas turbulent modulation can be attenuated by fine particles, reaching up to 25% for swirling particle dispersions. As for the particle dispersions of swirling gas-particle two-phase flow in a combustor, the isotropic or homogeneous particle dispersion characteristics were denoted by references [13,14,21,24] and isotropic dispersions by references [4,5,[29][30][31][32][33][34]. In addition, the studies regarding particle dispersions of binary-particle mixture turbulent flow in a liquid-fluidized bed are also carried out. Kim et al. (1983) [25] found that the particle dispersion increased with gas flow rate, liquid surface tension, and viscosity for liquid-solid axial expansion. Asif et al. (1993) [26] proposed a particle dispersion coefficient based on mass balance laws. Galvin et al. (2006) [27] established a generalized model with adjustable parameter of 0.7 to describe the dispersion coefficient for steady-state segregation and dispersion of a binary system of particles in a liquid-fluidized bed. Khan et al. (2020) [28] found that the energy dissipation rate enhances particle dispersion. However, these models failed to describe anisotropic particle dispersion behaviors due to the isotropic hypothesis. Zhou et al. (2017Zhou et al. ( , 2018 [29,30] and   [31] established and developed a set of second-order moment gas-particle turbulent models to successfully reveal the anisotropic dispersions of discrete particle and bubble phases, as well as the hydrodynamics of swirling two-phase flows [32][33][34]. The two-phase anisotropic behaviors, Reynolds stresses redistributions, and mixing and separation characteristics were indicated. To date, the effects of larger spans of particle Stokes numbers-a function of density and diameter-on hydrodynamics have not been investigated. The purpose of this work is to explore the gas turbulence modulations and particle dispersions using the expanded graphite particle (EGP) and the Copper particle (CP) with distinctly different sizes of particle diameters and particle Stokes numbers. A kinetic friction stress model coupled with a second-order moment model and particle temperature equation were proposed to predict the gas turbulence modulations and particle dispersions and analyze the effects of particle density and diameter. Experimental validation using glass particles was performed with acceptable error [37]. Effects of the larger spans of particle density and diameter on gas turbulence modulations and particle dispersions are investigated in detail for the first time.

Modeling and Governing Equations
A second-order moment two-fluid model to describe both gas and particle phases is developed. As for the momentum transfer between gas and particle interactions, a four-way coupling method that fully considers the relationship between gas and particle, and between particle and particle collisions, is adopted. Anisotropic characteristics of gas and particle dispersion were solved by the two-phase Reynolds stresses transport equations and particle-particle collision was modelled by the particle temperature equation. Conservation laws of mass and momentum are satisfied for each phase, no mass and heat transfer occurred between gas and particle phases. The governing equations in relation with correlations and closure transport equations are given below.

Continuity and Momentum Equations
The continuity equations for gas and particle phases (k = gas, particle) are where α k is the concentration of phase k, u ki is the velocity vector of phase k, and ρ k is the density of phase k (k = g,p). The momentum equation for the gas and particle phase are ∂ α g ρ g u gi ∂t ∂ α p ρ p u pi ∂t where g is the gravity acceleration, p is the thermodynamic pressure, β pg is the momentum transfer coefficient between gas and particle, and τ g and τ p are the viscous stress tensors of gas and particle phases.

Momentum Transfer between Gas and Particle Phases
The Huilin-Gidaspow bending drag model that combines the Wen-Yu model with the Ergun equation at solid volume fractions above 0.2, is employed to describe the momentum transfer between gas and particle phases (Lu et al., 2003 [39]). A blending function is used to smooth out the discontinuity between the two equations and is given as follows:

Reynolds Stress Transport Equations of Gas and Particle Phases
The gas Reynolds stress equation is ∂ α g ρ g u gi u gj ∂t + ∂ α g ρ g u gk · u gi u gj ∂x k = D g,ij + P g,ij + Π g,ij − ε g,ij + G g,gp,ij , (9) where the terms on the right-hand side are the diffusion term, shear production term, pressure-strain term, dissipation rate, and gas-particle interaction term, respectively. They are represented by The particle Reynolds stress equation is where the terms on the right-hand side are the diffusion term, shear production term, the pressure-strain term, dissipation rate, and gas-particle interaction term, respectively. They are expressed as follows: The transport equation is used to model and close the anisotropic interactions between gas and particle, shown in Equation (21), where the terms on the right-hand side are the diffusion term, shear production rate term, pressure-strain term, dissipation term, and gas-particle interaction term, respectively.

Particle Temperature Equation Coupled with Kinetic Frictional Stress Model
A kinetic frictional stress model is developed to calculate the particle fluctuating energy-particle temperature θ. It is defined as where τ p is particle stress and Γ p is the conductivity coefficient of granular temperature. The effective coefficient e eff of particle restitution is introduced to consider the Magnus lift force for particle rotation and particle-particle collision (Johnson et al., 1990 [40]).
The net fluctuation energy exchange between the gas and particles (Gidaspow et al., 1994 [2]) and rate of energy dissipation from the transfer of gas phase fluctuation (Koch et al., 1990 [41]) are denoted by the last two terms on the right-hand side in Equation (29).
Particle pressure is used to describe the particle normal forces due to particle-particle interactions. It is improved by the kinetic friction stress p f p , Particle viscosity can be a summation of the shear viscosity from the kinetic theory of granular flows and particle-particle frictional contact in a slow flow regime, it is where µ p f is the definition for frictional viscosity, I 2D is the second variant of the deviator of the rate-of-strain tensor, and ω is angle of internal friction with the value of 28.5.
The translational fluctuation energy dissipation rate is The bulk particle viscosity is The radial distribution function is (44)

Experimental Setting
The experimental combustor consists of the primary jetting and the secondary corner regions (see Figure 1). The diameters of the jetting, corner, and outer chamber are 32.0 mm, 64.0 mm, and 194.0 mm, respectively. The total length of the test-section, oriented vertically with gravity acting in the direction of the flow, is 960.0 mm. Predictions were carried out using the EGP and CP, and glass particle was used for the experimental validation. The primary central flow rate is 9.9 g/s, annular flow rate is 38.3 g/s, inlet Reynolds number is 26,200, and particle loading is 0.034. Swirling number is defined as in Equation (45) Particle Stokes number is a function of particle diameter and density, it is defined as St = τ p /τ f in Equations (46) and (47) and is used to describe the ability particle entrainment carried by the gas phase.
All parameters in the simulation and experimental measurements are shown in Table 1. The effects of larger spans of density, diameter, and St numbers on gas turbulence modulations are discussed in detail.

Numerical Algorithms
The aforementioned governing equations are solved by the finite volume approach based on the staggered grid scheme in the Cartesian coordinate system. Difference strategies using the quadratic upstream interpolation and central difference scheme are set to separate the convective and diffusion terms. Numerical algorithms are employed, namely, the semi-implicit pressure linked correction, the tridiagonal marching way, the line-by-line iteration, and the under-relaxation method. The residual mass for gas and particle phases is set to 1.0 × 10 −4 . The uniform and parabolic distributions for both initial velocity and isotropic and eddy viscosity assumptions of the normal and shear Reynolds stresses are adopted, as well as Johnson's wall boundary condition [40].

Experimental Validations for Gas and Particle Velocity
Glass particle is used to validate the proposed model and source codes. Three sets grid sizes using a coarse system with 6.0 mm by 1.9 mm, medium system with 2.0 mm by 0.48 mm, and fine system with 1.3 mm by 0.33 mm, take advantage of the grid independence test for axial particle velocity (see Figure 2). The medium system is adopted because of the economic computation time and acceptable errors. Swirling flow exhibits strongly anisotropic dispersions in the shear stress zones. Regarding the turbulent structures in the primary and secondary regions, the (x,r) center coordinates of cores are 117.4 mm and 61.1 mm, and 50.4 mm and 73.3 mm, respectively. The center reattachment position is 53.8 mm and 93.0 mm, and the maximum negative velocity at the recirculation region is located at x = 117.4 mm, (as shown in Figure 3 and listed in Table 2), which agreed with Sommerfeld et al. (1991) [37] within 5% error, and the LES simulations by Apte et al. (2003) and Oefelein et al. (2007) [13,14]. Tangential gas and particle velocity at the spanwise regions in comparison to the measurements are given in Figure 4. W-shaped profiles and Rankine vortex structure generated in shear flow zones were slightly overpredicted at x= 155.0 mm, indicating that they correlate well with experiments. It can be explained that the RANS algorithm fails to obtain the instantaneous message and the coherent turbulent structures.

Gas Turbulence Modulations
As shown in Figures 5 and 6, the distributions of gas streamlines and vorticities affected by particles Stokes number are indicated. The labels within contour maps represent the values where negative is the recirculation and positive is the regular flow. Two typical recirculating regions are observed, in which the primary region is located at the nearaxial region and the corner region is located at the near-wall region due to the sudden expansion sections of the combustor, turbulent diffusions, turbophoresis, and inertia force of particles. Since particle Stokes number is a function of particle density and diameter, it takes great effects the evolution of vortices. When St is tiny (St ≤ 1.0), the particles are more easily carried into the gas flow, indicating that they exhibit the excellent followability characteristics for multiphase turbulent flows. Under this condition, the turbulence can be considered as homogeneous flow. When St is far greater than 1, particles hardly respond to the gas fluctuations and only move through the vortex determined by its inertia. Table 3 shows the effect of Stokes numbers on flow behaviors. Ultralight density particles rapidly respond to gas flow in recirculation, resulting in strong entrainment in secondary recirculation zones for strong followability. Shedding vortices from the inlet region rolled up by the boundary layer with vortices from the shear layer experienced the processes of pairing, merging, and breakup. Large, heavy particles will penetrate the reversed flow regions instead of undergoing to the wall-normal region due to the large inertia and lagging velocity with gas inlet velocity. Therefore, distinct vortex structures could not be generated due to the perturbance of particle inertia.     Figure 7 shows the modulations of axial and tangential gas velocities by particle density and diameter with a larger span of St numbers. We can see the axial negative velocities are observed at both central and near-wall regions. Decaying and maximum shifting toward the wall for axial velocity phenomena are disclosed. Both axial and tangential velocities are changed significantly after x = 155.0 mm due to disappearances of the inlet effects; moreover, they are negligible for small particle St numbers even if they are at near-inlet regions. It is noted that turbulence attenuation caused by large, heavy particles is distinctive. Surrounding the inlet regions, lighter and smaller particles are ready to be pushed from the central flows, then moved toward outer regions as a result of centrifugal force. In sharp contrast, higher particle St numbers attentuated modulation due to their larger  Figure 7 shows the modulations of axial and tangential gas velocities by particle density and diameter with a larger span of St numbers. We can see the axial negative velocities are observed at both central and near-wall regions. Decaying and maximum shifting toward the wall for axial velocity phenomena are disclosed. Both axial and tangential velocities are changed significantly after x = 155.0 mm due to disappearances of the inlet effects; moreover, they are negligible for small particle St numbers even if they are at near-inlet regions. It is noted that turbulence attenuation caused by large, heavy particles is distinctive. Surrounding the inlet regions, lighter and smaller particles are ready to be pushed from the central flows, then moved toward outer regions as a result of centrifugal force. In sharp contrast, higher particle St numbers attentuated modulation due to their larger slip velocity and particle inertia. As for the developed flow position, axial velocity is dominantly increased and tangential velocity is slightly reduced even in the case of smaller St numbers. As for particle velocity, large, heavy particles have a greater axial velocity than light, smaller particles because they readily penetrate the central primary flow regions for higher inertia. As seen from the distribution of light particles with diameter of 1000.0 µm and heavy particles with diameter of 45.0 µm, the effect of diameter as a parameter is distinctively in comparison with density at the same order of St numbers. Hence, a adverse trend for tangential particle velocity is demonstrated in Figure 8. The axial gas velocity of 45.0 µm EGP is approximately a maximum of 5.0 times greater than 1000.0 µm Copper particles, and the axial particle velocity of 1000.0 µm Copper particles is approximately a maximum of 4.0 times greater than that of 45.0 µm EGP.   Figure 9 shows the effects of St number on the root-mean-four-sided (RMS) fluctuation velocities of gas. Significantly anisotropic behaviors are found at the shear stress regions, associated with the larger gradient that decreases toward central and jet boundaries. Then, it gradually weakens along the downstream flow direction, caused by the smaller particle mass flux at local positions. The contributions of St numbers are the reinforcement and reduction for gas fluctuation velocity. With the condition of the effect of initial particle inertia, jetting particles firstly went through the recirculation region, were slowed down gradually, ceased, and proceeded to backflow. Meanwhile, the tangential velocity became vigorously, while particles gradually shifted toward wall regions driven by centrifugal force. Lighter and smaller particles with smaller St numbers are favorable for cutting down the fluctuations in the central region due to being regulated and slowed down in the recirculating jet region. However, large, heavy particles with higher St numbers are favorable for intensifying gas turbulence fluctuations. As for the axial particle fluctuation velocity in Figure 10a, small, light particles are larger, vice versa for those of large, heavy particle in the tangential direction because of excellent followability and higher inertia, respectively (see Figure 10b). Particle diameter is a sensitive parameter compared with density and St number for particle fluctuation velocity. As for the gas RMS fluctuation velocity, 45.0 µm EGP are a maximum of 3.0 times greater than 1000.0 µm Copper particles, and 2.0 times greater for the particle RMS fluctuation velocity.

Effects of St Number on Turbulence Kinetic Energy (TKE)
As seen in Figure 11, the turbulence kinetic energy is enhanced by Copper particle at the central reversed regions and occurred at x > 155.0 mm by EGP. Moreover, the modulation degree by EGP are because of rapid energy absorption. Mohanarangam et al. (2007) [42] indicated that the particle movement, slip velocity of particles, and drag force responding to the gas turbulence are always dependent on the particle time constants or particle St numbers on the condition that there is a larger ratio of particle to fluid density when the particle diameter is smaller than the turbulent Kolmogorov scales. As soon as the turbulent Kolmogorov scale is larger than the particle diameter, a slight distortion to turbulent vortex may occur. Thus, even for the same St number consisting of different particle diameters and particle densities, it could regulate the kinetic energy via several pathways. Figure 11. Effects of particle density and diameter on gas turbulence kinetic energy. Figure 12 shows the effects of particle Stokes numbers on the Reynolds shear stresses, which exhibit the significantly anisotropic behaviors. At shear flow regions with a larger velocity gradient, the axial-axial normal stresses are larger than those of axial-radial and radial-tangential correlations. At the (x,r) section of 50.4 mm and 93.0 mm, the maximum value is generated, and it subsequently begins to decrease toward the central and boundary primary jetting regions. Values at the inlet region are larger than at the development flow regions because Copper particles readily penetrate central regions. However, EGPs with negative velocity can quickly respond to the reversed flow. Thus, modulation processes are determined by centrifugal forces, turbulent diffusions, and moment transferring between gases and particles. The degree of modulation in the axial-radial direction is larger than those of radial-tangential and axial-tangential directions. Both enhancement and reduction by smaller EGP and larger Copper particles are obtained, respectively. The reason is that the particle accumulation concentration of higher preference is easier to generate at the feebish vorticity regions. Reynolds shear stresses of larger Copper particles are less than those of smaller EGP particles (see Figure 13).

Effects of St Number on the Tensor Stress Invariants
Stress invariants of the Reynolds stress imply the geometrical and physical information of turbulence velocity fluctuations, which are denoted as I g,x , I g,y , and I g,z , using the components of normal and shear Reynolds stress (see Equations (50)-(52)). Now that invariants are closely related with stress components, it is more reasonable to observe the effects on invariants than those of any individual parameter. I g,x = λ g,x + λ g,y + λ g,z = u g,x u g,x + u g,y u g,y + u g,z u g,z , I g,y = λ g,x λ g,y + λ g,x λ g,z + λ g,z λ g,y = u g,y u g,y u g,z u g,y u g,y u g,z u g,z u g,z + u g,x u g,x u g,z u g,x u g,x u g,z u g,z u g,z + u g,x u g,x u g,y u g,x u g,y u g,y u g,y u g,y I g,z = λ g,x λ g,y λ g,z = u g,x u g,x u g,x u g,y u g,x u g,z u g,y u g,x u g,y u g,y u g,y u g,z u g,z u g,x u g,z u g,y u g,z u g,z .
(52) Figure 14 shows the effects of particle Stokes number on the distributions of I g,x , I g,y , and I g,z . We can see that modulations in every direction are intensified by the increased particle Stokes number. This indicates that turbulent flows, vortex rotations, and vortices evolutions are affected by the different density and diameter of particles, even with the same space and geometry configurations.
(52) Figure 14 shows the effects of particle Stokes number on the distributions of Ig,x, Ig,y, and Ig,z. We can see that modulations in every direction are intensified by the increased particle Stokes number. This indicates that turbulent flows, vortex rotations, and vortices evolutions are affected by the different density and diameter of particles, even with the same space and geometry configurations. In summary, larger particle Stokes number of large-heavy particles more easily resist the action of centrifugal force and penetrate the central flow regions, resulting in the decreasing particle dispersions. Swirling two-phase turbulence is characterized by the complicated and complex chaos status, indicating that mechanisms have not been understood as a result of multiphase diffusions, interactions of gas-particle and particle collisions, turbulent scales, and particle Stokes number variations. The primary parameters related with gas turbulence modulations and particle dispersions have not been obtained clearly.

Concluding Remarks
In this investigation, effects of large ranges of particle Stokes numbers on modulations and dispersions in swirling gas-particle flow was discussed. Gas turbulence modulations and particle dispersions were numerically simulated. Small, light EGP significantly enhanced the axial velocity and RMS fluctuation velocity compared with larger, heavy Copper particles, reaching up to 5.0 times magnitude. However, effects on gas Reynold shear stress were weakened. Particle dispersions exhibit the anisotropic behaviors as those caused by preferential particle accumulations at lower magnitude vortices. Greater understandings for the mechanism of gas turbulence modulations and particle dispersions, and how to identify the dominant factors from a group of Stokes number Figure 14. Effects of particle density and diameter on gas stress tensor invariants: (a) I g,x , (b) I g,y , and (c) I g,z .
In summary, larger particle Stokes number of large-heavy particles more easily resist the action of centrifugal force and penetrate the central flow regions, resulting in the decreasing particle dispersions. Swirling two-phase turbulence is characterized by the complicated and complex chaos status, indicating that mechanisms have not been understood as a result of multiphase diffusions, interactions of gas-particle and particle collisions, turbulent scales, and particle Stokes number variations. The primary parameters related with gas turbulence modulations and particle dispersions have not been obtained clearly.

Concluding Remarks
In this investigation, effects of large ranges of particle Stokes numbers on modulations and dispersions in swirling gas-particle flow was discussed. Gas turbulence modulations and particle dispersions were numerically simulated. Small, light EGP significantly enhanced the axial velocity and RMS fluctuation velocity compared with larger, heavy Copper particles, reaching up to 5.0 times magnitude. However, effects on gas Reynold shear stress were weakened. Particle dispersions exhibit the anisotropic behaviors as those caused by preferential particle accumulations at lower magnitude vortices. Greater understandings for the mechanism of gas turbulence modulations and particle dispersions, and how to identify the dominant factors from a group of Stokes number parameters, particle density or diameter, and particle Reynold numbers, need to be validated in further detail.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to ethical.