Review of CFD-DEM Modeling of Wet Fluidized Bed Granulation and Coating Processes

Abstract

Wet fluidized bed granulation and coating processes have been widely used in the pharmaceutical and food industries. The complex gas–solid flow coupled with heat and mass transfer in such processes made it hard to form complete control over the apparatuses. To serve better design, scaling-up, and optimization of granulators and coaters, the underlying micro-scale mechanisms must be clarified. Computational fluid dynamics coupled with the discrete element method (CFD-DEM) provides a useful tool to study in-depth the gas-solid hydrodynamics of the granulation and coating processes. This review firstly introduced the fundamental theory of CFD-DEM from governing equations, force calculation, and coupling schemes. Then the application of CFD-DEM in simulating wet fluidized bed granulation and coating was presented. Specifically, the research focus and the role of CFD-DEM in resolving issues were discussed. Finally, the outlook on the development of CFD-DEM in the context of granulation and coating was given.

1. Introduction

The coating process has been widely used in many industries such as the pharmaceutical, chemical, and food industries [1,2]. An appropriate coating can significantly improve the quality of the final particulate products by controlling the release of active ingredients, protecting the particles from external impurities, and masking the undesired taste. From the perspective of particles, a typical coating process primarily involves collision with the sprayed droplets (particle wetting) and then drying under an ambient environment or reinforced drying via a hot air flow [3]. The coating grows gradually under multiple wetting–drying cycles.

Fluidized bed granulators or coaters have been widely used in past decades to provide high-quality final products. However, their wide application largely depends on the personal experience of engineers. There is still a big gap between engineering application and granulation/coating theories. Representative unresolved problems in fluidized bed granulation/coating include gas–solid hydrodynamics under special conditions (such as fine or cohesive powders), heat and mass transfer (temperature and moisture distribution), liquid bridge force modeling, and microscopic particle growth mechanisms. Many researchers have attempted to unravel the veil of granulation and coating processes using analytical methods [4,5,6,7], experimental methods [8,9,10,11,12,13,14,15,16,17,18,19], numerical simulation [20,21,22,23,24,25,26,27,28,29,30,31,32,33], or a combination of experimental and numerical methods [34,35,36,37,38,39,40,41,42,43,44,45]. At the early stage, experimental research was largely conducted on a macro/global scale and confined to a specific apparatus (see Figure 1). This might help form a particular understanding of the specific particulate processes, whilst the lack of micro-dynamics information on particles made it difficult to truly understand the relationship between the micro-mechanisms and the macro-phenomena [46]. Although some experimental techniques, such as tracer particle techniques [47,48,49,50,51], have already been used to capture the microscopic motion characteristics of particles, the cost was fairly high. This sometimes puts researchers in a dilemma where it is very difficult, if not impossible, to completely disclose the detailed micro-dynamic information of particles at a low cost merely by experimental techniques.

Fortunately, numerical simulation techniques, especially those discrete ones, provide an effective tool for us to gain insight into the micro-dynamic information of particles [52]. Straightforwardly, the discrete element method (DEM) has been gaining popularity because it directly calculates the per-atom properties such as trajectory and velocity. This would benefit the microscopic mechanistic research of the complex gas–solid flows. In its early stage, DEM was not so popular due to the big computation cost. However, its simulation scale has been largely broadened with the rapid development of DEM software and computer hardware. Currently, with the aid of super-computers, DEM can be even used to simulate some industrial processes.

Figure 1. Representative configurations of wet fluidized bed granulators and coaters: (a) Top spray granulator, (b) bottom spray granulator, (c) Wurster granulator, (d) rotor granulator, and (e) spouted bed granulator. Reproduced from reference [53] permitted by Elsevier.

Figure 1. Representative configurations of wet fluidized bed granulators and coaters: (a) Top spray granulator, (b) bottom spray granulator, (c) Wurster granulator, (d) rotor granulator, and (e) spouted bed granulator. Reproduced from reference [53] permitted by Elsevier.

Two types of DEM approaches are most commonly seen in the open literature, i.e., the hard-sphere model and the soft-sphere model. The hard-sphere model assumes the particles to be ideally rigid and the particle collisions are instantaneous. It does not explicitly calculate the collisional force between particles. Note that the hard-sphere model only considers binary collisions and cannot handle multi-body collisions [54]. Obviously, this largely confines its application to rapid and dilute granular flows. By contrast, the soft-sphere model originally proposed by Cundall and Strack [55] allows for slight particle deformation based on which collision forces are explicitly calculated. Then the motion of particles can be calculated according to Newton’s laws of motion. The soft-sphere model is capable of dealing with multi-body collisions and has been extensively used in various industrial processes such as hopping [56,57,58,59,60], blending [61,62,63], and granulation [31,32,64,65,66,67,68]. Due to the large-scale application, in the following remainder of the article, the term ‘DEM’ refers to the soft-sphere approach by default if not noted elsewhere.

The coupling with CFD has further broadened the application of DEM since there are two or more phases in some industrial processes such as fluidization, pneumatic conveying, and hydro-cyclones. In general, the jargon ‘CFD-DEM’ indicates that the volume-averaged RANS CFD is used to describe the fluid flow which is treated as a continuum while DEM is used to track the motion of particles which is treated as a discrete phase. This treatment, albeit quite natural, is not the only choice. In terms of modeling scale, there can be many variants. For example, the simulation scale for fluid flow can span from discrete molecular dynamics simulation (MDS), the lattice Boltzmann method (LBM), the pseudo-particle method (PPM) to continual direct numerical simulation (DNS), large eddy simulation (LES), and the conventional RANS CFD. Likewise, particles can also be simulated as a continuum such as the popular two-fluid model (TFM). Hence, other combinations can also be seen in the literature, such as TFM, LES-DEM, DNS-DEM, and LB-DEM. Nevertheless, CFD-DEM is still one of the most popular simulation techniques for gas–solid flows because it could capture the micro-dynamics of particles while resolving fluid flow at a low computation cost [54,69,70].

This work serves to create an entry for researchers or engineers who are interested in using CFD-DEM to simulate wet fluidized bed granulation or coating processes. Firstly, the fundamental theory of CFD-DEM was introduced, including governing equations, force calculation, and coupling schemes. Then, the application of CFD-DEM in the granulation and coating processes was presented where the research focus and the role of CFD-DEM in solving these problems were discussed. Finally, conclusions and our outlook on the development of CFD-DEM in the context of granulation and coating were provided.

2. Governing Equations

Since CFD-DEM is a combination of the volume-averaged RANS CFD and DEM, the governing equations that describe the gas flow and particle motion can be discussed in sequence.

2.1. CFD

In the open literature, two types of governing equations can be seen, i.e., Model A and Model B as follows:

Model A:

$$\frac{\partial {\epsilon }_f}{\partial t}+\nabla \cdot \left({\epsilon }_f{\overrightarrow{v}}_f\right)=0$$

(1)

$$\frac{\partial \left({\epsilon }_f{\rho }_f{\overrightarrow{v}}_f\right)}{\partial t}+\nabla \cdot \left({\epsilon }_f{\rho }_f{\overrightarrow{v}}_f{\overrightarrow{v}}_f\right)=-{\epsilon }_f\nabla p+\nabla \cdot \left({\epsilon }_f{\overline{\overline{\tau }}}_f\right)-{\overrightarrow{F}}^A+{\epsilon }_f{\rho }_f\overrightarrow{g}$$

(2)

Model B:

$$\frac{\partial {\epsilon }_f}{\partial t}+\nabla \cdot \left({\epsilon }_f{\overrightarrow{v}}_f\right)=0$$

(3)

