# DSMC Prediction of Particle Behavior in Gas-Particle Two-Phase Impinging Streams

Authors: Min Du; Changsui Zhao; Bin Zhou; Yingli Hao
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Devices with impinging streams have been employed in various fields of chemical engineering, as a means of intensifying heat and mass transfer processes. The particle behavior in gas-particle two-phase impinging streams (GPISs), which is of essential importance for the research of transfer processes, was simulated by an Eulerian-Lagrangian approach in this paper. Collisional interaction of particles was taken into account by means of a modified direct simulation Monte Carlo (DSMC) method based on a Lagrangian approach and the modified Nanbu method. A quantitative agreement was obtained between the predicted results and the experimental data in the literature. The particle motion behavior and the distributions of particle concentration and particle collision positions were presented reasonably. The results indicate that the particle distribution in GPIS can be divided into three zones: particle-collision zone, particle-jetting zone, and particle-scattering zone. Particle collisions occur mainly in the particle-collision zone, which obviously results in a few particles penetrating into the opposite stream. The interparticle collision rate and the particle concentration reach their maximum values in the particle-collision zone, respectively. The maximum value of the particle concentration increases with the increasing inlet particle concentration according to a logarithmic function. The interparticle collision rate is directly proportional to the square of local particle concentration.


---

## Body

## 1. Introduction

The earliest emergence of impinging streams (ISs) can be traced back to the development and application of the Koppers-Totzek gasifier in 1953, and its scientific concept was proposed by Elperin [1] in the early 1960s. In an original gas-particle two-phase IS (GPIS), two gas-particle streams impinge against each other at high velocity. Particles from one stream penetrate into the other, which increases the relative velocity between the particles and the gas and prolongs the mean residence time of particles. This phenomenon yields the effect of intensifying the interphase heat and mass transfer in particulate systems [2]. Therefore, the GPIS devices have a broad application in many industrial processes, such as coal gasification, combustion, and drying [3]. The particle behavior in GPIS is of essential importance for research of the transfer processes; however, it is difficult to study by experimental means. This study will focus on the particle behavior in GPIS, such as particle motion and particle collision, from a viewpoint of mathematical modeling.There are two classical models in simulation of gas-particle flow, the Eulerian-Eulerian two-fluid model which treats particles as a continuous phase and the Eulerian-Lagrangian particle trajectory model which tracks individual particles. For the former model, it is difficult to give different characteristics of particle motion, interparticle collisions, and interactions between gas and particles in different zones in GPIS, because this model is based on all kinds of hypotheses. Therefore, the particle trajectory model, which needs fewer assumptions, is used in this study to calculate the gas turbulence and the motion of individual particle. On the other hand, an important character in GPIS is the severe interparticle interaction which plays a significant role in the characteristics of the particle motion and concentration distribution. Therefore, the effect of collisions between particles should be taken into consideration in addition to the particle-fluid interaction in the numerical study of GPIS.Many models have been reported for GPIS. In earlier studies, Kitron et al. [4] applied the particle Boltzmann transport equation [5, 6] including interparticle collisions to the study of GPIS. The Monte Carlo method was then used to solve this equation to obtain the velocity distribution of particle phase. However, the solution process of this method was too complex to be applied to the engineering practice. Guo et al. [7] and Ni et al. [8] used a Markov chain stochastic model to predict the residence time distribution of gas and particle in an opposed multiburner (OMB) gasifier. Other researchers [9–13] proposed various single-particle dynamics models to analyze the motion behavior of single particles in GPIS. These models neglected the interaction between particles and cannot give the particle collision effect.The direct simulation Monte Carlo (DSMC) method, based on a Lagrangian approach developed in rarefied gas dynamics for handling collisions of a large number of gas molecules [14, 15], is an effective method of dealing with particle motion and collisions in dense gas-particle two-phase flow [16, 17]. The authors of this paper first applied the DSMC method, called as the traditional DSMC method in this paper, to GPIS to deal with the interparticle collision with acceptable computational cost [18]. Li et al. [19] then established a 3D model of the impinging zone of an OMB gasifier with DSMC method. This model reveals the concentration and the mean velocity profiles of particles in the impinging zone. Nevertheless, the particle flow has its own characteristics in GPIS, such as nonuniform distributions of particle concentration and particle collision positions, and a period of short time required for a particle to pass through the impinging zone. According to these characteristics of multiple spatial scales and multiple time scales, the DSMC method is modified to be more suitable to the numerical study of GPIS as in our earlier work [20]. This modified DSMC method is not restricted to flow field cells in the calculation of particle motion and collisions, and the particle time step in this method is adaptive and calculated using local particle parameters, such as the particle concentration and the velocities of particles nearby. In this work, this modified DSMC method is validated quantitatively using experimental results in the literature. Then, this method is applied to the numerical study of the particle behavior in GPIS. The distributions of particle concentration and particle collision positions are analyzed based on the calculation results. More specifically, the effect of the inlet particle mass flow rate will also be given. These results lay the foundations for further heat and mass transfer study in GPIS.


## 2. Mathematical Model

