Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows

G Model PARTIC-730; No. of Pages 11

ARTICLE IN PRESS Particuology xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Particuology journal homepage: www.elsevier.com/locate/partic

Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows夽 Yongxiang He, Haibo Zhao ∗ , Haoming Wang, Chuguang Zheng State Key Laboratory of Coal Combustion, Huazhong University of Science and Technology, Wuhan 430074, China

a r t i c l e

i n f o

Article history: Received 13 December 2013 Received in revised form 28 April 2014 Accepted 12 May 2014 Keywords: Direct simulation Monte Carlo Differentially weighted method Gas–solid flow Particle–particle collision Four-way coupling

a b s t r a c t In gas–solid flows, particle–particle interaction (typical, particle collision) is highly significant, despite the small particles fractional volume. Widely distributed polydisperse particle population is a typical characteristic during dynamic evolution of particles (e.g., agglomeration and fragmentation) in spite of their initial monodisperse particle distribution. The conventional direct simulation Monte Carlo (DSMC) method for particle collision tracks equally weighted simulation particles, which results in high statistical noise for particle fields if there are insufficient simulation particles in less-populated regions. In this study, a new differentially weighted DSMC (DW-DSMC) method for collisions of particles with different number weight is proposed within the framework of the general Eulerian–Lagrangian models for hydrodynamics. Three schemes (mass, momentum and energy conservation) were developed to restore the numbers of simulation particle while keeping total mass, momentum or energy of the whole system unchanged respectively. A limiting case of high-inertia particle flow was numerically simulated to validate the DW-DSMC method in terms of computational precision and efficiency. The momentum conservation scheme which leads to little fluctuation around the mass and energy of the whole system performed best. Improved resolution in particle fields and dynamic behavior could be attained simultaneously using DW-DSMC, compared with the equally weighted DSMC. Meanwhile, computational cost can be largely reduced in contrast with direct numerical simulation. © 2014 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Published by Elsevier B.V. All rights reserved.

Introduction Gas–solid flows are frequently found in industrial processes such as fluidized beds, pneumatic conveying, boilers and furnaces, electrostatic precipitators, and particle separation in cyclones. These systems often involve complicated flow dynamics and interaction between flow constituents and their surroundings. A number of numerical studies of gas–solid flow have been conducted focusing on one-way coupling (Deutsch & Simonin, 1991) and two-way coupling (Squires & Eaton, 1990; Wang, Zhao, Guo, He, & Zheng, 2013). One-way coupling that the influence of the solid particle on the continuous phase is neglected, is reasonable when the particle fractional volume ˚v and mass loading ˚m are

夽 This paper is adapted from the presentation at the 4th UK-China International Particle Technology Forum, October 15–19, Shanghai, China, as recommended by Prof. Xiaoshu Cai and Dr. Jerry Heng, the co-chairs of the scientific committee. ∗ Corresponding author. Tel.: +86 27 87545526; fax: +86 27 87545526. E-mail addresses: [email protected], [email protected] (H. Zhao).

small (e.g., ˚v < 10−6 ). However, with an increase in fractional volume ˚v (e.g., 10−6 < ˚v < 10−3 ), the effect of the continuous phase on the dynamics of the dispersed phase and the feedback of the dispersed phase on the continuous phase dynamics should be considered simultaneously. The conservation equations of continuous phase include appropriate source terms resulting from the dispersed phase. The interaction between two phases is called twoway coupling. While the volume fraction of particles ˚v continues to increase (e.g., ˚v > 10−3 ), flows are referred as dense suspension. The particle–particle collision plays an important role on the profiles of continuous phase and dispersed phase, the term fourway coupling effect emerges. In fact, it is essential to consider particle–particle collision even though the fractional volume of the particles ˚v is small (e.g., ˚v < 10−3 ), because the turbulent transport effect and the preferential concentration effect lead to an increase in inter-particle collision rates by a factor of 30 (Wang, Wexler, & Zhou, 2000). Indeed, the average particle fractional volume is not the only measure of the importance of particle collision. Oesterle and Petitjean (1993) performed a developed horizontal channel flow calculation in a non-dilute gas–solid suspension flow.

http://dx.doi.org/10.1016/j.partic.2014.05.013 1674-2001/© 2014 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Published by Elsevier B.V. All rights reserved.

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model PARTIC-730; No. of Pages 11

ARTICLE IN PRESS Y. He et al. / Particuology xxx (2014) xxx–xxx

2

Nomenclature c d Ets Fother g l Mts N Nc Nf P Pts S t t Te u Vs w

particle velocity, m/s diameter, m initial total energy of the whole system external and inter-particle forces gravity acceleration, m/s2 grid length, m initial total mass of the whole system total particle numbers average collision rate per unit volume, m−3 total simulation particles the probability of one particle interacting with any other particle initial total momentum of the whole system source term time scale, s time step, s the integral Lagrange time scale, s fluid velocity, m/s; velocity volume of the whole system, m3 number weight

Greek letters ˇ collision cross-section, m3 /s  turbulent kinetic energy effective transport tensor  ε turbulent kinetic energy dissipation rate, m2 /s3 ı relative error time-accumulated relative error fluid density, kg/m3 ; density, kg/m3  p relaxation time scale of particle, s Kolmogorov time scale, s k ˚v particle volume fraction, m−3 Subscripts analytical value a e discarded particle i, j, k index of simulation particle coordinate direction k new condition after dynamic event old condition before dynamic event p particle phase property between two simulation particles r t detected value x, y, z coordinate directions indication of one of the decomposed particles 1 2 indication of the other decomposed particle ϕ the general transport variable Superscripts max maximum value minimum value min — ensemble average value * particle velocity after collision event  post-collision and pre-restoration; a new coordinate system  post-collision and post-restoration

It was found that particle–particle collision has a strong influence on the profile of the vertical particle concentration compared with dilute flow simulation, revealing that particle–particle collision plays a non-negligible role once the particle mass loading

