A Dynamic Coarse Grain Discrete Element Method for Gas-Solid Fluidized Beds by Considering Particle-Group Crushing and Polymerization

The discrete element method (DEM) coupled with computational fluid dynamics (CFD) is used extensively for the numerical simulation of gas-solid fluidized beds. In order to improve the efficiency of this approach, a coarse grain model of the DEM was proposed in the literature. In this model, a group of original particles are treated as a large-sized particle based on the initial particle distribution, and during the whole simulation process the number and components of these particle-groups remain unchanged. However, collisions between particles can lead to frequent crushing and polymerization of particle-groups. This fact has typically been ignored, so the purpose of this paper is to rationalize the coarse grain DEM-CFD model by considering the dynamic particle-group crushing and polymerization. In particular, the effective size of each particle-group is measured by a quantity called equivalent particle-group diameter, whose definition references the equivalent cluster diameter used by the energy-minimization multi-scale (EMMS) model. Then a particle-group crushing criterion is presented based on the mismatch between the equivalent diameter and actual diameter of a particle-group. As to the polymerization of two colliding particle-groups, their velocity difference after collision is chosen as a criterion. Moreover, considering the flow heterogeneity induced by the particle cluster formation, the EMMS drag force model is adopted in this work. Simulations are carried out by using a finite volume method (FVM) with non-staggered grids. For decoupling the Navier-Stokes equations, the semi-implicit method for pressure linked equations revised (SIMPLER) algorithm is used. The simulation results show that the proposed dynamic coarse grain DEM-CFD method has better performance than the original one.


Introduction
Inside a gas-solid fluidized bed, particles driven by a gas flow exhibit complex and intriguing flow patterns. The mutual interaction between the discrete particles and continuum gas brings desirable features for many industrial applications, which include rapid heat and mass transfer, good mixing of solids and fast chemical reaction [1][2][3]. Over the past few decades, many attempts have been made to model the complex gas-solid flow behavior in fluidized beds [4][5][6][7].
The discrete element method (DEM) which was first described by Cundall and Strack [8] is one of the most used simulation tools for dense granular flows [9,10]. In the DEM, each particle is tracked individually based on Newton's second law of motion. The particle-particle and particle-wall interactions are modeled by using springs, dashpots and friction sliders. To be able to model a gas-solid fluidized bed, the DEM is coupled with computational fluid dynamics (CFD) to set up a DEM-CFD important physical process. Consequently, its ability for describing the heterogeneity caused by the complex mesoscale structures of gas-solid systems may be not enough. The system is likely to be overly averaged due to the unchangeable statistical averaging done at the beginning of the simulation.
The purpose of this paper is to rationalize the coarse grain DEM-CFD model by considering the dynamic particle-group crushing and polymerization to overcome the unrealistic coarse grain treatment of particles. In particular, the effective size of each particle-group is measured by a quantity called the equivalent particle-group diameter, whose definition references the equivalent cluster diameter used by the energy-minimization multi-scale (EMMS) model [37,38]. Then a particle-group crushing criterion is presented based on the mismatch between the equivalent diameter and actual diameter of a particle-group. As to the polymerization of two colliding particle-groups, their velocity difference after collision is chosen as a criterion. Moreover, particle clusters can be induced by the complex gas-solid flow. Once a particle cluster is formed, the local fluidized behavior will be greatly changed. The flowing gas is more likely to bypass the particle cluster rather than go through it. Therefore, the drag force acting on the particle cluster will be less than the summation of those acting on the individual particles. That is the drag force drop observed in experiments [39]. Considering the flow heterogeneity induced by the particle cluster formation, the EMMS drag force model is used in this work to predict the drag force drop in case of particle clustering. Simulation results show that the proposed dynamic coarse grain DEM-CFD model is more accurate than the previous one.

Coarse Grain DEM Model for Solid Phase
In the DEM, the trajectory of particle i with mass m i and volume V i is tracked based on solving Newton's equation of transitional motion: Here v i is the translational velocity of particle i. p indicates the gas phase pressure. F drag,i , F g,i and F cont,i are the drag force, gravitational force and contact force acting on particle i, respectively.

