Discrete particle simulation of gas-solid flows (From dilute to dense flows)

Yutaka Tsuji Department of Mechanical Engineering, Faculty of Engineering Osaka University* Increasing computer power has made discrete particle simulation a practical means for predicting particle motion in industrial processes. Simulation techniques for various particle motions, from dilute to dense phase flows, are explained here. In simulations of this kind, the most important problem is how to model the interactions of particles with walls, with other particles and with the fluid. Calculations based on various models for these interactions are discussed. Particle segregation phenomena, particle mixing in a rotating vessel, dense phase pneumatic conveying, and fluidized bed are shown as examples of simulations which have been made in the author's laboratory.


Introduction
Since numerical simulation of single phase flows is considered to have attained maturity, attention is being paid to multiphase flows. Fields in which such numerical simulations will play an important role as a prediction method are widespread. For example, flows in manufacturing industries, flows in nature and flows in living bodies are often multiphase. In the simulation of multiphase flows, we must treat both the fluid phase and discrete phase motions. Owing to the great achievements in single phase flows, computation techniques for fluid motion have been established satisfactorily. Now techniques for particle motion are being developed and most techniques can be applied to powder technology.
Multiphase flows are classified into the following four categories; gas-liquid, liquid-solid, gas-solid and gas-liquid-solid. As far as formulation is concerned, the gas-solid flow is among the simplest, because mass transfer and surface phenomena which make the formulation complicated need not usually be taken into account. Simulation of gas-solid flows started from dilute phase flows where particles are dispersed in gas at a low concentration. In the recent several years, the author's research group in Osaka University has been interested in flows with high concentrations. The higher the particle concentration, the stronger the impact of simulation on practical applications. For instance, dense phase flows are closely * Suita, Osaka 565 Japan t Received July 8,1993 KONA No.ll (1993) related to powder technology. Before, it seemed a dream for engineers working in powder technology to design machines based on computer simulation, but the situation is gradually changing. Nowadays computers enable us to deal with the motion of individual particles in dense phase flows where particles are in contact with each other.

Eulerian vs. Lagrangian
According to the wellknown classification, numerical methods of multiphase flows are classified as Eulerian or Lagrangian depending on the treatment of the particulate phase. The Eulerian method regards the particulate phase as a continuum. This method is further divided into one-fluid (homogeneous) and two-fluid models. An advantage of the Eulerian method is to apply techniques developed for single phase flows to the particulate phase with fewer modifications. However, it should be noted that a sufficient number of particles must be included uniformly in a control volume for calculation in order to use the Eulerian method. Unfortunately, cases satisfying this condition are limited.
The Lagrangian method computes the motion of individual particles. Its basic ideas and formulations are simple compared with the Eulerian method. However, as the number of particles increases, large memories and long computation time are needed. If the Lagrangian calculation is to be performed rigorously for all the particles in a practical multi-phase flow, enormous amount of memory and computation time would be required. This is a substantial disadvantage of the Lagrangian method. Fortunately, progress in computer hardware development is still continuing, which is helpful to the Lagrangian method. Also, our group has made an attempt to apply the DSMC (Direct Simulation Monte Carlo) method, developed in molecular dynamics, to multi-phase flows to avoid the above demerit of the Lagrangian method.
The above classification has been applied to dispersed phase flows in general. The same classification holds true for dense phase flows or granulate flows where the particle concentration is so dense that the particles are in continuous contact with each other. When the Eulerian method is used for such dense phase flows, the most important assumption which has to be made is about the relation between stress and strain of the continuum consisting of particles. Many efforts have been made to obtain this relation, experimentally and theoretically. However, it is not easy to deduce relations for general use.
It has not been long since the Lagrangian method was applied to dense phase flows. Compared with the Eulerian method, the Lagrangian method requires less assumptions and is expected to be a powerful means for predicting phenomena in powder technology.