exceeds unity. In many cases, the particle loading is moderate to high and hence particle–particle collisions significantly influence the flow field. Therefore, the particle–particle collision and fourway coupling between the continuous and discrete phases should be processed for the purpose of accurate description and control in dispersed systems. Various models have been proposed to consider particle–particle collision in the general framework of Eulerian–Lagrangian models, where the continuous flow behavior is studied in Euler coordinates and the motion of dispersed phase is described by Lagrangian equation. Wassen and Frank (2001) classify these models in two principally different ways: the trajectory calculation (TC) method and simultaneous particle tracking (SPT) method. The inherent assumption of the TC method is that each trajectory represents a constant flow of particles with identical physical properties, and macroscopic properties of the dispersed phase can be acquired by averaging all trajectories that cross the local numerical cell. As a consequence, this method is limited to applications of steady flow. While the SPT method is inherently unsteady because each particle represents a certain number of real particles with identical physical properties, and the macroscopic properties of the particulate phase are achieved by averaging all particles in the same cell. To take particle–particle collision into account in the frame of the TC method, Oesterle and Petitjean (1993) proposed an iterative technique with which particle–particle collision is treated stochastically on the basis of the local concentration and velocities obtained from the previous iteration. In the frame of the SPT method, collisions can be computed either deterministically or stochastically. The most straightforward deterministic trajectory method, in which all particles have to be tracked simultaneously through the flow field and the occurrence of particle–particle collision can be judged by particle position and relative motion during one time step, is based on molecular dynamics. Sundaram and Collins (1996) used direct numerical simulation (DNS) method to deal with particle–particle collision (hard-sphere collisions, similar to the present study) to investigate the particle collision rate in isotropic turbulent flows. Two different techniques (proactive and retroactive) are used for particle collisions in the DNS method. The former (proactive) technique anticipates all the collisions that will occur within one time step and then enacts them in order. While the latter (retroactive) looks for particle–particle overlaps at the end of one time step, and then resolves these overlaps by enacting particle collisions retroactively in the order in which they occurred until all overlaps have been considered (Chen, Kontomaris, & McLaughlin, 1999; Sundaram & Collins, 1996). However, the DNS method is not able to handle multiple particle collisions that are of great importance for the modeling of quasi-static systems. The distinct element method (DEM) first proposed by Cundall and Strack (1979) for granular dynamic simulation provides new possibilities for discrete particle simulation to calculate the dense phase flows using a soft-sphere collision model to deal with particle collision dynamics. Tsuji, Kawaguchi, and Tanaka (1993) performed a two-dimensional gas-fluidized bed simulation where particle collisions are modeled by DEM and observed the phenomena of the formation of bubbles and slugs and the process of particle mixing which is consistent with experiment. Note that the inherent computation cost of the deterministic trajectory method while choosing colliding particle pairs is proportional to N2 , where N is the number of traced particles. For this reason, the deterministic calculation is of no practical relevance in some engineering applications such as high mass loading impinging streams (Xu, Zhao, & Zheng, 2014). In fact, in the majority of the applications of the SPT method, collisions were treated stochastically. The direct simulation Monte Carlo (DSMC) method developed by Bird (1976, 1994) is a physically

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model PARTIC-730; No. of Pages 11

ARTICLE IN PRESS Y. He et al. / Particuology xxx (2014) xxx–xxx

accurate method for the computation of non-equilibrium gas flows. This technique is most useful in circumstances where there are insufficient numbers of collisions in the flow to maintain the equilibrium forms of the distribution functions describing the various energy modes of the gas. Several stochastic models were developed to reduce the numerical effort to a linear proportionality of the total number of simulation particles. Bird (1976, 1989) proposed the event-driven time-counter method and time-driven no-timecounter (NTC) method for gas molecular collisions. In the frame of the NTC method, the number of collisions in each grid is calculated and then an acceptance–rejection method is used to judge whether a collision event will occur. In addition, other stochastic models (Illner & Neunzert, 1987; Ivanov & Rogasinsky, 1988; Koura, 1986; Nanbu, 1980) for selecting colliding pairs were also developed in the dilute molecular gas dynamics. These stochastic models take advantage of the analogy between the motion of discrete particles in gas–solid flows and the motion of molecules in dilute gas flow to consider particle–particle collision in gas–solid flows. Tanaka, Yonemura, and Tsuji (1995) simulated a two-dimensional upward gas–solid flow in a rectangular domain with a periodic boundary condition to study the effect of the particle property on the structure of the particle cluster. Tsuji, Kawaguchi, and Tanaka (1998) simulated the riser of a circulating fluidized bed to investigate the formation of particle aggregation compared with the two-fluid model. Furthermore, a stochastic particle–particle collision model for particle-laden flow was suggested by Sommerfeld (2001) and Oesterle and Petitjean (1993), in which a generated simulation particle, whose size and velocity are sampled from local probability density functions of particle diameters and velocities, collides with the selected particle during each time step. The post-velocity of the selected simulation particle is solved by the conservation equation for linear and angular momentum, while the properties of the simulation particle generated are of no further interest. The advantage of this model is that it does not require information on the position of surrounding particles. However, it needs to conserve the particle size and velocity distribution of each grid. The studies mentioned above focused on the collision of simulation particles with equal number weight. In many applications of interest, the traced particles possess only a small portion of multi-species or polydisperse particle population, and the equally weighted DSMC faces enormous difficulties or complex particle dynamics in terms of either huge computational cost or significant statistical noise or even both. To circumvent this difficulty, a scheme to apply different weight to simulation particles is required, depending on the species or polydisperse particle population. A conservative weighting scheme (Boyd, 1996) was developed to improve the resolution of the very small quantities in multispecies non-equilibrium gas flows. Nonetheless, the scheme does not conserve linear momentum and energy explicitly at each collision. Zhao, Kruis, and Zheng (2009) proposed a new differentially weighted time-driven Monte Carlo (MC) method, which could capture the coagulation dynamics in dispersed systems with low noise and track the size distribution over the full size range simultaneously. So far, the differentially weighted scheme has been successfully adopted for particle coagulation in the population balance (Zhao & Kruis, 2014; Zhao, Kruis, & Zheng, 2010; Zhao, Kruis, & Zheng, 2011; Zhao & Zheng, 2009a,b, 2013); however, there are few reports on particle collision in gas–solid two-phase flows. If differentially weighted (DW) simulation particles are tracked in gas–solid flows, the key issues include how to design a rule (model) for a collision event between two differentially weighted simulation particles, restoring the number of simulation particles after a simulation collision event, and coupling the DW-DSMC with the Eulerian–Lagrangian model. A splitting and restoring procedure was proposed for particle collision with different number weight.