Drag Force
For calculating F drag,i , a proper drag force model should be selected. In the literature, the Gidaspow type drag models [40][41][42] and the EMMS drag models [43][44][45][46][47][48] are most used. Generally, the Gidaspow type drag models neglect any heterogeneities induced by the complex gas-solid flow. They lie on the assumption that in each computational cell homogeneous conditions exist. Research work by Wang and Li [44] has demonstrated that the momentum transfer between solid phase and gas phase is over estimated by the Gidaspow type drag models. They are incapable of describing the drag force drop caused by mesoscale effect. By contrast, the EMMS drag models take the flow heterogeneity induced by the particle cluster formation into account though dividing the gas-solid flow into a dilute phase (gas dominated zone) and a dense phase (particle dominated zone). Because the heterogeneity is considered, the EMMS drag force models can predict the drag force drop in case of particle clustering which has been observed in experiments [39]. For the above reason, the following EMMS drag force model [44] is adopted in this study: Appl. Sci. 2020, 10, 1943 4 of 24 where ε i , β i and u i are the void fraction (volume fraction of gas), interphase momentum exchange coefficient and gas velocity associated with particle i, respectively. The β i is given by: where d p , ρ g , C d and H d are the particle diameter, gas density, drag coefficient and heterogeneity index, respectively. The drag coefficient C d for an isolated particle depends on the particle Reynolds number Re p and is given by:

Re p
Re p ≤ 1000 0. 44 Re p > 1000 (4) where Re p is expressed as: Here, µ g is the viscosity of the fluidizing gas phase. The heterogeneity index H d depends on both locally transient and globally averaged information of two-phase velocities and void fraction. It may be expressed as: where a, b and c are parameters associated with the void fraction ε as listed in Table 1 [49,50]. Table 1. Parameters a, b and c in the EMMS drag force model.
These parameters are obtained by researchers through comprehensive optimization. They have been shown to have the ability to produce good results over a wide range of applications, although they may not be optimal for many special cases. Because the parameter optimization is not an easy task, the widely applicable parameters given in Table 1 are preferred in this work. Generally, the value of H d is between 0 and 1. When H d trends to 1, it represents that the particles are locally homogenously distributed and the corresponding fluidizing state is stable. Whereas if H d trends to 0, severe particle clustering exists and the system shows strong heterogeneity.

Contact Force
In case of a collision between two particles (treated as soft spheres), the contact force can be evaluated using a linear spring dashpot model, supplemented with a frictional slider [8]. Consider two colliding partners i and j, the contact force F cont,ij is divided into a normal component F cont n ,ij and a tangential component F cont t ,ij : Here κ n is the normal spring constant (contact stiffness), δ n,ij is the normal overlap between particles i and j, n ij is the unit normal vector pointing from particle i to collision partner j, η n is the normal damping coefficient, v n,ij = (v ij · n ij )n ij is the relative normal velocity of particle i and collision partner j, κ t is the tangential spring constant, δ t,ij is the tangential overlap between particle i and collision partner j, t ij is the unit tangential vector evaluated by v t,ij / v t,ij , η t is the tangential damping coefficient, v t,ij = v ij − v n,ij is the relative tangential velocity of particle i and collision partner j, and µ p is the friction coefficient. The relative velocity at the contact point v ij is defined as follows: where r p and ω are the radius and angular velocity of the colliding bodies, respectively. The spring constant κ depends on the geometry and the elastic properties of the colliding bodies. We have: for the collisions between particles, and: for the collisions between particles and solid walls. Here E p and σ p are respectively the elasticity modulus and the Poisson ratio of the collisions between particles, and σ w and E w are those of the collisions between particles and solid walls. The equivalent shear modulus of the colliding bodies G p in Equations (10) and (11) is given by: The damping coefficient η can be obtained by: where m is the effective mass of the contact bodies evaluated by: for the collisions between particles, and: for the collisions between particles and solid walls. Now the contact force F cont,i acting on particle i can be expressed as the summation about all the particles and solid walls who collide with particle i, i.e.,:

Coarse Graining
The coarse grain DEM model treats a group of original particles with similar properties as a large-sized particle, generally termed as a coarse grained particle, a parcel or a particle-group.
As illustrated in Figure 1 [30], it assumes that all the original particles are spheres with the small diameter, and the size of each particle-group is l times larger than the original particle. Here l is called as coarse grain ratio. The sphericities of both original particle and coarse grained particle are considered as 1. According to the Newton's second law of motion, it is reasonable to assume that the translational motion of the particle-group is the average of that of the original particles governed by: The damping coefficient η can be obtained by: where m is the effective mass of the contact bodies evaluated by: for the collisions between particles, and: for the collisions between particles and solid walls. Now the contact force cont,i F acting on particle i can be expressed as the summation about all the particles and solid walls who collide with particle i , i.e.:

Coarse Graining
The coarse grain DEM model treats a group of original particles with similar properties as a large-sized particle, generally termed as a coarse grained particle, a parcel or a particle-group. (b) Coarse grain particle (particle-group).
As illustrated in Figure 1 [30], it assumes that all the original particles are spheres with the small diameter, and the size of each particle-group is l times larger than the original particle. Here l is called as coarse grain ratio. The sphericities of both original particle and coarse grained particle are The subscript CGM refers to the coarse grain model. Under the assumption that the total kinetic energy of the particle-group is consistent with that of the original particles, we have: where I is the moment of inertia. The bar over the variables and the subscript O refer to average and original particle, respectively. According Equations (17) and (18), the following formulas hold: where g is the gravitational acceleration.

