DEM Study of Wet Cohesive Particles in the Presence of Liquid Bridges in a Gas Fluidized Bed

A modified discrete element method (DEM) was constructed by compositing an additional liquid-bridge module into the traditional soft-sphere interaction model. Simulations of particles with and without liquid bridges are conducted in a bubbling fluidized bed.The geometry of the simulated bed is the same as the one in Müller’s experiment (Müller et al., 2008). A comparison between the dry and the wet particular systems is carried out on the bubble behavior, the bed fluctuation, and the mixing process. The bubble in the dry system possesses a regular round shape and falling of scattered particles exists while the bubble boundary of the wet particles becomes rough with branches of agglomerates stretching into it. Themixing of the dry system is quicker than that of the wet system. Several interparticle liquid contents are applied in this work to find their influence on the kinetic characteristic of the wet particle flow. With an increase of liquid content, the mixing process costs more time to be completed. Symmetrical profiles of the velocity and granular temperature are found for two low liquid contents (0.001% and 0.01%), while it is antisymmetrical for the highest liquid content (0.1%).


Introduction
Particular systems are very common in petroleum, food, chemical, and energy industries.These engineering processes involve highly complex phenomena of gas-solid twophase flow, which include particle mixing and segregation, formation of agglomerate, heat transfer, drying, and liquid transfer between particles in wet systems.The DEM method, originally proposed by Cundall and Strack [1], provides a robust technique for investigating the kinetic characteristic of particle flows.It is able to trace the trajectories and velocities of dispersed particles.The computational fluid dynamics (CFD) is commonly employed to account for the existence of the gas.Combination of the DEM and CFD methods has emerged as one of the most important tools in probing granular flows in spouted bed [2,3], gas-fluidized bed [4][5][6][7], vibrating bed [8,9], mixer [10,11], and hopper [12].
Most studies on particle fluidization are focused on dry systems.However, the system behavior in the presence of liquid is also very important from both industrial and scientific viewpoints.Due to the significantly varied characteristics on surface, the liquid bridge forces between particles and between a particle and a wall cannot be negligible compared to the inertia forces, especially at a scale of submillimeter or micrometer.Consequently, it causes the flow behavior of the wet particles greatly different from that of the dry particles and even defluidization phenomenon.Some challenging works have been reported by Seville and Clift [13], McLaughlin and Rhodes [14], and Wormsbecker and Pugsley [15] on the transition of particle flow from Geldart's [16] Group B to Group A and Group C due to an addition of small quantities of a nonvolatile liquid.In general, experimental work on wet particular system is still at a qualitative stage because of the great technical difficulty to directly measure forces in cohesive particle fluidized beds.Accordingly, the DEM method coupling with the CFD method becomes a predictive tool to investigate the detailed phenomena in fluidized beds of wet cohesive particles.Preliminary simulations applying the DEM-CFD method are conducted in centrifugal tumbling granulator by Muguruma et al. [17], bladed mixer by Lekhal et al. [18] and Radl et al. [19], and hopper by Anand et al. [20].Nevertheless, seldom literatures so far have been reported on the numerical simulation of wet cohesive particles in fluidized beds because the simulation of dry noncohesive particles has been the major issue of the last few decades.
As the mechanics of cohesive gas-particle systems is still poorly understood, the objective of the present work is to develop a numerical DEM-CFD method for the fluidization of wet cohesive particles in a gas-fluidized bed.At first, a modified soft-sphere interaction model was developed by compositing an additional module to account for the presence of liquid bridge.Then, a comparison between the dry and the wet particular systems is carried out and the influence of the interparticle liquid content is investigated on the kinetic characteristic of wet cohesive particle flow.

DEM-CFD Model
2.1.Governing Equations of the Gas Phase.The gas is treated as a continuous phase and the local mean Navier-Stokes equation, proposed by Anderson and Jackson [21], is used to describe its behavior.The local mean variables are derived from their point ones and averaged over the regions that are large compared to the particle size but small with respect to the characteristic dimension of the whole system.The resulting mass and momentum balances for the gas phase are as follows: where   is the void fraction of the gas in each computing cell, -;   is the gas density, kg⋅m −3 ; u  is the gas velocity, m⋅s −1 ;  is the pressure, Pa; F  is the rate of exchange of momentum between the gas phase and the particles, kg⋅m −2 ⋅s −2 ; g is the gravity, m⋅s −2 .The viscous stress tensor   is modeled by considering the gas phase as a Newtonian fluid: where   is the gas shear viscosity, Pa⋅s; I is the identity matrix, -.