3

The first stage was to split a simulation particle with larger weight into two particles. Mass, momentum and energy conservation schemes were then adopted to restore the number of simulation particles. The purpose of this study is to validate the DW-DSMC model with different schemes for particle number restoration and to set up a general simulation strategy. Model description Simulation strategy The equations governing continuous phase operation are conveniently expressed in the Eulerian frame of reference. By focusing on a small finite control volume in the suspension, the conservation equations for the continuous phase operation can then be written as (Lun & Liu, 1997): ∂(uj ϕ) ∂ ∂(ϕ) = + ∂t ∂xj ∂xj



ϕ

∂ϕ ∂xj



+ Sϕ ,

(1)

where  is the fluid density;  ␸ is an effective transport tensor; S␸ is a ϕ-dependent source term that is the feedback influence of the dispersed phase on the continuous phase; and ϕ is a general transport variable that could be 1 (for mass conservation), uj (fluid velocity for momentum conservation), and E (total energy for energy conservation). Considerable research has focused on S␸ and many models have been developed such as the Reynolds stress model (Zhou, 1993). The position and velocity of every single particle can be described as follows using Newton’s second law (Lun & Liu, 1997): dxpi dt dupi dt

= upi , = gi +

(2) 1 (u − upi ) + Fother , p i

(3)

where xpi and upi are the particle position and velocity, respectively; the first term (gi ) on the right-hand side of Eq. (3) is the gravitational acceleration; the second term is the drag force;  p = p dp 2 /(18 ) is the particle relaxation time (where p is the particle density, dp is the particle diameter, and is the fluid viscosity); and the last term (Fother ) expresses the external and inter-particle forces contributing to the particle transport. Our particular focus here is on the coupling strategy of the differentially weighted Monte Carlo and multiphase flow model, and thus the detailed models for multiphase flows are not discussed. As is known, most numerical strategies solving the gas–solid two-phase flow in the framework of Eulerian–Lagrangian model are time-driven. Therefore, it is possible to embed the timedriven differentially weighted Monte Carlo method with the Eulerian–Lagrangian model into a single framework. By setting an appropriate time increment t, it is thought that the three processes, fluid flow, particle motion, and particle collisions, are uncoupled. As a result, the flow fields characterized by Eq. (1), the particle fields characterized by Eqs. (2) and (3), and particle collisions (see Section “Treatment of particle collisions”) are independently solved, and four-way coupling is realized by the exchange of the environmental variables for the continuous and discrete phases within t. The numerical simulation flowchart for the four-way coupling model is shown in Fig. 1. The probability of particle–particle collision A probabilistic collision rule was developed here for colliding particles with different number weights, analogous to the probabilistic coagulation rule (Zhao et al., 2009; Zhao & Zheng, 2009a,b).

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model

ARTICLE IN PRESS

PARTIC-730; No. of Pages 11

Y. He et al. / Particuology xxx (2014) xxx–xxx

4

Determination of the time step As noted above, the most important factor in the simulation strategy for coupling the differentially weighted Monte Carlo for particle collision with the Eulerian–Lagrangian model for hydrodynamics is the appropriate time step within which the fluid transport, the particle transport, and the particle collision are uncoupled and then separately simulated. Therefore, the time step t needs to be smaller than the particle relaxation time scale  p and the flow integral time scale Te to satisfy the general two-way coupling condition. In reality, to narrow the scope of potential colliding particle pairs, the maximum displacement in each direction of a simulation particle is constrained to be smaller:

 t ≤ min



lk

vmax pk

,

k = 1, 2, 3

(7)

where vpk max and lk are the maximum velocity and grid length along the kth direction, respectively. Furthermore, to couple the particle–particle collision, each particle is restricted to participate in, at most, one collision event. The waiting time between two successive collision events of the simulation particle i can be obtained: t ≤ tcollision =

1 . Vs Ci

(8)

To sum up, the time step t is restricted as follows: t ≤ min(p , Te , min(lk /vmax ), min(1/Vs Ci )), pi,k ∀i,∀k

∀i

(9)

where vpi,k max is the maximum velocity of particle i in the kth direction. Selection of colliding particle pairs Fig. 1. Flowchart of simulation strategy for four-way coupling.

For the simulation particle i and j with weight wi and wj , diameter dpi and dpj , mass mpi and mpj , and velocity upi and upj , respectively, each real particle from the simulation particle i undergoes a collision event with a probability of min(wi , wj )/wi , and each real particle from j does so with a probability of min(wi , wj )/wj . As a result, only min(wi , wj ) real particles from i or j participate in real collisions on average. The collision rate of simulation particle i and j can be described as (Zhao et al., 2009; Zhao & Zheng, 2009a,b): 2

ˇij = Pij =

(dpi + dpj ) cr 4

,

(4)

ˇij wj t

2wi wj

Vs

(wi + wj ) min(wi , wj )

=

ˇij wj t 2 max(wi , wj ) Vs

(wi + wj )

, (5)

where Vs is the volume of the whole system, ˇij is the collision cross-section of particle i and j in unit time, cr is the relative velocity of particle i and j, and cr = |upi − upj |. Thus, the total collision rate of simulation particle i with any other simulation particles is:

Nf Ci

=

P j=1, j = / i ij Vs

=

1 Vs2

Nf  j=1,j = / i

ˇij wj

2 max(wi , wj ) (wi + wj )

,

(6)

where Nf is the total number of simulation particles. Note that 2 Pij = (dpi + dpj ) cr wt/(4Vs ), if wi = wj = w. Thus, the equally weighted collision rule is only a special case of the differentially weighted probabilistic collision rule.

