Lattice Boltzmann based discrete simulation for gassolid fluidization
Abstract
Discrete particle simulation, a combined approach of computational fluid dynamics and discrete methods such as DEM (Discrete Element Method), DSMC (Direct Simulation Monte Carlo), SPH (Smoothed Particle Hydrodynamics), PIC (ParticleInCell), etc., is becoming a practical tool for exploring labscale gassolid systems owing to the fast development of parallel computation. However, gassolid coupling and the corresponding fluid flow solver remain immature. In this work, we propose a modified lattice Boltzmann approach to consider the effect of both the local solid volume fraction and the local relative velocity between particles and fluid, which is different from the traditional volumeaveraged NavierStokes equations. A timedriven hard sphere algorithm is combined to simulate the motion of individual particles, in which particles interact with each other via hardsphere collisions, the collision detection and motion of particles are performed at constant time intervals. The EMMS (energy minimization multiscale) drag is coupled with the lattice Boltzmann based discrete particle simulation to improve the accuracy. Two typical fluidization processes, namely, a single bubble injection at incipient fluidization and particle clustering in a fast fluidized bed riser, are simulated with this approach, with the results showing a good agreement with published correlations and experimental data. The capability of the approach to capture more detailed and intrinsic characteristics of particlefluid systems is demonstrated. The method can also be used straightforward with other solid phase solvers.
keywords:
Discrete particle simulation, Lattice Boltzmann method, Computational fluid dynamics, Fluidization, Multiphase flow, Simulation1 Introduction
Gassolid fluidization systems are widely encountered in both physical and chemical processes for many industries, for instance, fluid catalytic cracking (FCC), circulating fluidized bed combustion (CFBC), coal gasification, and sulfide roasting. Earlier studies of these systems mainly focused on experimental investigations including measurement of macroscopic hydrodynamic behavior and development of some corresponding correlations. In recent decades, to quantitatively understand the complex hydrodynamics of gassolid fluidization, the computational fluid dynamics approach is adopted in many cases, and a lot of numerical methods in the hydrodynamic modeling and simulation of gassolid fluidization have been proposed, such as twofluid model (TFM) (Anderson and Jackson, 1967; Ishii, 1975), quadraturebased moment methods (QBMM) (Fox, 2008, 2009a,b; Desjardins et al., 2008;), discrete particle simulation (DPS) (Tsuji et al., 1993; Hoomans et al., 1996; Xu and Yu, 1997), and direct numerical simulation (DNS) (Hu et al., 1992; Ma et al., 2006; Wang et al., 2010; Xiong et al., 2012).
Among these numerical methods, the most frequently used TFM treats the gas and solid phases as two interpenetrating continua, and locally averaged quantities such as volume fractions, velocities, species concentrations, and temperatures of gas and solid phases appear as dependent field variables (Anderson and Jackson, 1967; Ishii, 1975). To derive TFM using ensemble averaging techniques, terms such as effective stresses and the interphase interaction have to be introduced, which require constitutive equations for closure. Only under very limited conditions, those constitutive equations can be obtained rigorously from the kinetic theory of granular flow (Gidaspow, 1994), otherwise we have to resort to empirical models. The accuracy and effectiveness of TFM are, therefore, still unsatisfactory in many circumstances. The recently developed QBMM permits to solve population balance equation (PBE) in commercial CFD codes at relatively low computational cost. However, its application to the context of multiphase flows is to be explored (Mazzei, 2011). Comparably, DNS not only fully resolves the motion of each individual solid particle and fluid flow, but also directly calculates the hydrodynamic force acting on each individual solid particle from the stress on the fluidsolid boundary. Due to its capability in detailed solution around each particle, DNS has been regarded as the most accurate method for the simulation of gassolid flow. Unfortunately, DNS is too costly for predicting the hydrodynamics in large industrial scale fluidized beds even at low Reynolds numbers, let alone the high Reynolds number cases where their grid size and time step are limited by the Kolmogorov length scale and the turbulence time scale (Xiong et al., 2012).
For numerical modeling of gassolid fluidized beds mentioned above, TFM is computationally more economic but inaccurate, while DNS is computationally more accurate but expensive. So it is natural to ask whether there exists a better alternative combining the advantages of the two methods for modeling the gassolid flows. As a particlescale approach, DPS is somehow in between these two ends and seems to give a good balance among accuracy, cost and efficiency. Specifically, DPS resolves the continuum fluid flow at the scale of computational cells in CFD, describes the motion of individual particles by the wellestablished Newton’s equations of motion, and models particleparticle interactions through different collision models such as the hardsphere model and the softsphere model, which has been proven to be effective in modeling various particle flow systems (Deen et al., 2007; Zhu et al., 2007, 2008), such as slugging fluidized bed (Xu et al., 2007), spouted bed (Zhao et al., 2008), pneumatic conveying (Kuang et al., 2008), bubbling fluidized bed (Geng and Che, 2011), soundassisted fluidized bed (Wang et al., 2011), and cyclone separator (Chu et al., 2011). However, in all these mentioned work, the fluid motion with suspended solids is commonly governed by the volumeaveraged NavierStokes equations or their simplified forms (Tsuji et al., 1993; Hoomans et al., 1996; Xu and Yu, 1997; Mikami et al., 1998), and those equations are solved based on implicit schemes no matter by Fluent (Chu and Yu, 2008a, b; Chu et al., 2009a, b, 2011; Wu et al., 2010; Zhao et al., 2010), OpenFOAM (Su et al., 2011; Goniva et al., 2012), MFIX (Darabi et al., 2011; Garg et al., 2012; Li and Guenther, 2012; Li et al., 2012a,b; Gopalakrishnan and Tafti, 2013), or inhouse codes (Ouyang and Li, 1999a, b; Zhou et al., 2004a, b; Zhao et al., 2009; Wang et al., 2009; Wu et al., 2009). With implicit methods the discretized equations are solved simultaneously which inevitably requires some kind of global data dependence and hence global communication. Therefore, most algorithms involved suffer from relatively lower scalability and parallel efficiency, which becomes a grand challenge for fast simulation of largescale industrial systems.
As a smoothed alternative to lattice gas automata (LGA), lattice Boltzmann method (LBM) (McNamara and Zanetti, 1988) is an efficient secondorder flow solver capable of solving various systems for hydrodynamics owing to its explicit solution of particle distribution function, algorithmic simplicity, natural parallelism, and flexibility in boundary treatment (Chen and Doolen, 2003). Therefore, LBM becomes an increasingly popular approach to simulation of complex flows (Aidun and Clausen, 2010) and can be easily incorporated into DPS. Filippova and Hanel (1997) proposed a combination of latticeBGK model and Lagrangian approach, and performed threedimensional simulation of gasparticle flow through filters with oneway coupling, where the fluid affected the particles but the particles did not affect the fluid. Chen et al. (2004) simulated particleladen flow over a backwardfacing step with twoway coupling, where a modified latticeBGK model was developed for the fluid flow and a Lagrangian approach for particles. But they did not consider the effect of solid volume fraction on gas flows. Sungkorn et al. (2011) proposed a gasliquid LagrangianLBM to simulate turbulent gasliquid bubbly flows with a relatively low gas holdup. Specifically, they solved the continuous liquid phase by singlephase lattice Boltzmann equation (LBE) incorporated with large eddy simulation (LES) (Smagorinsky, 1963), and evolved the dispersed gas phase (i.e. the individual bubbles) by Lagrangian trajectories, but did not include the gas volume fraction in the conservation equations and its effect on drag force.
In this paper, we proposed a modified LBE to model the fluid flow and developed the corresponding fluidsolid interaction model in the framework of DPS. The effects of both the local solid volume fraction and the local relative velocity between particles and fluid are considered. The equations of motion governing individual particles are solved with timedriven hardsphere (TDHS) model. It is noteworthy that the computational strategy herein has ever been implemented in direct simulation of particlefluid systems (Wang et al., 2010) where the modified LBE was used with particle size much larger than the cell spacing. In the present work, the partial saturation concept has been extended to model the objects much smaller than the cell spacing (i.e. porous media), and both the linear and nonlinear drag effects of the solid phase (media) have been considered in the lattice Boltzmann based discrete particle simulation for the first time.
2 Numerical approach
The objective of this research is to develop a lattice Boltzmann based numerical method for discrete simulation of gassolid fluidization systems. For illustration, we used twodimensional ninevelocity (D2Q9) lattice Boltzmann model as an example, and the solid particles distributed in the lattice cell are described by the timedriven hard sphere model. A schematic diagram of this method is shown as Fig. 1.
2.1 LBM for gas phase hydrodynamics
The particle distribution function is the number of particles per unit volume with velocity at time t and position , which is often used in nonequilibrium statistical mechanics. The evolution of the distribution function f is described by the Boltzmann equation (Chapman and Cowling, 1970):
(1) 
The lefthand side of this equation is an advection term, while the righthand side contains the collision operator which describes the interaction between particles. External forces are neglected in this derivation. After a finite difference discretization of the Boltzmann equation in time, space and velocity, LBE is simplified as (McNamara and Zanetti, 1988),
(2) 
where the subscript i represents the orthogonal and diagonal directions in the Cartesian coordinates, is the particle distribution function in direction i, is the discrete velocity in direction i, x and t are the lattice simulation position and time, respectively, and is collision term in direction i. Assuming that the collision frequency is a function of spatial coordinates and time, and independent of the molecular velocity, the collision term can be approximated by the socalled BGK (Bhatnagar et al., 1954) model:
(3) 
Here is the relaxation time related to kinematic viscosity and the discrete time step , is the equilibrium distribution function expressed as (Qian et al., 1992)
(4) 
where is defined as for and for , and . The weights are given by , for , for , and is the speed of sound and equals to in lattice unit, and u are density and velocity, respectively.
The macroscopic quantities such as mass density and momentum density can then be obtained by evaluating the hydrodynamic moments of the distribution function , namely,
(5) 
By the ChapmanEnskog procedure (Sterling and Chen, 1996; Guo et al., 2002), the incompressible NavierStokes equations can be obtained from LBE in the limit of small Mach number. Therefore, the standard LBE is just suitable for singlephase flow. To model incompressible gas flow through suspended particles, an immersed moving boundary method for a medium with porosity needs to be introduced in LBM framework.
The key idea of the immersed moving boundary method is that the effect of boundary is modeled using a source term in the momentum equations. Here, the particlefluid coupling is implemented by immersed moving boundary method (Noble and Torczynski, 1998; Cook et al., 2004), which introduces a term that depends on the percentage of the cell saturated with fluid to modify the collision operator and represent the effect of both the solid volume fraction and the gassolid slip velocity on hydrodynamics. Although most of applications of the immersed moving boundary are associated with the objects whose size is much larger than the cell spacing (Noble and Torczynski, 1998; Cook et al., 2004; Feng et al., 2007; Wang et al., 2010; Zhou et al., 2011; Xiong et al., 2012; Wei et al., 2013), the method is theoretically not confined to this scale. Therefore, the partial saturation concept is extended to model the objects with much smaller size than the cell spacing, and the corresponding LBE modified for partially saturated cells (Noble and Torczynski, 1998) reads
(6) 
where is an additional collision term that accounts for the bounce back of the nonequilibrium part of the distribution function, and is given by (Noble and Torczynski, 1998)
(7) 
where denotes the component of the distribution function opposite to and is the average velocity of solid phase at the computational lattice. can be approximately calculated as
(8) 
where is the velocity of solid particle k and is the number of solid particles in fluid cell (see Fig.1).
The weighting function in Eq. (6) depends on the relaxation time and the solid volume fraction of each cell, which is given, among other forms, by (Noble and Torczynski, 1998)
(9) 
In our case, since the gas velocity presents a spatial average over areas fully, partially or not occupied by the solid phase, a different expression of may be developed. But as a rough estimation, we will at first neglect the actual difference of gas velocity at different locations within each lattice and continue to use the expression of Eq. (9). Obviously, and represent pure fluid and pure solid states, respectively.
The solid volume fraction is related to the porosity of lattice cell , namely . , an important parameter that influences both the gasphase and solidphase motions, can be calculated based on the area occupied by the particles in the lattice cell. The twodimensional porosity of solid phase is given by
(10) 
where r is the radius of solid particle and h is lattice spacing between nodes. As the dragforce correlation is based on 3D systems, we also use a 3D porosity transformed from 2D porosity for this purpose. According to (Ouyang and Li, 1999a),
(11) 
To capture turbulent structures in gassolid fluidization systems which are usually associated with large Reynolds numbers or become turbulent in nature, mesoscale turbulence models are incorporated into the volumeaveraged NavierStokes equations in DPS (Zhou et al., 2004a, b; Gui et al. 2008; Liu and Lu, 2009), here LES incorporated into LBM is used (Krafczyk et al., 2003; Yu et al., 2005). The effect of the small subgrid eddies on the largescale flow structures are modeled through an additional turbulent viscosity .
In Smagorinskybased subgrid scale (SGS) model (Smagorinsky, 1963), depends on the strain rate:
(12) 
Here is the Smagorinsky constant and strain rate tensor can be obtained directly by computing the momentum fluxes , which are secondorder moments of the nonequilibrium distribution function (Yu et al., 2005):
(13) 
where is the kth component of the lattice velocity . Based on an eddy relaxation time assumption for subgrid scale stress, the effect of the flow structure at the unresolved scale is modeled through an effective collision relaxation time , which is the eddy relaxation time with respect to the eddy viscosity . The eddy viscosity is then incorporated into the LBE by using instead of in Eq. (6). Accordingly,
(14) 
where c is the lattice speed given by . From Eq. (13) and Eq. (14), a quadratic equation is obtained, which yields
(15) 
It should be pointed out that Eq. (15) is suitable for both laminar and turbulent flows. On the one hand, under high gas velocity, Q is large and hence , which expresses the effect of motion at unresolved scale. On the other hand, under low gas velocity, Q vanishes and , describing the laminar flow.
2.2 TDHS for particle dynamics
Each individual solid particle, treated as a pointvolume particle, has four properties: mass (m), radius (r), position (P), and velocity (v). According to Newton’s second law of motion, the equation of motion for solid particle k is
(16) 
where g is the gravitational acceleration, is the drag force, is the collision force, is the particle’s volume, and is the pressure gradient.
Timedriven hardsphere model (Hopkins and Louge, 1991) is used for particleparticle collisions; while particles interact as binary and quasiinstantaneous hardspheres, and both the collisions are detected and the particle state are updated at constant time intervals. Namely, when two particles are in contact with each other and if the internal product of and is negative, which means that the particles are moving closer to each other, they will collide (see Fig. 2). For the collision between particle 1 and particle 2, the postcollisions velocities can be obtained by the following formulae:
(17) 
(18) 
where “*” means postcollision velocities, e is the restitution coefficient that represents the ratio of speeds after (post) and before (pre) collision, and the velocities in the righthand side of the equations are precollision velocities. The walls involved in simulations are composed of infinitely heavy particles with zero velocity, and the collisions between particles and the wall also satisfy Eqs. (17) and (18). In the next time step, the particles move to new positions with their new velocities.
It should also be noted that particleparticle interactions are important in simulating gassolid fluidization, especially for dense gasfluidized beds. In timedriven hardsphere model mentioned above, only the binaryparticle collision mechanism and the normal component of contact force are considered, which limits the model to low solid holdup conditions (less than 3040% by volume). For dense gasfluidized beds, discrete element method (DEM) (Cundall and Strack, 1979) might be expected to reflect both the multiparticle collision mechanism and the tangential component of contact force.
2.3 Interphase momentum transfer
The drag correlations are of great importance in numerical models that predict the flow behavior of gassolid fluidized beds. The interphase momentum transfer coefficient based on the traditional Ergun and Wen & Yu correlations (Gidaspow, 1994; Gidaspow and Jiradilok., 2009) is only valid for homogeneous particle suspensions, and insufficient to capture the heterogeneous structures of gassolid fluidization. A structuredependent drag based on the energy minimization multiscale model (EMMS) is used here, that is (Yang et al., 2003)
(19) 
According to Yang et al. (2003), the drag correction factor
(20) 
The drag coefficient for isolated rigid spherical particle can be calculated by the Schiller and Nauman (1935) equation
(21) 
Here the particle Reynolds number is defined as follows:
(22) 
It has been demonstrated that the drag correction factor , varying in value between 0.0152 and 4.2876, plays an important role in simulation of gassolid fluidization (Yang et al., 2003), and hence is employed in a similar way to reflect the heterogeneous structures. The modified LBE is thus written as
(23) 
Accordingly, the volumeaveraged gas velocity must be redefined as (Guo and Zhao, 2002; Guo et al., 2002)
(24) 
It should be stressed here that the gas velocity calculated from the modified LBE is the superficial gas velocity and needs to be converted into the actual gas velocity.
The hydrodynamic force acting on the fluid in a computational cell, , can be calculated by summing up the momentum transfer that occurs in all discrete directions as
(25) 
Then a preestimated drag force acting on solid particlek, , can be calculated according to the following equation (Hoomans et al., 1996; Li and Kuipers, 2003):
(26) 
Here the local interphase momentum transfer coefficient , the local gas velocity and the local voidage are all calculated by bilinear interpolation using the values of four surrounding grid nodes. In each computational cell (see Fig. 3), the preestimated total drag force acting on all the particles, , can be obtained as
(27) 
which is not necessarily equal to the hydrodynamic force acting on the fluid. To satisfy Newton’s third law, the preestimated drag force acting on solid particle k, , is then adjusted to
(28) 
where and are the x and ycomponents of the drag force acting on solid particle , respectively. The total drag force acting on all the particles, , can be obtained
(29) 
With this drag force correction, the value of the total drag force acting on all the particles automatically equals to the value of the hydrodynamic force acting on the fluid in each computational cell and the Newton’s third law, , is also automatically satisfied in Eqs. (24) and (29).
2.4 Computational algorithm and implementation
In traditional DPS, the drag force acting on an individual solid particle is first calculated by the empirical formula Eq. (26), then the total drag force acting on all the particles in each computational cell is calculated by Eq. (27), and the total drag force, as a source term of force, is added into the momentum equation for the fluid phase, which is finally solved to obtain the fluid phase velocity field.
The procedure in this work is slightly different from the case that the fluid phase velocity field is first obtained from the evolution of the modified LBE, the total drag force acting on all the particles is simultaneously calculated by Eq. (25), and then the total drag force is distributed to each particle in the computational cell with the corrections described in Eqs. (26)(28). In fact, the preestimated drag force acting on the solid particles plays the role of a weight function for this distribution. The stepbystep procedure for the lattice Boltzmann based discrete simulation is outlined in Fig.4, which mainly includes the following steps:

Input the initial data such as the size and geometry of the simulation domain (for both the solid and fluid phases), the specified boundary conditions (i.e. particle distribution, initial flow field, velocity profile, etc.). The wall boundary condition given in this paper is the bounceback scheme for fluid and specular reflections for solid particles (Wang et al., 2006).

Perform the statistical calculations of the average velocity and the porosity of solid phase at the fluid lattice nodes according to Eq. (8) and Eq. (11), respectively.

Perform the consecutive propagation and collision process over a discrete fluid lattice cell, and then evolve the fluid flow according to the modified LBE (Eq. (23)), finally obtain the gas flow field and pressure filed.

The local gas velocities, the local pressure gradient, the local porosity, and the local interphase momentum transfer coefficient at the center of the particle are calculated by bilinear interpolation using the values of the surrounding computational lattices. Compute from Eq. (25) and from Eqs. (26) and (27) at all lattice nodes. Substitute these values of and into Eq. (28), obtain the drag force acting on the each particle .

Integrate the particle motion equation (Eq. (16)) numerically to calculate the velocities and the trajectories of individual particles.

Detect the instantaneous collisions of solid particles according to new particle positions and velocities. If the instantaneous collisions occur, the postcollision velocities of individual particles are determined according to Eqs. (17) and (18).

Evolve the next time and go back to step 2 unless the set time step is reached.
3 Simulations and results
To validate the proposed model, two typical fluidization processes, namely, a single bubble injected into a fluidized bed at incipient fluidization condition and particle clustering in the riser of a circulating fluidization bed, are simulated and compared with published correlations and experimental data. To simplify the simulation, the gas weight is neglected, but the effective gravity acting on the particles is corrected as .
3.1 A single bubble injected into a fluidized bed at incipient fluidization
In this example, a pulse of gas is injected at the bottom of a fluidized bed at incipient fluidization, the gas velocity on both sides of the central orifice is set as the minimum fluidization velocity (background gas velocity) , and the injection lasts 0.006 seconds. The parameters used for the simulation are summarized in Table 1.
Item  Dimensional  Dimensionless 
Bed width  0.00675  25 
Bed height  0.027  100 
Orifice width  0.00081  3 
Solid particle diameter  0.2  
Cell or lattice size  1  
Time step  1  
Minimum fluidization velocity  0.0025  
Injection velocity  0.8  0.04 
Gas density  1.1795  1 
Solid density  930.0  788.5 
Gas viscosity  
Coefficient of restitution  0.9  0.9 
Particle number  15000  15000 
Figure 5 presents the snapshots of bubble formation at different times. We observe the growth of the spherical cap bubble from an initially small perturbation in the bulk of the bed (see Fig. 5(a)(e)) to its subsequent detachment and rising (see Fig. 5(f)(h)). Eventually, the bubble wake is followed when the kidneybubble detaches from the bottom of the bed (see Fig. 5(j)(n)). Finally, the bubbles arrive at the surface, and it lifts some particles from wake into the bubble and causes the breakup of the bubble (see Fig. 5(o)(p)). This phenomenon is in qualitative agreement with the previous simulation and experimental results of bubbling (Rowe and Yacono, 1976; Nieuwland et al., 1996). It should be noted that the rising bubble shape is more sphericalcap than the elliptic bubble shape in the experimental and numerical results of Bokkers et al. (2004) due to the smaller solid particle diameter.
Figure 6 shows that the gas phase leaves the roof of bubble, and the downward moving particles near the wall drag the gas to the bottom of the bubble where it reenters the bubble region, resulting in a pair of symmetrical vortices which are observed in the neighborhood of the rising bubble. In addition, the simulated fluidizing gas streamlines passing through a rising bubble are close to those predicted by the twophase model of a fluidized bed (Davidson and Harrison, 1963). These results illustrate that the modified LBE is capable of capturing accurately the detailed flow structure of the gas phase.
The volumeaveraged equivalent bubble diameter is defined as the diameter of a circle with equal area to the region where the gas voidage is higher than 0.85 (Nieuwland et al., 1996). The computed maximum bubble diameter at different injection velocities is shown in Fig. 7, together with the correlations of Cai et al. (1994). The results show that increases with the increasing of injection velocities (or equivalently overall superficial gas velocity), and the deviation is below 5%, suggesting that the simulation results are reasonable.
3.2 Riser behavior of fluidization systems
The riser hydrodynamics of circulating fluidized beds (CFB) have been investigated extensively in both experiments and simulations in past decades. Our simulation refers to a CFB established by Li and Kwauk (1994), for which a large number of experimental data have been collected. The 2D riser configuration has been reported enough to mimic the real system (Cammarata et al., 2003; Xie et al., 2008), and the comparison of planar 2D simulation results with cylindrical 3D experimental data could also be found widely in the literature to study gassolid fluidized beds (Wang et al., 2008; Lan et al., 2009; Lu et al., 2011; Chen et al., 2012; Li et al., 2012). Table 2 lists the parameters used in this present simulation.
Item  Dimensional  Dimensionless 
Bed width  0.0162  60 
Bed height  0.194  720 
Initial solids volume fraction  0.09007  0.09007 
Solid particle diameter  0.2  
Cell or lattice size  1  
Time step  1  
Gas velocity at the inlet  1.52  0.0338 
Gas density  1.1795  1 
Solid density  930.0  788.5 
Gas viscosity  
Coefficient of restitution  0.9  0.9 
Particle number  123921  123921 
Figure 8 shows some snapshots from this simulation. Apparently, radial and axial heterogeneous structures are formed spontaneously and gradually from the homogeneous initial state and appear all the time hereafter.
Figure 9 shows the evolution of the simulated solids flux with time and its comparison with experimental data of Li and Kwauk (1994). Obviously, the simulated solids fluxes are in good agreement with the experimental value .
Figure 10 shows the comparison of the axial voidage profiles between the simulated results and the experimental data of Li and Kwauk (1994). The simulated profiles are calculated based on timeaveraged from 4.0s6.6s. The predicted sigmoid distribution of voidage in the axial direction is in reasonable agreement with experimental results, and the deviation may be partly ascribed to the unrealistic setup of inlet and outlet boundary conditions in the 2D simulation.
Figure 11 shows the comparison of the radial voidage profiles between the simulated results and the experimental correlations proposed by Tung et al. (1988). The simulated profiles are also calculated based on timeaveraged from 4.0s6.6s, which are in good agreement with empirical correlations in the core region, but overpredict the voidage in the annulus region. The coexistence of a dense annulus and a dilute core can be recognized, and the radial profile in the top region is larger than that in the bottom region. This implies that we have successfully captured a region with the coexistence of a dilute section at the top and a dense section at the bottom.
4 Discussion and conclusions
In this paper, a new numerical modeling method is presented for discrete particle simulation of gassolid fluidization; while the lattice Boltzmann method is applied to model the gas flow at a scale larger than the particles, the timedriven hardsphere model is employed to describe the motion of individual particles and the EMMS drag is used to correct gassolid interaction. The bubble formation of a single jet into the fluidized bed and the heterogeneous flow structures in the riser of a circulating fluidized bed are simulated successfully by the proposed approach.
Although the simulated results of the proposed approach seem reasonable, the current investigation is still very preliminary. The following aspects need to be clarified further:

The corresponding macroscopic equations of the modified LBE. At present, the immersed boundary method is introduced to the framework of LBM by an additional collision term to ensure that the LBE varies smoothly between the nodes occupied by pure fluid or by pure solid. This is a convenient approximation but lacks rigorous theoretical derivation. The exact form of the LBE presenting the volumeaveraged NS equation for the gas flow in DPS is not yet to be established.

Gassolid coupling. The drag force between gassolid may be still a source of controversy in DPS due to the uncertainty and the variety of different options. Although the present work only considers noslip discrete force () with the EMMS drag correction factor () as a source term (Eq. (23)), it is very flexible to incorporate other source terms or corrections into the framework of LBE. In addition to , which accounts for the effect of nonuniform particle distribution, the gas velocity distribution below the lattice scale, which is fully resolved in DNS but neglected in this work, may have to be accounted for in a simplified way.

Fully threedimensional simulation. For simplicity of implementation, the current simulation lies in the restriction of twodimensional fluid and threedimensional individual spherical particles, and then it can be regarded as a quasithreedimensional simulation. Further, its extension to threedimensions is straightforward in principle, namely using the threedimensional lattice Boltzmann models (i.e. D3Q15, D3Q19 and D3Q27) (d’Humieres et al., 2002) for gas flows. For threedimensional physical problems of interest, despite the relatively high numerical efficiency of LBM and TDHS, excessive computational cost is still required. Naturally, we will resort to parallel computation due to the advantage of the natural parallelisms inherent in LBM and TDHS. Therefore, implementing threedimensional massive parallel computation of the proposed model is highly desirable in future.

