CFD-DEM Simulation of Spouted Bed Dynamics under High Temperature with an Adhesive Model

: Particle adhesion is of great importance to coating processes due to its effect on fluidization. Currently, Computational Fluid Dynamics-Discrete Element Method (CFD-DEM) has become a powerful tool for the study of multiphase flows. Various contact force models have also been proposed. However, particle dynamics in high temperature will be changed with particle surface properties changing. In view of this, an adhesion model is developed based on approaching-loading-unloading-detaching idea and particle surface change under high temperature in this paper. Analyses of the adhesion model are given through two particle collision process and validated by experiment. Effects of inlet gas velocity and adhesion intensity on spouted bed dynamics are investigated. It is concluded that fluidization cycle will be accelerated by adhesion, and intensity of fluidization will be marginally enhanced by slight adhesion. Within a certain range, increasing inlet gas velocity will lead to strong intensity of particle motion. A parameter sensitivity comparison of linear spring-damping model and Hertz-Mindlin Model is given, which shows in case of small overlaps, forces calculated by both models have little distinction, diametrically opposed to that of large overlaps.


Introduction
Spouted beds (SBs) have drawn much attention in recent years for a broad range of applications, such as coating [1], drying [2], biomass [3] or waste pyrolysis [4,5] and gasification [6]. In fact, these processes are conducted at high temperatures and particle adhesion can give way to severe operational problems, even causing bed de-fluidization [7]. Accordingly, it is necessary to study the adhesion process to better control the adverse effects of it.
Due to the high fluid-solid contact efficiency, especially for coarse particles in drying and coating processes [8][9][10], the spouted bed has unique characteristics as a prototype reactor for research. In addition to their ability to handle coarse particles, spouted beds also possess certain structural and flow characteristics that are very desirable in some chemical reaction systems because of their perfect particle periodic dynamics. For instance, a central spout formation, where the particles move upwards in a dilute phase, and a semi-packed bed annular space between the spout zone and the vessel wall, where particles move slowly downwards and radially inwards, are typical periodic dynamics phenomena. They are quite beneficial for particle coating.
CFD modelling has become a powerful tool for the study of multiphase flows with the development of computational power and the advance of numerical algorithms. There are two kinds of methods that are called the Eulerian-Eulerian approach (two-fluid model, TFM) and the Eulerian-Lagrangian approach (CFD-DEM) in general. In TFM, the kinetic ∂ ερ ∂t + ∇ • ερ u = 0 (1) ∂ ερ u ∂t + ∇ • ερ uu = −ε∇p − ∇ • ετ − F + ερ g (2) τ g =-λ g -2 3 μ (∇•u)I-μ((∇u)+(∇u) T ) For an individual particle, its motion equation, i.e., Newton's laws, are listed as (5)-(6). A multi-particle fluid drag force Di Felice model [29] is used as equations (7)- (10). The drag scaling factor is also applied considering the nonuniformity effect of fluid drag on the particles inside and outside an agglomerate [30]. This particle contact force shown as equation (11) can be calculated by linear spring model and Hertz-Mindlin model. A detailed comparison is given in the Section 3. In this paper, the flow field information is solved by commercial software Fluent 16.1 while the particle motion model is solved by User Defined Function (UDF) using in-house code. m i v =F C +F f +m i g (5) χ=3.7-0.65exp(-(1.5-log 10 Re) 2 2 ) (10)