Several stochastic models such as the no-time-counter (Bird, 1989) and Nanbu’s methods (Nanbu, 1980) have been proposed as described in the introduction section. In this study, modified Nanbu’s method is employed to deal with particle collision and it can be described as follows: Step 1: generate a uniform random number R ∈ (0, 1) and the other particle probably involved in the collision event of a particle i is j = int(R × N) + 1 (int is the integer formulation, N is total particle number in the local grid); Step 2: calculate the probability of particle i and j according to Eq. (5); Step 3: particles i and j will collide with each other if R > (j/N − Pij ), which will be introduced in Section “Treatment of particle collisions”; otherwise, return to Step 1 until all simulation particles are checked. In fact, the modified Nanbu’s method is similar to the acceptance–rejection method. The computational cost of this technique is linearly proportional to the numbers of simulation particles of the whole system, because the collision partner j is determined by a mathematical formula. The time step t is more restricted because the modified Nanbu’s method divides the unit length into equal portions of N with each of 1/N shown in Fig. 2. With respect to time step, the following condition needs to be satisfied: Pij =

ˇij wj t 2 max(wi , wj ) Vs

(wi + wj )

<

1 . N

(10)

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model PARTIC-730; No. of Pages 11

ARTICLE IN PRESS Y. He et al. / Particuology xxx (2014) xxx–xxx

5

Fig. 2. Implementation of the modified Nanbu’s method.

Treatment of particle collisions Collision dynamics between two equally weighted simulation particles For the collision event between two equally weighted simulation particles (the common number weight w = wi = wj ), the number of collision events between two real particles is equal to the common number weight w. Because the simulation particle is an indicator of the real particles, the collision dynamics between two equally weighted simulation particles is fully identical to that of two represented real particles. A new coordinate system (x , y , z ) is established at the collision point of the two particles, where the x -direction is along the centerline from particle i and j; the y -direction is situated in the plane determined by the x -direction and the vector (upj − upi ), and is perpendicular to the x’-direction; and the z -direction is always perpendicular to both the x -direction and y -direction. The post-collision velocities of particles i and j are u∗pi and u∗pj , respectively. Using the law of conservation of momentum, the relationships between the pre- and post-collision velocities in the coordinate system (x , y , z ) can be written as (Lun & Liu, 1997):

⎧ Jx ∗ ⎪  = wpi,z  ⎨ u∗pi,x = upi,x − m ; v∗pi,y = vpi,y ; wpi,z pi

⎪ ∗ ⎩ u∗pj,x = upj,x + Jx ; v∗pj,y = vpj,y ; wpj,z  = wpj,z 

(11)

mpj

where the impulse Jx = (upj,x − upi,x )(2mpi mpj )/(mpi + mpj ). Rotational movement of particle is ignored in this study. However, the post-collision properties, i.e. linear and angular velocities, of simulation particle j and 1 can be obtained by solving the conservation equation for linear and angular momentum, i.e., the impulse equations in connection with Couloumb’s law of friction (Tsuji, Tanaka, & Yonemura, 1998), if particle rotation is considered. Therefore, the angular velocity of simulation particle 1 and j should be updated after the collision event. In conclusion, the DW-DSMC method can be extended to systems with the consideration of particle rotation. The post-collision velocities of the two particles in the original coordinate system (x, y, z) are obtained through transforming the coordinate system (x , y , z ) into (x, y, z) since the impact position of the two particles at the moment of collision first needs to be determined. The impact point of the two particles can only be located if two of the three angles are known. A stochastic process is used to determine this impact point. It is assumed that the probability density of finding the impact point on the surface of one particle is uniform. Thus the new coordinate system (x , y , z ) is transformed into the original coordinate system (x, y, z) and then the post-collision velocities of the two particles in (x, y, z) can be obtained. Collision dynamics between two differentially weighted simulation particles This study developed a simple binary collision model for two differentially weighted simulation particles i and j. As described in Section “The probability of particle–particle collision”, once a

collision between simulation particles i and j is considered, the real particles of simulation particle i interact with the real particles of simulation particle j. It is thought that only min(wi ,wj ) real particles of simulation particles i and j undergo a real collision event. In other words, the number of real collision events between real particles is wj provided that wi > wj . Therefore, the simulation particle i is split into two simulation particles: 1 with number weight wj and 2 with number weight (wi –wj ). The other environmental parameters of simulation particles 1 and 2 are inherited from the parent particle i. Simulation particle 1 interacts with simulation particle j using the hard-sphere model, which belongs to the variant of collision dynamics between two equally weighted simulation particles (see Section “Collision dynamics between two equally weighted simulation particles”), while simulation particle 2 does not participate in the dynamic process as illustrated in Fig. 3. Parametric variation during the splitting process can be described as:

⎧ (w ) = (w ) ; j old ⎨ 1 new ⎩

(m1 )new = (mi )old ;

(up1 )new = u∗pm

(w2 )new = (wi )old − (wj )old ;

(m2 )new = (mi )old ;

(up2 )new = (upi )old

(wj )new = (wj )old ;

(mj )new = (mj )old ;

(upj )new = u∗pj

(12)

where subscripts “old” and “new” mean pre- and post-collision, respectively. It is noticeable that each collision event between two differentially weighted simulation particles results in a net generation of simulation particles. As a result, the increasing simulation particles have to be traced simultaneously, resulting in high CPU cost. Some measures need to be adopted to restore the number of simulation particles. Lin, Lee, and Matsoukas (2002) and Smith and Matsoukas (1998) developed a Monte Carlo algorithm called constant number Monte Carlo in which the number of simulation particles remains constant throughout the simulation, regardless of whether the actual growth process leads to a net loss (coagulation) or gain (fragmentation). In the light of the constant number Monte Carlo, we proposed mass, momentum and energy conservation schemes to restore the number of simulation particles. These schemes are described as follows: Mass conservation scheme: Once a collision event between two differentially weighted simulation particles has occurred, a simulation particle randomly chosen from the subsystem (here this indicates the grid in which a collision event occurs) is discarded. The number weight of the remaining simulation particles in the subsystem is recalculated to maintain the mass unchanged. Momentum conservation scheme: Discard a simulation particle randomly as in the mass conservation scheme. The total momentum is conserved to alter the number weight of the remaining simulation particles. Energy conservation scheme: Different from mass conservation scheme, conserved energy of the subsystem is used as an indicator to update the weight of the simulation particles. Table 1 shows the variation of the subsystem properties at different stages while dealing with the collision event, where Mts , Pts , and Ets mean the total mass, momentum, and energy before the collision event (pre-collision), respectively; and we , mpe , and vpe are the weight, mass, and velocity of the discarded particle.

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model