CFD Model for Gas Phase
In a gas-solid fluidized beds, it is difficult to directly calculate the mutual interaction between the continuum gas and discrete particles. An alternative way to solve this problem is the local volume average technique which has been widely used in the continuum approach of modelling fluidization process [42]. According to this volume average technique, the governing equations of the gas phase are the fluid continuity and Navier-Stokes equations expressed in terms of the local mean variables over the computational cell: where ε, ρ g , u, p and g have the same meaning as mentioned above. S p and τ g are the volumetric gas-particle interaction force and viscous stress tensor, respectively. According to Newton's third law of motion, the value of S p at a given CFD computational cell is defined as the volume average of the drag forces acting on all the particles located in this cell: where V cell and N c are the volume of the computational cell and the number of particles located in this cell, respectively. The viscous stress tensor τ g is given by: By definition, the void fraction ε in a computational cell is the ratio of the void volume (gas volume) to the volume of the cell. If V i is the volume of particle i inside a computational cell, then the void fraction in this cell is: The void fraction ε i and gas velocity u i at the center of mass of particle i can be calculated by interpolation from the flow variables stored at the grid nodes. As illustrated in Figure 2, particle i divides the cell, which it belongs to, into four parts. Each part has the area S ij ( j = 1, 2, 3, 4). If the Appl. Sci. 2020, 10, 1943 8 of 24 void fraction and gas velocity valued at the four vertexes of this cell are ε ij and u ij ( j = 1, 2, 3, 4), then ε i and u i can be given by: in the two-dimensional cases.
The void fraction i ε and gas velocity i u at the center of mass of particle i can be calculated by interpolation from the flow variables stored at the grid nodes. As illustrated in Figure 2, particle i divides the cell, which it belongs to, into four parts. Each part has the area ( 1, 2,3, 4) ij S j = . If the void fraction and gas velocity valued at the four vertexes of this cell are ij ε and then i ε and i u can be given by: (31) in the two-dimensional cases. Figure 2. Illustration of the interpolation from gas to particles (two dimensional cases).

Particle-Group Crushing and Polymerization Model
Because of the complex particle-particle and gas-particle interactions inside a gas-solid fluidized beds, the particle-groups of the coarse grain DEM-CFD model can hardly to maintain their initial configurations during the whole simulation process. Frequent crushing and polymerization of the particle-groups should be taken into consideration. However, the current coarse grain DEM-CFD model ignores this important physical process. Consequently, its ability for describing the heterogeneity caused by the complex mesoscale structures of gas-solid systems may be not enough. The system is likely to be overly averaged due to the unchangeable statistical averaging done at the beginning of the simulation. In order to rationalize the coarse grain DEM-CFD model by overcoming