CFD-DEM with the Adhesion Model
Due to temperature rise, chemical reaction or coating in fluidized bed, a new surface of a particle may be generated by melting or corroding locally, resulting in the change of surface properties as shown in Figure 1. The liquid layer generated by melting under high temperature will decrease the elasticity and increase the plasticity of particles, thus making them adhesive. For simulating this phenomenon accordingly, an adhesion model can be developed based on the CFD-DEM frame. The adhesion model in this paper is derived from Luding's idea [30] and integrated with force relationships in JKR, as well as a derivation of phase transition due to high temperature. There are four stages in linear hysteresis model: approaching, loading, unloading, and detaching. At approaching stage, van der Waals force, which is the dominant force in nanoparticle simulation, can be negligible because of the large particle spouted here. At loading stage, the normal force is gradually increasing with normal overlap rising until the particle velocity decreases to 0. At the unloading stage, contact force and adhesion force work together, where particle velocity increases inversely and the overlap keeps down to 0. In other words, particles are rebounding. At detaching stage, adhesion force still exists and particles continue to set apart until liquid bridge rupture when the connection is suddenly broken, and the two particles are completely separated. In the JKR model, the breaking force is 5/9 of "pull off" force, while the adhesion force is 8/9 of that. Therefore, in Equation (20), when x 2 = 1, the breaking force should be 5/8 of adhesion force.
Normal force at approaching stage is calculated by equation (12): Normal force at loading stage is calculated by equation (13): F n (δ)=-k load n δ n (13) Normal force at unloading stage is calculated by equations (14)- (19): F n (δ)= -k unload n δ n -δ p , F n ≥f ij,min (1-x 1 )( -F adhesion )+x 1 f ij,min , F n <f ij,min (14) x 1 = δ n δ ij,min (15) k unload δ ij,min = f ij,min k unload n +δ p (18) Normal force at detaching stage is calculated by equations (20)- (21): The adhesion phenomenon of concern in this paper is closely relevant to the detaching stage, so this part is further modeled. It should be noted that adhesion between particles is controlled by a type of capillary force owing to a quite thin melting liquid layer as shown in Figure 1. Hence, adhesion force can be analogous to the capillary force model, in which the key parameter of liquid content is the same as liquid layer in our model. The adhesion model is relevant with temperature because of surface melting. The following is the modeling process of adhesion force with temperature.
Previous studies classified liquid bridge force into capillary term and viscous term [31,32], where lubrication effect (produced by viscous term) is quite small thus being ignored in this paper [33]. For original capillary force model, it is not solved directly in DEM [32]. Therefore, it is usually employed the form of contact angle and the rupture distance [31,33]. As aforementioned, the adhesion force can be shown as: where Γ is the surface tension, θ 0 is the contact angle between liquid and grain surface. D is the gap between the spherical surfaces. Liquid content V is an important parameter in the capillary force model. In the melting process, assuming the thickness of liquid layer is quite thin, so the diameter of particles can be regard as constant as well as particle density. The heat transfer between gas and particles is fast, so the heat equilibrium is obtained: where C 1 and C 2 are the specific heat capacity of particles and gas, respectively. ∆T p and ∆T g are the temperature change of particles and gas. h is latent heat of melting and q m is the mass flow rate of gas. M = NρV, N is the number of total particles. ρ is the particle density. V is the volume of liquid layer of for a single particle. Compared to the change of particle temperature, the gas temperature decreases very little, so it can be considered constant when simulating fluidization phenomenon. Here, we assume ∆T = k ∆T . The particles are in solid-liquid coexistence state, so the internal temperature of particles remains the same as the surface temperature. Finally, the liquid content can be obtained as follows: where Q d is the energy loss due to heat dissipation. It should be noticed that the focus on the paper is the effect of adhesion on fluidization, while the rupture distance dc plays a crucial role in adhesion phenomenon. Therefore, dc should be selected properly. Consequently, a coefficient k1 is applied to the traditional capillary force model in the CFD-DEM simulation: where k1 is changed according to the adhesion state.

Geometry and Parameters
A typical 2D spouted bed is simulated here, as shown in Figure 2. Geometric parameters and all the conditions are shown in Tables 1 and 2. The apparatus is composed of a cuboid and a conical body with a reverse cone angle of 60°. Total height is 500 mm. The bottom width of the cone is 15 mm while that of the cuboid is 150 mm. The jet inlet width is 9 mm. Overall thickness of the apparatus is 15 mm, which is 7.5 times the diameter of the particle, thus eliminating both the three-dimensional effect on bed thickness and wall friction effect. The apparatus is made of transparent plexiglass and the process can be recorded with a high-speed camera. The high-speed camera has a resolution of 1000 fps. To validate the simulation, a typical working condition is tested as follows. Particles are loaded into the spouted bed slowly. Particle loading is 100 mm. Compressed air is used as the source gas. Gas inlet flow rate is controlled and measured by the rotameter, which is adjusted to 26 m/s (inlet gas velocity) and the fluidization state is recorded by the highspeed camera when stabilization spouted state is achieved.