ARTICLE IN PRESS

PARTIC-730; No. of Pages 11

Y. He et al. / Particuology xxx (2014) xxx–xxx

6

Fig. 3. Schematic of collision process for differentially weighted particles.

wi , mpi , and vpi are the weight, mass, and velocity of the simulation particle i before the collision event, wi , mpi , and vpi are the weight, mass, and velocity of the simulation particle i after the collision dynamics and before the restoration, and wi , mpi , and vpi are the weight, mass, and velocity of the simulation particle i after the collision dynamics and after the restoration. Here we simply adjust the particle number weight to keep the system conservation parameter (mass, or momentum, or energy) constant before and after the restoration. Therefore, mpi = mpi = mpi ; vpi = vpi , while vpi = / vpi because of the collision dynamics. Variation in the individual weights of the remaining simulation particles in the subsystem with the mass conservation scheme is described as follows:



 Mts =

i=1

 = Mts



Ns +1

Ns

wi mpi =



wi =

wi mpi =





Ns



wi

wi mpi

=

Nc = wi =

i=1

+ we mpe  wi ,  Mts

1+

we mpe v2Pe



 2Ets

wi

(15)

The average inter-particle collision rate, Nc , is a key parameter that describes the numbers of particle collisions and can be calculated as follows:

i=1

,

(14)

Model validation

wi mpi + we mpe  Mts

 + w m v Pts e pe pe  wi  Pts

for energy conservation scheme,

Ns

i=1

= Mts =



Therefore, the weight variation using the momentum and energy conservation schemes can be described as:for momentum conservation scheme,

(13)

Nct , Vs t

(16)

where Nct is the detected collision events within one time step. The analytical average collision rate can be obtained in some limited conditions such as the zero-inertia particle motion ( p   k , where  k is the turbulent Kolmogorov time scale) in simple homogeneous shear flow (Wang, Wexler, & Zhou, 1998) and the free molecular motion of high-inertia particles ( p  Te ) (Abrahamson, 1975). In this study, the high-inertia particle flow is simulated to validate our model. As the particle motion is not affected by the gas flow, analogous to that of molecular gas, the particle velocity is constant

where Mts is the total mass after dealing with the collision dynamics and before restoring the simulation particle number  means the total mass after (post-collision and pre-restoration), Mts discarding a randomly selected simulation particle and restoring the sample (post-collision and post-restoration). Table 1 Variation in the subsystem properties at different stages dealing with particle collision. System parameter

Pre-collision

Post-collision and pre-restoration

Post-collision and post-restoration

Total particle numbers

Ns

Ns + 1

Ns

Total particle mass

Mts =

Ns 

wi mpi

 Mts =

i=1

Total particle momentum

Pts =

Ns 

Ets =

1 2

wi mpi + we mpe

 Mts =

i=1

wi mpi vpi

 Pts

=

i=1 Ns

Total particle energy

Ns 

 i=1

Ns 

 Ets =

1 2

wi mpi

i=1

wi mi pi

 Pts

v + we mpe vpe

=

i=1 Ns

wi mpi v2pi

Ns 

 i=1

Ns 

wi mpi vpi

i=1 Ns

wi mpi v2pi +

1 w m 2 2 e pe pe

v

 Ets =

1 2



wi mpi v 2pi

i=1

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model

ARTICLE IN PRESS

PARTIC-730; No. of Pages 11

Y. He et al. / Particuology xxx (2014) xxx–xxx

within a time step increment t. In addition, the high-inertia particle flows are simulated using the DNS and equally weighted DSMC method for particle collisions as references. The technique used to identify collisions in DNS is the retroactive method. To reduce the computational cost, an optimization technique is implemented on the time-sequenced descending order collision list (Chen, Liu, & Zheng, 2004; Sundaram & Collins, 1996). The particles are first moved backward in time to the point of collision of the first overlapping particle pair, the elastic collision between the pair is affected, and the particles are then moved forward in time. Overlap computations are redone for collided particles and new overlaps are added to the schedule list while nonexistent overlaps are deleted, after which the entire list is re-sorted in descending order. These instructions are carried out until there are no remaining overlaps to consider. Case 1: monodisperse high-inertia particle flow Abrahamson (1975) developed a turbulent collision kernel with regard to high-inertia particles. As described above, the particle motion is similar to that of molecular gas. Thus, when the initial particle velocities follow a Maxwell distribution, the mean particle collision rate can be written as follows by applying molecular kinetic theory:

Nca = 0.5n2 dp2

0.5

16 v2p 3

,

(17)

where Nca is the analytical average collision rate per unit volume and unit time, n the particle number concentration, and vp 2

the mean square value of the particle fluctuation velocity, which

7

Table 2 Initial parameters of monodisperse high-inertia particle flow. Initial parameter

Heavy particles

Computational region Boundary condition Number of cells Initial particle velocity distribution Total particle numbers, N Total fictitious particles, Nf Particle diameter, dp Time step, t (s) Time evolution, T (s)

0 < x < 2,0 < y < 2,0 < z < 2 Periodical boundary in three direction 8×8×8 Maxwell distribution and vp 2 ≈ 3 105 8 × 104 0.1 0.0001 0.01