$$\frac{\partial \left({\epsilon }_f{\rho }_f{\overrightarrow{v}}_f\right)}{\partial t}+\nabla \cdot \left({\epsilon }_f{\rho }_f{\overrightarrow{v}}_f{\overrightarrow{v}}_f\right)=-\nabla p+\nabla \cdot \left({\epsilon }_f{\overline{\overline{\tau }}}_f\right)-{\overrightarrow{F}}^B+{\epsilon }_f{\rho }_f\overrightarrow{g}$$

(4)

where ${\epsilon }_f$ is the volume fraction of fluid phase in a computation cell, $p$ is the fluid pressure, and ${\overline{\overline{\tau }}}_f$ is the viscous stress tensor of fluid; for incompressible flows, ${\overline{\overline{\tau }}}_f$ can be computed by:

$${\overline{\overline{\tau }}}_f={\mu }_f\left(\nabla {\overrightarrow{v}}_f+{\left(\nabla {\overrightarrow{v}}_f\right)}^T\right)$$

(5)

${\overrightarrow{F}}^A$ and ${\overrightarrow{F}}^B$ are the volume-averaged fluid–solid interaction forces, which can include different terms under different situations. Note that these two terms can be mutually converted by:

$${\overrightarrow{F}}^B={\overrightarrow{F}}^A/{\epsilon }_f-{\rho }_f{\epsilon }_s\overrightarrow{g}$$

(6)

To date, there has not been any agreement on which model is better. Feng and Yu [71] and Kafui et al. [72] proposed that there was little difference between the two models when they were applied to model the mono-sized particles. Subsequently, Feng and Yu [73] investigated the fluidization of binary particle mixtures and found the two models displayed a significant discrepancy. By comparing them with the experimental results, they suggested Model B should be favored when modeling a binary particle mixture. To resolve this important issue, Zhou et al. [70] sourced the origins of the models and analyzed their applicability to three representative gas–solid systems, i.e., fluidization, pneumatic conveying, and hydrocyclones. They proposed that there should be three sets of governing equations, as shown in Equations (7)–(12), rather than two, to describe gas–solid flows.

Set I:

$${\rho }_f{\epsilon }_f\left(\frac{\partial {\overrightarrow{v}}_f}{\partial t}+\nabla \cdot \left({\overrightarrow{v}}_f{\overrightarrow{v}}_f\right)\right)=\nabla \cdot {\overline{\overline{\tau }}}_f-n{\overrightarrow{f}}_i+{\rho }_f{\epsilon }_f\overrightarrow{g} \left(fluid phase\right)$$

(7)

$${\rho }_s{\epsilon }_s\left(\frac{\partial {\overrightarrow{v}}_s}{\partial t}+\nabla \cdot \left({\overrightarrow{v}}_s{\overrightarrow{v}}_s\right)\right)=n \overrightarrow{\mathsf{\Phi }}-\nabla \cdot S+n{\overrightarrow{f}}_i+{\rho }_s{\epsilon }_s\overrightarrow{g} \left(solid phase\right)$$

(8)

Set II:

$${\rho }_f{\epsilon }_f\left(\frac{\partial {\overrightarrow{v}}_f}{\partial t}+\nabla \cdot \left({\overrightarrow{v}}_f{\overrightarrow{v}}_f\right)\right)={\epsilon }_f\nabla \cdot {\overline{\overline{\tau }}}_f-n{\overrightarrow{f}}_i^{\prime }+{\rho }_f{\epsilon }_f\overrightarrow{g} \left(fluid phase\right)$$

(9)

$${\rho }_s{\epsilon }_s\left(\frac{\partial {\overrightarrow{v}}_s}{\partial t}+\nabla \cdot \left({\overrightarrow{v}}_s{\overrightarrow{v}}_s\right)\right)={\epsilon }_s\nabla \cdot {\overline{\overline{\tau }}}_f+n{\overrightarrow{f}}_i^{\prime }+{\rho }_s{\epsilon }_s\overrightarrow{g}+\nabla \cdot {\overline{\overline{\tau }}}_s \left(solid phase\right)$$

(10)

Set III:

$${\rho }_f{\epsilon }_f\left(\frac{\partial {\overrightarrow{v}}_f}{\partial t}+\nabla \cdot \left({\overrightarrow{v}}_f{\overrightarrow{v}}_f\right)\right)=\nabla \cdot {\overline{\overline{\tau }}}_f-\left[n{\overrightarrow{f}}_i^{\prime }/{\epsilon }_f-{\rho }_f{\epsilon }_s\overrightarrow{g}\right]+{\rho }_f{\epsilon }_f\overrightarrow{g} \left(fluid phase\right)$$

(11)

$${\rho }_s{\epsilon }_s\left(\frac{\partial {\overrightarrow{v}}_s}{\partial t}+\nabla \cdot \left({\overrightarrow{v}}_s{\overrightarrow{v}}_s\right)\right)=\nabla \cdot {\overline{\overline{\tau }}}_s+\left[n{\overrightarrow{f}}_i^{\prime }/{\epsilon }_f-{\rho }_f{\epsilon }_s\overrightarrow{g}\right]+{\rho }_s{\epsilon }_s\overrightarrow{g} \left(solid phase\right)$$

(12)

where $\overrightarrow{\mathsf{\Phi }}$ is the local averaged particle–particle interaction force, S is the ‘Reynolds stress tensor’ of solid phase, ${\overrightarrow{f}}_i$ is the local averaged force on particle i exerted by the surrounding fluid, n is the number of particles per unit volume, ${\overline{\overline{\tau }}}_s$ is the stress tensor of the solid phase, and ${\overrightarrow{f}}_i^{\prime }$ is the local averaged particle–fluid interaction force arising from the velocity fluctuations as the fluid passes around the particles or through the interstices among the particles.

The authors recognized that set I was the original form of the governing equations in the continuum approach, and sets II and III, corresponding to Model A and Model B, respectively, were all derived from set I. They also found that sets II and III were essentially the same except for some minor differences caused by the different mathematical treatments of a few terms in the original model. Moreover, set III was found to be a simplified form of set I and should only be applied in steady and uniform flow. Finally, sets I and II were recommended to model complex particle-fluid flow.

2.2. DEM

A given particle generally has two types of motion, i.e., transitional and rotational. In a gas–solid system, particles may collide with the neighboring particles or walls that are in contact with them. Furthermore, the particles would interact with the surrounding fluid and momentum, and energy can be exchanged between the two phases. When computing the motion of particles, most commercial or open-source DEM software only considers the interaction with the closest particles or fluid by default. This surely makes sense, but the underlying assumption that should be pointed out is that the effect of disturbance waves derived from some distant particles or fluid is ignored. Thus, in practical DEM calculations, the timestep should be controlled under a critical value within which the disturbance waves cannot spread over a particle or CFD cell scale. Given all the immediate neighboring forces, the motion of particles (particle i as an example) can be described by Newton’s laws of motion as follows:

$$m_i\frac{d{\overrightarrow{v}}_i}{dt}={\overrightarrow{F}}_{con,i}+{\overrightarrow{F}}_{ncon,i}+{\overrightarrow{F}}_{f,i}+m_i\overrightarrow{g}$$

(13)

$$I_i\frac{d{\overrightarrow{\omega }}_i}{dt}={\overrightarrow{T}}_i$$

(14)

where $m_i$ is the mass of particle i, ${\overrightarrow{F}}_{con,i}$ and ${\overrightarrow{F}}_{ncon,i}$ are the total contact force and total non-contact force acting on particle i, respectively, ${\overrightarrow{F}}_{f,i}$ is the interaction force with the surrounding fluid, $I_i$ is the moment of inertia of particle i, ${\overrightarrow{\omega }}_i$ is the angular velocity of particle i, and ${\overrightarrow{T}}_i$ is the total torque on particle i.