Particle-Group Crushing and Polymerization Model
Because of the complex particle-particle and gas-particle interactions inside a gas-solid fluidized beds, the particle-groups of the coarse grain DEM-CFD model can hardly to maintain their initial configurations during the whole simulation process. Frequent crushing and polymerization of the particle-groups should be taken into consideration. However, the current coarse grain DEM-CFD model ignores this important physical process. Consequently, its ability for describing the heterogeneity caused by the complex mesoscale structures of gas-solid systems may be not enough. The system is likely to be overly averaged due to the unchangeable statistical averaging done at the beginning of the simulation. In order to rationalize the coarse grain DEM-CFD model by overcoming the unrealistic coarse grain treatment of particles, a dynamic particle-group crushing and polymerization model is present in this section.
For capturing the crushing of a particle-group, the effective size of each particle-group should be measured. Following the definition of the equivalent cluster diameter used by the EMMS model [37,38], we use d cl to represent the equivalent diameter of a particle-group. In a gas-solid fluidized bed, d cl is generally inversely proportional to the energy input from the gas E input : where K is a coefficient to be determined and E input equals to the difference between the total energy input E st and the energy input for the state of generalized minimum fluidization (E st ) m f : For the EMMS model, the energy consumption is divided into three parts [38]. That is: Appl. Sci. 2020, 10, 1943 9 of 24 where the subscripts c, f and b refer to the dense phase, dilute phase and inter phase, respectively. ρ p , n, F, f and U are the particle density, estimated particle number of each phase, drag force acting on single particle or cluster in suspension, volume fraction of dense phase and apparent gas velocity, respectively. Since the drag force and apparent gas velocity have the same direction, only their norms are considered in Equation (34). At the state of generalized minimum fluidization, the energy input (E st ) m f can be calculated by: where U m f , U p , ε m f and g are the minimum fluidization velocity, apparent particle velocity, void fraction at the minimum fluidization and gravitational constant, respectively. According to the relationship d cl = d p at the limiting condition ε = ε max , the coefficient K is determined by: For the DEM model, because the motion of each particle can be directly tracked, there is no need to divide the system into dense, dilute and inter phases. The number of particles in any local area can be exactly calculated. Therefore, the expression of E st for a single particle-group i may be simplified as: where n i is the number of original particles inside the particle-group i. Note that Equation (32) only takes the influence of the external gas on the equivalent diameter of a particle-group into consideration. In fact, the mutual interactions between particle-groups also have considerable influence on the equivalent diameter of a particle-group. For this reason, the energy input from other particle-groups E cc is considered in our coarse grain DEM-CFD model. In this situation, we have: where F cont, CGM ij , v CGM ij m CGM i and N p are the contact force exerted by particle-group j to particle-group i, relative velocity of particle-groups i and j, mass of particle-group i and number of particle-groups colliding with particle-group i. Finally, the equivalent particle-group diameter d cl for the coarse grain DEM-CFD model can be given by: From Equations (35)-(39), if ε max is given, then d cl can be calculated. The analysis of Matsen [51] indicated that all the particles would be completely dispersed in gas-solid fluidized beds at a critical void fraction 0.9997. So it is reasonable to set ε max as this critical value 0.9997. Now, a particle-group crushing criterion can be presented based on the mismatch between the equivalent diameter (d cl ) i and actual diameter (d CGM ) i of a particle-group i. Specifically, if: holds, then particle-group i will break into several child particle-groups. The number of the child particle-groups n new is determined by: where · refers to the operation of rounded down. For simplicity, it assumes that there are n new − 1 child particle-groups having the same diameter (d cl ) i and mass m k , and the n new th child particle-group consisted by the remainder mass has the following mass and diameter: The velocity of the child particle-groups is inherited from their parent particle-group i. Then, as to the polymerization of two colliding particle-groups, their velocity difference after collision is chosen as a criterion. Considering the collision between particle-groups i and j, if their velocities after collision have the following relationship: then the two particle-groups are regarded as aggregated. After polymerization, the mass and diameter of the new particle-group are given by: The velocity of the new particle-group is set as either v CGM i or v CGM j .

Numerical Method
For two dimensional cases, the CFD model for gas phase has the following component form: where the subscripts x and y refers to the components in x and y directions, respectively. The finite volume method (FVM) with non-staggered grids is used to discretize Equation (47). Figure 3 illustrates the control volume for a grid node P. By using this control volume, the discrete forms of Equations (47b) and (47c) at the grid node P can be uniformly written as: where the components of the source term B are given by a a a a t ε Δ Δ = + + + + Δ for the central scheme. Note that the momentum Equation (48) can be solved only when the pressure field is given or is somehow estimated. The semi-implicit method for pressure linked equations revised (SIMPLER) is used for solving the pressure and thus the velocity. In order to derive an equation for pressure, the continuity Equation (47a) at the grid node P is discretized as: by using the FVM under the assumption that the void fraction ε maintains constant within a time step t Δ . To avoid the problem of volume locking on non-staggered grids, the pressure is removed from Equation (48) to obtain a pseudo velocity u  : It can be noted that P u  is composed of the neighbor velocities and contains no pressure. In turn, the gas velocity at the interface of the control volume can be written as: Here u 0 x and u 0 y are the gas velocities at the previous time step. The coefficients are: for the upwind scheme and a P = a E + a W + a N + a S + ε P ∆x∆y ∆t for the central scheme.
Note that the momentum Equation (48) can be solved only when the pressure field is given or is somehow estimated. The semi-implicit method for pressure linked equations revised (SIMPLER) is used for solving the pressure and thus the velocity. In order to derive an equation for pressure, the continuity Equation (47a) at the grid node P is discretized as: by using the FVM under the assumption that the void fraction ε maintains constant within a time step ∆t. To avoid the problem of volume locking on non-staggered grids, the pressure is removed from Equation (48) to obtain a pseudo velocity u: It can be noted that u P is composed of the neighbor velocities and contains no pressure. In turn, the gas velocity at the interface of the control volume can be written as: By inserting Equation (53) into Equation (50), we get: ε e a e + ε w a w ∆y + ε n a n + ε s a s ∆x p P = ε e a e p E + ε w a w p W ∆y + ε n a n p N + ε s a s p S ∆x Now if a guessed velocity u is given, the pseudo velocity u can be calculated from Equation (51) and hence a guessed pressure p * can be solved from Equation (54). By employing p * , a velocity can be solved from Equation (48). Because the pressure p * is not correct in general, the resulting velocity will not satisfy the continuity equation. Such an imperfect velocity field based on the guessed pressure p * is denoted by u * . Suppose the correct pressure is p = p * + p , where p is called the pressure correction. Next, the corresponding velocity correction u can be introduced in a similar manner u = u * + u . At this point, the velocity at the interface of the control volume, taking (u x ) e for example, can be expressed as: By inserting the interfacial velocities in form of Equation (55) into the discrete continuity Equation (50), we obtain an equation for p : ε e a e + ε w a w ∆y + ε n a n + ε s a s ∆x p P = ε e a e p E + ε w a w p W ∆y + ε n a n p N + ε s a s p S ∆x Now the numerical method for solving the CFD model for gas phase can be summarized as follows: (1) Calculate the pseudo velocity u from Equation (51)  Finally, the flow chart for the whole coarse grain DEM-CFD method with considering the particle-group crushing and polymerization is shown in Figure 4. Finally, the flow chart for the whole coarse grain DEM-CFD method with considering the particle-group crushing and polymerization is shown in Figure 4.