Turbulence. Twodimensional LES just provides a starting point for modeling high Re number gas flows and probably makes no sense here since the turbulence is inherently threedimensional in fluidized beds.
In closing, LBM preserves the advantages of explicit methods in terms of fast speed and parallelism. While, timedriven hardsphere model, which can be regarded as a combination of molecular dynamics (MD) and direct simulation Monte Carlo (DSMC) (Bird, 1994), inherits the virtue of computational efficiency in DSMC and physical picture in MD. Therefore, it can be expected that the proposed approach, combining LBM and TDHS, should be a fast numerical method, especially for highly parallel computing. In addition, we would like to point out that the proposed LBM based fluid flow solver can be combined with other solid phase solvers, such as DEM, DSMC (Liu and Lu, 2009), smoothed particle hydrodynamics (SPH) (Gingold and Monaghan, 1977; Monaghan, 1992; Xiong et al., 2011), and particleincell (PIC) (Snider, 2001). We would like to extend the present method to complex cases involving nonspherical particles in the future work.
Notation
{deflist}[A]\defitem\deftermlattice speed, m/s \defitemc\deftermunit velocity of lattice, m/s\defitem\deftermdrag coefficient, kg/(s) \defitem\deftermspeed of sound, m/s \defitem\deftermSmagorinsky constant number \defitem\deftermsolid particle diameter, m \defitem\deftermcoefficient of restitution \defitem\deftermdistribution functions \defitem\deftermforce, m/ \defitem\deftermparticle equilibrium distribution function \defitem\deftermgravitational acceleration, m/ \defitem\deftermcell or lattice size, m \defitem\deftermbed height, m \defitem\deftermmass of solid particle, kg \defitem\deftermnumber of particles in a cell \defitem\deftermtotal number of solid particles in the simulation \defitem\deftermposition of solid particle, m \defitem\deftermmomentum fluxes \defitem\deftermradius of solid particle, m \defitem\deftermorifice width, m \defitem\deftermparticle Reynolds number \defitem\deftermstrain rate tensor, \defitem\deftermtime, s \defitem\deftermtime step, s \defitem\deftermgas velocity, m/s \defitem\deftermparticle velocity, m/s \defitem\defterminjection velocity, m/s \defitem\deftermgas velocity at the inlet, m/s \defitem\deftermminimum fluidization velocity, m/s \defitem\deftermvolume of particle, \defitem\deftermbed width, m \defitem\deftermposition of lattice, m \defitem\deftermgrid size, m
Greek letters
\defitem\deftermdrag coefficient, kg/(s) \defitem\deftermvoidage \defitem\deftermgas volume fractions \defitem\deftermsolid volume fractions \defitem\deftermtwo dimensional porosity \defitem\deftermthree dimensional porosity \defitem\deftermweighting function \defitem\deftermgas viscosity, kg/(ms) \defitem\deftermkinematic viscosity, /s \defitem\defterminitial solids volume fraction \defitem\deftermgas density, kg/ \defitem\deftermsolid density, kg/ \defitem\deftermlattice weight or drag correction factor \defitem\deftermcollision term
Subscripts
\defitem\deftermthe numbers of solid particles \defitem\deftermcollion \defitem\deftermdrag \defitem\deftermturbulent eddy or average equivalent \defitem\deftermgas phase \defitem\deftermthe ith direction of lattice discrete velocity \defitem\deftermthe kth solid particle \defitem\deftermparticle \defitem\deftermpreestimated \defitem\deftermsolid phase \defitem\deftermtotal \defitem\deftermxdirection \defitem\deftermydirection
Abbreviations
\defitemBGK\defterm BhatnagarGrossKrook \defitemCFB\defterm Circulating Fluidized Beds \defitemCFBC\defterm Circulating Fluidized Bed Combustion \defitemCFD\defterm Computational Fluid Dynamics \defitemDEM\defterm Discrete Element Method \defitemDPS\defterm Discrete Particle Simulation \defitemDNS\defterm Direct Numerical Simulation \defitemDSMC\defterm Direct Simulation Monte Carlo \defitemD2Q9\defterm 2Dimensional 9Velocity Lattice Boltzmann Model \defitemD3Q15\defterm 3Dimensional 15Velocity Lattice Boltzmann Model \defitemD3Q19\defterm 3Dimensional 19Velocity Lattice Boltzmann Model \defitemD3Q27\defterm 3Dimensional 27Velocity Lattice Boltzmann Model \defitemEMMS\defterm Energy Minimization MultiScale \defitemFCC\defterm Fluid Catalytic Cracking \defitemLGA\defterm Lattice Gas Automata \defitemLBE\defterm Lattice Boltzmann Equation \defitemLBM\defterm Lattice Boltzmann Method \defitemLES\defterm Large Eddy Simulation \defitemMD\defterm Molecular Dynamics \defitemPBE\defterm Population Balance Equation \defitemPIC\defterm ParticleInCell \defitemQBMM\defterm QuadratureBased Moment Methods \defitemSGS\defterm SubGrid Scale \defitemSPH\defterm Smoothed Particle Hydrodynamics \defitemTDHS\defterm TimeDriven HardSphere \defitemTFM\defterm TwoFluid Model
Acknowledgement
This work is financially supported by the National Natural Science Foundation of China under Grants No. 21106155, the Specialized Fund from the Youth Innovation Promotion Association of the Chinese Academy of Science and the Chinese Academy of Sciences under Grant No. XDA07080303. We would like to thank Prof. Junwu Wang for his valuable suggestions and Dr. Qingang Xiong for useful discussions on this work.
References
References
 Anderson (T.B.) Anderson, T.B., Jackson, R., 1967. Fluid mechanical description of fluidized beds. Industrial Engineering & Chemistry Fundamentals 6, 527539.
 Aidun (C.K.) Aidun, C.K., Clausen, J.R., 2010. LatticeBoltzmann method for complex flows. Annual Review of Fluid Mechanics 42, 439472.
 Bhatnagar (P.L.) Bhatnagar, P.L., Gross, E.P., Krook, M., 1954. A model for collision processes in gases. I. Small amplitude processes in charged and neutral onecomponent systems. Physical Review 94, 511525.
 Bird (G.A.) Bird, G.A., 1994. Molecular gas dynamics and the direct simulation of gas flows. Oxford University Press, Oxford.
 Bokkers (G.A.) Bokkers, G.A., Annaland, M.V.S., Kuipers, J.A.M., 2004. Mixing and segregation in a bidisperse gassolid fluidised bed: a numerical and experimental study. Powder Technology 140, 176186.
 Cai (P.) Cai, P., Schiavetti, M., De Michele, G., Grazzini, G., Miccio, M., 1994. Quantitative estimation of bubble size in PFBC. Powder Technology 80, 99109.
 Cammarata (L.) Cammarata, L., Lettieri, P., Micale, G.D.M., Colman, D., 2003. 2D and 3D CFD simulations of bubbling fluidized beds using EulerianEulerian Models. International Journal of Chemical Reactor Engineering 1, 114.
 Chapman (S.) Chapman, S., Cowling, T.G., 1970. The mathematical theory of nonuniform gases. Cambridge University Press, London.
 Chen (C.) Chen, C., Li, F., Qi, H., 2012. Modeling of the flue gas desulfurization in a CFB riser using the Eulerian approach with heterogeneous drag coefficient. Chemical Engineering Science 69, 659668.
 Chen (S.) Chen, S., Doolen, G.D., 2003. Lattice Boltzmann method for fluid flows. Annual Review of Fluid Mechanics 30, 329364.
 Chen (S.) Chen, S., Shi, B., Liu, Z., He, Z., Guo, Z.L., Zheng, C., 2004. LatticeBoltzmann simulation of particleladen flow over a backwardfacing step. Chinese Physics 13, 16571664.
 Chu (K.) Chu, K., Yu, A., 2008a. Numerical simulation of complex particlefluid flows. Powder Technology 179, 104114.
 Chu (K.) Chu, K., Yu, A., 2008b. Numerical simulation of the gassolid flow in threedimensional pneumatic conveying bends. Industrial & Engineering Chemistry Research 47, 70587071.
 Chu (K.) Chu, K., Wang, B., Yu, A., Vince, A., Barnettc, G.D., Barnettc, P.J., 2009a. CFDDEM study of the effect of particle density distribution on the multiphase flow and performance of dense medium cyclone. Minerals Engineering 22, 893909.
 Chu (K.) Chu, K., Wang, B., Yu, A., Vince, A., 2009b. CFDDEM modelling of multiphase flow in dense medium cyclones. Powder Technology 193, 235247.
 Chu (K.) Chu, K., Wang, B., Xu, D., Chen, Y., Yu, A., 2011. CFDDEM simulation of the gassolid flow in a cyclone separator. Chemical Engineering Science 66, 834847.
 Cook (B.K.) Cook, B.K., Noble, D.R., Williams, J.R. 2004. A direct simulation method for particlefluid systems. Engineering Computations 21, 151168.
 Cundall (P.) Cundall, P., Strack, O., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 4765.
 Davidson (J.F.) Davidson, J.F., Harrison, D., 1963. Fluidised particles. Cambridge University Press, London.
 Darabi (P.) Darabi, P., Pougatch, K., Salcudean, M., Grecov, D., 2011. DEM investigations of fluidized beds in the presence of liquid coating. Powder Technology 214, 365374.
 Deen (N.G.) Deen, N.G., Van Sint Annaland, M., Van der Holf, M.A., Kuipers, J.A.M., 2007. Review of discrete particle modeling of fluidized beds, Chemical Engineering Science 62, 2844.
 Desjardins (O.) Desjardins, O., Fox, R., Villedieu, P., 2008. A quadraturebased moment method for dilute fluidparticle flows. Journal of Computational Physics 227, 25142539.
 eng (Y.) Feng, Y., Han, K., Owen, D. 2007. Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows: Computational issues. International Journal for Numerical Methods in Engineering 72, 11111134.
 Filippova (O.) Filippova, O., Hanel, D., 1997. LatticeBoltzmann simulation of gasparticle flow in filters. Computers & Fluids 26, 697712.
 d’Humieres (D.) d’Humieres, D., Ginzburg, I., Krafczyk, M., Lallemand, P., Luo, L., 2002. Multiplerelaxationtime lattice Boltzmann models in three dimensions. Philosophical transactions of the Royal Society of London A 360, 437451.
 Fox (R.) Fox, R., 2008. A quadraturebased thirdorder moment method for dilute gasparticle flows. Journal of Computational Physics 227, 63136350.
 Fox (R.) Fox, R.O., 2009a. Higherorder quadraturebased moment methods for kinetic equations. Journal of Computational Physics 228,77717791.
 Fox (R.) Fox, R.O., 2009b. Optimal moments sets for multivariate direct quadrature method of moments. Industrial & Engineering Chemistry Research 48,96869696.
 Fox (R.) Garg, R., Galvin, J., Li, T., Pannala, S., 2012. Opensource MFIXDEM software for gassolids flows: Part IVerification studies. Powder Technology 220, 122137.
 Anderson (T.B.) Geng, Y., Che, D., 2011. An extended DEMCFD model for char combustion in a bubbling fluidized bed combustor of inert sand. Chemical Engineering Science 66, 207219.
 Anderson (T.B.) Gidaspow, D., 1994. Multiphase flow and fluidization: continuum and kinetic theory descriptions. Academic Press, New York.
 Anderson (T.B.) Gidaspow. D., Jiradilok.V., 2009. Computational techniques: the multiphase CFD approach to fluidization and green energy technologies, Nova Science Publishers, New York.
 Anderson (T.B.) Gingold, R., Monaghan, J., 1977. Smoothed particle hydrodynamics: theory and application to nonspherical stars. Monthly Notices of the Royal Astronomical Society 181, 375389.
 Anderson (T.B.) Goniva, C., Kloss, C., Deen, N.G., Kuipers, J.A.M., Pirker, S., 2012. Influence of rolling friction on single spout fluidized bed simulation. Particuology 10, 582591.
 Anderson (T.B.) Gopalakrishnan, P., Tafti, D., 2013. Development of parallel DEM for the open source code MFIX. Powder Technology 35, 3341.
 Anderson (T.B.) Gui, N., Fan, J., Luo, K., 2008. DEMLES study of 3D bubbling fluidized bed with immersed tubes. Chemical Engineering Science 63, 36543663.
 Anderson (T.B.) Guo, Z.L., Zhao, T.S., 2002. Lattice Boltzmann model for incompressible flows through porous media. Physical Review E 66, 036304.
 Anderson (T.B.) Guo, Z.L., Zheng, C.G., Shi, B.C., 2002. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Physical Review E 65, 046308.
 Anderson (T.B.) Hoomans, B.P.B., Kuipers, J.A.M., Briels, W.J., Van Swaaij, W.P.M., 1996. Discrete particle simulation of bubble and slug formation in a twodimensional gasfluidised bed: A hardsphere approach. Chemical Engineering Science 51, 99118.
 Anderson (T.B.) Hopkins, M.A., Louge, M.Y., 1991. Inelastic microstructure in rapid granular flows of smooth disks. Physics of Fluids A: Fluid Dynamics 3, 4757.
 Anderson (T.B.) Hu, H., Joseph, D. D., Crochet, M.J., 1992. Direct simulation of fluid particle motion. Theorectical and Computational Fluid Dynamics 3,285306.
 Anderson (T.B.) Ishii, M., 1975. Thermofluid dynamic theory of twophase flow. Eyrolles, Paris.
 Anderson (T.B.) Krafczyk, M., Tolke, J., Luo, L.S., 2003. Largeeddy simulations with a multiplerelaxationtime LBE model. International Journal of Modern Physics B 17, 3340.
 Anderson (T.B.) Kuang, S., Chu, K., Yu, A., Zou, Z., Feng, Y., 2008. Computational investigation of horizontal slug flow in pneumatic conveying. Industrial & Engineering Chemistry Research 47, 470480.
 Anderson (T.B.) Lan, X., Xu, C., Wang, G., Wu, L., Gao, J., 2009, CFD modeling of gassolid flow and cracking reaction in twostage riser FCC reactors. Chemical Engineering Science 64, 38473858.
 Anderson (T.B.) Li, F., Song, F., Benyahia, S., Wang, W., Li, J., 2012. MPPIC simulation of CFB riser with EMMSbased drag model. Chemical Engineering Science 82, 104113.
 Anderson (T.B.) Li, J., Kuipers, J., 2003. Gasparticle interactions in dense gasfluidized beds. Chemical Engineering Science 58, 711718.
 Anderson (T.B.) Li, J., Kwauk, M., 1994. Particlefluid twophase flow: the energyminimization multiscale method. Metallurgical Industry Press, Beijing.
 Anderson (T.B.) Li, T., Guenther, C., 2012. MFIXDEM simulations of change of volumetric flow in fluidized beds due to chemical reactions. Powder Technology 220, 7078.
 Anderson (T.B.) Li, T., Garg, R., Galvin, J., Pannala, S., 2012a. Opensource MFIXDEM software for gassolids flows: Part IIValidation studies. Powder Technology 220, 138150.
 Anderson (T.B.) Li, T., Gopalakrishnan, P., Garg, R., Shahnam, M., 2012b. CFDDEM study of effect of bed thickness for bubbling fluidized beds. Particuology 10, 532541.
 Anderson (T.B.) Liu, H., Lu, H., 2009. Numerical study on the cluster flow behavior in the riser of circulating fluidized beds. Chemical Engineering Journal 150, 374384.
 Anderson (T.B.) Lu, B., Wang, W., Li, J., 2011. Eulerian simulation of gassolid flows with particles of Geldart groups A,B and D using EMMSbased mesoscale model. Chemical Engineering Science 66, 46244635.
 Anderson (T.B.) Ma, J., Ge, W., Wang, X., Wang, J., Li, J., 2006. Highresolution simulation of gassolid suspension using macroscale particle methods. Chemical Engineering Science 61, 70967106.
 Anderson (T.B.) Mazzei, L., 2011. Limitations of quadraturebased moment methods for modeling inhomogeneous polydisperse fluidized powders. Chemical Engineering Science 66, 36283640.
 Anderson (T.B.) McNamara, G.R., Zanetti, G., 1988. Use of the Boltzmann equation to simulate latticegas automata. Physical Review Letters 61, 23322335.
 Anderson (T.B.) Mikami, T., Kamiya, H., Horio, M., 1998. Numerical simulation of cohesive powder behavior in a fluidized bed. Chemical Engineering Science 53,19271940.
 Anderson (T.B.) Monaghan, J., 1992. Smoothed Particle Hydrodynamics. Annual Review of Astronomy and Astrophysics 30, 543574.
 Anderson (T.B.) Nieuwland, J.J., Veenendaal, M.L., Kuipers, J.A.M., Van Swaaij, W.P.M., 1996. Bubble formation at a single orifice in gasfluidised beds. Chemical Engineering Science 51, 40874102.
 Anderson (T.B.) Noble, D.R., Torczynski, J.R., 1998. A latticeBoltzmann method for partially saturated computational cells. International Journal of Modern Physics C 9, 11891201.
 Anderson (T.B.) Ouyang, J., Li, J., 1999a. Particlemotionresolved discrete model for simulating gassolid fluidization. Chemical Engineering Science 54, 20772083.
 Anderson (T.B.) Ouyang, J., Li, J., 1999b. Discrete simulations of heterogeneous structure and dynamic behavior in gassolid fluidization. Chemical Engineering Science 54, 54275440.
 Anderson (T.B.) Qian, Y., d’Humieres, D., Lallemand, P., 1992. Lattice BGK models for NavierStokes equation. Europhysics Letters 17, 479484.
 Anderson (T.B.) Rowe, P.N., Yacono, C.X.R., 1976. The bubbling behaviour of fine powders when fluidised. Chemical Engineering Science 31, 11791192.
 Anderson (T.B.) Schiller, L., Naumann, A., 1935. A drag coefficient correlation. Vdi Zeitung 77, 318320.
 Anderson (T.B.) Smagorinsky, J., 1963. General circulation experiments with the primitive equations. Monthly Weather Review 91, 99164.
 Anderson (T.B.) Snider, D., 2001. An incompressible threedimensional multiphase particleincell model for dense particle flows. Journal of Computational Physics 170, 523549.
 Anderson (T.B.) Sterling J.D., Chen S., 1996. Stability analysis of lattice Boltzmann methods. Journal of Computational Physics 123, 196206.
 Anderson (T.B.) Su, J., Gu, Z., Xu, X., 2011. Discrete element simulation of particle flow in arbitrarily complex geometries. Chemical Engineering Science 66, 60696088.
 Anderson (T.B.) Sungkorn, R., Derksen, J.J., Khinast, J.G., 2011. Modeling of turbulent gasliquid bubbly flows using stochastic Lagrangian model and latticeBoltzmann scheme. Chemical Engineering Science 66, 27452757.
 Anderson (T.B.) Tsuji, Y., Kawaguchi, T., Tanaka, T., 1993. Discrete particle simulation of twodimensional fluidizedbed. Powder Technology 77, 7987.
 Anderson (T.B.) Tung, Y., Li, J., Kwauk, M., 1988. Radial voidage profiles in a fast fluidized bed, in: Kwauk, M., Kunii, D. (Eds.), Fluidization’88: Science and Technology. Science Press, Beijing, pp. 139145.
 Anderson (T.B.) Wang, J., Ge, W., Li, J., 2008. Eulerian simulation of heterogeneous gassolid flows in CFB risers: EMMSbased subgrid scale model with a revised cluster description.Chemical Engineering Science 63, 15531571.
 Anderson (T.B.) Wang, L., Ge, W., Li, J., 2006. A new wall boundary condition in particle methods. Computer Physics Communications 174, 386390.
 Anderson (T.B.) Wang, L., Zhou, G., Wang, X., Xiong, Q., Ge, W., 2010. Direct numerical simulation of particlefluid systems by combining timedriven hardsphere model and lattice Boltzmann method. Particuology 8, 379382.
 Anderson (T.B.) Wang, S., Lu, H., Li, X., Wang, J., Zhao, Y., Ding, Y., 2009. Discrete particle simulations for flow of binary particle mixture in a bubbling fluidized bed with a transport energy weighted averaging scheme. Chemical Engineering Science 64, 17071718.
 Anderson (T.B.) Wang, S., Li, X., Lu, H., Liu, G., Wang, J., Xu, P., 2011. Simulation of cohesive particle motion in a soundassisted fluidized bed. Powder Technology 207, 6577.
 Anderson (T.B.) Wei, M., Wang, L., Li, J., 2013. Unified stability condition for particulate and aggregative fluidizationExploring energy dissipation with direct numerical simulation. Particuology 11, 232241.
 Anderson (T.B.) Wu, C., Berrouk, A.S., Nandakumar, K., 2009. Threedimensional discrete particle model for gassolid fluidized beds on unstructured mesh. Chemical Engineering Journal 152, 514529.
 Anderson (T.B.) Wu, C., Cheng, Y., Ding, Y., Jin, Y., 2010. CFDDEM simulation of gassolid reacting flows in fluid catalytic cracking (FCC) process. Chemical Engineering Science 65, 542549.
 Anderson (T.B.) Xie, N., Battaglia, F., Pannala, S., 2008. Effects of using two versus threedimensional computational modeling of fluidized beds: part I, hydrodynamics. Powder Technology 182, 113.
 Anderson (T.B.) Xiong, Q., Deng, L., Wang, W., Ge, W., 2011. SPH method for twofluid modeling of particlefluid fluidization. Chemical Engineering Science 66, 18591865.
 Anderson (T.B.) Xiong, Q., Li, B., Zhou, G., Fang, X., Xu, J., Wang, J., He, X., Wang, X., Wang, L., Ge, W., Li, J., 2012. Largescale DNS of gassolid flows on Mole8.5. Chemical Engineering Science 71, 422430.
 Anderson (T.B.) Xu, B., Yu, A., 1997. Numerical simulation of the gassolid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chemical Engineering Science 52, 27852809.
 Anderson (T.B.) Xu, M., Ge, W., Li, J., 2007. A discrete particle model for particlefluid flow with considerations of subgrid structures. Chemical Engineering Science 62, 23022308.
 Anderson (T.B.) Yang, N., Wang, W., Ge, W., Li, J., 2003. CFD simulation of concurrentup gassolid flow in circulating fluidized beds with structuredependent drag coefficient. Chemical Engineering Journal 96, 7180.
 Anderson (T.B.) Yu, H., Girimaji, S.S., Luo, L., 2005. DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method. Journal of Computational Physics 209, 599616.
 Anderson (T.B.) Zhao, X., Li, S., Liu, G., Yao, Q., Marshall, J.S., 2008. DEM simulation of the particle dynamics in twodimensional spouted beds. Powder Technology 184, 205213.
 Anderson (T.B.) Zhao, Y., Jiang, M., Liu, Y., Zheng, J., 2009. Particlescale simulation of the flow and heat transfer behaviors in fluidized bed with immersed tube. AIChE Journal 55, 31093124.
 Anderson (T.B.) Zhao, Y., Ding, Y., Wu, C., Cheng, Y., 2010. Numerical simulation of hydrodynamics in downers using a CFDDEM coupled approach. Powder Technology 199, 212.
 Anderson (T.B.) Zhou, H., Flamant, G., Gauthier, D., 2004a. DEMLES of coal combustion in a bubbling fluidized bed. Part I: gasparticle turbulent flow structure. Chemical Engineering Science 59, 41934203.
 Anderson (T.B.) Zhou, H., Flamant, G., Gauthier, D., 2004b. DEMLES simulation of coal combustion in a bubbling fluidized bed Part II: coal combustion at the particle level. Chemical Engineering Science 59, 42054215.
 Anderson (T.B.) Zhu, H., Zhou, Z., Yang, R., Yu, A., 2007. Discrete particle simulation of particulate systems: theoretical developments, Chemical Engineering Science 62, 33783396.
 Anderson (T.B.) Zhu, H., Zhou, Z., Yang, R., Yu, A., 2008. Discrete particle simulation of particulate systems: A review of major applications and findings. Chemical Engineering Science 63, 57285770.