Note that Equations (13) and (14) are only a generalized framework to calculate the motion of particles in DEM. The detailed formulas for the calculation of the forces and torques are discussed in Section 3.1.

3. Force Calculation

In this section, the calculation of the particle–particle and particle–wall contact forces is discussed first, followed by a discussion of non-contact forces. To be more specific, the non-contact forces in wet fluidized beds primarily involve the liquid bridge force. Finally, the interaction forces between the particle phase and the fluid phase are covered in detail.

3.1. Contact Forces

If a particle system is not very dilute, there is a big chance that the particles can contact/collide with each other. The contact area is not a single point. Rather, it is usually a finite area due to particle deformation, which is represented by overlaps in DEM. To calculate the contact force explicitly, it is generally decomposed into two components: One is the normal force, which is normal to the contact plane; another is the tangential force, which lies in the contact plane. Then the estimation of the two components can be performed by the corresponding models, respectively.

Typically, a spring-dashpot system (Figure 2) is used to represent the particle collision for the sake of contact force calculation. The spring stands for elastic deformation and the dashpot represents the viscous dissipation of the mechanical energy. The most basic and intuitive model, called the linear spring-dashpot (LSD) model, assumes that the normal force is proportional to the normal displacement. This model, considered simple and computation-saving, was first used by Cundall and Strack [55] as they introduced the soft-sphere model in 1979. Due to its simplicity and easy implementation, the LSD model has been widely applied. However, the LSD model was challenged due to some theoretical flaws. For example, the contact force at the start and end of collisions is non-zero in the LSD model, which is physically unrealistic and should be avoided in DEM implementation [74]. Moreover, the contact time and the coefficient of restitution have little to do with the impact velocity in this model, but this point was defied by some experimental results [75,76]. Partly due to these theoretical flaws, certain results [77] given by the LSD model are less accurate when compared with those non-linear contact force models. However, this does not necessarily mean the LSD model should be abandoned when performing DEM simulation because there are still cases [78,79,80,81] found in open literature where the results of the LSD model comply very well with the experimental results. Hence, to be on the safe side, we advise the LSD model should be used with caution in DEM simulation. Some efforts have been made to remedy the flaws of the LSD model, such as Kruggel-Emden et al. [75] and Maio et al. [82], while these modified LSD models have their own limitations and thus have not been widely used yet.

Nonlinear viscoelastic models arise from the fact that the normal displacements and the normal contact force do not obey a linear relationship, which has been verified by various experiments and finite element analysis. Currently, the most popular force-displacement model is the Hertz model [83] and its modifications [84,85,86,87]. Often, the Hertz model is coupled with the Mindlin and Deresiewicz model [88], which accounts for the tangential contact force. According to these two models, the contact force not only depends on the instantaneous change rate of the displacements but also on the loading/unloading history of the particles. The Hertz model in conjunction with the Mindlin and Deresiewicz model has been widely thought to be accurate and reliable. Based on these two models, some modified models have also been proposed either to increase the computation accuracy, such as Maw et al. [89], Stevens et al. [76], and Thornton et al. [90], or to mitigate thecomputation burden, such as Thornton et al. [91], Langston et al. [58], Zhou et al. [92], and Zhu et al. [93]. Essentially, all these nonlinear viscoelastic models use the Hertz model to account for the elastic effect while applying different models to represent the viscous effect, and thus particle velocities after a collision are different. Force models on the normal and tangential directions have been frequently used miscellaneously (such as Schafer et al. [94] and Latzel et al. [95]), but this is sometimes questionable since most of the models only consider one or two aspects and neglect the interactive effects of many factors. Nonetheless, nonlinear viscoelastic models have successfully avoided the aforementioned flaws of the LSD model. In addition to viscoelastic models, there are also certain models that take the plastic effect into account such as Thornton et al. [96] and Vu-Quoc et al. [87], but these models are not as commonly used as viscoelastic models.

For realistic particle contact, the contact force occurs within a finite area rather than the mass center of the particles, and this would impose a torque on the particles and induce their rotation. Determining the torques is thus becoming the key to accurately describing the rotational motion of particles. The torque is usually composed of two components: The tangential contact force and the asymmetric normal contact force. The contribution of tangential force to the torque is intuitive and relatively easy to estimate. In contrast, the distribution of the asymmetric normal contact force is fairly difficult to determine because it depends on many factors such as particle morphology, material properties of the particles, and instantaneous impact status. DEM packages usually ignore this contribution (also called the rolling friction torque) by default, but in some cases, this may not make sense. For example, in the simulation of sand heaping/piling, the rolling friction torque plays an important role in the transition between static and dynamic status. The neglect of the rolling friction torque may lead to inaccuracies or even errors in the estimation of some important parameters such as the repose angle. Hence, the determination of the rolling friction torque needs further research in the future.

3.2. Non-Contact Forces

As mentioned before, the non-contact forces in wet fluidized bed granulators and coaters primarily involve the liquid bridge force and van der Waals force. However, the effect of the van der Waals force is only significant when the particle size is small enough, and hence we focus on the discussion of the liquid bridge force in this section.

When two particles with free liquids on the surface come into ‘contact’, a liquid bridge may form instantly or at a finite rate. This bridge often applies an attraction force to the particles, and sometime this force may play a dominant role in the motion of the particles. Therefore, it is worth investigating the liquid bridge force.

Basically, the liquid bridge force is composed of two components: Capillary force and viscous force. Calculation of the two components would be given in sequence. As shown in Figure 3, the capillary force can be estimated theoretically as follows:

$${\overrightarrow{F}}_{ij}^{cap}=\left(2\pi R\sigma \sin \varphi \sin \left(\varphi +\theta \right)-\pi R^2\Delta p {\sin }^2\varphi \right){\overrightarrow{n}}_{ij}$$

(15)

where R is the harmonic mean radius, and it can be calculated as:

$$R=2{\left(\frac{1}{R_i}+\frac{1}{R_j}\right)}^{-1}$$

(16)

$\sigma$ is the surface tension of the liquid, ${\overrightarrow{n}}_{ij}$ is the unit vector pointing from the center of particle i to the center of particle j, and $\Delta p$ stands for the contribution of the pressure gradient and can be calculated by:

$$\Delta p=\sigma \left(\frac{1}{R_i}+\frac{1}{R_j}\right)$$

(17)

In principle, an implicit equation can be obtained by coupling Equations (15) and (17) through which the capillary force, angle $\varphi$, angle $\theta$, and the mean curvature of the bridge $\mathsf{\Gamma }$ can be related. Fisher [97] first studied the capillary force between two particles via a toroidal model. Following his work, Hotta et al. [98] introduced an explicit approximation (called the ‘gorge’ method), which is written as:

$${\overrightarrow{F}}_{ij}^{cap}=\pi \sigma {\rho }_2\left(1+\frac{{\rho }_2}{{\rho }_1}\right){\overrightarrow{n}}_{ij}$$

(18)

$${\rho }_1=\frac{l/2+R\left(1-\cos \varphi \right)}{\cos \left(\theta +\varphi \right)}$$

(19)

$${\rho }_2=R\sin \varphi -\left(1-\sin \left(\theta +\varphi \right)\right){\rho }_1$$

(20)

The ‘gorge’ method has been reported to improve the calculation accuracy where the error could be less than 10%. This is acceptable for most engineering applications. However, Mazzone et al. [99] found that this method might underestimate the capillary force at large separation distances. Later, Lian et al. [100] advised that the ‘gorge’ method could be improved by applying scaling coefficients. The models mentioned above largely rely on the implicit analytical method. However, for the sake of DEM implementation, explicit methods, albeit sometimes with lower accuracy, are much more preferred.