Simulation Results
Two simulation cases are performed to test the dynamic coarse grain DEM-CFD method. In the first simulation case, a two-dimensional fluidized bed with 90 mm width and 500 mm height is considered. The geometry is discretized as 18 × 100 cells by setting the grid size to be 5 mm. Gas is injected from the base, where the superficial velocity is set to be 0.9 m/s. The number of solid particles

Simulation Results
Two simulation cases are performed to test the dynamic coarse grain DEM-CFD method. In the first simulation case, a two-dimensional fluidized bed with 90 mm width and 500 mm height is considered. The geometry is discretized as 18 × 100 cells by setting the grid size to be 5 mm. Gas is injected from the base, where the superficial velocity is set to be 0.9 m/s. The number of solid particles is 4080. All the particles have equal size and are assumed to be spheres with sphericity equal to 1. The physical properties of the particles and gas are illustrated in Table 2. The time step size used in the simulations is 0.05 ms. Because the simulation conditions are very similar with those adopted in the experiments made by Van Wachem et al. [18], we can compare our simulation results with the experimental ones. Before we start the simulations, some comparative studies are made to demonstrate that the EMMS drag force model fits our need. Both the EMMS drag force model and the famous Wen and Yu drag force model are used to the present simulation case. The original DEM-CFD method is used to generate simulation results. Both qualitative and quantitative comparisons show that the EMMS results are more consistent with the experimental ones. These facts as well as the research work done by many other researchers give us the reason to believe that the EMMS drag force model fits our need.
Although the particle-group crushing and polymerization model is proposed for particle-groups, we are still curious whether the simulation results will be improved if it is applied to original particles. For this purpose, we first directly use our dynamic coarse grain DEM-CFD method to the original particles under the assumption that the particles are divisible. In this case, each particle is treated as a particle-group. Figure 5 shows the results obtained by the current method together with those obtained by experiments and traditional coarse grain DEM-CFD method. For facilitating the comparisons, we call our dynamic coarse grain DEM-CFD method involving the particle-group crushing and polymerization model as PCPM and the traditional coarse grain DEM-CFD method as TM for short. The results show that the particles are in the fluidized state. Bubbles are generated periodically at the bottom of the bed and gradually rise and swell. The bubble nose is getting thinner and thinner. When the maximum bed height is reached, the top particles begin to fall back and the bubbles eventually disappear. No particles can flee from the computational domain. They will bounce when they collide with the wall. The particle-wall interactions are modeled by using springs, dashpots and friction sliders as introduced in Section 2. This consideration prevents the particles from going through the solid wall. On the other hand, if the gas velocity is too high, the particles may escape from the outlet of the fluidized bed. But for the simulation cases we considered, the gas velocity is moderate and the gravitational force prevents the particles escape from the outlet. By comparison, the results obtained by the PCPM are more consistent with the experimental results than those obtained by TM. That is the PCPM has the ability to improve the simulation results of the original DEM-CFD method (TM with a coarse grain ratio equal to 1). The transverse and longitudinal sizes and positions of the large bubbles predicted by the PCPM are much closer to the experimental results. Therefore, it appears that the PCPM can simulate the macroscopic flow properties of the fluidized bed. The particles at the top layer have a longer suspension time, which indicates that the particles at the top of the bed are crushing and separated, thus reducing the drag force required for suspension of these particles. These results indicate that the PCPM can better capture the phenomena of bubbles popping out of the bed and particle-group crushing than the TM.  Figure 7 shows the change of bed height over observation time. The bed height predicted by the PCPM is slightly higher than that of the TM, as shown in this figure. This phenomenon is mainly caused by the consideration of particle-group crushing. That is, when the bubbles move to the top of the bed during the bed moving, the particles in the top layer will be blown higher due to the reduction The above is just an apparently qualitative comparison. We now compare the simulation results of the two methods and the experimental results quantitatively. Figure 6 shows the pressure drop at the 45 mm horizontal section above the fluidized bed distributor. It can be seen that both the PCPM and TM predicted the periodic fluctuation of the pressure. However, in terms of amplitude, the result of PCPM is closer to the experimental result. The result of TM shows rather flat pressure valleys which indicates a long duration of low pressure, while the PCPM improves this low pressure effect.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 25 variations of the total particle number and particle-group crushing amount during the simulation of PCPM are shown in Figures 8 and 9. These results indicate that the collisions between particles cause violent crushing and polymerization of particle-groups during the simulation, and thus the total particle number presents a dynamic fluctuation.   Figure 7 shows the change of bed height over observation time. The bed height predicted by the PCPM is slightly higher than that of the TM, as shown in this figure. This phenomenon is mainly caused by the consideration of particle-group crushing. That is, when the bubbles move to the top of the bed during the bed moving, the particles in the top layer will be blown higher due to the reduction of their sizes caused by particle-group crushing. Although there are some differences between the two simulation results and experimental result, both the simulations are acceptable. The dynamic variations of the total particle number and particle-group crushing amount during the simulation of PCPM are shown in Figures 8 and 9. These results indicate that the collisions between particles cause violent crushing and polymerization of particle-groups during the simulation, and thus the total particle number presents a dynamic fluctuation.
Next, we redo the above simulations at a coarse grain ratio l = 2. The particle-group distribution predicted by the PCPM and TM are shown in Figure 10. The results of both PCPM and TM show a qualitative agreement with the experimental results about the distribution of the particles. Hence, the macroscopic flow properties of the fluidized bed obtained from the coarse grain models can be regarded to be equivalent to the original particle system to some extent. Still, there are several differences between the results of the two methods. First, the particle distribution predicted by the TM is relatively homogeneous, while that predicted by the PCPM shows strong aggregations. Second, the top layer of the bed predicted by the TM changes mildly when the particles rise. By contrast, the PCPM can capture the drastic movement of the top layer particles when the bubbles come out of the bed. Finally, the PCPM is capable of showing the presence of multiple bubbles in the bed compared to the TM. These results indicate that the PCPM is better that the TM in describing the heterogeneity of the particle distribution in the fluidized bed. Moreover, as shown in Figure 11, the crushing and polymerization of particle-groups in PCPM make the pressure fluctuation more violent. This property alleviates the problem of too flat pressure drop predicted by the TM.   Next, we redo the above simulations at a coarse grain ratio 2 l = . The particle-group distribution predicted by the PCPM and TM are shown in Figure 10. The results of both PCPM and TM show a qualitative agreement with the experimental results about the distribution of the particles. Hence, the macroscopic flow properties of the fluidized bed obtained from the coarse grain models can be regarded to be equivalent to the original particle system to some extent. Still, there are several differences between the results of the two methods. First, the particle distribution predicted by the TM is relatively homogeneous, while that predicted by the PCPM shows strong aggregations. Second, the top layer of the bed predicted by the TM changes mildly when the particles rise. By contrast, the  Next, we redo the above simulations at a coarse grain ratio 2 l = . The particle-group distribution predicted by the PCPM and TM are shown in Figure 10. The results of both PCPM and TM show a qualitative agreement with the experimental results about the distribution of the particles. Hence, the macroscopic flow properties of the fluidized bed obtained from the coarse grain models can be regarded to be equivalent to the original particle system to some extent. Still, there are several differences between the results of the two methods. First, the particle distribution predicted by the TM is relatively homogeneous, while that predicted by the PCPM shows strong aggregations. Second, the top layer of the bed predicted by the TM changes mildly when the particles rise. By contrast, the PCPM can capture the drastic movement of the top layer particles when the bubbles come out of the bed. Finally, the PCPM is capable of showing the presence of multiple bubbles in the bed compared to the TM. These results indicate that the PCPM is better that the TM in describing the heterogeneity of the particle distribution in the fluidized bed. Moreover, as shown in Figure 11, the crushing and polymerization of particle-groups in PCPM make the pressure fluctuation more violent. This property alleviates the problem of too flat pressure drop predicted by the TM.  For further validating the performance of PCPM, we increase the initial particle number to 9600 and adjust the superficial gas velocity to 1.9 m/s accordingly. Figure 12 shows the evolution of particle distribution at the very early stage. As can be seen, the system exhibits more complex fluidization phenomenon when the number of particles is increased. Influenced by the wake vortexes of the large bubble in the middle region of the bed, airflow escapes from the small bubbles below to the large bubble. At the same time, the particles aggregating at the top of the small bubbles are collapsing. This process causes the disappearance of small bubbles and the growth of large bubble. Because of the thick layer of aggregated particles at the top of the bed, the large bubble cannot pass through the bed  For further validating the performance of PCPM, we increase the initial particle number to 9600 and adjust the superficial gas velocity to 1.9 m/s accordingly. Figure 12 shows the evolution of particle distribution at the very early stage. As can be seen, the system exhibits more complex fluidization phenomenon when the number of particles is increased. Influenced by the wake vortexes of the large bubble in the middle region of the bed, airflow escapes from the small bubbles below to the large bubble. At the same time, the particles aggregating at the top of the small bubbles are collapsing. This process causes the disappearance of small bubbles and the growth of large bubble. Because of the thick layer of aggregated particles at the top of the bed, the large bubble cannot pass through the bed For further validating the performance of PCPM, we increase the initial particle number to 9600 and adjust the superficial gas velocity to 1.9 m/s accordingly. Figure 12 shows the evolution of particle distribution at the very early stage. As can be seen, the system exhibits more complex fluidization phenomenon when the number of particles is increased. Influenced by the wake vortexes of the large bubble in the middle region of the bed, airflow escapes from the small bubbles below to the large bubble. At the same time, the particles aggregating at the top of the small bubbles are collapsing. This process causes the disappearance of small bubbles and the growth of large bubble. Because of the thick layer of aggregated particles at the top of the bed, the large bubble cannot pass through the bed quickly. Squeezed by the top layer particles, the bubble almost transforms into rectangular in shape and becomes large enough to completely fill the cross-section of the bed. Until 1.68 s, cracks appear in the top layer of the bed and the gas trapped in the large bubble escapes quickly.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 20 of 25  Figure 13 shows the evolution of particle distribution at the later stage. The phenomenon of bubble merging along the transverse direction can be observed from these figures. In this process, the bubbles on the left and right sides near the walls grow rapidly due to the vertical bubble merging. Soon two slender bubbles are formed on the left and right sides. The particle groups between the two bubbles gradually thin and disappear under the influence of air flow on both sides, and finally the two bubbles merge horizontally. In the case shown in the figures, since the two bubbles have reached the top layer of the bed by the time they finally merge, so the complete merging morphology is not seen. These results reflect the complexity of the fluidization pattern change in dense particle flow systems. It suggests that the heterogeneous mesoscale structures have great influence on gas-solid flows. In addition to the traditional spherical bubbles, there are also bubbles with complex structures, which will play a leading role in the dynamic behavior of the whole bed.   Figure 13 shows the evolution of particle distribution at the later stage. The phenomenon of bubble merging along the transverse direction can be observed from these figures. In this process, the bubbles on the left and right sides near the walls grow rapidly due to the vertical bubble merging. Soon two slender bubbles are formed on the left and right sides. The particle groups between the two bubbles gradually thin and disappear under the influence of air flow on both sides, and finally the two bubbles merge horizontally. In the case shown in the figures, since the two bubbles have reached the top layer of the bed by the time they finally merge, so the complete merging morphology is not seen. These results reflect the complexity of the fluidization pattern change in dense particle flow systems. It suggests that the heterogeneous mesoscale structures have great influence on gas-solid flows. In addition to the traditional spherical bubbles, there are also bubbles with complex structures, which will play a leading role in the dynamic behavior of the whole bed.
At last, Figure 14 shows some results obtained at coarse grain ratios l = 2 and l = 3. As can be seen, in the initial stage of simulation (t = 0.4s), the bed experienced uniform expansion under the effect of gas velocity. The bed morphology simulated by the coarse grain system with l = 2 is relatively close to that of the original system. However, when the coarse grain ratio increases to l = 3, the bed shows a morphology between the scattered bed and the bubbling bed. Although the bubbles generated are very small, the whole bed is still in a fluidized state. In the middle stage of simulation (t = 6.44s), the system with coarse grain ratio l = 2 can still keep a similar state to the original system, while the result of l = 3 is seriously distorted. Accordingly, the coarse grain ratio of the simulation system needs to be carefully selected. Excessive coarse grain will make the coarse grain system severely deviate from the original system, and proper coarse grain will greatly reduce the simulation time. For the current simulation, when the coarse grain ratio is l = 2, the simulation time of the coarse grain system is about 10% of the original system. two bubbles merge horizontally. In the case shown in the figures, since the two bubbles have reached the top layer of the bed by the time they finally merge, so the complete merging morphology is not seen. These results reflect the complexity of the fluidization pattern change in dense particle flow systems. It suggests that the heterogeneous mesoscale structures have great influence on gas-solid flows. In addition to the traditional spherical bubbles, there are also bubbles with complex structures, which will play a leading role in the dynamic behavior of the whole bed. At last, Figure 14 shows some results obtained at coarse grain ratios 2 l = and 3 l = . As can be seen, in the initial stage of simulation ( 0.4s t = ), the bed experienced uniform expansion under the effect of gas velocity. The bed morphology simulated by the coarse grain system with 2 l = is relatively close to that of the original system. However, when the coarse grain ratio increases to 3 l = , the bed shows a morphology between the scattered bed and the bubbling bed. Although the bubbles generated are very small, the whole bed is still in a fluidized state. In the middle stage of simulation ( 6.44s t = ), the system with coarse grain ratio 2 l = can still keep a similar state to the original system, while the result of 3 l = is seriously distorted. Accordingly, the coarse grain ratio of the simulation system needs to be carefully selected. Excessive coarse grain will make the coarse grain system severely deviate from the original system, and proper coarse grain will greatly reduce the simulation time. For the current simulation, when the coarse grain ratio is 2 l = , the simulation time of the coarse grain system is about 10% of the original system.