Particle Dynamics.
In the DEM approach, the trajectories of individual particles are determined from the Newton's second law of motion.The time evolution of these trajectories then gives the global description of the whole granular system.The translational motion of a wet cohesive particle is governed by the pressure gradient force, the drag force, the gravitational force, the contact force, and the liquid bridge force: where   is the particle mass, kg⋅m −3 ; r is the position of particle center, m;  is the time, s;   is the particle volume,  m 3 ; v  is the particle velocity, m⋅s −1 ; F  is the contact force, N; F lb is the liquid bridge force, N. The angular velocity   under the action of a torque is given as follows: where   is the moment of inertia, kg⋅m 2 ; T  is the torque, N⋅m.The tangential forces due to both contact and liquid bridge interactions contribute to the torque T  .

Contact Force.
As shown in the dashed box of Figure 1, a linear spring-dashpot (LSD) soft-sphere model is applied to compute the force F  when a particle comes into contact with another or a wall.In the LSD model, a position-dependent spring obeying Hooke's law is used for describing the elastic deformation (energy conservation) and a velocity-dependent dashpot accounts for the damping effect (energy dissipation).
The normal component F  and the tangential component F  of the contact force F  are expressed as follows: where  is the spring stiffness, N⋅m −1 ;  is the damping coefficient, N⋅s⋅m −1 ;  is the elastic deformation, m; v  is the relative collision velocity, m⋅s −1 ;   is the coefficient of sliding friction, -.The subscripts  and  represent the normal and tangential directions, respectively.

Liquid Bridge Force.
The forces, arising in the wet cohesive particle flows, could be modeled using the concept of liquid bridge.To take account of the presence of liquid bridge, an additional module is added into the traditional soft-sphere interaction model and the whole model for wet cohesive particle is given in Figure 1.The liquid bridge force  lb consists of two parts, that is, the static capillary force  cp and the dynamic viscous force  V .The following assumption of liquid bridge is made for a wet particular system.(i) The total liquid of the wet system is uniformly distributed to each particle, which is covered by a liquid layer of the same thickness.The liquid volume on the surface of each particle is equivalent to six bridges (twelve half-bridges at maximum in a 3D coordinate system).
(ii) When the particles are directly contacting and an elastic formation (overlap) exits, within the region AB in Figure 2, both the liquid bridge force and the contact forces act on the particles.The whole softsphere model shown in Figure 1 is used to calculate all the interaction forces.
(iii) When the particles are separated (both approaching and separating) and the distance between them is less than the critical rupture distance  cr , within the region BC in Figure 2, a steady liquid bridge exists and the wet particles are only subject to this force in the absence of contact forces.
(iv) When the interparticle separation distance exceeds the critical rupture distance for a given liquid bridge volume, within the region CD in Figure 2, the liquid bridge collapses and the liquid returns back to the particle surface.Since the particles possess the same properties, the transfer of liquid between particles is neglected.This means, in the whole simulation, a particle carries the same content of liquid as in the initial state.
(v) The same manner is applied to the interaction between a particle and a wall, except that the volume of liquid bridge only comes from the liquid carried by the particle and the liquid returns back to the particle while the wall remains dry after collision.
(1) Capillary Force.(I) Particle-Particle Interaction.A total energy theory could be used to calculate the capillary force of a fixed liquid volume.When two wet particles are interacting as shown in Figure 3, the total energy  tot of the nonvolatile liquid bridge (constant volume), based on the work of Israelachvili [22], is expressed as follows: where  is the surface tension, N⋅m −1 ;  is the particle radius, m;  is the half-filling angle, rad;  is the contact angle, rad.The normal capillary force could be obtained from the derivation of the following: where  is the distance between the two particles, m.From the geometry of liquid bridge in Figure 3, the liquid volume between the two particles is approximately given by Under a constraint of constant volume for the liquid bridge,  lb / = 0, substituted by (8), From the geometry, the immersion height  has the following relationship: In the case of a small volume liquid bridge (i.e.,  ≪ 1), sin(/2) ≈ /2, and then Combining ( 7), (9), and ( 11), the capillary force can be expressed as The above expression of the capillary force comes from the pressure difference arising from the neck curvature and the influence of the liquid bridge surface tension is not considered.To take into account this effect, the whole normal capillary force is According to Lambert et al. [23] and Liu et al. [24], although an overestimate exists at a relatively high liquid volume, (13) gives a satisfactory prediction on the experimental data at a low liquid volume and is suitable for the conditions in the present work.Moreover, the immersion height  is a function of the interparticle distance  and the liquid bridge volume  lb .The expression proposed by Rabinovich et al. [25] is given by Combining ( 11) and ( 14), we get the relationship between the half-filling angle and the liquid bridge volume  lb : (II) Particle/Wall Interaction.Similar to the particle/particle interaction, the normal capillary force between a particle and a wall, shown in Figure 4 is given by When a direct contact occurs between two particles or between a particle and a wall, an overlap exits (within the region AB in Figure 2) and the interparticle distance  gets an unreasonable negative value in the above theory.In fact, there is an elastic deformation at the contact point.Liquid exists in the gap due to the particle roughness and the distance  is still positive.Here, we assume that the capillary force remains a constant value with a minimum distance  of 1 × 10 −5 m (particle roughness).
(2) Viscous Force.Besides the static capillary force, dynamically squeezing out and stretching the interparticle liquid result in a viscous resistance on the motion of particles.This effect becomes as significant as the capillary force if the relative collision velocity increases [26].Based on the lubrication theory [27] and the work of Adams and Edmonson [28], the normal viscous force is expressed as follows: where  lb is the liquid bridge viscosity, Pa⋅s.Lian et al. [29] suggest a solution of the viscous force in the tangential direction from the work by Goldman et al. [30]:   The same as the situation of the capillary force, a minimum value of 1 × 10 −5 m (particle roughness) is employed for the interparticle distance  to keep the viscous force bounded.
The liquid bridges remain in place until the particles reach a critical rupture distance, where the solid contact will cease.Lian et al. [31] found the following relationship between the critical rupture distance and the liquid bridge volume: where Vlb is the dimensionless liquid bridge volume:

Gas-Particle Interaction.
As the primary mode of the interphase momentum transfer, the drag between gas and particles couples the discrete simulation to the continuous fluid flow.The rate of the interphase exchange of momentum F  , computed by the sum of the drag forces acting on all individual particles in a computing cell, is expressed as a function of the product of the interphase drag coefficient and the relative velocities between the gas phase and the particles: where  cell is the computing cell volume, m 3 ;   is the number of particles in a computing cell, -;  is the interphase momentum exchange coefficient, kg⋅m −3 ⋅s −1 .It is extremely difficult to calculate the drag force theoretically.In this work, an expression from the correlations of Ergun [32] and Wen and Yu [33] is applied: where   is the particle diameter, m.The drag coefficient   is calculated from the correlations based on the particle Reynolds number: Re (1 + 0.15 Re 0.687 ) Re < 1000, 0.44 Re ≥ 1000. ( And the particle Reynolds number is defined as 2.4.Boundary Conditions.In this work, simulations of the particle flow with and without liquid bridges are conducted in a bubbling fluidized bed.The geometry of the cuboid bed is the same as that of the experimental bed by Müller et al. [34] and the dimensions are 44 mm × 10 mm × 120 mm.The total number of solid particles is 9240 with an initial distribution in the bed (0 s) shown in Figure 5.The monosize particles have a diameter of 1.2 mm and are all spherical, which is a little different from the kidney-like-shape particles in experiment.
The gas entering through the bottom distributor is uniform.For the gas phase, a no-slip boundary condition is applied on the walls.The finite volume approach is applied with discretization on a staggered grid for the hydrodynamic model.At first, all simulations of both the dry and the wet particular systems are carried out for 15 s and the time-averaged results are obtained from the last 10 s.Then, in order to calculate the minimum fluidization velocity  mf , another 10 seconds are conducted for the wet particle flows with the superficial gas velocity gradually decreasing to zero.The parameters and properties, the same as those in Müller's experiment, used in the simulation are summarized in Table 1.

Comparison between Simulation and Experiment.
The motion of dry noncohesive particles is first simulated.Figure 6 shows the 3D spatial distribution of the dry particles and the bubble predicted in the bed.The superficial gas velocity is set to 0.9 m/s.As we can see from the figure, a typical bubble forms at the bottom distributor, grows up in the middle, and collapses when reaching the open surface.A single bubble repeatedly occurs during the whole procedure of the simulation and it is comparable to the dimensions of the fluidized bed.In Figure 7, the transverse profile of the vertical particle velocity   is time averaged at three different heights, (a) 10 mm, (b) 20 mm, and (c) 30 mm.A comparison between the DEM simulation and the experiment using magnetic resonance (MR) by Müller is conducted.Combined with Figure 6, it indicates that the particles at the center are forced to move up and fall down near the side walls to occupy the void space arising from the up motion of the large bubble.The numerical simulation is found to agree well with the measured data except an underestimation near the two side walls, where the particle velocity increases in the experiment.This is because, in order to guarantee the validity of the flow field, the size of cells in DEM simulation is required to be several times greater than the particle diameter.This coarse numerical grid together with the no-slip boundary condition amplifies the wall shear effect and underestimates the gas velocity and then the particle velocity near the side walls.
Moreover, a small overestimation of the central particle velocity exists at 10 mm high while an underestimation presents at 30 mm high.One of the main reasons is the shape effect of particles on the bubble hydrodynamics, because the motion of a bubble strongly influences the kinetic characteristics of particles in the bubbling fluidized bed.The kidneylike-shape particles in the experiment may present a stronger resistance to the particle movement than the spherical particles in the dense phase zone.This leads to a lower particle velocity in the experiment than that in the simulation near the distributor.On the other hand, the stronger resistance from the kidney-like-shape particles generates a higher pressure in the bubble.Once the bubble collapses when reaching the open surface (nearly 30 mm under the present conditions), the kidney-like-shape particles will be sprayed out and get a higher velocity from the gas phase.Another reason for the higher particle velocity at 30 mm in the experiment may be the effect of the nonspherical particles on the drag force.At 30 mm, where the bubble reaches the open surface, the particles are mainly subjected to the drag from the gas and the interaction between the particles becomes less important.
The granular temperature, marking the fluctuating intensity of particles, is one of the most important parameters to measure the kinetic behavior of a particulate system.Jung et al. [35,36] proposed a granular temperature that is caused by the formations and motions of bubbles and gives rise to the normal Reynolds stresses.This granular temperature is related to the second moment of the velocity fluctuation: where  and  represent the directions , , and ;   is the number frames that are averaged in time, -.The timeaveraged mean velocity V is defined as And the cell-averaged mean velocity V  is defined as Finally the time-averaged granular temperature is given by Figure 14: Influence of liquid volume on pressure drop and bed height.The DEM prediction of the granular temperature   in the vertical direction is obtained at different heights, as shown in Figure 8. Two superficial gas velocities  = 0.6 m/s and  = 0.9 m/s are simulated in the fluidized bed.As we can see, the choice of the gas velocity has a great influence.A higher gas velocity increases the granular temperature, whose value nearly doubles at  = 0.9 m/s to that at  = 0.6 m/s.For both gas velocities, the maximum appears at the bed center where the bubble moves up.This indicates a more powerful interphase exchange of moment.A comparison is conducted between the DEM simulations and the experiment of Müller.
The numerical simulation gives an acceptable prediction of the MR measured data, since it is extremely hard to accurately calculate the granular temperature in fluidized beds [37,38].

Bed Behavior.
The bubble behavior of the wet cohesive particles is illustrated in Figure 9.The content of liquid added to the granular system is 0.01% by volume relative to the particles.Combined with Figure 6, the bubble behavior of the wet cohesive particles is found to be similar to the dry noncohesive particles in the fluidized bed.A single large bubble moves gradually from the distributor to the open surface, in the process of which the particles are pushed up at the head and entrained at the tail of the rising bubble.However, there are great differences between the dry and wet systems.The bubble possesses a regular round shape and falling of scattered particles exists in the bubble for the dry system.On the other hand, individual particles are gathered by agglomerates in and around the bubble for the wet system.Due to the strong viscous effect and the enhanced sliding friction arising from the liquid bridges, particles are not easy to move tangentially and the bubble boundary becomes rough.Branches of agglomerates, stretching into the bubble, contribute to the irregular characteristic of the bubble.
The computed fluctuation of the pressure drop and the bed height of the fluidized bed is shown for 15 s in Figure 10.As shown in the figure, we found that the fluctuation of the pressure drop in the wet system is much stronger than that in the dry system and the interparticle liquid contributes to a small increase of the average pressure drop in the wet system.This is because the rising bubble has to overcome an extra resistance between particles, which is from an additional cohesive force by the interparticle liquid bridges.The same reason is for the similar variance of the bed height besides the support of the wall with liquid bridges.The power spectral magnitude of the dynamic pressure as a function of frequency is shown in Figure 11.It is found that the main frequency for the dry and wet systems is 3.24 Hz and 4.80 Hz, respectively.Moreover, the magnitude of the fluctuations in the wet system is higher than in the dry system due to more violent pressure oscillations.

Mixing Index.
The mixing of particles in a fluidized bed is important as it greatly influences the mass, momentum, and energy exchange between particles and between a particle and the gas.The mixing extent could be quantified by a mixing index, which is based on the variance  2 of the sample particle composition: where  is the number of sample cells, -;   and  are the local and average concentration of particles, respectively.Lacey's mixing index  [39] is employed in this work and is defined as follows: where  2 0 is the variance of a complete segregated system, -;  2   is the variance of a randomly mixed system, -.For a binary mixture system,  2 0 = (1 − ) and  2  = (1 − )/;  is the proportion of the sample component in the whole mixture and  = 0.5 if the portions of the two components are equal.Lacey's mixing index ranges from 0 to 1 corresponding to completely random mixed and segregated states, respectively.
To investigate the mixing characteristic of a monosize particular system, the particles are divided into two groups with the same portions and traced by different colors (grey for lower half particles and black for upper half particles x (mm) Figure 17: Influence of liquid content on particle velocity ( = 0.9 m/s).
at the initial state).The bin size of the sample cells, used in calculating the mixing index, is the same as the size of simulation grid and has a dimension of 3.67 mm × 3.33 mm × 5.0 mm.The snapshots of the particle mixing evolution are illustrated in Figure 12 for a wet system with a liquid volume of 0.01%.As we can see, the particles are gradually mixed uniformly under a repeatedly occurrence and motion of a bubble.The quantitative comparison of the mixing index is shown in Figure 13.Both the dry and wet systems finally reach the same dynamic equilibrium state ( = 0.91).The dry system results in a larger slope of mixing to the wet system.It takes the dry system only 2.5 s to achieve the maximum Lacey mixing index while the wet system costs 3.3 s.

Influence of Liquid Content on Wet Particle Behaviors.
The influence of the liquid content on the pressure drop and bed height is plotted in Figure 14.Besides the first 15 seconds of operation of the wet cohesive bed, additional 10 seconds are conducted with the superficial gas velocity decreasing linearly from 0.9 m⋅s −1 to zero.From the figure, the standard deviation of the pressure fluctuation increases obviously to 31.1 Pa, 38.4 Pa, and 44.3 Pa with a relative liquid content of 0.001%, 0.01%, and 0.1%, respectively.This means a great enhancement of the pressure fluctuation is obtained by adding liquid into wet particle systems.This is because in a dry system, the contact force dominates the behavior of particles and, in a wet system with liquid, the capillary  force (attractive) and viscous force (repulsive) become more important.
The minimum fluidization velocity  mf is one of the most convenient indices to characteristics of the particle system in a fluidized bed.To determine  mf by the computer simulation, the superficial gas velocity was decreased gradually from 0.9 m⋅s −1 to zero during 10 s.The evolutions of the bed pressure drop and bed height are plotted in Figure 14.The particles are defluidized at 18.0 s, 18.5 s, and 19.3 s for liquid contents of 0.001%, 0.01%, and 0.1%, respectively.The corresponding minimum fluidization velocities  mf are determined as 0.63 m⋅s −1 , 0.585 m⋅s −1 , and 0.513 m⋅s −1 .When the gas velocity decreases below their minimum fluidizing velocity, the higher bed height with an increase of the liquid volume indicates a higher bed voidage.This is because the wet particles are supported by the wall with liquid bridges.
Figure 15 gives the results of liquid content on the mixing index.The final dynamic equilibrium state is independent from the liquid content.The increase of the liquid content significantly delays the mixing process of particles.It costs 2.6 s, 3.1 s, and 4.3 s for most particles reaching the maximum Lacey mixing index for liquid contents of 0.001%, 0.01%, and 0.1%, respectively.A great delay shows for the wet cohesive system and it is found that the bed behavior under this condition is significantly different from the other two (see the following part).
The spatial distributions of the time-averaged vertical velocities of the wet systems are plotted in Figure 16.It shows, from (a) to (c), an addition of liquid of (a) 0.001%, (b) 0.01%, and (c) 0.1%, respectively.The particle velocity near the distributor is found to be extremely low compared to the superficial gas velocity.In the simulation, a core of high positive velocity is generated at the center of the bed, aside which two cores of high negative velocities exist, for the liquid of 0.001% and 0.01%.This means a single bubble arises in the bed center to push the particles up and then transversely to the side walls.A different phenomenon occurs at a liquid content of 0.1%.A positive core and a negative core occupy either side of the walls, respectively.This is because the interparticle liquid bridge forces increase to such a great value that the particles move together as a group and the gas is pushed away to the corner to pass.
The profiles of the vertical velocities are shown in Figure 17 at different heights.A symmetrical profile presents both for the two low liquid contents, which is according to the spatial distributions in Figure 16.Moreover, an antisymmetrical profile occurs for 0.1% liquid because the bubble turns up at the right wall rather than at the center.The transverse profiles indicate that the particles move up along with the rising bubble and the particles around the bubble flow back to occupy the void space.A higher force between a particle and a wall arises if more liquid content is added.The particles become hard to move and their velocities decrease.A less movement of the near-wall particles decreases the supplement of the particles pushed away by a bubble and consequently reduces the velocity of particles there.This is according to the velocity variance in Figure 17 that an increase of liquid content decreases the scope of velocities.
The influence of the liquid content on the granular temperature   in the vertical direction is obtained at different heights, shown in Figure 18.Similar to the variances of the velocities, a symmetrical profile for 0.001% and 0.01% and an antisymmetrical profile for 0.1% also occur.For 0.001% and 0.01%, two kinds of granular temperature peaks are presented in the fluidized bed.The first kind is along the center of a bubble trajectory, where the interphase exchange between the particles and the gas is most acute.The second kind with two peaks, symmetrically distributing at both sides of the first kind, exists at the junction of the upward and downward particles where the collision between particles is most acute.There is only the first kind of peak for 0.1% because the particles are agglomerated as a group in the absence of collision between particles.An increase of the liquid decreases the granular temperature.

Conclusions
To investigate the hydrodynamics of the wet cohesive particles in a gas-fluidized bed, a modified DEM-CFD numerical method was developed by compositing an additional liquidbridge module into the traditional soft-sphere interaction model.Simulations of the dry noncohesive and wet cohesive particles are conducted in a bubbling fluidized bed.A comparison between the dry and wet particular systems is carried out on the bubble behavior, the bed fluctuation, and the mixing process.Different liquid contents are applied in this work to find their influence on the kinetic characteristic of the wet particular system.
Good agreements were found between the DEM simulation of the dry particles and the MR measurement of Müller.The bubble in the dry system is similar to the one in the wet system except that the bubble possesses a regular round shape and falling of scattered particles exists in the bubble for the dry system while the bubble boundary becomes rough and branches of agglomerates stretch into the bubble.The fluctuations of the pressure drop and bed height of the wet system are much larger than those of the dry system.The mixing of the dry system is found to be quicker than that of the wet system and the mixing process costs more time to be completed with an increase of the liquid content.An increase of the liquid content decreases the scope of the particle velocities.Symmetrical profiles of the velocity and granular temperature are shown for the two low liquid contents (0.001% and 0.01%) while it is antisymmetrical for a higher liquid content of 0.1%.

Figure 3 :Figure 4 :
Figure 3: Schematic of liquid bridge between two particles.

Figure 5 :
Figure 5: Schematic of the bubbling fluidized bed and initial particle positions.

Figure 15 :
Figure 15: Influence of liquid volume on mixing index.

Figure 18 :
Figure 18: Influence of liquid content on granular temperature.

Table 1 :
Parameters applied in the simulation.