Figure 3. A schematic of meniscus liquid bridge between two particles. Adapted from [101] permitted by John Wiley and Sons.

Figure 3. A schematic of meniscus liquid bridge between two particles. Adapted from [101] permitted by John Wiley and Sons.

Many explicit approximate approaches have been proposed, such as Simons et al. [102] and Willett et al. [103]. The explicit models could be easily implemented in DEM and provide applicable results. On the other hand, the computation errors increased dramatically as the volume ratio of liquid to particles increased. This brought much uncertainty when a relatively large amount of liquid was injected into the system.

It is noteworthy that Mikami et al. [104] systematically investigated the prediction of the capillary force through regression analysis. In their research, the capillary force was related to several important parameters of the liquid bridge, and the equations are shown below:

$${\overrightarrow{F}}^{cap}=\pi R\sigma \left(\exp \left(A{l}^{*}+B\right)+C\right)\overrightarrow{n}$$

(21)

$l^{*}$ is the dimensionless separation distance defined by:

$$l^{*}=\frac{l}{R}$$

(22)

$l$ is the separation distance. A, B, and C are the functions of the dimensionless liquid bridge volume $V^{*}$ and the contact angle $\theta$. $V^{*}$ is defined as:

$$V^{*}=\frac{V_{lb}}{R^3}$$

(23)

$V_{lb}$ is the liquid bridge volume. For particle–particle collision, A, B, and C are calculated as:

$$A=-1.1V^{*-0.53}$$

(24)

$$B=\left(-0.34\ln V^{*}-0.96\right){\theta }^2-0.019\ln V^{*}+0.48$$

(25)

$$C=0.0042\ln V^{*}+0.078$$

(26)

For particle–wall collision, they are calculated as:

$$A=-1.9V^{*-0.51}$$

(27)

$$B=\left(-0.016\ln V^{*}-0.76\right){\theta }^2-0.12\ln V^{*}+0.12$$

(28)

$$C=0.013\ln V^{*}+0.18$$

(29)

Based on their work, Willett et al. [103] and Shi and McCarthy [105] investigated the prediction of the capillary force between particles of unequal size. Willett et al. [103] suggested that the Derjaguin method [106] could be applied to alleviate the computation burden in the case of different-sized particles. However, this approximation may display a non-negligible error when the separation distance is small for contact or large for rupture.

The research by Maugis [107] should also be mentioned where the liquid bridge was assumed to be cylindrical and a simple equation was proposed to calculate the capillary force:

$${\overrightarrow{F}}_{ij}^{cap}=2\pi R\sigma X_v\cos \left(\theta \right){\overrightarrow{n}}_{ij}$$

(30)

$$X_v=1-\frac{1}{\sqrt{1+2V_{lb}/\left(\pi R{l}^2\right)}}$$

(31)

Butt and Kappl [108] reviewed the estimation of normal capillary forces where the effects of particle shape and surface roughness were discussed. Detailed information can be found in their work.

The viscous force arises from the liquid viscosity. However, the viscous force is ‘dynamic’ and would also be influenced by the relative velocity of the colliding particles. The calculation of the viscous force is usually based on the lubrication theory [109] where the equation is given by:

$${\overrightarrow{F}}_{ij}^{vis,n}=\frac{3}{2}\pi \mu R^2\frac{X_V^2}{l}{\overrightarrow{v}}_{ij}^n$$

(32)

where ${\overrightarrow{v}}_{ij}^n$ is the normal component of the relative velocity of the particles. Note that in the formula above, the viscous force would extend to infinity as the separation distance approaches zero. To avoid this issue, a ‘cut-off’ (minimum) distance is usually applied in DEM implementation.

Some studies have been conducted to investigate the relative effect of the capillary force and the viscous force. Based on the capillary number (Ca), Ennis et al. [110] concluded that the viscous effect is insignificant when Ca is smaller than 10⁻³. Pitois et al. [109] calculated the capillary force and the viscous force using Equations (30) and (32), respectively, and then compared the results with experimental data. It was found that the calculated results agreed well with the experimental data. Likewise, they also found that the viscous force was prominent only at small separation distances, and the capillary force would play a major role in large separation distances.

Most simulation work omitted the tangential viscous force by default based on the assumption that it was trivial compared with the normal component. Goldman et al. [111] proposed that the tangential component could be estimated by:

$${\overrightarrow{F}}_{ij}^{vis,t}=6\pi \mu R\left(\frac{8}{15}\ln \left(\frac{R}{l}\right)+0.9588\right){\overrightarrow{t}}_{ij}$$

(33)

Goldman et al. [111] and Pitois et al. [109] showed that at a small separation distance, the normal component is indeed much larger than the tangential one, while the latter is larger than the former when the separation distance surpasses ~12% of the particle radius. This seems to imply that the tangential viscous force cannot always be omitted.

The liquid bridge would rupture and the liquid would redistribute to the particles or walls if the separation distance was beyond a critical value, which is called the rupture distance. This value is used to judge whether the liquid bridge still exists in DEM implementation and is thus an important parameter. Lian et al. [100] committed a theoretical study and quantified the rupture distance using the following formula:

$$\frac{l_{rup}}{R}=\left(1+\frac{\theta }{2}\right){\left(\frac{V_{lb}}{R^3}\right)}^{1/3}$$

(34)

This equation shows that the rupture distance increases with the liquid bridge volume. Pitois et al. [112] argued that the rupture distance should also depend on the viscosity of the liquid and the relative velocity of the colliding particles. They conducted experiments and proposed a modified version as follows:

$$\frac{l_{rup}}{R}=\left(1+C{a}^{1/2}\right)\left(1+\frac{\theta }{2}\right){\left(\frac{V_{lb}}{R^3}\right)}^{1/3}$$

(35)

In this model, they introduced the Ca number to represent the effect of the viscosity of the liquid and the relative velocity of the colliding particles. Equations (34) and (35) are consistent when the viscosity of the liquid or the relative velocity approaches zero. Again, Mikami et al. [104] investigated the rupture distance through regression analysis:

For particle–particle:

$$l_{rup}^{*}=\left(0.62\theta +0.99\right)V^{*0.34}$$

(36)

For particle–wall:

$$lrup\ast =0.22\theta +0.95V\ast 0.32l_{rup}^{\ast }=\left(0.22\theta +0.95\right)V^{*0.32}$$

(37)

The $l_{rup}^{*}$ is the dimensionless rupture distance defined by:

$$l_{rup}^{*}=\frac{l_{rup}}{R}$$

(38)

Willett et al. [103] investigated the rupture distance when the spheres had unequal size and the formula was given as:

$$l_{rup}^{*}=\left(1+\frac{\phi }{4}\left(\frac{R_2}{R_1}+1\right)\right)\left(V_{lb}^{*1/3}+\left(\frac{R_2}{2R_1}-\frac{2}{5}\right)V_{lb}^{*2/3}\right)$$

(39)

In practical DEM implementation, the liquid bridge volume is sometimes assumed to be a fixed percent (e.g., 1%) of the volume of particles for simplicity. To obtain higher accuracy, some researchers also proposed formulas to calculate the bridge volume. Pietsch and Rumpf [113] suggested the following equation:

$$\begin{array}{c}{V}_{lb}=2\pi \left(\left({\rho }_1^2+{\left({\rho }_1+{\rho }_2\right)}^2\right){\rho }_1\cos \left(\varphi +\theta \right)-\frac{1}{3}{\rho }_1^3{\cos }^3\left(\varphi +\theta \right)\right.\\ -{\rho }_1^2\left({\rho }_1+{\rho }_2\right)\cos \left(\varphi +\theta \right)\sin \left(\varphi +\theta \right)\left(\frac{\pi }{2}-\varphi -\theta \right)\\ \left.-\frac{1}{24}\left(2+\cos \varphi \right){\left(1-\cos \varphi \right)}^2\right)\end{array}$$

(40)

Kuwagi et al. [114] assumed a spherical interface and fixed the contact angle to be zero and gave the following formula:

$$V_{lb}=2\pi \left(\left(C^2+r_0^2\right)a-C\left(a\sqrt{r_0^2-a^2}+\alpha r_0^2\right)-\frac{a^3-b^2\left(3-b\right)}{3}\right)$$

(41)

$$\alpha =\frac{\pi }{2}-\varphi$$

(42)

$$C=\left(R_i+\frac{l}{2}\right)\tan \varphi$$

(43)

$$r_0=\frac{2R_i+l}{2\cos \varphi }-R_i$$

(44)

$$a=R_i\left(1-\cos \varphi \right)+\frac{l}{2}$$

(45)

$$b=R_i\left(1-\cos \varphi \right)$$

(46)

Some approximate formulas were also proposed. Weigert and Ripperger [115] gave the following equation and claimed the error was less than 4%:

$$V_{lb}=0.96R^3\left(1+3l^{*}\right)\left(1+1.1\sin \theta \right){\sin }^4\varphi$$

(47)

Rabinovich et al. [116] proposed the following equation:

$$V_{lb}=\pi R^2{\varphi }^{2}l+\frac{1}{2}\pi R^3{\varphi }^4$$

(48)

Liquid transfer is another important problem during liquid bridge formation and rupture. Concerning this topic, the key issue is how fast and how much the liquid can transfer between the particles and the liquid bridge. Shi and McCarthy [105] first systematically investigated the liquid transfer problem. Assuming the liquid bridge could form instantaneously, and ignoring the evaporation of the liquid, they found that the liquid redistribution ratio depends on the particle size and the contact angle. For the liquid bridge formation between two particles of equal size, the liquid would redistribute evenly to the two particles after rupture. Wu et al. [117] employed direct numerical simulation (DNS) to study the formation rate of the liquid bridge and proposed that a portion of the bridge could form instantaneously followed by a filling process at a finite rate, which depends on surface tension and the viscosity of the liquid.

3.3. Gas-Solid Interaction Forces

In a gas–solid fluidization system, the particle phase would interact with the gas phase if there was a relative velocity between them. The interactive forces may involve drag, lift force, pressure gradient force, viscous force, buoyancy, virtual mass force, Basset force, and so on. However, for a gas–solid fluidized bed, drag force usually acts as the driving force and predominates the motion of particles, and hence most studies investigating the gas–solid fluidization process did not consider other interactive forces except drag force. This generally makes sense. Therefore, we focus on the discussion of drag force models in this section.

For a single spherical particle, the calculation of the drag coefficient has been well-established for many years. However, for a multi-particle system, the calculation of the drag coefficient has become much more complicated due to the effect of other particles. Experimental methods relating to the bed pressure drop or the bed expansion height have been frequently used to correlate the drag coefficient to local voidage, particle Reynolds number, and the properties of both gas and particles. Derived from such complex gas–solid flow systems, the drag models correlated by experiments are sometimes confined to a certain situation. More recently, numerical methods such as DNS or LBM have been reported more and more frequently to derive the drag formulas from a micro scale. Such micro-scale simulation was often conducted in a rather simple system with many assumptions due to the limitation of the computation resource. This makes the generalization of such models questionable. From the perspective of practical implementation of CFD-DEM, gas–solid drag models are all supposed to undergo validation processes in specific environments, whether they were obtained from experimental or numerical studies. As a reference, some representative gas–solid drag models are listed in

Table 1. Drag models for spherical particles.