Irregular bouncing
In dilute phase pneumatic conveying, the gas flow gives momentum and kinetic energy to the particles. That is, particles are accelerated by gas through the drag force. If a particle in a horizontal pipe did not collide with the wall, the particle velocity would be equal to the gas velocity. However the time-averaged particle velocity is always smaller than the gas velocity, because the particles lose momentum due to collisions with the wall. Since the additional pressure drop due to the particles corresponds to the fluid drag which is related to the velocity difference between the fluid and the particles, the energy loss caused by collision must be estimated precisely in the simulation.
Although particle-to-wall collision is an important phenomenon in gas-solid pipe flows, the Eulerian method has disadvantages concerning collision because the effects of collision can not be considered in a direct way. In contrast, the Lagrangian method is suitable for dealing with the collision problem. When the coefficients of restitution and dynamic friction are given as particle physical properties, the collision problem can be solved using the classical 58 dynamics of a rigid body. However, the problem is not as simple as expected. That is, if the simulation were based on the simple assumption that spherical particles collide with a flat plate with a coefficient of restitution smaller than one, particles colliding repeatedly with the wall in a long horizontal pipe would ultimately slide on the wall. To avoid such an unrealistic result, we must take the irregularity of collisions into consideration. If the irregularity is neglected, particles can not be suspended in the gas. Previously, Magnus or shear lift forces and fluid turbulence were often considered as mechanisms suspending particles in a horizontal pipe. It is true that these fluid dynamic forces act on particles and that the effect of these forces are dominant for small particles, but the mechanism suspending large particles in horizontal pipes is the irregular bouncing of the particles against the wall. The irregularity is caused by particle shape deviating from a true sphere and roughness of the wall.
Regarding the irregular bouncing, several models have been proposed so far (Matsumoto et al. 1970a(Matsumoto et al. , 1970b, Tsuji eta!. 1987, Tsuji et a!. 1989, Frank et a!. 1991). In general the relation between particle velocities before and after collision are derived using the impulsive equations with the coefficient of restitution and friction coefficient given. We (Tsuji et a!., 1991) attributed the irregularity to the non-sphericity of particle shape, and demonstrated a calculation method for three dimensional collision of a nonspherical particle. This method is applicable to particles of arbitrary shape.
We calculated the trajectories of nearly spherical particles in a two-dimensional channel. Fig. 1 shows the shape of the particle. Calculated results of particle trajectories are shown in Fig. 2, where the scale of the longitudinal distance is greatly reduced so that the motion of the particles in the whole channel can be Restitution coefficient e = 0.93, Friction coefficient J1 r = 0.28 seen. It was found that even a small deviation from the sphere leads to a considerable irregular bouncing motion.
3.2 Inter-particle collision As the particle concentration becomes higher, interparticle collision occurs more often and the loss of kinetic energy of the particles due to this collision cannot be neglected. Fortunately, as long as the particulate phase is dispersed, it is sufficient to consider only simple collision and not multiple collision. As in the particle-to-wall collision, the impulsive equations are solved to obtain the translation and angular velocities after collision assuming that the coefficients of restitution e and friction p, f are known (Tanaka et a!., 1991).

DSMC (direct simulation monte carlo) method
The trajectories of individual particles are obtained by numerically integrating the equations of particle motion with a time interval M. This integration process is straightforward. As described in the above section 3.2, the particle velocities after binary collision are calculated from the values before collision and the physical properties such as the coefficients of restitution and friction. However it is not easy to find pairs of particles colliding with each other when a large number of particles exist in the flow field. If we attempt to find such pairs from the trajectories of individual particles, the computation time would become extremely long as the number of particles increases. The DSMC method, proposed by Bird (1976) to solve the Boltzmann equation of rarefied gas flow has been developed to resolve this difficulty.
In dispersed gas-solid flow involving large particles, the particle inertia forces are dominant, and collisions against a particle or wall are assumed to occur instantaneously. Therefore the motion of the solid phase is similar to that of molecules in rarefied gas except that the kinetic energy of the particle fluctuat-KONA No.ll (1993) ing motion tends to decay due to the inelastic and frictional characteristics of collision.
The DSMC method has already been applied to calculate particle motion in gas-solid flows by Kitron et a!. (1990), Tanaka et a!. (1991) and Tanaka et al. (1993). Tanaka eta!. (1991) applied the DSMC method to fully developed flows in a vertical pipe and obtained results which agree satisfactorily with whose predicted by the deterministic method . The outline of the DSMC method is as follows. (i) Each simulated particle represents a large number of "physical" particles, (ii) The time internal for numerical integration is taken sufficiently small compared with the mean free time of the particles, and thus inter-particle collisions are uncoupled from free particle motion. (iii) The flow field is divided into small cells over which change in flow properties is small. Particles are allowed to collide through a Monte Carlo procedure. (iv) The calculation of particle motion is carried out by repeating the following procedure: First, the motions of all simulated particles in the time interval M are calculated using the equations of motion without regard for interparticle collisions. If the calculated path of a particle crosses a solid wall, the velocity is replaced with the post-rebound velocity. Secondly, inter-particle collision during the time interval .6t is searched by means of the Monte Carlo procedure. The collision counterpart and geometry of collision are also chosen with the Monte Carlo method.
If a particle collides with another particle, the post -collision velocities of the colliding pair are calculated using the impulsive equations as described in Section 3.2. The particle velocities are replaced by the post-collision velocities, but without changing their positions. Several schemes to find out whether a particle collides with another particle within a time interval have been proposed. We used the modified Nanbu method (Illner and Neunzert, 1987).The probability of inter-particle collision for a particular particle is evaluated from the particle properties included in the cell. The collision probability Pi of particle i during a time interval .6 t is given by, (1) where N is the number of simulated particles in the cell and Pij , the probability of collision between the particle i and particle j in the cell during the time interval .6t.
In evaluating Pij, it must be noted that each simulated particle represents many physical particles. Pij is not a collision probability between a single particle i and a single particle j but that for a physical particle in a field over which physical particles represented by simulated particles i and j are distributed homogeneously. Assuming the particles are spheres having a uniform diameter d 5 , pii is given by, (2) Where n is the number density of the real flow at the corresponding position, and gii is the magnitude of relative velocity between both particles.
According to the modified Nanbu method, the occurrence of inter-particle collision and the collision counterpart are decided using a random number RND obtained from a uniform distribution which ranges from zero to unity. A ''candidate'' collision counterpart j is selected first using the following equation.
where [ [RND x N]] is defined as the integer part of RND x N. Particle i is then assumed to collide with particle k during the previous time interval if: (4) In this procedure, Pii for any particle combination in a cell must not exceed 1/N. This condition is satisfied by selecting an appropriate time interval, since Pij is proportional to .1 t as shown in Eq. (2).
If the particle i collides with the particle k, the velocity of particle i (though not that of k) is replaced by the velocity after the collision. The relative 0 position of the collision 1s also given using random numbers.

Cundall & Strack model
As the particle concentration increases, almost all the particles are in contact with neighboring ones. Since each particle has elasticity, a particle assembly forms a complicated vibration system having multiple degrees of freedom. Moreover the contact points vary with time. The influence of a particle on other particles far from it propagates in a disturbance wave. It is very difficult to consider the interaction between one particle and remote ones. If the time interval in numerical calculation is chosen sufficiently small, it can be assumed that during a single time interval disturbances do not propagate from any particle further than its immediate neighbors. In other words, the instantaneous motion of each particle is determined by the contact forces acting between that particle and particles with which it is in contact. This assumption, proposed by Cundall and Strack (1979), makes it possible to save memory space. Based on this assumption, the motion of each particle in dense phase flows can be obtained by integrating the equations of motion step by step in which the contact forces between a particle and its immediate neighbors are taken into account. However, unlike dispersed phase flows, a method such as DSMC can not be used and trajectory calculation must be performed for all the physical particles.
In this simulation, the contact forces are more important then the inter-particle collision. Cundall and Strack (1979) proposed the model shown in where o nii displacement of particle caused by the normal force ---> v rij velocity vector of particle i relative to particle j ---> n ii unit vector drawn from the center of particle i to that of particle j where kt and o tii are, respectively, the stiffness and displacement in the tangential direction. In the above equations, the suffixes n and t represent the components corresponding to the normal and tangential directions, respectively. ~sii is the slip velocity of the contact point, which is given by where r 5 is the radius of the sphere. If the following relation is satisfied ---> I f Ctij I > JAr I f Cnij I , (8) particle i slides and the tangential force is given by instead of Eq. (6). Eq. (9) is the Coulomb-type The same relations as the above equations are derived for the contact with the wall if particle j is replaced by the wall. In general, several particles are in contact with particle i at the same time. Therefore the total force and torque acting on particle i is obtained by taking the sum of the above forces with respect to j.

Determination of parameters in Cundall &
Strack model The next step after modeling the contact forces is to determine the values of the stiffness k, damping coefficient 17 and friction coefficient JAr • Among these parameters, the friction coefficient JAr is measurable and regarded as a parameter given empirically.
Fortunately the stiffness can be calculated using the Hertzian contact theory when the physical properties such as the Young's modulus and Poisson ratio are known. According to the Hertizian contact theory, the relation between the normal force P n and displacement 0 n is given by where r 5 Es radius Young's modulus of the particle Poisson ratio of the particle (12) The above equation signifies that the force varies with the 3/2 power of the displacement. Therefore if the above results are applied to the model of the contact forces given in the foregoing section, Eq. (5) can be replaced by (13) As for the relation between the tangential force Pt and displacement, theories developed by Mindlin (1949) and Mindlin and Deresiewicz (1953) are used (Tsuji eta!., 1992).
Regarding the damping coefficient, Cundall and Strack (1979) proposed two expressions, given by the following equations, (14) which are derived from the condition of the critical damping of a single degree-of-freedom system consisting of a mass, spring and dash-pot.
We (Tsuji et al., 1992) proposed another method for determining the damping coefficient. In our method, the damping coefficient is related to the coefficient of restitution which is regarded as one of particle's physical properties or can be measured in a simple experiment.

Equation of particle motion
Individual particles have two types of motion; translational and rotational motions. The translational motion is caused by the contact force, fluid force and gravitational force. The rotational motion in dispersed phase flows is caused by collision with the wall and particles. The rotational motion in the fluid is affected by viscous dissipation. In dense phase flows where particles are in contact with neighboring ones, the rotational motion is caused by the tangential components of the contact forces. The Regarding the fluid forces, drag, lift (due to rotation and velocity shear) and torque are considered. Various theoretical or empirical formulae are available for such forces acting on particles freely suspended in the fluid. In dense phase flows, the lift and torque due to viscosity can be neglected, and only the drag is taken into account. Such drag can be estimated by considering the pressure gradient caused by the fluid passing through the particle bed, as might be found from established empirical formulae such as the Ergun equation (Ergun, 1952).

Equation of fluid motion
It is almost impossible even for modern supercomputers to solve the problem of the instantaneous flow field of small scale between moving particles in multiphase flows. Furthermore, the scale of the phenomena of interest to us is much larger than the size of individual particles. Therefore it is reasonable to perform calculations based on locally averaged quantities according to Anderson and Jackson (1967). As is usual in many numerical calculations of flow fields, the finite difference method can be used to obtain the flow field in the simulation of multiphase flows. Therefore, the flow domain is divided into cells, the size of which is smaller than the macroscopic motion but larger than the particle size. All the quantities such as pressure p and velocity u are averaged in the cell using a weight function. The void fraction 8 of each cell can be defined by the number of particles existing in the cell.
The equation of continuity is given by (20) The equation of fluid motion is where p is the fluid density. The last term in Eq.
(21) denotes the mutual interaction between particles and fluid. Eq. (21) signifies that the fluid is assumed to be inviscid. The interaction term depends on the fluid properties, the void fraction and the relative velocity between particles. For details, refer to Tsuji eta!., (1993).

Calculation procedure
As an example of the calculation procedure, flow charts in the simulation of a fluidized bed are given in Fig. 4. 9. Examples of simulation 9.1 Without fluid Our group in Osaka University has developed computer programs for fluidized beds and pneumatic conveying, in which the fluid forces play important roles. However, simulations neglecting fluid flow also have many applications in powder technology. In fact we have attempted to predict various particulate flows such as particle discharge from hoppers and particle mixing in rotating vessels. In this section, such examples are shown first. Fig. 5 illustrates the mixing of two types of particles in a rotating square duct. This simulation was made three-dimensionally. From this simulation, the time required for mixing and the optimal rotating speed can be predicated. Fig. 6 illustrates a particle flow from a round hopper. Initially, two types of particles are uniformly mixed in the upper hopper. The dark color particles are heavy particles and the light color particles are light ones. As the particles are discharged, particle segregation is observed, that is, the heavy particles tend to occupy the lower and center parts of the lower hopper. Fig. 7 also illustrates particle flows from round hoppers. The effect of a cone near the gate on the flow pattern can be clearly observed. The case without the cone shows a funnel-type flow pattern while the case with the cone shows a mass flow. Fig. 8 (Tanaka et a!., 1993) illustrates a dispersed gas-solid flow in a vertical channel, where the DSMC method was used. Note that the vertical scale of the channel is compressed in these figures. When the superficial gas velocity is large or the mean solid mass flux is small, the particle concentration has a nearly homogeneous distribution. As the gas velocity decreases and the solid loading increases, of particles in dense clusters can even become negative, meaning that they drop down along the wall. Fig. 9 represents a flow pattern obtained in a channel wider than that of Fig. 8, giving rise to clusters forming not only along the walls but also in the central part of the channel. The close-up of a typical V -shaped cluster formed in the central part of the channel is also shown in Fig. 10. The scale of the close-up is not compressed. Fig. 11 (Tsuji et al., 1992) illustrates a plug flow based on simple models, the formation of bubbles and slugs is observed in the same way as in physical experiments. Fig. 12 represents an example of such simulations. In this case, the particle motion was treated two-dimensionally. The air is issued from three nozzles provided at the bottom. The author has not reviewed sufficiently the studies made by other researchers in this review. It should be noted that discrete particle simulations without fluids have been made actively in the field of soil and rock mechanics. power increases, this method will be all the more attractive and useful. As far as particle motion is concerned, we have already succeeded in obtaining qualitatively good correlations. It is expected that from now, research in discrete particle simulation will be developed along the two directions; seeking quantitative correlations or including various influential factors such as effects of humidity, cohesion forces, particle shape and so on. Recently, a commercial software "P-TAK" based on part of our program has been developed for the prediction of particle flows in hoppers and fluidized beds by Yabushita eta!., (1993). gravity acceleration vector moment of inertia of particle stiffness particle mass number of particles unit vector in the normal direction elastic force pressure particle mass flux Reynolds number position vector of particle's gravity center particle radius time unit vector in the tangential direction torque bulk gas velocity gas velocity particle velocity vector displacement void fraction (porosity) damping coefficient friction coefficient air density Poisson ratio angular velocity vector of particle normal component tangential component KONA No.ll (1993)