is dependent on the particle fields. Table 2 presents the initial parameters of the monodisperse high-inertia particle flow case. The particle volume fraction ( v = Nf dp3 /482 ) is 16.7% to ensure multiple collision events during each time step. The total number of real particles for the whole system is 105 and the total number of simulation particles 8 × 104 , hence the average number weight of each simulation particle using an equally weighted DSMC is 1.25 (=N/Nf ). In this case, to validate the differentially weighted scheme, different number weights are applied artificially to the simulation particles. The weight ratio of one half of the simulation particles to the other half is 1.5. Simulation particles are initially distributed uniformly in the computational region using a stochastic game. A constant time step is used throughout the simulation because of the collision dynamic event and hard-sphere elastic collision model, reducing computational time for evolution in a time step. Fig. 4(a) shows the time evolution of the average collision rate per unit volume predicted by the DW-DSMC, DSMC, and DNS, compared with analytical results according to Eq. (17). The results

Fig. 4. Comparison of statistical parameters between different numerical methods: (a) the average collision rate per unit volume, (b) the relative error in the average collision rate per unit volume, and (c) the time-accumulated relative error of average collision rate per unit volume.

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model PARTIC-730; No. of Pages 11

ARTICLE IN PRESS Y. He et al. / Particuology xxx (2014) xxx–xxx

8

obtained using the stochastic model for particle collisions fluctuate up and down around analytical data, and are of the same order of accuracy. It can be inferred that the variation in the average particle collision rate using DNS is larger than that using the other two numerical methods. We further calculated the relative error ıNC in the average collision rate per unit volume, ıNC = |Nc − Nc,a |/Nc,a . We found that the relative error using the stochastic model for particle collisions fluctuated within 7.5% while the relative error using DNS is in the range of 0.048–18.5% as shown in Fig. 4(b). Furthermore, the time-accumulated relative error for the average collision rate NC was calculated, NC (t) =

t

Nc dt − tNc,a |/tNc,a . Fig. 4(c) shows the evolution of the time-

| 0

accumulated relative error for the average collision rate per unit volume. It was inferred that the time-accumulated relative error NC (t) using each numerical method reached an asymptotic value and it was generally constrained within 3.7% for a relatively long evolution time (the 20th time step, t = 0.002 s). It was obvious that the time-accumulated relative errors using the momentum and energy conservation schemes were slightly larger than those when using the mass conservation scheme. Comparison of the proposed schemes for the restoration of the number, mass, momentum, and energy of the simulation particles for the total system at the end of each time step, expressed as Mi , Pi , Ei , respectively, was also observed. Fig. 5(a) and (c) presents the time evolution of the normalized total system mass, momentum, and energy. It was found that not only the momentum of the

total system, but also the mass and energy were conserved during each time step when using the momentum conservation scheme. Meanwhile, the mass conservation scheme sustained high statistical noise levels for the energy of the whole system, although it fitted the momentum of the whole system well as shown in Fig. 5(b). It was obvious that the mass and momentum of the whole system increased with time after 30 time steps using the energy conservation scheme, as shown in Fig. 5(a) and (b). Therefore, the momentum conservation scheme for the restoration of simulation particles was adopted in this study. Quantitative comparisons of the detailed information of particle fields such as particle number density and particle turbulent kinetic energy are also discussed. Fig. 6(a) and (b) shows the information on the detailed particle fields at a specified time-point (the 50th time step t = 0.05 s) at a specific two-dimensional section (in the central section of the z-direction (z = ␲)). Irrespective of the particle number density or particle turbulent kinetic energy, the results obtained using DW-DSMC were in good agreement with those using DNS and DSMC. Therefore, we concluded that DW-DSMC was able to reasonably predict not only the statistical parameters such as the average collision rate per unit volume (shown in Fig. 4(a)), but also the temporal and spatial evolution of particle fields, which were determined by the ensemble averaging procedure over all simulation particles in the subsystem, such as the number density at a specific time-point (shown in Fig. 6(a)). The numerical differences between the results from DNS, DSMC, and DW-DSMC can be ascribed to the following two factors. (1) Initialization of particle number weight and velocity. The actually traced simulation particles for each simulation method were

Fig. 5. Total system properties as a function of time for monodisperse high-inertia particle flow normalized by the value of the total system initial properties for: (a) mass, (b) momentum, and (c) energy.

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model

ARTICLE IN PRESS

PARTIC-730; No. of Pages 11

Y. He et al. / Particuology xxx (2014) xxx–xxx

9

Table 3 Total number of simulation particles and particle number weight for each species using different simulation methods. Numerical method

dp 0.05

DNS DSMC DW-DSMC

0.1

Nf, 0.05

w0.05

Nf, 0.1

w0.1

1000 100 1000

1 10 1

9.9 × 104 9.9 × 103 9.9 × 103

1 10 10

not the same. Thus, statistical noise would be present for the initial velocity distribution at the start of each simulation despite the same seed for random number generator used for all numerical simulations. (2) Strategies adopted for particle collision. In the frame of the SPT method, collision using DNS that is deterministic is ascertained based on particle trajectory crossing, while it depends on the probability between selected particle pairs in DSMC and DW-DSMC. Furthermore, the collision dynamics between two simulation particles are the dominant difference between DSMC and DW-DSMC in the monodisperse high-inertia particle flow. Case 2: bidisperse high-inertia particle flow In this case, we conducted a bidisperse numerical experiment to illustrate the advantage of DW-DSMC. It considered a collection of 105 real particles and the initial particle velocity distribution followed a Maxwell distribution. The total system consisted of 1% of the particles that had a diameter of 0.05 and 99% of the particles that had a diameter of 0.1. Table 3 shows the weight and number of initial simulation particles with diameters of 0.05 and 0.1 using different numerical methods, where Nf,0.05 and Nf,0.1 mean the total simulation particles with diameters of 0.05 and 0.1, respectively. Other initial simulation requirements were the same as those in the monodisperse case. Note that the system initial experimental conditions for simulation particles with a diameter of 0.1 using