Drag Models	Equations
Gidaspow [118] (a combination of Ergun [119] and Wen-Yu [120] models)	${\overrightarrow{F}}_D=\beta \left({\overrightarrow{v}}_f-{\overrightarrow{v}}_p\right)/{\rho }_f, \beta =\left\{\begin{array}{c}150\frac{{\left(1-{\epsilon }_f\right)}^2}{{\epsilon }_f}\frac{{\mu }_f}{{\left({\phi }_p{d}_p\right)}^2}+1.75\left(1-{\epsilon }_f\right)\frac{{\rho }_f}{{\phi }_p{d}_p}\left|{\overrightarrow{v}}_f-{\overrightarrow{v}}_p\right|, {\epsilon }_f\le 0.8\\ \frac{3}{4}{C}_D\frac{{\rho }_f\left(1-{\epsilon }_f\right)}{d_p}{\epsilon }_f^{-2.7}\left|{\overrightarrow{v}}_f-{\overrightarrow{v}}_p\right|, {\epsilon }_f>0.8\end{array}\right.$
Di Felice [121]	${\overrightarrow{F}}_D=f\left({\epsilon }_f\right){\overrightarrow{F}}_{D0},f\left({\epsilon }_f\right)={\epsilon }_f^{-\left(\alpha +1\right)},\alpha =3.7-0.65\exp \left(\begin{array}{c}-{\left(1.5-{\mathrm{logRe}}_p\right)}^2/2\end{array}\right), {\overrightarrow{F}}_{D0}=\frac{1}{2}{\rho }_f{C}_D{A}_{\perp }\left({\overrightarrow{v}}_f-{\overrightarrow{v}}_p\right)\left|{\overrightarrow{v}}_f-{\overrightarrow{v}}_p\right|$
Koch and Hill [122]	${\overrightarrow{F}}_D=A\frac{{\left(1-{\epsilon }_f\right)}^2}{{\epsilon }_f}+B\left(1-{\epsilon }_f\right){\mathrm{Re}}_p,A=\left\{\begin{array}{c}180, {\epsilon }_f\le 0.6\\ \left(\frac{18{\epsilon }_f^3}{1-{\epsilon }_f}\right)\frac{1+\frac{3}{\sqrt{2}}\sqrt{1-{\epsilon }_f}+\frac{135}{64}\left(1-{\epsilon }_f\right)\ln \left(1-{\epsilon }_f\right)+16.14\left(1-{\epsilon }_f\right)}{1+0.681\left(1-{\epsilon }_f\right)-8.48{\left(1-{\epsilon }_f\right)}^2+8.16{\left(1-{\epsilon }_f\right)}^3},{\epsilon }_f>0.6\end{array}\right.$$B=\left(0.6057+1.908\left(1-{\epsilon }_f\right)\right){\epsilon }_f^3+0.209{\epsilon }_f^{-2},{\mathrm{Re}}_p=\frac{{\epsilon }_f{\rho }_f\left|{\overrightarrow{v}}_f-{\overrightarrow{v}}_p\right|d_p}{{\mu }_f}$
Beetstra [123]	$\beta =\frac{18{\mu }_f{\epsilon }_f{\epsilon }_s}{d_p^2}f\left({\epsilon }_s,{\mathrm{Re}}_p\right),f\left({\epsilon }_s,{\mathrm{Re}}_p\right)=\frac{10{\epsilon }_s}{{\epsilon }_f^2}+{\epsilon }_f^2\left(1+1.5{\epsilon }_s^{0.5}\right)+\frac{0.413{\mathrm{Re}}_p}{24{\epsilon }_f^2}\left(\frac{{\epsilon }_f^{-1}+3{\epsilon }_f{\epsilon }_s+8.4{\mathrm{Re}}_p^{-0.343}}{1+{10}^{3{\epsilon }_s}{\mathrm{Re}}_p^{-\left(1+2{\epsilon }_s\right)}}\right),{\mathrm{Re}}_p=\frac{{\epsilon }_f{\rho }_f\left|{\overrightarrow{v}}_f-{\overrightarrow{v}}_p\right|d_p}{{\mu }_f}$
Rong [124]	${\overrightarrow{F}}_D=\frac{C_D}{24}{\mathrm{Re}}_p{\epsilon }_f^{-\chi },C_D={\left(0.63+\frac{4.8}{\sqrt{{\mathrm{Re}}_p}}\right)}^2,\chi =2.65\left(1+{\epsilon }_f\right)-\left(5.3-3.5{\epsilon }_f\right){\epsilon }_f^{-2}\exp \left(-\frac{{\left(1.5-{\mathrm{logRe}}_p\right)}^2}{2}\right),{\mathrm{Re}}_p=\frac{{\epsilon }_f{\rho }_f\left|{\overrightarrow{v}}_f-{\overrightarrow{v}}_p\right|d_p}{{\mu }_f}$
Tang [125]	$\beta =\frac{18{\mu }_f{\epsilon }_f{\epsilon }_s}{d_p^2}f\left({\epsilon }_s,{\mathrm{Re}}_p\right),f\left({\epsilon }_s,{\mathrm{Re}}_p\right)=\frac{10{\epsilon }_s}{{\epsilon }_f^2}+{\epsilon }_f^2\left(1+1.5{\epsilon }_s^{0.5}\right)+\left(0.11{\epsilon }_s\left(1+{\epsilon }_s\right)-\frac{0.00456}{{\epsilon }_f^4}+\left(0.169{\epsilon }_f+\frac{0.0644}{{\epsilon }_f^4}\right){\mathrm{Re}}_p^{-0.343}\right){\mathrm{Re}}_p$${\mathrm{Re}}_p=\frac{{\epsilon }_f{\rho }_f\left|{\overrightarrow{v}}_f-{\overrightarrow{v}}_p\right|d_p}{{\mu }_f}$

.

4. Coupling Schemes

In a gas–solid fluidization system, there might be mass, energy, and momentum exchanges between the two phases, but it is unnecessary to always include all of them in practical implementation. For example, we may only pay special attention to the momentum exchange in a pure hydrodynamics simulation. We herein briefly introduce the coupling scheme of momentum as an example, and the scheme can be similarly extended to the coupling of mass and energy.

According to Norouzi et al. [101], the coupling scheme can be categorized into four types based on the concentration of the solid phase. In a dilute gas–solid system, e.g., a solid volume fraction <10⁻⁶, the motion of particles is largely controlled by the gas flow, whereas the effect of particle motion on the gas flow can be neglected. Hence, it is sensible to apply a one-way coupling scheme that only considers the effect of gas flow on the motion of the particles. As the concentration of particles increases to a certain value, e.g., a solid volume fraction >10⁻⁶, the mutual influence between the two phases should be implemented, which is called two-way coupling. If the effect of the nearby disturbance wave on the motion of particles was also considered, then the scheme was termed three-way coupling. However, we have not noticed the frequent implementation of such coupling in the simulation of gas–solid fluidization. Most often, a so-called four-way coupling scheme was used where particle–particle/particle–wall collisions were incorporated in addition to the interphase effect above. This usually occurs in a dense particle system where the solid volume fraction is typically larger than 10⁻². In fact, there is not any agreement on when to adopt a certain coupling scheme, and this was sometimes determined by the specific research goals.

The concrete strategy of momentum exchange varies by software. We herein introduce the coupling scheme of the open-source software CFDEM [126] as an example. As shown in Figure 4, the CFD solver served as the master and would iterate to convergence at the first coupling time specified by the users. Then the local voidage and slip velocity (relative velocity between the two phases) were estimated and passed to the coupling module for the calculation of the drag, which was subsequently passed to the DEM solver. The DEM solver was then able to update the motion status of the particles at the same time, such as the position and velocity. The new motion status of the particles was then transferred back to the coupling module to update the local voidage and the drag, which would then be used by the CFD solver to iterate the gas status to the next coupling time. More discussion on the interphase momentum exchange is detailed by Feng and Yu [73].

Sometimes, the terminology of ‘loose coupling’ or ‘tight coupling’ can be seen in the literature. Theoretically, a higher coupling frequency is a better approximation of reality and should thus be preferred, but this would increase the computation cost. Hence, the coupling frequency during practical implementation is somewhat more of an empirical problem, which is a trade-off between the expected approximation and affordable computation cost.

5. Research Focus and the Role of CFD-DEM

In a typical fluidized bed granulator or coater, the gas flow provided by a distributor is used to fluidize the particles. The binder or coating liquids can be introduced into the system via different locations, such as bottom spray, top spray, or tangential spray. One of the most prominent advantages of the fluidized bed granulation or coating processes is the flexible control capability of the operation conditions. For example, the inlet gas flow rate can be increased to intensify particle collisions for granulation or be decreased to mitigate particle collisions for coating. This often enables granulation or coating to be implemented in the same apparatus, thus saving time and cost.

Granulation in the wet fluidized bed often refers to the agglomeration process where primary fine particles are stuck together by the binder liquid to form larger particles. This is often used to enlarge the particle diameter or make porous structures. Another kind of granulation, usually called layer granulation, is sometimes also needed. In this process, the liquid dissolved/suspended with the desired materials is deposited onto the surface of the particles, which grow larger layer by layer. Compared with the agglomeration process, this is often used to make relatively dense particle products. Quite similar to layered granulation, in the coating process, coating liquids are also sprayed onto the surface of the particles without substantial deposition of other materials. The particle size after coating usually shows an insignificant change, which can be neglected. It is often hard to distinguish the three concepts above in practice because they usually occur simultaneously. The point is to make one of them predominant, which needs a profound understanding of the underlying mechanisms.

The basic research goal of granulation or coating processes is to form complete capability of the design, scaling-up, and optimization of the processes. This goal can be achieved from studies on multiple scales. Essentially, macroscopic phenomena, microscopic mechanisms, and the relationship between them need to be systematically investigated.

Early studies usually started from the hydrodynamics of granulation or coating processes. This is rational because it would help create a basic command of the operation principles of the apparatuses. In 2011, L. Fries et al. [53] used CFD-DEM to investigate the hydrodynamics of a fluidized bed spray granulator (see Figure 5). In their study, two different configurations of the granulator were studied and found that the Wurster coater could maintain a narrower residence time distribution of the particles in the spray zone compared with the top-spray granulator. This indicated more homogeneous particle wetting could be achieved in the Wurster coater. The effect of operation parameters and geometric details was also evaluated. For example, jet velocity above 100 m/s coupled with a gap distance below 15 times the particle diameter could lead to stable particle circulation. This work is one of the earliest studies that investigated the hydrodynamics of fluidized bed granulators through CFD-DEM. Later, L. Fries et al. [127] studied the collision dynamics in a fluidized bed granulator via the same method. Different from the prior work, this study focused on the particle–particle interactions, such as the average particle velocity, angular velocity, and collision frequency, through which the agglomeration quality can be evaluated. Three different geometric configurations were compared, and the pros and cons were assessed for a specific application. The top-spray granulator was found to be suitable for large-scale granulation processes while the wetting intensity and growth rate are relatively low. The spouted bed showed the most intensive gas–solid–liquid contact, and compact agglomerates were produced. The Wurster coater achieved the fastest growth rate and the most homogeneous particle wetting. The two works above, to some extent, set up a framework for hydrodynamics study in fluidized bed granulators or coaters by CFD-DEM.

The research above focused on mono-dispersive particle systems, while binary- or poly-dispersive particle systems were also investigated by some researchers. Li et al. [23] studied the residence time and cycle time distributions of particles in a Wurster fluidized bed. The effect of particle size was assessed, and large particles were found to spend longer in the spray zone and tube zone in each cycle and could move closer to the spray nozzle on average than small particles. This might induce an inter-particle non-uniformity in the coating process. However, the non-uniformity could be partly compensated by the fact that large particles traveled through the spray zone less frequently. Jiang et al. [128] conducted a similar study to investigate the residence time, droplet deposition, and collision velocity of a binary particle mixture in a Wurster fluidized bed coater. The droplets were inserted into the system as solid-like particles. The CFD-DEM simulation data were found to be in good agreement with the positron emission particle tracking (PEPT) results. Nonetheless, there was a distraction between the simulated and experimental results for the appearance probability of nonideal particle circulation cycles. Other representative studies that investigate the hydrodynamics of particle mixtures include Jiang et al. [48], Feng et al. [129], and Olaofe et al. [130].

Special treatments were also considered in the investigation of the hydrodynamics of particles. Limtrakul et al. [131] employed a vibrated fluidized bed to enhance the fluidization and mixing of cohesive fine powders. Using CFD-DEM simulation, they found that good mixing could be obtained in such an apparatus, which was also confirmed by the experimental data. Furthermore, the simulation results suggested that high vibration amplitude and frequency coupled with high superficial gas velocity could lead to better fluidization and mixing quality. Namdarkedenji et al. [132] studied the effect of inlet flow pulsation on fluidization characteristics through CFD-DEM. In detail, three types of pulsation, i.e., sinusoidal, rectangular, and relocating, were applied, and some key indices were paid special attention to, such as minimum fluidization velocity, pressure drop, and bubble formation. The simulation results showed that the pulsation amplitude imposed a more significant impact on the minimum fluidization velocity than the pulsation frequency. The inlet pulsation increased the bed expansion and average particle velocity whilst undermining the pressure drop and the interaction forces between particles. In addition, high vibration frequencies were found to induce more regular bubbles. Essentially, special treatments were usually applied to enhance the fluidization of certain types of particles, which were hard to fluidize normally, such as cohesive or fine powders, magnetic particles, or charged particles. Similar research includes Wang et al. [133], Hao et al. [134], and Pei et al. [135].

Heat and mass transfer in granulation and coating processes has also attracted much attention in the past decade. Initial hydrodynamics research often focused on the fluidization of dry particles due to its simplicity. However, to conduct a more realistic study, the injection of binder or coating liquids coupled with heat and mass transfer is necessary. As an extension of Van Buijtenen et al. [136], Sutkar et al. [137] successfully simulated the granulation process inside a spouted bed with the consideration of heat and mass transfer and liquid injection. In their work, the inclusion of heat and mass transfer was verified by experimental results from a combination of particle image velocimetry (PIV) and infrared thermography (IRT), and good agreement was achieved. In addition to flow pattern, pressure drop, and particle velocity, distributions of the moisture and gas temperature were also analyzed. More recently, Che et al. [45] focused on the gas–solid heat and mass transfer in a Wurster coater using CFD-DEM. In this study, much attention was paid to interphase heat transfer, coating liquid spraying, and evaporation. The temperature and humidity distributions of the gas phase were simulated and then verified by the experimental measurements. Three different drying modes were observed in the annular region. Likewise, Madlmeir et al. [24] investigated two important problems in the Wurster coating process, i.e., the evaporation of multicomponent spray liquids and the quantification of losses due to spray drying. Specifically, the effects of the inlet air flow rate, inlet air temperature, and spray flow rate on the coating yield were assessed. The simulation results, validated by the temperature and coating yield experiments, showed that appropriate operation parameters could reduce the spray-drying losses from 28.2% to 6.1%. Furthermore, optimization of the input parameters was able to maximize the coating yield and reduce the coating time by 75%. Modeling works that integrated heat and mass transfer displayed a huge potential for practical application. More works investigating the heat and mass transfer in granulation or coating processes can be found in [138,139,140].

In recent years, the coupling of CFD-DEM with the population balance model (PBM) has been gaining attraction. The coupling of CFD-DEM with PBM arose from the need to reconcile the discrepancy between the scales of particle dynamics and its growth or agglomeration. In 2014, Sen et al. [141] applied the hybrid CFD-DEM-PBM to study a fluidized bed granulation process where CFD-DEM served to calculate the particle dynamics and PBM was used to quantify the change of particle size. This model monitored the evolution of certain key process variables, such as the average particle diameter, particle size distribution, and particle liquid content, which was then qualitatively verified by the experimental results. Such a hybrid model showed the advantage of clarifying the mechanisms of the granulation process on multiple scales. Heinrich et al. [142] performed a multi-scale analysis of the coating process inside a Wurster coater using CFD-DEM-PBM. In their study, CFD-DEM played a major role and was used to determine the dynamics of particles. The coating processes were analyzed on four time and length scales. Specifically, PBM served to describe the particle growth on the macro scale. It was found that the gap distance had a significant effect on the particle dynamics while the atomization flow rates in the studied range did not. Currently, CFD-DEM usually plays a major role in such hybrid models, and PBM often serves as a supplementary tool to monitor particle growth. This type of one-way coupling stratagem can be used to consider the influence of particle dynamics on particle growth, and the reverse effect was usually neglected. Hence, the two-way coupling between CFD-DEM and PBM should be expected in the future. More information relating to this topic can be found in [143,144,145].

Artificial neural networks (ANNs) have also been used to assist granulation or coating simulation. This was displayed in two aspects. For one thing, CFD-DEM was only suitable for small time-scale simulation and ANNs were used to perform time extrapolation based on the datasets from CFD-DEM. Such studies include but are not limited to [21,145,146,147]. Furthermore, ANNs were used to predict the interactions between gas and solid phases, such as drag models. Traditional drag models usually correlate the gas–solid drag to several key parameters such as the particle Reynolds number and local void fraction. Much information relating to the interactions might be lost in such simplifications. ANNs were capable of preserving more complexity of gas–solid flows. For example, Yan et al. [148] trained a neural network based on a radial basis function and then employed it to predict the drag coefficient of non-spherical particles; Jiang et al. [149] used a neural network to filter the drag corrections on a sub-grid scale in their TFM simulations.

The use of CFD-DEM in fluidized bed granulation or coating processes mainly serves hydrodynamics research, which would contribute to the design, scaling-up, and optimization of the apparatuses or processes. To date, the fluidization of mono-sized particles or binary-/poly-sized particles has both been studied. Generally, the CFD-DEM simulation technique itself can be regarded as mature. However, for some special particles, such as fine powders or cohesive particles, simulating fluidization of them requires special treatments. Vibrated or pulsed fluidized beds, or further assistance by external magnetic/electric fields or microwaves, are worthy of consideration for such hard-fluidized particles. Heat and mass transfer are becoming increasingly popular in fluidized bed granulation and coating research. Binder or coating liquid injection and evaporation and intra-/interphase heat transfer make the granulation and coating processes even more complicated. To make the simulation results more realistically applicable, sub-models for heat and mass transfer are suggested to be incorporated in CFD-DEM simulation. As a powerful tool for hydrodynamics simulation, it is usually hard for CFD-DEM to directly simulate particle growth, which is, however, a rather important process. The coupling with PBM can help to overcome the time scale discrepancy between dynamics and growth/agglomeration for particles. Although the hybrid model CFD-DEM-PBM has been applied by some researchers, some problems remain to be resolved, such as the two-way information exchange problem, hence the coupling of CFD-DEM with PBM is expected to attract more attention in the future. ANNs have been used to assist the simulation of granulation and coating processes primarily from two aspects: Time extrapolation of the datasets and prediction of interphase interactions. Some representative applications of CFD-DEM in simulating wet fluidized bed granulators and coaters can be found in

Table 2. Representative modeling studies of wet fluidized bed granulators and coaters using CFD-DEM.

Ref.	Publication Year	Variables/Problems Considered	Highlights
Li, H., Liu, D., Ma, J., et al. [40]	2022	Heat and mass transfer; Liquid injection rate; Circulation pattern; Particle cohesion	A cohesive contact model was integrated into CFD-DEM to model the hydrodynamics, heat and mass transfer of a Wurster coater.
Madlmeir, S., Forgber, T., Trogrlic, M., et al. [31]	2021	Coating mass and thickness distributions over the entire process time; Inter-particle coating variability; Airflow rate; Spray rate	CFD-DEM was used to estimate the input parameters for a novel Monte-Carlo simulation approach.
Jiang, Z., Rieck, C., Bück, A. et al. [32]	2020	Inter- and intra-particle coating uniformity; Cycle time distribution, residence time distribution, coating coverage, uniformity of porosity and layer thickness distributions; Particle agglomeration	The effect of multiple factors on the Wurster coating process was considered, including gas flow, particle motion, droplet deposition, and the drying and solidifying of droplets on particle surfaces. The Monte Carlo approach was integrated into CFD-DEM.
Farivar, F., Zhang, H., Tian, Z.F., et al. [150]	2020	Particle shape; Residence time distributions (RTD) of particles; Positions of spray droplet-particle contacts; Final product’s particles size distribution; Coefficient of variation (COV) for the coating mass received by the particles	The coating process of non-spherical particles was studied by CFD-DEM-DDM.
Che, H.Q., Liu, D., Tian, W.B., et al. [151]	2020	Gas-solid flow regimes in a Wurster coater	Various gas-solid flow regimes inside a Wurster coater were identified on different fluidization conditions.
Tamrakar, A. and Ramachandran, R. [144]	2019	The coupling of CFD-DEM and Population Balance Model (PBM)	A coupled CFD–DEM–PBM framework for simulation of fluidized bed wet granulation was developed and validated.
Kieckhefen, P., Lichtenegger, T., Pietsch, et al. [21]	2019	Particle residence time distribution in a fictitious spray zone; Particle surface coverage distribution	A Lagrangian recurrence CFD method was used to simulate spouted beds successfully.
Breuninger, P., Weis, D., Behrendt, I., et al. [152]	2019	Particle and gas dynamics; Collision dynamics	The spouting behavior of fine and cohesive powders inside a Wurster spouted bed was investigated.
Vollmari, K. and Kruggel-Emden, H. [153]	2018	Effect of operational parameters on the particle residence time	CFD-DEM was used to simulate a continuously operated dual-chamber fluidized bed.
Boyce, C.M., Ozel, A., Kolehmainen, J., et al. [154]	2017	Growth and breakup of wet agglomerates; How liquid spreads when agglomerates interact with dry fluidized particles	A new way to map agglomerate growth and breakup behavior based on Bo and Ca was identified.
Breinlinger, T., Hashibon, A., and Kraft, T. [155]	2015	The role of surface tension on the evolution of the granule morphology	The influence of surface tension on granule morphology was simulated on a micro(granule) scale.
Barrasso, D. and Ramachandran, R. [156]	2015	Particle collision frequencies and liquid distribution; Particle size distribution	Development of a numerical framework suitable for the complex sub-processes in wet granulation.
Sen, M., Barrasso, D., Singh, R. et al. [141]	2014	Evolution of important process variables (average particle diameter, particle size distribution (PSD) and particle liquid content) over time; Distributions of collision frequencies, particle velocity and particle liquid content in different sections of the fluid bed granulator (FBG)	Development of a multi-scale hybrid CFD-DEM-PBM model for a fluidized bed granulation process.
Hilton, J.E., Ying, D.Y., and Cleary, P.W. [68]	2013	Intra- and inter-particle coating qualities; Effects of varying geometry and system operating conditions	A new method based on a spherical harmonic formulation for mapping the coating coverage over each particle was used in the CFD-DEM simulation.
Fries, L., Antonyuk, S., Heinrich, S. et al. [127]	2013	Collision dynamics of particles; Homogeneity of particle wetting; Granulator configurations; Agglomeration probability, breakage and growth rate and agglomerate strength	Evaluation of the effect of granulator configurations (top-spray granulator, spouted bed, Wurster coater) on the product properties.
Dosta, M., Antonyuk, S., Heinrich, S. [143]	2013	Breakage of agglomerates; Particle dynamics	A multi-scale strategy was proposed for the simulation of a fluidized granulator where the breakage of agglomerates was considered.
Li, L., Remmelgas, J., and van Wachem [23]	2011	Particle cycle and residence time distributions in different regions	The particle cycle and residence time distributions in a laboratory-scale Wurster fluidized bed coater were analyzed. The effect of particle size on coating uniformity was investigated.
Fries, L., Antonyuk, S., Heinrich, S. et al. [53]	2011	Granulator configurations; Process parameters (air flow rate); Residence time distribution; Homogeneity of wetting	A numerical model of a fluidized bed granulator coupling the gas and particle dynamics and considering particle wetting was developed.
Suzzi, D., Radl, S., and Khinast, J.G. [157]	2010	Coating uniformity; Tablets shape; Evolution of the liquid film on the surface of the tablets	The effect of droplets on the evolution of film surrounding the particles were investigated, where droplets were simulated by a Discrete Droplets Method and liquid evaporation and particle/wall interactions are considered.

.

6. Conclusions and Outlook

This work started with the introduction of the fundamental theory of CFD-DEM as a basic entry for beginners who are interested in simulating granulation or coating processes using CFD-DEM. Later, we discussed the concrete application of CFD-DEM, especially highlighting the research focus and role of CFD-DEM in resolving these issues. Using CFD-DEM to simulate wet fluidized bed granulation or coating processes has displayed a booming trend in the past two decades. This should be owed to its prominent capability of disclosing the microscale dynamics information of particles, which is difficult to obtain from experimental measurements. Moreover, the rapid development of both computer hardware and large-scale parallelization techniques has shed much light on the popularity of CFD-DEM. Nonetheless, the fact that CFD-DEM implementation is computationally intensive has not changed yet. Nowadays, the application of CFD-DEM is still confined to relatively small time and length scales if applied without the aid of super-computers. From the perspective of computation burden mitigation, the development of more efficient computing architectures or software can be expected to further broaden the application of CFD-DEM. For example, large-scale clusters, GPU-based, or CPU-GPU hybrid solvers could be quite promising.

For the application of CFD-DEM in granulation or coating, coupling with other models/methods such as PBM and ANNs should be another hot spot in the near future. Wet fluidized bed granulation or coating processes innately involve multiple scales, indicating a multi-scale framework spanning from micro to macro scales should be set up to elucidate fully the complexity of such processes. At present, coupling with PBM has been attempted in some studies as mentioned in Section 5, but most of them are not very mature. For example, only one-way information flow, mostly from DEM to PBM, was considered in many cases, and the inverse influence was neglected. Nevertheless, the integration of CFD-DEM into a multi-scale framework should benefit the granulation and coating research. ANNs have been used to perform time extrapolation of datasets and predict interphase interactions. However, the application of ANNs in assisting the simulation of wet fluidized bed granulation and coating processes is still immature. We believe ANNs would bring strong wings to CFD-DEM simulation in the future.