Comparison of Contact Force Models
A detailed comparison of linear spring-damping model and Hertz-Mindlin model is given as follows: Normal force models were tested for particles with a radius of 0.001 m, and the results were shown in Figure 3 and Table 3. When normal overlap increases from 0 to 0.001 m, normal force calculated by linear model changes by about 5 N, while that of Hertz-Mindlin model (short for "H-M") changes by about 13 N, which indicates normal force with H-M has a larger range. In the scale of small overlap (<0.3 r), the difference between two models is not obvious, while in the scale of large overlap (>0.6 r), the force (absolute value) calculated by the linear model is apparently lower than that of H-M. Figure 3 also illustrates that as relative normal velocity of particles increases, the critical point, at which the two models are equal, will shift to the upper right and the section between two critical points will broaden slowly. When the relative normal velocity is 90 m/s, the overlap at critical point is 0.3 r, and when the velocity decreases to 10 m/s, the overlap is about 0.24 r.
β= lne n π 2 +ln 2 e n S n =2E * R * δ n Input parameters Normal spring constant (k n ) Normal coefficient of restitution (e n ) Particle radius R 1 , R 2 Particle Poisson ratio ν 1 , ν 2 Shear modulus G 1 , G 2 Final expression F n =-k n δ n -2 (m * k n ) lne n π 2 +ln 2 e n v ij,n F n =-4 3 E * R * 1 2 ⁄ δ n  According to the Table 4, the two models were compared and analyzed with different tangential velocities as well as normal overlaps, as shown in Figure 4. Figure 4a suggests that tangential forces of both models show a linear relationship with tangential overlap, and the results are extremely consistent. However, with relative tangential velocity rising up, the difference between two models gradually gets obvious, and the result of H-M is slightly lower than that of linear model. It can be noted from Table 4 that the tangential force of the H-M is not only related to relative tangential velocity, but also connected with normal overlap, while it has nothing to do with the normal overlap in linear model. Therefore, Figure 4b explores the influence of various normal overlaps. As can be seen from the Figure 4b, when normal overlap is about 0.1 r, the results of the two models are in the best agreement. When overlap is minimal (<0.05 r), the results (absolute value) of H-M are lower than those of linear model, but the discrepancy is acceptable. However, when overlap is large (>0.6 r), the results (absolute value) of H-M are much higher than those of linear model, and this difference should be carefully considered. Table 4. Comparison of tangential force models.

Item Linear Spring-Damping Model Hertz-Mindlin Model
Tangential force η t = 2 (m * k t )lne t π 2 +ln 2 e t S t =8G * R * δ n Input parameters Tangential spring constant (k t ) Tangential coefficient of restitution (e t ) Particle radius R 1 , R 2 Particle Poisson ratio ν 1 , ν 2 Shear modulus G 1 , G 2 Final expression F t =-k t δ t -2 (m * k t ) lne t π 2 +ln 2 e t v ij,t F t =-8G * R * δ n δ t -4 5 3 G * where L i is the distance from the point of contact to the center of the sphere, F t is tangential force, µ r is the coefficient of rolling friction F n , is normal force, ω is angular velocity unit vector, R * is the equivalent radius. Directions of the moments by these models are different. F n is in the same direction as L i , while v = R * ω (linear velocity) is in the same plane as F t . The angle among F t and v is assumed to be θ. According to the right-hand rule and its geometric relationship, there is also an angle between torque M1 and M2, which is the same angle as θ.
In summary, both linear spring damping model and H-M can be expressed as "spring + damping", but their coefficients are different. The former is constant, and the latter is relevant to the overlap. Forces calculated by both models have little difference when the overlap is small, but forces calculated by H-M are larger than those by linear mode when the overlap is big. Therefore, it is suggested that a proper k of linear model should be selected according to the overlap range in the particle collision process. The value of k n with the highest probability of overlaps in the real collision process is appropriate. Torque of two models has an angle, which is equal to the angle between tangential force and linear velocity. As the particles used in this paper are glass beads, which have large stiffness and will not produce large overlap, a linear model will be adopted in the following simulation for its flexibility and simplicity for calculation.

Results and Discussion
In this part, the CFD-DEM model with and without the adhesion model are used to study the spout behaviors, which was validated by experiments at first. Subsequently, the effects of adhesion intensity and inlet gas velocity in spouted bed are investigated.