Fig. 6. Distribution of the detailed information of particle fields using different numerical methods undergoing particle collision in the monodisperse high-inertia flow at a specific time (t = 0.05 s) at the specific two-dimensional section (in the central section in the z-direction (z = ␲)), for: (a) particle number density and (b) particle turbulent kinetic energy.

DSMC and DW-DSMC, and 0.05 using DNS and DW-DSMC, were the same (see Table 3). Fig. 7(a) and (b) presents the distribution of the number density and turbulent energy of particles with a diameter of 0.1 at a specified time-point (the 50th time step t = 0.05 s) at a specific

Fig. 7. Distribution of the detailed information of particle fields using different numerical methods undergoing particle collision in the bidisperse high-inertia flow at a specific time (t = 0.05 s) at a specific two-dimensional section (in the central section in the z-direction (z = ␲)), for: (a) number density of particles with a diameter of 0.1, (b) turbulent energy of particles with a diameter of 0.1, (c) number density of particles with a diameter of 0.05, and (d) turbulent energy of particles with a diameter of 0.05.

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model

ARTICLE IN PRESS

PARTIC-730; No. of Pages 11

Y. He et al. / Particuology xxx (2014) xxx–xxx

10

Table 4 Computational cost for high-inertia particle flows using different numerical methods. Computational cost

DSMC

DNS

DW-DSMC (momentum conservation scheme)

Case 1 (CPU t/s) Case 2 (CPU t/s)

77.695 3.0313

815.828 926.924

99.348 3.350

two-dimensional section (in the central section in the z-direction (z = ␲)). It was obvious that there was no any difference between the results using DSMC and DW-DSMC, which were in the vicinity of those using DNS, because of the statistical noise as illustrated above. Fig. 7(c) and (d) presents the distribution of the number density and turbulent energy of particles with a diameter of 0.05. The results obtained using DW-DSMC were in excellent agreement with those using DNS. However, considerable statistical noise was observed for both the distribution of the particle number density and particle turbulent energy using DSMC to deal with particle collisions, because only a small number of simulation particles with a diameter of 0.05 lay in a grid. In conclusion, the results using DWDSMC were correct enough and captured the collision dynamics with low noise in regions where there were insufficient particle numbers. Table 4 shows the computational cost using different numerical methods for the monodisperse and bidisperse high-inertia particle flows. All the simulations in this study were performed on the same desktop PC equipped with a Pentium (R) Dual-Core E6700 @3.2 GHz CPU, and 1.96 GB of memory. The reduction of computational cost in the DSMC and DW-DSMC simulations were the result of lower particle numbers and the use of stochastic processes to capture collision events. Compared with DSMC, the increased computational cost for the DW-DSMC method was because of the restoration of the simulation particle numbers.

Conclusions In this study, a new differentially weighted DSMC (DW-DSMC) for particle collisions with different weights is proposed, and a simulation strategy of coupling DW-DSMC with the general Eulerian–Lagrangian models for hydrodynamics within the same framework is developed. The mass, momentum and energy conservation schemes are proposed for simulation particle number restoration using a constant number Monte Carlo approach. Limiting monodisperse and bidisperse cases for high-inertia particle flows were numerically simulated to validate the DW-DSMC method. From this work, we conclude the following: (1) The monodisperse high-inertia particle flow simulation confirms that the DW-DSMC method is able to correctly predict both the statistical parameters and temporal–spatial evolution of particle fields. (2) The momentum conservation scheme for DW-DSMC, in which the total momentum is conserved while discarding a simulation particle randomly, shows little fluctuation in the mass and energy of the whole system. Meanwhile, the mass and energy conservation schemes cannot conserve these properties simultaneously. (3) Improved resolution can be obtained simultaneously using DW-DSMC in the bidisperse particle flow case, in comparison with the equally weighted DSMC in which considerable statistical noise in the particle properties (dp = 0.05) (e.g., number density, turbulent energy) exists.

(4) Computational cost can be largely reduced using a stochastic model (DSMC, DW-DSMC) for particle collisions compared with the deterministic trajectory model (DNS). In short, this study focused on the development of a differentially weighted method, different schemes for simulation particle restoration and a simulation strategy for four-way coupling, but paid little attention to the complex two-phase turbulent models and other particle systems involved with particle rotation and polydisperse size distribution, which is easily incorporated with the DW-DSMC and require further study in future. The stochastic essence dealing with particle collision makes DW-DSMC method unable to obtain particle force information in the process of collision dynamics. Therefore, this method is limited to account for binary collision. As illustrated by Bird (1994), the statistical fluctuations associated with variable particle weight appear to be larger than normally to be expected from equally weighted DSMC simulation because of the random walk effect. In reality, adequate particle properties can be gained without consideration of variable weight when the ratio of traced particles to numbers of particles in the whole system exceeds 0.1. Therefore, the variable particle weighted scheme should only be used where absolutely necessary, for example, ratio of traced particles to numbers of particles in the whole system should be 0.1 or less (Boyd, 1996). It is also recommended that the total simulation particle numbers at a small portion of particle number concentration should be equal to those with higher portion by modulating the weight of the traced particles. This approach can weaken the random walk effect. Furthermore, investigation into the operability of the numerical method in simulating the complex particulate processes, and a new method for differentially weighted particle collision that could retain the system properties mathematically, are underway. Acknowledgements The authors were supported by the National Natural Science Foundation of China (51276077 and 51390494) and the National Key Basic Research and Development Program (2010CB227004). References Abrahamson, J. (1975). Collision rates of small particles in a vigorously turbulent fluid. Chemical Engineering Science, 30, 1371–1379. Bird, G. A. (1976). . Molecular gas dynamics. NASA STI/Recon technical report A (Vol. 76). Bird, G. A. (1989). Perception of numerical methods in rarefied gas dynamics. Progress in Astronautics and Aeronautics, 117, 211–226. Bird, G. A. (1994). Molecular gas dynamics and the direct simulation Monte Carlo of gas flows. Oxford: Clarendon Press. Boyd, I. D. (1996). Conservative species weighting scheme for the direct simulation Monte Carlo method. Journal of Thermophysics and Heat Transfer, 10, 579–585. Chen, M., Kontomaris, K., & McLaughlin, J. (1999). Direct numerical simulation of droplet collisions in a turbulent channel flow. Part I: collision algorithm. International Journal of Multiphase Flow, 24, 1079–1103. Chen, Y., Liu, Z., & Zheng, C. (2004). Direct numerical simulation method of interparticle collision. Chinese Journal of Computational Physics, 21(5), 421–426. Cundall, P. A., & Strack, O. D. (1979). A discrete numerical model for granular assemblies. Geotechnique, 29, 47–65. Deutsch, E., & Simonin, O. (1991). Large eddy simulation applied to the modelling of particulate transport coefficients in turbulent two-phase flows. In 8th Symposium on turbulent shear flows, Vol. 1 (pp. 10–11). Illner, R., & Neunzert, H. (1987). On simulation methods for the Boltzmann equation. Transport Theory and Statistical Physics, 16(2–3), 141–154. Ivanov, M., & Rogasinsky, S. (1988). Analysis of numerical techniques of the direct simulation Monte Carlo method in the rarefied gas dynamics. Russian Journal of Numerical Analysis and Mathematical Modelling, 3(6), 453–465. Koura, K. (1986). Null-collision technique in the direct-simulation Monte Carlo method. Physics of Fluids, 29, 3509–3511. Lin, Y., Lee, K., & Matsoukas, T. (2002). Solution of the population balance equation using constant-number Monte Carlo. Chemical Engineering Science, 57, 2241–2252.

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013