Real processes in practical engineering GPIS are very complex generally, which include multiphase flow, heat and mass transfer, and chemical reactions sometimes. To simplify the issue, the following assumptions are made in the present modeling:(1) heat and mass transfer is ignored in flow process; (2) no chemical reactions are considered; and (3) particles are treated as rigid hard spheres with the same diameter.


### 2.1. Gas-Particle Two-Phase Flow

An Eulerian-Lagrangian approach is adopted to describe the gas-particle two-phase flow in GPIS as mentioned above. The interactions between gas phase and particle phase are determined by means of Newtonian third law. Therefore, four-way coupling between discrete phase and continuous phase is carried out. The key features of the mathematical model are briefly descried as follows.


#### 2.1.1. Governing Equations for Gas Phase

The gas flow field in GPIS is generally a turbulent flow field with backflows, and it is calculated using the realizablek-ε model [21] in this simulation. In this model, the conservation equations for mass and momentum can be written as follows:
(1)∂ρg∂t+∇·(ρgu)=0,(2)∂(ρgu)∂t+∇·(ρguu)=-∇p+∇·(τ)+ρgg+Fg,
where Fg is the external body forces that arise from interaction with the particle phase, which results in the momentum exchanges between the gas and particle phases. And Fg is calculated from the reverse action of particles. To close (2), the modeled transport equations for the turbulence kinetic energy, k, and its dissipation rate, ε, are
(3)∂(ρgk)∂t+∇·(ρgku)=∇·[(μ+μtσk)∇k]+Gk+Gb-ρgε-YM,∂(ρgε)∂t+∇·(ρgεu)=∇·[(μ+μtσε)∇ε]+ρgC1Sε-ρgC2ε2k+νε+C1εC3εGbεk,
where C1=max[0.43,η/(η+5)]; η=Sk/ε;μt=ρgCμk2/ε;S=(2SijSij)1/2; and Gk=μtS2 represents the generation of turbulence kinetic energy due to the mean velocity gradients. Gb is the generation of turbulence kinetic energy due to buoyancy. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. The model constants are C1ε=1.44,  C2=1.9,  σk=1.0,  and  σε=1.2 [22].


#### 2.1.2. Particle Motion Equations