Validation by Literature and Experimental Results
The fluidization behavior of particles in the 2D spouted bed with different densities are studied by our group in previous studies [13]. It was found that the movement of glass beads had a bed oscillation cycle period of 165 ms in this literature. Therefore, our selfcoded model is verified in two ways: comparing with previous simulation results by commercial software and experimental results.
The results calculated by the linear spring-damping model (self-coded), Hertz-Mindlin model (EDEM 12.6, DEM Solutions Ltd., Edinburgh, England; commercial software) and the experimental results are given from left to right as shown in Figure 6. The calculation parameters by Hertz-Mindlin model were the same as literature [13]. It can be seen clearly that the overall is relatively consistent, such as the shape of the inlet as well as the spouted flow, the location of the spouted neck and the approximate spouted height.
The flow pattern of fluidization is almost the same, which can be divided into spouted zone, fountain zone and annular zone respectively. However, the details are slightly different, such as the dispersion of particles in the fountain zone, and the depression of particle accumulation in the annular zone. These differences may be attributed to the different contact force model models (self-coded: spring-damping model, EDEM: Hertz-Mindlin model).

Periodicity
It has been pointed out that typical characteristics of spouted beds are the periodicity of the vertical component of particle velocity as well as pressure drop and their periodicity are in agreement in literature [34]. The flow patterns of particle fluidization, where the color represents particle velocity are given in Figure 7. During spouting, a spout neck (red circle) gradually moves upwards from inlet to fountain area and finally disappears on the top of fountain. The spout neck has a periodic motion as well, with a period of 154 ms, which is close to the bed oscillation cycle time in Zhao (160 ms) [34] but a little smaller than that of our previous work (165 ms) [13]. The slight difference may be attribute to the particle-tofluid grid mapping relationship, resulting in a little difference of gas-solid drag force.

Particle Trajectory
The movement of a single particle can be obtained from CFD-DEM simulation results. Figure 8 shows a trajectory of a single particle. It can be clearly seen that the particle has passed through the spout zone, fountain zone and annulus zone. The residence times in each zone are also given in Figure 8. Also, periodicity can be found. The particle residence times are about 138, 220 and 1700 ms in the spout, fountain, and annulus zone, respectively. This result is related to the velocity characteristics of particles in each zone. It can be found that the flow pattern periodicity is close to the particle residence time in the spout zone. However, this periodicity of a single particle is not the same as that of flow patterns, because the flow pattern cycle is the statistical result of large number of particles, rather than a single particle.

Adhesion Model Validation
The adhesion model is based on the Luding's idea, which is a typical linear hysteresis model. Here, a head-on collision process of two particles with our model is obtained as shown in Figure 9. Here, dc is 2 µm and F adhesion is 8.2 × 10 −4 N. The left is the Hooke term (fn, "a" series), and the right is contact force term (Fn, including damping, "b" series) with spring, partially latching spring and hysteresis models. The key difference between these models lies in the Hooke term. In Figure 9(a1), the loading stiffness is the same as unloading stiffness, while those are different in Figure 9(a2) due to collision loss. The model in Figure 9(a3) is complex, and the whole process is divided into four stages with different stiffness as mentioned above. The stiffness of each stage has its characteristic parameters, like energy loss, adhesion force and rupture distance. It can be found that the framework of our force model is consistent with that of Luding's [16] and Liu's model [25], but it has its own feature for adhesion. As seen from the Figure 9, the minimum force is not huge, and the rupture distance is quite long, which means particles are more likely to keep adhesion.
For validation of the adhesion model, the experimental results from the literature [33] is shown as Figure 10 where the density of particles is 2380 kg/m 3 , weight of particles is 86.1 g, number of particles is 4698, liquid content is 1.313 × 10 −11 m 3 , superficial gas velocity is 1.08 m/s (corresponding to inlet gas velocity is 18.24 m/s). The geometric structure used here is shown in Figure 2. Figure 10 illustrates the longitudinal profile of vertical particle velocity along the spout axis by both experiment [33] and simulation. For the deduced part, C 1 is 0.966 kJ/(kg·K), C 2 is 1.142 kJ/(kg·K), q m is 15.48 m 3 /h, k temp is 10,000, ΔTp = 1000 K, h = 300 kJ/kg, Q d is a half of energy input, thus getting V = 1.790 × 10 −11 m 3 . All the curves have three similar stages under similar liquid contents (experiment: 1.313 × 10 −11 m 3 ; simulation: 1.313 × 10 −11 m 3 (Li et al. [35]), 1.790 × 10 −11 m 3 (our work)): rapid acceleration (stage I), slow acceleration (stage II) and fountain (stage III). The following figures show the flow patterns of our adhesion model, which has higher spouted height compared with the results without adhesion model. The cycle time is about 110 ms, which implies that adhesion may accelerate cycle period.