Conclusions
In this work, a dynamic coarse grain DEM-CFD method for gas-solid fluidized beds is presented by considering the particle-group crushing and polymerization. Compared with the original coarse grain DEM-CFD method, the coarse graining system is updated at every simulation step. It allows a dynamic adjustment of the number and components of the particle-groups. Big particle-groups can be crushed into several small particle-groups or original particles. Conversely, several small particlegroups or original particles can be polymerized into a big particle-group. That is, the number and size of the particle-groups are not invariable but change dynamically during the simulation. The size of each particle-group relies on its equivalent diameter, which is measured by an EMMS type model. If the actual diameter of a particle-group is greater than its equivalent diameter, then crushing of this particle-group takes place. As to the polymerization of two colliding particle-groups, their velocity difference after collision is chosen as a criterion. If the velocities of the two particle-groups after collision are fairly close, then the two particle-groups can be regarded as one new particle-group. Due to the diversity of particle-group sizes, the resulting particle-group heterogeneity as well as the flow heterogeneity induced by the particle cluster formation should be considered. For this reason, the EMMS drag force model is adopted in this work. A non-staggered FVM and the SIMPLER algorithm are used to perform the simulations. Flow patterns and pressure drop predicted by the current dynamic coarse grain DEM-CFD method agree better with the experimental results than the original method. For fluidized beds involving more particles, the new method can reveal the complex fluidization pattern evolution in dense particle flow systems. The simulation results also suggest that