In GPIS, the spherical particles are entrained by gas flow. The translational motion of a particle in the gas phase is governed by Newton’s second law of motion and can be written as(4)mpdvpdt=Fd+F,
where F is an additional force term used in the simulation, including virtual mass force, pressure gradient force, and Saffman’s lift force. In addition, Fd represents the drag force of gas phase acting on the particle which is the most important force acting on the particle calculated using the following equation:
(5)Fd=0.125πDpμCdRe(u-v).
The drag coefficient, Cd, is determined by the following equation:
(6)Cd={24Re(1+16Re2/3)Re≤1000,0.424Re>1000,
where Re is the relative Reynolds number, which is defined as
(7)Re≡ρgDp|u-v|μ.


### 2.2. Interparticle Collision

The modified DSMC method, based on a gridless approach and local time stepping with automatic adaption, is applied for modeling the interparticle collisions in GPIS here. The basic ideas of this method are as follows:(1)
real particles are replaced by smaller number of sampled particles. And the trajectories of sampled particles are calculated. In this way, the computer memory and computation time can be saved significantly. This item is just the same as the traditional DSMC method;(2)
sampled particle motion is decoupled into the movement and collision processes as the traditional DSMC method. Sampled particle movement obeys the single particle motion model, and the collision process follows the interparticle collision dynamics. In the modified DSMC method, the particle time step, that is, the time required for a particle in each movement, is adaptive, which is different from the traditional DSMC method;(3)
collision pairs are found through the theory of collision probability instead of using trajectories as the traditional DSMC method. Nevertheless, the calculation of the collision probability between particles and the searching of collision pairs is based on the spherical region generated from the local particle parameters rather than the flow field cells.To update the particle time step, the local collision mean-free path, that is, the average distance a particle travels between collisions with other moving particles, is calculated firstly using the following equation:(8)L=vpυ,
where vp is the velocity magnitude of the particle being tracked. υ is the collision frequency, that is, the number of collisions occurred in unit time, which is obtained according to the local distribution of particles. To satisfy the “principle of uncoupling” in the modified DSMC method, the particle time step Δtp is given by
(9)Δtp=min[L3vx,L3vy,L3vz],
where vx, vy, and vz are components of particle velocity vp. Δtg is the time step for gas phase.In the modified DSMC method, the collision pairs are found through the spherical region, that is, searching scope, with the particle tracked being considered as center andr as its radius given below
(10)r=max[vpΔtp,vr,maxΔtp],
where vr,max is the maximum value of relative velocity between the particle being tracked and the sampled particles nearby. Then, the probability of collision between particle i being tracked and particle j during a time step Δtp is calculated by
(11)Pij=0.25π(Di+Dj)2|G(0)|np,jΔtpV,
where Di and Dj are the diameters of particle i and particle j, respectively. G(0) is the relative velocity vector between particles i and j. np,j is the real particle number represented by the sampled particle j. V is the volume of the searching scope and equal to 4πr3/3.The modified Nanbu method [23] is used to search a candidate collision partner j, which probably collide with particle i in the searching scope during a time step Δtp, and decide whether the collision occurs. First, a random number R with a uniform distribution from zero to unity is extracted from a generator. The candidate collision partner j is selected according to
(12)j=int[R×N]+1,j≠i,
where int[R×N] is defined as the integer part of R×N and N is the total number of all sampled particles in the searching scope. If the relation
(13)R>jN-Pij
is satisfied, particle i would collide with particle j during the time step Δtp, and the velocities of particles i and j are replaced with the postcollision velocities, but without changing their positions. The post-velocities of the two particles are as follows:when(n·G(0)/|Gct(0)|)<1/(ξ(1+e))(14)vi=vi(0)-(n+ξt)(n·G(0))(1+e)mjmi+mj,vj=vj(0)+(n+ξt)(n·G(0))(1+e)mimi+mj,when(n·G(0)/|Gct(0)|)≥1/(ξ(1+e))(15)vi=vi(0)-[(1+e)(n·G(0))n+|Gct(0)|t]mjmi+mj,vj=vj(0)+[(1+e)(n·G(0))n+|Gct(0)|t]mimi+mj,where n is the normal unit vector directed from particle i to particle j at the moment on contact. t is the unit vector in the tangential direction. Gct(0) is the tangential component of the relative velocity.


### 2.3. Numerical Solution

The numerical solutions of the gas-particle two-phase flow, including the governing equations for gas continuous phase, motion equations for sampled particles, and the momentum exchanges between two phases, are obtained using the CFD code FLUENT. The calculations of particle time step and interparticle collision process are performed using self-written program in VC language, which is taken as self-defined function (UDF) embedded in the calculation of gas-particle two-phase flow in GPIS. To obtain the rapid convergence of the calculation of gas-solid two-phase flow, the gas flow field is calculated to be convergent first, and then the particles are introduced. When a particle reaches the outlet, the particle tracking is terminated. If the number of particles in GPIS does not change greatly, the calculation is assumed to be convergent.


## 2.1. Gas-Particle Two-Phase Flow

An Eulerian-Lagrangian approach is adopted to describe the gas-particle two-phase flow in GPIS as mentioned above. The interactions between gas phase and particle phase are determined by means of Newtonian third law. Therefore, four-way coupling between discrete phase and continuous phase is carried out. The key features of the mathematical model are briefly descried as follows.


### 2.1.1. Governing Equations for Gas Phase

The gas flow field in GPIS is generally a turbulent flow field with backflows, and it is calculated using the realizablek-ε model [21] in this simulation. In this model, the conservation equations for mass and momentum can be written as follows:
(1)∂ρg∂t+∇·(ρgu)=0,(2)∂(ρgu)∂t+∇·(ρguu)=-∇p+∇·(τ)+ρgg+Fg,
where Fg is the external body forces that arise from interaction with the particle phase, which results in the momentum exchanges between the gas and particle phases. And Fg is calculated from the reverse action of particles. To close (2), the modeled transport equations for the turbulence kinetic energy, k, and its dissipation rate, ε, are
(3)∂(ρgk)∂t+∇·(ρgku)=∇·[(μ+μtσk)∇k]+Gk+Gb-ρgε-YM,∂(ρgε)∂t+∇·(ρgεu)=∇·[(μ+μtσε)∇ε]+ρgC1Sε-ρgC2ε2k+νε+C1εC3εGbεk,
where C1=max[0.43,η/(η+5)]; η=Sk/ε;μt=ρgCμk2/ε;S=(2SijSij)1/2; and Gk=μtS2 represents the generation of turbulence kinetic energy due to the mean velocity gradients. Gb is the generation of turbulence kinetic energy due to buoyancy. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. The model constants are C1ε=1.44,  C2=1.9,  σk=1.0,  and  σε=1.2 [22].


### 2.1.2. Particle Motion Equations

In GPIS, the spherical particles are entrained by gas flow. The translational motion of a particle in the gas phase is governed by Newton’s second law of motion and can be written as(4)mpdvpdt=Fd+F,
where F is an additional force term used in the simulation, including virtual mass force, pressure gradient force, and Saffman’s lift force. In addition, Fd represents the drag force of gas phase acting on the particle which is the most important force acting on the particle calculated using the following equation:
(5)Fd=0.125πDpμCdRe(u-v).
The drag coefficient, Cd, is determined by the following equation:
(6)Cd={24Re(1+16Re2/3)Re≤1000,0.424Re>1000,
where Re is the relative Reynolds number, which is defined as
(7)Re≡ρgDp|u-v|μ.


## 2.1.1. Governing Equations for Gas Phase

The gas flow field in GPIS is generally a turbulent flow field with backflows, and it is calculated using the realizablek-ε model [21] in this simulation. In this model, the conservation equations for mass and momentum can be written as follows:
(1)∂ρg∂t+∇·(ρgu)=0,(2)∂(ρgu)∂t+∇·(ρguu)=-∇p+∇·(τ)+ρgg+Fg,
where Fg is the external body forces that arise from interaction with the particle phase, which results in the momentum exchanges between the gas and particle phases. And Fg is calculated from the reverse action of particles. To close (2), the modeled transport equations for the turbulence kinetic energy, k, and its dissipation rate, ε, are
(3)∂(ρgk)∂t+∇·(ρgku)=∇·[(μ+μtσk)∇k]+Gk+Gb-ρgε-YM,∂(ρgε)∂t+∇·(ρgεu)=∇·[(μ+μtσε)∇ε]+ρgC1Sε-ρgC2ε2k+νε+C1εC3εGbεk,
where C1=max[0.43,η/(η+5)]; η=Sk/ε;μt=ρgCμk2/ε;S=(2SijSij)1/2; and Gk=μtS2 represents the generation of turbulence kinetic energy due to the mean velocity gradients. Gb is the generation of turbulence kinetic energy due to buoyancy. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. The model constants are C1ε=1.44,  C2=1.9,  σk=1.0,  and  σε=1.2 [22].


## 2.1.2. Particle Motion Equations

In GPIS, the spherical particles are entrained by gas flow. The translational motion of a particle in the gas phase is governed by Newton’s second law of motion and can be written as(4)mpdvpdt=Fd+F,
where F is an additional force term used in the simulation, including virtual mass force, pressure gradient force, and Saffman’s lift force. In addition, Fd represents the drag force of gas phase acting on the particle which is the most important force acting on the particle calculated using the following equation:
(5)Fd=0.125πDpμCdRe(u-v).
The drag coefficient, Cd, is determined by the following equation:
(6)Cd={24Re(1+16Re2/3)Re≤1000,0.424Re>1000,
where Re is the relative Reynolds number, which is defined as
(7)Re≡ρgDp|u-v|μ.


## 2.2. Interparticle Collision

The modified DSMC method, based on a gridless approach and local time stepping with automatic adaption, is applied for modeling the interparticle collisions in GPIS here. The basic ideas of this method are as follows:(1)
real particles are replaced by smaller number of sampled particles. And the trajectories of sampled particles are calculated. In this way, the computer memory and computation time can be saved significantly. This item is just the same as the traditional DSMC method;(2)
sampled particle motion is decoupled into the movement and collision processes as the traditional DSMC method. Sampled particle movement obeys the single particle motion model, and the collision process follows the interparticle collision dynamics. In the modified DSMC method, the particle time step, that is, the time required for a particle in each movement, is adaptive, which is different from the traditional DSMC method;(3)
collision pairs are found through the theory of collision probability instead of using trajectories as the traditional DSMC method. Nevertheless, the calculation of the collision probability between particles and the searching of collision pairs is based on the spherical region generated from the local particle parameters rather than the flow field cells.To update the particle time step, the local collision mean-free path, that is, the average distance a particle travels between collisions with other moving particles, is calculated firstly using the following equation:(8)L=vpυ,
where vp is the velocity magnitude of the particle being tracked. υ is the collision frequency, that is, the number of collisions occurred in unit time, which is obtained according to the local distribution of particles. To satisfy the “principle of uncoupling” in the modified DSMC method, the particle time step Δtp is given by
(9)Δtp=min[L3vx,L3vy,L3vz],
where vx, vy, and vz are components of particle velocity vp. Δtg is the time step for gas phase.In the modified DSMC method, the collision pairs are found through the spherical region, that is, searching scope, with the particle tracked being considered as center andr as its radius given below
(10)r=max[vpΔtp,vr,maxΔtp],
where vr,max is the maximum value of relative velocity between the particle being tracked and the sampled particles nearby. Then, the probability of collision between particle i being tracked and particle j during a time step Δtp is calculated by
(11)Pij=0.25π(Di+Dj)2|G(0)|np,jΔtpV,
where Di and Dj are the diameters of particle i and particle j, respectively. G(0) is the relative velocity vector between particles i and j. np,j is the real particle number represented by the sampled particle j. V is the volume of the searching scope and equal to 4πr3/3.The modified Nanbu method [23] is used to search a candidate collision partner j, which probably collide with particle i in the searching scope during a time step Δtp, and decide whether the collision occurs. First, a random number R with a uniform distribution from zero to unity is extracted from a generator. The candidate collision partner j is selected according to
(12)j=int[R×N]+1,j≠i,
where int[R×N] is defined as the integer part of R×N and N is the total number of all sampled particles in the searching scope. If the relation
(13)R>jN-Pij
is satisfied, particle i would collide with particle j during the time step Δtp, and the velocities of particles i and j are replaced with the postcollision velocities, but without changing their positions. The post-velocities of the two particles are as follows:when(n·G(0)/|Gct(0)|)<1/(ξ(1+e))(14)vi=vi(0)-(n+ξt)(n·G(0))(1+e)mjmi+mj,vj=vj(0)+(n+ξt)(n·G(0))(1+e)mimi+mj,when(n·G(0)/|Gct(0)|)≥1/(ξ(1+e))(15)vi=vi(0)-[(1+e)(n·G(0))n+|Gct(0)|t]mjmi+mj,vj=vj(0)+[(1+e)(n·G(0))n+|Gct(0)|t]mimi+mj,where n is the normal unit vector directed from particle i to particle j at the moment on contact. t is the unit vector in the tangential direction. Gct(0) is the tangential component of the relative velocity.


## 2.3. Numerical Solution

The numerical solutions of the gas-particle two-phase flow, including the governing equations for gas continuous phase, motion equations for sampled particles, and the momentum exchanges between two phases, are obtained using the CFD code FLUENT. The calculations of particle time step and interparticle collision process are performed using self-written program in VC language, which is taken as self-defined function (UDF) embedded in the calculation of gas-particle two-phase flow in GPIS. To obtain the rapid convergence of the calculation of gas-solid two-phase flow, the gas flow field is calculated to be convergent first, and then the particles are introduced. When a particle reaches the outlet, the particle tracking is terminated. If the number of particles in GPIS does not change greatly, the calculation is assumed to be convergent.


## 3. Model Validation

To validate the established model, two experimental cases (Case 1 and Case 2) in the literature [24] are calculated first. Figure 1 shows the schematic diagram of the experimental apparatus with two opposed nozzles. The inner diameter of either nozzle was 8 mm. The distance between two nozzles was 140 mm. Particles were conveyed by air into the nozzles. The parameters of the gas and particle phases injecting from nozzles for two cases are also summarized in Table 1.Table 1
Conditions and parameters used in simulation.
Parameter
Value
Gas phase
Gas type
Air
Inlet velocity
100 m/s
Apparatus
Nozzle diameter,d
8 mm
Distance between nozzles,l
140 mm
Particle spray angle
14°
Mass ratio of gas particle
1
Particle phase
Particle shape
Sphere
Density
2500 kg/m3
Restitution coefficient,e
0.9
Friction coefficient,ξ
0.2
Inlet velocity
25 m/s (Case 1) 55 m/s (Case 2)
Average diameter,Dp
77μm (Case 1) 33 μm (Case 2)Figure 1
Schematic diagram of the experimental apparatus.The coordinate system is shown in Figure1. The origin of the coordinate system is placed on the midpoint of the line connecting the two nozzles. The axial direction x-axis is on the horizontal line through nozzle axes and defined as positive rightward. The x-axis is called the axial line of the GPIS device. Figure 2 shows the meshed geometry for the calculation domain, and the domain is meshed with about 80000 hexahedral cells. The mesh refinement in the zone between two nozzles is carried out to raise the accuracy of gas phase solutions because of the big variation gradients of various variables in this zone. The minimum mesh size is about 1.3 mm.Figure 2
Grid model and the drawing of partial enlargement.Figure3 shows the comparison between experimental and simulated dimensionless particle concentration profiles along the axial line, in which C0 is the particle concentration at the outlet of nozzle. It can be seen that the particle concentration decreases sharply with the increasing distance from the outlet of nozzle along the axial line because of dispersion, and it increases obviously near the center because of the particle collisions. The predicted results of the particle concentration are in reasonable agreement with the experimental data, which indicates the validity of the present method.Comparison between experimental and simulated dimensionless particle concentration profiles along the axial line.
(a)
Case 1
(b)
Case 2


## 4. Results and Discussion

Figure4 shows the sketch of a laboratory-scale GPIS device with two opposed nozzles [25]. The impinging chamber is made of stainless steel of 0.5 m in diameter. The distance between the two nozzles (l) is set to 200 mm; both nozzles have a diameter of 42 mm. The outlet diameter of the impinging chamber is relatively small to collect particles. The following part of work is to model the particle behavior in the impinging zone. To reduce the computational cost, the calculation domain is located near the nozzles framed by the black dash-line rectangle as shown in Figure 4. The calculation domain is meshed with hexahedral cells. The mesh refinement in the zone near the height of nozzles is carried out to raise the accuracy of gas phase solutions. In order to obtain grid-independent solutions, several mesh sizes, that is, 62740, 142394, and 274836, were tested. The results of the particle concentration along the axial line computed using 142394 and 274836 meshes were closer to each other; the maximum differences in the predicted results were indeed less than 10%. On the other hand, finer meshes required excessive computational time; therefore, mesh size of 142394 was selected in this study (see Figure 4).Figure 4
Sketch of a laboratory-scale GPIS device and grid meshing of calculation domain.The parameters used in the simulation are summarized in Table2 with reference to the parameters used practically. The boundary conditions are “velocity-inlet” and “outflow” for inlet and outlet, respectively. All walls are treated as nonslip boundaries in the lateral direction for the gas phase with standard wall function. To obtain the rapid convergence of the calculation of gas-solid two-phase flow, the gas flow field is calculated to be convergent first, and then the particles are introduced. When a particle reaches the outlet, the particle tracking is terminated. If the number of particles in GPIS does not change greatly, the calculation is assumed to be convergent.Table 2
Conditions and parameters used in simulation.
Parameter
Value
Size of computational domain
Distance between nozzles,l
400 mm
Nozzle diameter,d
42 mm
Calculating parameter [19]
Restitution coefficient,e
0.9
Friction coefficient,ξ
0.2
Gas time step,Δtg
1
×
10
-
4 s
Particle spray anglea
0°
Initial condition for inletb
Particle inlet velocity
20 m/s
Particle density
2450 kg/m3
Particle average diameter
200μm
Gas type
Air
Gas inlet velocity
25 m/s
Gas density
1.2 kg/m3
Inlet particle mass flow rate (single nozzle), Wp
0.01, 0.02, 0.027, 0.035, 0.04, 0.045, 0.05, 0.054, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.14, 0.16, 0.18, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 kg/s
Initial condition for outlet
Pressure
101325 Pa
aTo reduce the effect of dispersion on the particle behavior, the particle spray angle is set to 0°.
bThe gas-solid two-phase flow is fully developed at the outlet of nozzle.The coordinate system is shown in Figure4. The origin of the coordinate system is placed on the midpoint of the line connecting two nozzles. The rightward x-axis on the nozzle axes is called the axial line of the GPIS device. In addition, coordinate surface yz is called impinging plane.


### 4.1. Particle Motion Behavior

Figure5 shows the particle position distribution at different moments for the case with Wp=0.2 kg/s in which the colormap denotes the particle residence time, and Figure 6 shows the velocity distribution of particles between two nozzles. From those two figures, it is found that the particle velocity remains about the same before the time point 0.01 s at which two particle flows start to meet each other near the origin of the coordinate system. This is mainly because that the collision between particles is caused by difference in particle velocities. Before this time point, particles move towards the impinging plane along the axial direction. The velocity difference between neighboring particles is small, which results in the small interparticle collision probability. In addition, the change in particle velocities is not obvious due to the small impacting impulse if two particles collide with each other.Particle position distributions at different time moments,Wp=0.2 kg/s.
(a)
0.006 s
(b)
0.01 s
(c)
0.02 sParticle velocity distribution between two nozzles,Wp=0.2 kg/s, t=0.8 s.
(a)
Particle velocity component on the axial line
(b)
Particle velocity component in the impinging planeWhen two gas-particle flows with a great number of particles meet each other at 0.01 s in the area around the origin of the coordinate system, particles from the opposite streams collide with each other violently. This area is called impingement zone generally. The collisions in this zone increase the velocity components in the impinging plane and decrease the axial velocity components of particles. Therefore, particles cannot pass through this area and penetrate into the opposite stream easily. These particles spread out from this area after several collisions as the particle position distribution at 0.02 s shown in Figure5(c).


### 4.2. Particle Concentration Distribution

Figure7 shows the residence time distribution of particles between two nozzles. It is found that the violent particle collisions increase the particle residence time in the impingement zone. In other words, the violent particle collision makes the particles accumulate in this zone, which results in the highest particle concentration near the origin as shown in Figure 8. And the high particle concentration strengthens the particle collision in turn. After particles leave the impingement zone in all directions, the particle concentration decreases drastically because of the increasing volume occupied by the particles.Figure 7
Particle residence time distribution between two nozzles,Wp=0.2 kg/s, t=0.8 s.Particle concentration distribution,Wp=0.2 kg/s, t=0.8 s.
(a)
Horizontal plane through nozzles (z=0)
(b)
Impinging plane (x=0)Figure9 shows the particle concentration profiles along the axial line for different inlet particle mass flow rates. It can be seen that particle concentration along the axial line increases with the increasing inlet particle mass flow rate reasonably. For a certain case, the particles injected into the device move towards the impinging plane along the axial line and the collisions with neighboring particles obviously can not change the particle paths as mentioned above. Therefore, the particle concentration does not change significantly along the axial line in the areas between either the nozzle or the impingement zone. In the impingement zone, the particle concentration increases sharply and reaches a maximum value near the origin due to the particle collisions. The relation between the maximal particle concentration Cmax in the impingement zone and the inlet particle concentration C0 is drawn in Figure 10. It can be seen that the maximum value Cmax of particle concentration increases sharply with the increasing inlet particle concentration C0 for lower C0. The relative particle concentration, ratio of Cmax and C0, also increases sharply at first. With the further increase in C0, the increase in Cmax with C0 tends to slow down, and the ratio of Cmax and C0 decreases obviously. The change of Cmax with C0 can be fitted by a logarithmic function with a correlation coefficient 0.9860 as follows:
(16)Cmax=aln(C0+b)-c,
where a=46.317,  b=1.6264, and  c=22.527. The ratio of  Cmax with  C0 reaches the maximum value as  C0=7.2 kg/m3.Figure 9
Particle concentration profiles along the axial line for different inlet particle mass flow rates.Figure 10
Maximal particle concentration in impingement zone versus inlet particle concentration.


### 4.3. Particle Collision Position Distribution

Figure11 shows the distribution of particle collision positions in 0.01 s for the case with Wp=0.09 kg/s, where one red point represents one collision occurring there. In the 0.01 s, there are 184059 collisions between sampled particles in total. According to the distribution of particle collision positions, the particle distribution can be divided into three zones: particle-collision zone, particle-jetting zone, and particle-scattering zone. The statistical result shows that particle collisions occur mainly in the particle-collision zone in which the number of collisions accounts for about 90% of the total. The particle-collision zone is an ellipsoid whose center is located at the origin of coordinates. The shorter axis on x-axis and the longer axis on y-axis of the ellipsoid are about 0.075l and 0.1l for the case with Wp=0.09 kg/s, respectively. In this zone, particles from the opposite streams collide with each other violently, which causes two results. First, particle collisions make the particles accumulate in this zone, which results in the highest particle concentration shown in Figure 8. Second, particle collisions make the particles spread out from this zone to the whole device. Between the particle-collision zone and two nozzles, the two long and narrow zones are called particle-jetting zone. Although the particle concentration in this zone is high, the smaller collision probability between particles caused by small velocity difference results in the smaller number of collisions in these two zones. The broad zone outside the particle-collision zone and the particle-jetting zones is called the particle-scattering zone. Particles spreading out from the particle-collision zone scatter quickly in this zone. The collision probability between particles in this zone is small because of the low particle concentration shown in Figure 8. Therefore, there are only a few collisions occurring near the particle-collision zone and few collisions occur in the area away from the particle-collision zone.Figure 11
Statistical result of particle collision positions in 0.01 s,Wp=0.09 kg/s.Figure12 shows the profiles of interparticle collision rate Nc along the axial line for different inlet particle mass flow rates, where the interparticle collision rate Nc is defined as the number of interparticle collisions per unit time per unit volume. It is found that the profiles of interparticle collision rate along the axial line for different inlet particle mass flow rates are similar to those of particle concentration shown in Figure 9, and the maximum value of the interparticle collision rate is reached in the particle-collision zone due to the violent collisions. This phenomenon indicates that the interparticle collision rate has direct relation with the particle concentration, and the violent particle collisions occur in the particle-collision zone with high particle concentration. Figure 13 shows the relation between the average interparticle collision rate at the nozzle outlet and the local particle concentration for different cases. The fitting curve is as follows:
(17)Nc,0=a×C02±0.00316,
where a is a constant and equal to 1.48×108. The correlation coefficient of this fitting curve is 0.9966. From this formula, it can be seen that the interparticle collision rate is directly proportional to the square of particle concentration, which is in accordance with the statement in the literature [26].Figure 12
Profiles of interparticle collision rateNc along the axial line for different inlet particle mass flow rates.Figure 13
Interparticle collision rate at the nozzle outlet versus local particle concentration.


## 4.1. Particle Motion Behavior

Figure5 shows the particle position distribution at different moments for the case with Wp=0.2 kg/s in which the colormap denotes the particle residence time, and Figure 6 shows the velocity distribution of particles between two nozzles. From those two figures, it is found that the particle velocity remains about the same before the time point 0.01 s at which two particle flows start to meet each other near the origin of the coordinate system. This is mainly because that the collision between particles is caused by difference in particle velocities. Before this time point, particles move towards the impinging plane along the axial direction. The velocity difference between neighboring particles is small, which results in the small interparticle collision probability. In addition, the change in particle velocities is not obvious due to the small impacting impulse if two particles collide with each other.Particle position distributions at different time moments,Wp=0.2 kg/s.
(a)
0.006 s
(b)
0.01 s
(c)
0.02 sParticle velocity distribution between two nozzles,Wp=0.2 kg/s, t=0.8 s.
(a)
Particle velocity component on the axial line
(b)
Particle velocity component in the impinging planeWhen two gas-particle flows with a great number of particles meet each other at 0.01 s in the area around the origin of the coordinate system, particles from the opposite streams collide with each other violently. This area is called impingement zone generally. The collisions in this zone increase the velocity components in the impinging plane and decrease the axial velocity components of particles. Therefore, particles cannot pass through this area and penetrate into the opposite stream easily. These particles spread out from this area after several collisions as the particle position distribution at 0.02 s shown in Figure5(c).


## 4.2. Particle Concentration Distribution

Figure7 shows the residence time distribution of particles between two nozzles. It is found that the violent particle collisions increase the particle residence time in the impingement zone. In other words, the violent particle collision makes the particles accumulate in this zone, which results in the highest particle concentration near the origin as shown in Figure 8. And the high particle concentration strengthens the particle collision in turn. After particles leave the impingement zone in all directions, the particle concentration decreases drastically because of the increasing volume occupied by the particles.Figure 7
Particle residence time distribution between two nozzles,Wp=0.2 kg/s, t=0.8 s.Particle concentration distribution,Wp=0.2 kg/s, t=0.8 s.
(a)
Horizontal plane through nozzles (z=0)
(b)
Impinging plane (x=0)Figure9 shows the particle concentration profiles along the axial line for different inlet particle mass flow rates. It can be seen that particle concentration along the axial line increases with the increasing inlet particle mass flow rate reasonably. For a certain case, the particles injected into the device move towards the impinging plane along the axial line and the collisions with neighboring particles obviously can not change the particle paths as mentioned above. Therefore, the particle concentration does not change significantly along the axial line in the areas between either the nozzle or the impingement zone. In the impingement zone, the particle concentration increases sharply and reaches a maximum value near the origin due to the particle collisions. The relation between the maximal particle concentration Cmax in the impingement zone and the inlet particle concentration C0 is drawn in Figure 10. It can be seen that the maximum value Cmax of particle concentration increases sharply with the increasing inlet particle concentration C0 for lower C0. The relative particle concentration, ratio of Cmax and C0, also increases sharply at first. With the further increase in C0, the increase in Cmax with C0 tends to slow down, and the ratio of Cmax and C0 decreases obviously. The change of Cmax with C0 can be fitted by a logarithmic function with a correlation coefficient 0.9860 as follows:
(16)Cmax=aln(C0+b)-c,
where a=46.317,  b=1.6264, and  c=22.527. The ratio of  Cmax with  C0 reaches the maximum value as  C0=7.2 kg/m3.Figure 9
Particle concentration profiles along the axial line for different inlet particle mass flow rates.Figure 10
Maximal particle concentration in impingement zone versus inlet particle concentration.


## 4.3. Particle Collision Position Distribution

Figure11 shows the distribution of particle collision positions in 0.01 s for the case with Wp=0.09 kg/s, where one red point represents one collision occurring there. In the 0.01 s, there are 184059 collisions between sampled particles in total. According to the distribution of particle collision positions, the particle distribution can be divided into three zones: particle-collision zone, particle-jetting zone, and particle-scattering zone. The statistical result shows that particle collisions occur mainly in the particle-collision zone in which the number of collisions accounts for about 90% of the total. The particle-collision zone is an ellipsoid whose center is located at the origin of coordinates. The shorter axis on x-axis and the longer axis on y-axis of the ellipsoid are about 0.075l and 0.1l for the case with Wp=0.09 kg/s, respectively. In this zone, particles from the opposite streams collide with each other violently, which causes two results. First, particle collisions make the particles accumulate in this zone, which results in the highest particle concentration shown in Figure 8. Second, particle collisions make the particles spread out from this zone to the whole device. Between the particle-collision zone and two nozzles, the two long and narrow zones are called particle-jetting zone. Although the particle concentration in this zone is high, the smaller collision probability between particles caused by small velocity difference results in the smaller number of collisions in these two zones. The broad zone outside the particle-collision zone and the particle-jetting zones is called the particle-scattering zone. Particles spreading out from the particle-collision zone scatter quickly in this zone. The collision probability between particles in this zone is small because of the low particle concentration shown in Figure 8. Therefore, there are only a few collisions occurring near the particle-collision zone and few collisions occur in the area away from the particle-collision zone.Figure 11
Statistical result of particle collision positions in 0.01 s,Wp=0.09 kg/s.Figure12 shows the profiles of interparticle collision rate Nc along the axial line for different inlet particle mass flow rates, where the interparticle collision rate Nc is defined as the number of interparticle collisions per unit time per unit volume. It is found that the profiles of interparticle collision rate along the axial line for different inlet particle mass flow rates are similar to those of particle concentration shown in Figure 9, and the maximum value of the interparticle collision rate is reached in the particle-collision zone due to the violent collisions. This phenomenon indicates that the interparticle collision rate has direct relation with the particle concentration, and the violent particle collisions occur in the particle-collision zone with high particle concentration. Figure 13 shows the relation between the average interparticle collision rate at the nozzle outlet and the local particle concentration for different cases. The fitting curve is as follows:
(17)Nc,0=a×C02±0.00316,
where a is a constant and equal to 1.48×108. The correlation coefficient of this fitting curve is 0.9966. From this formula, it can be seen that the interparticle collision rate is directly proportional to the square of particle concentration, which is in accordance with the statement in the literature [26].Figure 12
Profiles of interparticle collision rateNc along the axial line for different inlet particle mass flow rates.Figure 13
Interparticle collision rate at the nozzle outlet versus local particle concentration.


## 5. Conclusions

The modified direct simulation Monte Carlo (DSMC) method is applied to predict the particle behavior in GPIS in this work. The main results are as follows:(1)
according to the distribution of particle collision positions, the particle distribution can be divided into three zones: particle-collision zone, particle-jetting zone, and particle-scattering zone. Particle collisions occur mainly in the particle-collision zone, and the distribution of particle collision positions is very similar to the particle concentration distribution;(2)
for GPIS, a few particles penetrate into the opposite stream due to the violent particle collisions and almost all particles spread out from the particle-collision zone;(3)
the particle concentration and the interparticle collision rate along the axial line have a similar variation trend. They increase sharply and reach their maximum values near the origin in the particle-collision zone. The maximum value of the particle concentration increases with the increasing inlet particle concentration according to a logarithmic function. The interparticle collision rate is directly proportional to the square of local particle concentration.



---

*Source: 254082-2013-11-28.xml*