Effect of Adhesion Intensity
The interaction force between particles will be changed by adhesion, thereby affecting the fluidization. In order to study the effect of adhesion on the fluidization, simulation results were obtained under different adhesion intensity using our model, as shown in Figure 11. It should be noticed that the temperature controls the liquid content, thus controlling the adhesion force. It can be noted from Figure 11 that there are 7, 7.5, 7.5 cycles in 1 s in simulation results with no adhesion, slight adhesion and strong adhesion respectively, that means the bed oscillation cycle periods are in turn 143, 133 and 133 ms. This indicates that the adhesion will accelerate fluidization, and enhance intensity of fluidization. It is because adhesion can promote particle agglomeration. Under the same gas velocity, agglomerates strengthen the interaction on gas flow. Therefore, a larger contact force will be obtained, resulting in a higher spouted height and larger upward acceleration, so the periodic time is shortening.
There is a significant difference between spouted height within and without adhesion as shown in Figure 12a. The same tendency is reflected in the particle velocities from Figure  12b which further supports this argument. Particle fluidization with adhesion have obvious and regular fluctuations, while that without adhesion is relatively stable. It is suggested that slight adhesion will increase particle velocity as well as movement fluctuations. It is due to the enhancement of granular interaction on gas. As a result, particles are more likely to rise up. Besides, it is also seen from Figure 12b that adhesion will cause a phase difference in the vertical component of the velocity. Compared with Figure 8, Figure 13 shows that there are still three stages, and the residence time of jetting region gets shorter as well as the annulus region. Residence times of each region are about 122, 174 and 674 ms in turn.

Effect of Inlet Velocity
In the spouted bed, the fluidization state is determined by not only the adhesion force (affecting the interaction force), but also the drag force. The effect of inlet gas velocity on fluidization was studied here.
It can be noted that in the adhesive state, the pressure drop cycle was more obvious by increasing inlet gas velocity, as shown in Figure 14. The reason is that increasing inlet gas velocity will enhance the drag force of particles. Therefore, particles behavior is much closer to the free spouting state where drag force is dominant. But in fact, there is still energy loss due to adhesion, even though the gas kinetic energy increase. The bed oscillation cycle periods at three cases are 111, 143 and 105 ms, respectively. Within a certain range, a higher gas velocity will lead to a higher spouted height, thus prolonging cycle time. In this case, both flow pattern and fountain width increase with gas velocity increasing. It is found that the inlet gas velocity is positively relevant to the spout height at steady state, as shown in Figure 15a. It also suggests that a large gas velocity will cause a high initial holding pressure, which leads to a high spouted height of particles. However, the spouted height changes relatively smoothly at 26 m/s. For the three gas velocities, the final average spouted heights are about 0.175, 0.25 and 0.325 m, respectively. Figure 15b shows total (the above) and vertical component (the below) of particle velocity at 26, 35.2, 45 m/s. It is obvious that high gas velocity will lead to high particle velocity at the steady state. High fluctuations in vertical velocity also occur at 35.2 m/s, indicating a strong intensity of particle motion. The cases of other two inlet gas velocities are relatively mild. Therefore, there is not a single trend in the effect of gas velocity on particles. The inlet gas velocity has a close impact on particle collision and mixing, suggesting that the drag force plays a significant role in fluidization of particles. It can be seen from the Figure 16 that the residence times of each region are about 184, 277 and 357 ms, in turn. This indicates that in a certain scale, increasing inlet gas velocities can prolong the residence time in jetting region, making the spouted height higher, while the time in the annulus region is greatly reduced. It suggests that increasing inlet gas velocities makes the mixing of particles in different zones more effective.