Conclusions
In this work, a dynamic coarse grain DEM-CFD method for gas-solid fluidized beds is presented by considering the particle-group crushing and polymerization. Compared with the original coarse grain DEM-CFD method, the coarse graining system is updated at every simulation step. It allows a dynamic adjustment of the number and components of the particle-groups. Big particle-groups can be crushed into several small particle-groups or original particles. Conversely, several small particle-groups or original particles can be polymerized into a big particle-group. That is, the number and size of the particle-groups are not invariable but change dynamically during the simulation. The size of each particle-group relies on its equivalent diameter, which is measured by an EMMS type model. If the actual diameter of a particle-group is greater than its equivalent diameter, then crushing of this particle-group takes place. As to the polymerization of two colliding particle-groups, their velocity difference after collision is chosen as a criterion. If the velocities of the two particle-groups after collision are fairly close, then the two particle-groups can be regarded as one new particle-group. Due to the diversity of particle-group sizes, the resulting particle-group heterogeneity as well as the flow heterogeneity induced by the particle cluster formation should be considered. For this reason, the EMMS drag force model is adopted in this work. A non-staggered FVM and the SIMPLER algorithm are used to perform the simulations. Flow patterns and pressure drop predicted by the current dynamic coarse grain DEM-CFD method agree better with the experimental results than the original method. For fluidized beds involving more particles, the new method can reveal the complex fluidization pattern evolution in dense particle flow systems. The simulation results also suggest that excessive coarse grain will make the coarse grain system severely deviate from the original system and proper coarse grain will greatly reduce the simulation time. In conclusion, the highlights of this work can be summarized as follows: • A dynamic particle-group crushing and polymerization procedure is considered in the coarse grain DEM-CFD model.

•
The dynamically coarse graining is based on an EMMS type model.

•
Heterogeneities are considered in the coarse graining systems by using the EMMS drag force model.

•
The method is validated by comparing the simulation results with the experimental ones.