G Model PARTIC-730; No. of Pages 11

ARTICLE IN PRESS Y. He et al. / Particuology xxx (2014) xxx–xxx

Lun, C. K. K., & Liu, H. S. (1997). Numerical simulation of dilute turbulent gas–solid flows in horizontal channels. International Journal of Multiphase Flow, 23, 575–605. Nanbu, K. (1980). Direct simulation scheme derived from the Boltzmann equation. I. Monocomponent gases. Journal of the Physical Society of Japan, 49, 2042–2049. Oesterle, B., & Petitjean, A. (1993). Simulation of particle-to-particle interactions in gas solid flows. International Journal of Multiphase Flow, 19, 199–211. Smith, M., & Matsoukas, T. (1998). Constant-number Monte Carlo simulation of population balances. Chemical Engineering Science, 53, 1777–1786. Sommerfeld, M. (2001). Validation of a stochastic Lagrangian modelling approach for inter-particle collisions in homogeneous isotropic turbulence. International Journal of Multiphase Flow, 27, 1829–1858. Squires, K. D., & Eaton, J. K. (1990). Particle response and turbulence modification in isotropic turbulence. Physics of Fluids A: Fluid Dynamics, 2, 1191–1203. Sundaram, S., & Collins, L. R. (1996). Numerical considerations in simulating a turbulent suspension of finite-volume particles. Journal of Computational Physics, 124, 337–350. Tanaka, T., Yonemura, S., & Tsuji, Y. (1995). Effects of particle properties on the structure of clusters. ASME-Publications-FED, 228, 297–302. Tsuji, Y., Kawaguchi, T., & Tanaka, T. (1993). Discrete particle simulation of twodimensional fluidized bed. Powder Technology, 77, 79–87. Tsuji, Y., Tanaka, T., & Yonemura, S. (1998). Cluster patterns in circulating fluidized beds predicted by numerical simulation (discrete particle model versus twofluid model). Powder Technology, 95, 254–264. Wang, L. P., Wexler, A. S., & Zhou, Y. (1998). On the collision rate of small particles in isotropic turbulence. I. Zero-inertia case. Physics of Fluids, 10, 266–276. Wang, L. P., Wexler, A. S., & Zhou, Y. (2000). Statistical mechanical description and modelling of turbulent collision of inertial particles. Journal of Fluid Mechanics, 415, 117–153.

11

Wang, H., Zhao, H., Guo, Z., He, Y., & Zheng, C. (2013). Lattice Boltzmann method for simulations of gas–particle flows over a backward-facing step. Journal of Computational Physics, 239, 57–71. Wassen, E., & Frank, T. (2001). Simulation of cluster formation in gas–solid flow induced by particle–particle collisions. International Journal of Multiphase Flow, 27, 437–458. Xu, H., Zhao, H., & Zheng, C. (2014). Simulation and investigation of periodic deflecting oscillation of gas-solid planar opposed jets. Chemical Engineering and Processing: Process Intensification, 76, 6–15. Zhao, H., & Kruis, F. E. (2014). Dependence of steady-state compositional mixing degree on feeding conditions in two-component aggregation. Industrial & Engineering Chemistry Research, 53, 6047–6055. Zhao, H., Kruis, F. E., & Zheng, C. (2009). Reducing statistical noise and extending the size spectrum by applying weighted simulation particles in Monte Carlo simulation of coagulation. Aerosol Science and Technology, 43, 781–793. Zhao, H., Kruis, F. E., & Zheng, C. (2010). A differentially weighted Monte Carlo method for two-component coagulation. Journal of Computational Physics, 229, 6931–6945. Zhao, H., Kruis, F. E., & Zheng, C. (2011). Monte Carlo simulation for aggregative mixing of nanoparticles in two-component systems. Industrial & Engineering Chemistry Research, 50, 10652–10664. Zhao, H., & Zheng, C. (2009a). A new event-driven constant-volume method for solution of the time evolution of particle size distribution. Journal of Computational Physics, 228, 1412–1428. Zhao, H., & Zheng, C. (2009b). Correcting the multi-Monte Carlo method for particle coagulation. Powder Technology, 193, 120–123. Zhao, H., & Zheng, C. (2013). A population balance – Monte Carlo method for particle coagulation in spatially inhomogeneous systems. Computers & Fluids, 71, 196–207. Zhou, L. X. (1993). Theory and numerical modeling of turbulent gas particle flows and combustion. Beijing: Science Press.

Please cite this article in press as: He, Y., et al. Differentially weighted direct simulation Monte Carlo method for particle collision in gas–solid flows. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.05.013