Conclusions
An adhesion model is developed based on approaching-loading-unloading-detaching idea and particle surface change under high temperature in this paper. Analyses of the adhesion model are given through two particle collision process and validated by experiment. Effects of inlet gas velocity and adhesion intensity on spouted bed dynamics are investigated. A parameter sensitivity comparison of linear spring-damping model and Hertz-Mindlin Model is also conducted. Some conclusions can be put forth as follows: 1. Comparing the linear spring model with Hertz-Mindlin model, it is found that in case of small overlap, the forces calculated by both models have little distinction. In case of large overlap, the forces calculated by Hertz-Mindlin are larger than those of the linear model, and the distinction cannot be ignored. Thus, it is suggested that a proper k n of linear model should be selected according to the maximum probability of overlap in the particle collision process. There is an angle between torques of the two models, which is the same angle between tangential force and linear velocity.
2. An adhesion model with temperature is developed according to Luding's idea and integrated with force relationships in JKR. It is validated by experiment results. 3. It is concluded that fluidization cycle will be accelerated by adhesion, and intensity of fluidization will be enhanced. The presence or absence of adhesion is of great difference. 4. In a certain adhesive state, the pressure drop cycle was more obvious and periodicity got longer with increasing inlet gas velocity. The fountain width and fountain height will increase as well. It also leads to strong intensity of particle motion. 5. The effect of the particle surface adhesion on the spouted state will be further studied by the adhesion model developed based on the energy loss and balance in the future.

Acknowledgments:
The authors thanks for the support of coding discussion of Liu Daoyin in Southeast University, China.

Conflicts of Interest:
The authors declare no conflict of interest.
Nomenclature a contact area of radius: m 2 C1 specific heat capacity of particles, J/(kg·K) C2 specific heat capacity of gas, J/(kg·K) CD drag coefficient of a particle D the gap between the spherical surfaces, m D0 diameter of nozzle, m Dc width of spouted bed, m dc rupture distance, m dp particle diameter, m E* equivalent Young's modulus, Pa en normal coefficient of restitution between particles en_w normal coefficient of restitution between a particle and a wall et tangential coefficient of restitution between particles et_w tangential coefficient of restitution between a particle and a wall Fadhesion the adhesion force at overlap of zero, N Ff interaction force that the particle suffers from the surrounding gas, N f ij,min the maximum magnitude of the cohesive force, N Fn normal contact force, N F damping term of normal force, N F spring term of normal force, N Fp source term, N Ft tangential contact force, N F t d damping term of tangential force, N G shear modulus, Pa G* equivalent shear modulus, Pa h latent heat of melting, J/kg H height of spouted bed, m H0 particle loading, m I the identity tensor Ii moment of inertia, kg·m 2 kn normal spring constant between particles, N/m kn_w normal spring constant between a particle and a wall, N/m kload loading stiffness, N/m kunload unloading stiffness, N/m kt tangential spring constant between particles, N/m kt_w tangential spring constant between a particle and a wall, N/m ktemp the ratio of particle temperature rising to gas temperature drop L thickness of spouted bed, m Li the distance from the point of contact to the center of the sphere, m m* equivalent mass, kg M1 the torque calculated by linear spring damping model, N·m M2 the torque calculated by Hertz-Mindlin model, N·m mi mass of ith particle, kg Mij torque that the jth particle exerts on the ith particle, N·m Mlayer mass of liquid layer for total particles, kg N number of total particles p pressure, Pa qm mass flow rate of gas, kg/s Qd energy loss due to heat dissipation, J R radius of particles, m R* equivalent radius, m Re particle Reynolds number tCFD CFD time step, s tDEM DEM time step, s u gas velocity, m/s V volume of liquid layer, m 3 Vcell volume of the CFD cell, m 3 vf velocity of inlet gas, m/s vij,n normal relevant velocity of i th and j th particles, m/s vij,t tangential relevant velocity of i th and j th particles, m/s Fc contact force, N ncj number of contact particles Vp volume of a single particle, m 3 vp particle velocity, m/s ΔTg temperature change of gas, K ΔTp temperature change of particles, K Greek letters β drag force coefficient, kg/(s·m 3 ) γ inverted cone angle Γ surface tension, N/m δ ij,min the particle overlap corresponding to , , m δ max maximum overlap during collision, m δ n normal overlap of two particles, m δ p the permanent plastic deformation, m δ t tangential overlap of two particles, m ε local voidage η n coefficient of particle normal dissipation, N·s/m η t coefficient of particle tangential dissipation, N·s/m θ the angle between tangential force and linear velocity θ 0 contact angle between liquid and grain surface λ g bulk viscosity of gas flow, Pa·s ν Poisson's ratio µ cinematic viscosity of gas flow, Pa·s µ r the coefficient of rolling friction µ s coefficient of sliding friction between particles µ s_w coefficient of sliding friction between a particle and a wall ξ drag scaling factor ρ particle density, kg/m 3 ρ g gas density, kg/m 3 τ g stress tensor of gas φ p plasticity index ω angular velocity vector, 1/s