Effect of Drag Models on Residence Time Distributions of Particles in a Wurster Fluidized Bed: a DEM-CFD Study

Abstract

Fluidized bed coating has been used to coat pellets or tablets with functional substances for a number of purposes. In this coating process, particle wetting, drying and film formation are coupled to particle motion. It is therefore of interest to study particle motion in such fluidized beds and to use the results to develop a model for predicting the quality of the final product. In this paper, we present results from DEM-CFD simulations, i.e. discrete element method and computational fluid dynamics simulations of particle motion in a laboratory-scale Wurster fluidized bed that was also employed in positron emission particle tracking (PEPT) experiments. As the drag force is the dominant interaction between the gas flow and the particle motion in this type of fluidized bed, the effect of drag models on the particle motion is investigated. More specifically, the particle velocity and residence time distributions of particles in different regions calculated from five different drag models are presented. It is found that the Gidaspow and Tang drag models predict both particle cycle and residence times well. The HKL and Beetstra drag models somewhat overestimate the particle velocity in the Wurster tube and therefore predict a reduced number of recirculations and a significantly shorter cycle time.

1. Introduction

The Wurster coating process (Wurster D.E., 1959), which takes place in a fluidized bed, has been utilized for pellet coating for several decades. In this process, particles are circulated in the fluidized bed and are coated in a spray zone for different functional needs (Teunou E. and Poncelet D., 2002; Ström D. et al., 2005). At the bottom of the fluidized bed, a fluidization air flow is supplied through a distributor plate. The distributor plate usually consists of a number of orifices that have been distributed so as to provide a specific flow distribution. An atomization air flow and a liquid solution in the form of droplets are introduced by one or more two-fluid nozzles at the bottom of the fluidized bed.

This type of fluidized bed can be divided into different regions (Christensen F.N. and Bertelsen P., 1997; Karlsson S. et al., 2009), as illustrated in Fig. 1. In an ideal coating cycle, a particle begins to receive coating in the spray zone, is then dried in the Wurster tube and the fountain region, and moves back towards the horizontal transport region through the downbed region. After travelling through these different regions, the particle is ready to begin the next coating cycle when it moves into the spray zone again. The sequence of coating and drying that occurs as the particles travel through these different regions is usually regarded as a particle coating cycle (Li L. et al., 2015a). It has been pointed out that the cycle time distribution (CTD) and the residence time distribution (RTD) in different regions are critical factors in determining the quality of the final product such as the film thickness and the film thickness variability (Cheng X.X. and Turton R., 2000; Shelukar S. et al., 2000; Li L. et al., 2015a). For example, the particles must spend sufficient time in the regions for drying before returning to the horizontal transport region; otherwise, the risk of agglomeration increases due to collisions between wet particles.

Fig. 1

A schematic of the Wurster process: (1) spray zone, (2) Wurster tube, (3) fountain region, (4) downbed region, and (5) horizontal transport region.

The CTD and RTD of particles in similar fluidized beds have been studied using both experimental (Mann V. and Crosby E., 1975; Shelukar S. et al., 2000; San José M.J. et al., 2013; Li L. et al., 2015a) and modeling (Fries L. et al., 2011; Yang S. et al., 2014; Li L. et al., 2015b) methods. It was found that particles spend most of the cycle time in the horizontal transport region and that different residence times in this region can cause a large variability in the cycle time, and therefore a large variability in the film thickness. Due to limitations of experimental techniques, it is difficult to measure how particles behave in detail in the spray zone. For example, in the positron emission particle tracking (PEPT) measurements (Parker D.J. et al., 1993), the γ-rays tend to be absorbed by the metal in the lower region of the fluidized bed, which results in poor data quality in the spray region (Li L. et al., 2015a). However, the particle behavior in the spray zone is important for coating. Thus advanced modeling methods, such as the discrete element method (DEM), in combination with computational fluid dynamics (CFD), offer an attractive way of addressing this kind of issue that cannot be handled experimentally.

The discrete element method was initially proposed for granular assemblies by Cundall P.A. and Strack O.D. (1979), and was introduced into particle flow in fluidized beds by Tsuji Y. et al. (1993). In this approach, every particle is modeled as an individual element and the gas phase is treated as a continuum. Particle-particle and particle-wall interactions can be accounted for in a straightforward fashion and it is also simple to include particles of different sizes. Owing to this benefit, a number of reviews and applications of this method have been presented in the literature (Kafui K.D. et al., 2002; Deen N.G. et al., 2007; Van der Hoef M.A. et al., 2008; Link J.M. et al., 2008; Zhu H.P. et al., 2008; Ketterhagen W.R. et al., 2009; Turton R., 2010; Darabi P. et al., 2011; Fries L. et al., 2013).

On the other hand, since every particle is modeled, the discrete element method is not practical for systems involving a considerable number of particles. For many particle coating processes in pharmaceutical development, the systems often have a batch size ranging from a few grams to a few hundred grams, which corresponds to from tens of thousands to a few million particles. Nowadays, simulating tens of thousands particles can be done at a reasonable computational cost, while modeling systems with a few million particles is challenging, though still feasible (Jajcevic D. et al., 2013). Hence, DEM-CFD modeling offers a great opportunity to assist process and product development.

In DEM-CFD simulations, particles and gas flow are coupled through an interphase momentum transfer term via a suitable drag model. It is well known that in fluidized beds, particles interact strongly with the gas flow, and consequently the drag model implemented affects the predicted particle motion, as well as the CTD and RTD. In the assessment of drag models by Du W. et al. (2006) based on the two-fluid model (TFM), the Gidaspow drag model (Gidaspow D., 1994) gave the best agreement with experimental observations in modeling of spouted beds. With the recent increase in computational power, it has become feasible to derive the drag model by fully solving the flow between particles (Kriebitzsch S.H.L. et al., 2013). New drag models have been proposed and developed in recent years based upon simulations using the Lattice-Boltzmann method (LBM) or immersed boundary method (IBM) (Hill R.J. et al., 2001a; 2001b; Van der Hoef M.A. et al., 2005; Benyahia S. et al., 2006; Beetstra R. et al., 2007; Tang Y. et al., 2015). In a combined experimental and simulation study, Link J.M. et al. (2005) evaluated the Gidaspow and Hill drag models and found that they could reproduce several important regimes in spout-fluid beds. At this moment, however, it is still unclear how the drag model affects the residence time distributions of particles in different regions of fluidized beds.

In the present study, the objective is to investigate the effect of drag models on the residence time distributions of particles in different regions of a Wurster fluidized bed. The CTD and RTD are obtained using DEM-CFD simulations with different drag models and are compared to recent PEPT experimental data (Li L. et al., 2015a). The performance of the selected drag models for fluidized bed coating applications is discussed.

2. Mathematical model

2.1 Equations of particle motion

For every particle, the particle motion is described using Newton’s second law   

$$m_i\frac{\text{d}{\overrightarrow{v}}_i}{\text{d}t}=\beta \frac{V_i}{{\varepsilon }_{\text{s}}}({\overrightarrow{u}}_{\text{g}}-{\overrightarrow{v}}_i)+m_i\overrightarrow{g}-V_i\nabla P+{\overrightarrow{F}}_{\text{c},\text{i}}$$(1)

where m_i is the particle mass, ${\overrightarrow{v}}_i$ is the particle velocity, ${\overrightarrow{u}}_g$ is the gas velocity, $\overrightarrow{g}$ is the gravitational acceleration, V_i is the particle volume, ∇P is the gas pressure gradient, ε_s is the particle volume fraction, ${\overrightarrow{F}}_{\text{c},\text{i}}$ is the contact force during particle-particle and particle-wall collisions and β is the interphase momentum transfer coefficient. It is noted that other forces such as a cohesion force (Thornton C. and Yin K.K., 1991) can be incorporated into Eqn. (1) in a straightforward manner.

The angular momentum of each particle is described using   

$$I_{\text{i}}\frac{\text{d}{\overrightarrow{\omega }}_i}{\text{d}t}=\sum {\overrightarrow{T}}_i$$(2)

where I_i is the moment of inertia of the particle, ${\overrightarrow{\omega }}_i$ is the rotational velocity and ${\overrightarrow{T}}_i$ is the total torque acting on the particle.

2.2 Equations of gas flow

As discussed by Li L. et al. (2015b), turbulence can be neglected in this study. Hence, the continuity and momentum equations are written as follows   

$$\frac{\partial }{\partial t}\left({\varepsilon }_{\text{g}}{\rho }_{\text{g}}\right)+\nabla \cdot \left({\varepsilon }_{\text{g}}{\rho }_g{\overrightarrow{u}}_{\text{g}}\right)=0$$(3)

  

$$\begin{array}{l}\frac{\partial }{\partial t}\left({\varepsilon }_{\text{g}}{\rho }_{\text{g}}{\overrightarrow{u}}_{\text{g}}\right)+\nabla \cdot \left({\varepsilon }_{\text{g}}{\rho }_{\text{g}}{\overrightarrow{u}}_{\text{g}}{\overrightarrow{u}}_{\text{g}}\right)=-{\varepsilon }_{\text{g}}\nabla P\\ -\nabla \cdot \left({\varepsilon }_{\text{g}}{\rho }_{\text{g}}{\mu }_{\text{g}}\left(\nabla {\overrightarrow{u}}_{\text{g}}+{\left(\nabla {\overrightarrow{u}}_{\text{g}}\right)}^T\right)\right)-\sum _{i=1}^n\beta \left({\overrightarrow{u}}_{\text{g}}-{\overrightarrow{v}}_i\right)\end{array}$$(4)

where ε_g = 1 − ε_s is the gas volume fraction, ρ_g is the density of air and μ_g is the dynamic viscosity of air, respectively.

2.3 Soft sphere model

In accounting for interparticle forces, i.e. ${\overrightarrow{F}}_{\mathrm{c,i}}$ (see Eqn. (1)), there are two common models available (Crowe C.T. et al., 2012), the hard sphere model and the soft sphere model. In the current study, the latter is selected since it is more appropriate for describing the behavior in a dense fluidized bed where particles can remain in contact for a long time. In this model, particle-particle collisions are described via a spring-dashpot system (Tsuji Y., 1993), as sketched in Fig. 2. In the present work, the contact forces in the normal and tangential directions are calculated using the Hertz-Mindlin theory (Hertz H., 1882).

Fig. 2

A schematic of the soft sphere model.

2.4 Drag models

In Eqn. (1), the interphase momentum transfer coefficient can be specified using different drag models. In this study, five drag models are examined: the Gidaspow, Min, HKL, Beetstra and Tang drag models. Their respective expressions are presented in Table 1 and a brief summary is given below.

Table 1 Drag models

Drag model	Expression
Gidaspow	${\beta }_{\text{Gidaspow}}=\{\begin{array}{ll}{\beta }_{\text{Ergun}}=\left(150\frac{{{\varepsilon }_{\text{s}}}^2}{{\varepsilon }_{\text{g}}}+1.75{\varepsilon }_{\text{s}}R{e}_{\text{p}}\right)\frac{{\mu }_{\text{g}}}{{d_{\text{p}}}^2}\hfill & {\varepsilon }_{\text{s}}>0.2\hfill \\ {\beta }_{\text{Wen}-\text{Yu}}=\frac{3}{4}R{e}_{\text{p}}{C}_{\text{D}}\frac{{\mu }_{\text{g}}{\varepsilon }_{\text{s}}}{{d_{\text{p}}}^2}{{\varepsilon }_{\text{g}}}^{-2.65}\hfill & {\varepsilon }_{\text{s}}\le 0.2\hfill \end{array}$
	$R{e}_{\text{p}}={\varepsilon }_{\text{g}}{\rho }_{\text{g}}\left|{\overrightarrow{u}}_{\text{g}}-{\overrightarrow{v}}_i\right|\frac{d_{\text{p}}}{{\mu }_{\text{g}}}$
	$C_{\text{D}}=\{\begin{array}{ll}24\left(\frac{1+0.15{R{e}_{\text{p}}}^{0.687}}{R{e}_{\text{p}}}\right)\hfill & R{e}_{\text{p}}<1000\hfill \\ 0.44\hfill & R{e}_{\text{p}}\ge 1000\hfill \end{array}$
Min	min(βErgun, βWen–Yu)
HKL	${\beta }_{\text{HKL}}=\frac{18{\mu }_{\text{g}}{\varepsilon }_{\text{g}}{{\varepsilon }_{\text{g}}}^2}{{d_{\text{p}}}^2}F({\varepsilon }_{\text{s}},R{e}_{\text{p}})$
	F(εs, Rep) = F2 + 0.5F3Rep
	$F_2=\{\begin{array}{ll}\omega \left[10\frac{{\varepsilon }_{\text{s}}}{{{\varepsilon }_{\text{g}}}^3}\right]+(1-\omega )\left[\frac{1+\frac{3}{\sqrt{2}}\sqrt{{\varepsilon }_{\text{s}}}+\frac{135}{64}{\varepsilon }_{\text{s}}\text{ln}({\varepsilon }_{\text{s}})+17.89{\varepsilon }_{\text{s}}}{1+0.681{\varepsilon }_{\text{s}}-11.03{{\varepsilon }_{\text{s}}}^2+15.41{{\varepsilon }_{\text{s}}}^3}\right]\hfill & {\varepsilon }_{\text{s}}<0.4\hfill \\ 10\frac{{\varepsilon }_{\text{s}}}{{{\varepsilon }_{\text{g}}}^3}\hfill & {\varepsilon }_{\text{s}}\ge 0.4\hfill \end{array}$
	$F_3=\{\begin{array}{ll}0.9351\hspace{0.17em}\hspace{0.17em}{\varepsilon }_{\text{s}}+0.03667\hfill & {\varepsilon }_{\text{s}}<0.0953\hfill \\ 0.0673+0.212\hspace{0.17em}\hspace{0.17em}{\varepsilon }_{\text{s}}+0.0232\hspace{0.17em}{{\varepsilon }_{\text{g}}}^{-5}\hfill & {\varepsilon }_{\text{s}}\ge 0.0953\hfill \end{array}$
	$w=\text{exp}\left(\frac{-10\left(0.4-{\varepsilon }_{\text{s}}\right)}{{\varepsilon }_{\text{s}}}\right)$
Beetstra	${\beta }_{\text{Beetstra}}=\frac{18{\mu }_{\text{g}}{\varepsilon }_{\text{s}}{\varepsilon }_{\text{g}}}{{d_{\text{p}}}^2}F({\varepsilon }_{\text{s}},R{e}_{\text{p}})$
	$F({\varepsilon }_{\text{s}},R{e}_{\text{p}})=\frac{10\hspace{0.17em}{\varepsilon }_{\text{s}}}{{{\varepsilon }_{\text{g}}}^2}+{{\varepsilon }_{\text{g}}}^2\left(1+1.5\sqrt{{\varepsilon }_{\text{s}}}\right)+\frac{0.413R{e}_{\text{p}}}{24{{\varepsilon }_{\text{g}}}^2}\left[\frac{{{\varepsilon }_{\text{g}}}^{-1}+3\hspace{0.17em}\hspace{0.17em}{\varepsilon }_{\text{s}}{\varepsilon }_{\text{g}}+8.4{R{e}_{\text{p}}}^{-0.343}}{1+{10}^{3\hspace{0.17em}{\varepsilon }_{\text{s}}}{R{e}_{\text{p}}}^{-(1+4{\varepsilon }_{\text{s}}/2)}}\right]$
Tang	${\beta }_{\text{Tang}}=\frac{18{\mu }_{\text{g}}{\varepsilon }_{\text{s}}{\varepsilon }_{\text{g}}}{{d_{\text{p}}}^2}F({\varepsilon }_{\text{s}},R{e}_{\text{p}})$
	$F({\varepsilon }_{\text{s}},R{e}_{\text{p}})=\frac{10\hspace{0.17em}{\varepsilon }_{\text{s}}}{{{\varepsilon }_{\text{g}}}^2}+{{\varepsilon }_{\text{g}}}^2\left(1+1.5\sqrt{{\varepsilon }_{\text{s}}}\right)+\left[0.11\hspace{0.17em}\hspace{0.17em}{\varepsilon }_{\text{s}}(1+{\varepsilon }_{\text{s}})-\frac{0.00456}{{{\varepsilon }_{\text{g}}}^4}+(0.169\hspace{0.17em}{\varepsilon }_{\text{g}}+\frac{0.0644}{{{\varepsilon }_{\text{g}}}^4}){R{e}_{\text{p}}}^{-0.343}\right]R{e}_{\text{p}}$

As mentioned in the introduction, Ergun S. (1952) obtained an expression for the average fluid particle force based on a classic study of pressure drop through packed beds. Wen C. and Yu Y. (1966) also conducted a series of fluidization experiments and presented their correlation. Later, Gidaspow D. (1994) suggested a model based on the Ergun equation for the dense regime and the Wen-Yu correlation for the dilute regime. The modification made by Gidaspow D. was found to be the best in modeling of spouted beds (Du W. et al., 2006) and has been widely used in engineering practice (Deen N.G. et al., 2007).

In other cases, the minimum value given by the Ergun equation and the Wen-Yu correlation is adopted (the Min drag model), for example, in the work of Link J.M. et al. (2005).

Using accurate numerical data from LBM simulations, Hill R.J. et al. (2001b) derived a new drag model for moderate Reynolds-number flows in ordered and random arrays of spheres. In this model, a dimensionless drag force, F(ε_s,Re_p), is defined and solved via a force balance. That is, the total force exerted on a particle by the gas flow, ${\overrightarrow{F}}_{\mathrm{g\to s}}$, is the sum of the drag force, ${\overrightarrow{F}}_d$, and the buoyancy-type force due to the gas pressure gradient, ${\overrightarrow{F}}_b$. Van der Hoef M.A. et al. (2005) pointed out that the drag force and the total force differ by a factor of ε_g, i.e. ${\overrightarrow{F}}_d={\varepsilon }_{\text{g}}{\overrightarrow{F}}_{\mathrm{g\to s}}$, and that the dimensionless drag force can be converted to the interphase momentum transfer coefficient in Eqn. (1) via the following expression:   

$$\beta =\frac{18{\mu }_{\text{g}}{\varepsilon }_{\text{s}}{\varepsilon }_{\text{g}}}{{d_{\text{p}}}^2}F({\varepsilon }_{\text{s}},R{e}_{\text{p}})$$(5)

Later on, this drag model was extended by Benyahia S. et al. (2006) to cover a wider range of Reynolds numbers and particle volume fraction. The model proposed is complicated and in Table 1, a simplification suggested by Benyahia S. et al. (2006) is provided¹.

Beetstra R. et al. (2007) used the same method to derive a drag model where a more complex functional form was proposed instead of a linear scaling with Re to describe the inertial effects. Moreover, their extensive LBM simulations covered a wider range of parameters and were performed for Reynolds numbers up to 1000.

Quite recently, Tang Y. et al. (2015) proposed a new drag model based on simulation results obtained using fully resolved IBM simulations. In order to obtain highly accurate (essentially grid-independent) results, they resolved the flow past fixed assemblies of monodisperse spheres. It was stated that the grid resolution effects are consistently better resolved compared to the previous LBM simulations (Hill R.J. et al., 2001b; Beetstra R. et al., 2007).

For these different drag models, the interphase momentum transfer coefficient, β, is plotted in Fig. 3 as a function of the particle volume fraction, ε_s. Since the characteristic slip velocity and the particle volume fraction are to a large extent different in different regions of the present fluidized bed, β is calculated for a slip velocity of 5 m/s and a particle volume fraction up to 0.2 for the dilute region, e.g. the Wurster tube. For the dense region, e.g. the horizontal transport region, β is calculated for a slip velocity of 1 m/s and a particle volume fraction ranging from 0.2 to 0.6. It can be seen that in the dilute region, β decreases in the following order: the HKL, Beetstra, Tang, Gidaspow and Min drag models. More specifically, β for ε_s < 0.2 due to the HKL and Beetstra drag models is nearly twice that for the other three. For the moderately dense region with 0.2 ≤ ε_s < 0.4, there is only a slight difference in the interphase momentum transfer coefficient between the different drag models. For the very dense region with ε_s ≥ 0.4, the Beetstra, Tang and HKL drag models give a significantly larger interphase momentum transfer coefficient than the Gidaspow and Min drag models.

Fig. 3

Comparison of drag models, the relative velocity between gas flow and particles is (a) 5 m/s and (b) 1 m/s, respectively.

1  It must be noted that this extension was developed following the original derivation by Hill R. J. et al, where the total force exerted on the particle was referred to as the drag force. Thus another factor of the gas volume fraction is required when converting the dimensionless drag force to the interphase momentum transfer coefficient using Eqn. (5).

3. Method

3.1 Experimental

Recent PEPT experimental data presented by Li L. et al. (2015a) are used in this study to evaluate different drag models by comparing the calculated particle velocity, cycle and residence time distributions with the measured values.

In the PEPT experiments, a tracer particle labeled with Fluorine-18 is incorporated into the fluidized bed. As the tracer particle moves around in the vessel, the radioisotope decays and releases back-to-back γ-rays. These γ-rays are detected via two large position-sensitive detectors located on both sides of the fluidized bed. Each pair of γ-rays can be used to calculate the location of the tracer particle by means of three-dimensional triangulation. More details about the technique and the algorithms for data processing can be found in the literature (Parker D.J. et al., 1993; 1997; 2002; Forster R. et al., 2000).

The fluidized bed used in the PEPT experiments was based on the STREA-1^TM laboratory fluidized bed from Aeromatic-Fielder. The fluidized bed is 380 mm high, with top and bottom diameters of 250 mm and 114 mm, respectively; detailed dimensions are provided elsewhere (Li L. et al., 2015a). A nozzle with a diameter of 5 mm was employed in order to avoid the numerical difficulties in simulating sonic or even supersonic flow while retaining the feature of high jet velocity due to the atomization air flow. The nozzle was only a big circular orifice and no liquid solution was introduced. The fluidization air flow was supplied through a bowl-shaped distributor, consisting of a central base and an outer annular region. The central base of the distributor is fully open, while the outer annular region of the distributor has a number of orifices that correspond to approximately 4 % of the area of the outer annulus region. A wire mesh screen was put over the distributor to prevent particles from falling down below the distributor.

Microcrystalline cellulose (MCC) pellets, a common material in the pharmaceutical industry, were employed in the PEPT experiments. The particle size was measured using QICPIC Particle Size Analysis (Sympatec GmbH), and the particle density was determined using mercury intrusion porosimetry (Micromeritics AutoPore III 9410). The material parameters are given in Table 2; they were selected to match those of MCC. A photograph of a sample of the MCC pellets is provided in Fig. 4.

Table 2 The material parameters used in the numerical simulations (Roberts R. et al., 1994; Bharadwaj R. et al., 2010; Darelius A. et al., 2007; Li L. et al., 2015b)

Variable	Unit	Value
Air density	kg/m3	1.204
Air viscosity	Pa·s	1.837 × 10−5
Particle diameter	μm	1749
Particle density	kg/m3	1420
Particle sphericity	—	0.85–0.95
Poisson’s ratio of particles	—	0.30
Young’s modulus	Pa	1.0 × 106
Particle-particle friction coefficient	—	0.53
Particle-particle restitution coefficient	—	0.83
Poisson’s ratio of walls	—	0.33
Particle-wall friction coefficient	—	0.20
Particle-wall restitution coefficient	—	0.80
Fluid time step	s	5.0 × 10−5
Particle time step	s	5.0 × 10−6

Fig. 4

A photograph of a sample of the MCC pellets.

The shape of particles may affect the drag force acting on the particles. In this case, the so-called sphericity factor can be used to account for the effect of the shape of particles (see, e.g. (Zastawny M. et al., 2012)). In order to evaluate the drag force more accurately, correlations are required according to, e.g. Crowe C.T. (2006), using the sphericity of particles provided in Table 2. In addition, the shape of particles may affect the porosity in the dense downbed region and therefore the air flow passing through this region. However, this is beyond the scope of the current study.

3.2 Numerical

3.2.1. General

A fully coupled multi-phase flow solver, MultiFlow (www.multiflow.org), was employed in this study to solve the equations of particle motion and gas flow. The fluid mesh, which is refined near the nozzle and in the Wurster tube, contains approximately 65000 hexahedral cells. Simulations of single-phase gas flow were performed to ensure grid independence.

In order to specify the boundary conditions in the DEM-CFD simulations, the flow rates of the atomization and fluidization air flow, as well as the flow distribution at the distributor are required. The global flow rates were measured using a rotameter in the PEPT experiments (Li L. et al., 2015a). However, it was not possible to measure the flow distribution in the experiments. Instead, a single- phase CFD model including the air supplying chamber, the distributor, the wire mesh screen, and the fluidized bed was developed in Ansys Fluent to determine the flow distribution at the distributor. The results (not shown) obtained using this single-phase CFD model show that the air flow passing through the central base of the distributor is between 53 % and 64 % of the total fluidization air flow. Since the particle velocity profile is an indication of the flow distribution, DEM-CFD simulations with different flow distributions at the distributor were performed to compare the particle velocity at different heights above the distributor with the corresponding experimental data. These additional simulations were used to verify that a flow distribution for which approximately 55 % of the fluidization flow enters normal to the boundary in the central base of the distributor and 45 % enters normal to the boundary in the outer annular region (Li L. et al., 2015b) gave adequate agreement with the experimental data.

The atomization nozzle is modeled as a flow normal to the inlet boundary, while the outlet of the domain is considered to be a pressure outlet. The walls of the fluidized bed and the Wurster tube are set to be no-slip for the gas flow.

In order to obtain a reasonable compromise between data quality and computational cost, the CTD and RTD in different regions were evaluated for different simulation times. It was found that the CTD becomes fairly steady after a simulation time of 25 s. It was also found that a simulation time of 10 s was appropriate to determine the RTD in different regions, except in the horizontal transport region.

3.2.2. Simulation parameters

The numerical parameters and the operating conditions are presented in Tables 2 and 3, respectively. The partition gap is set to 15 mm and the Wurster tube length is 150 mm. The fluidization and atomization air flow rates are 73.3 m³/h and 3.50 m³/h, respectively. A base case with a batch size of 200 g was selected in this study to investigate the drag models. In addition to the base case, which uses the Gidaspow drag model (Run #1), DEM-CFD simulations with the same operating conditions are performed for the other drag models in Runs #2–5 and for larger batch sizes in Runs #6–9. A typical simulation was run in parallel using 32–40 CPU cores and took about 12 hours wall clock time per simulated second.

Table 3 Set-up of simulations

Simulation #	Drag model	Batch size (g)
1	Gidaspow	200
2	Min	200
3	HKL	200
4	Beetstra	200
5	Tang	200
6	HKL	400
7	Tang	400
8	HKL	600
9	Tang	600

3.2.3. Post-processing

From the point of view of a typical coating cycle, we define the start of a particle coating cycle as the point when a particle enters the spray zone from the horizontal transport region or the lower region in the Wurster tube. The particle coating cycle then ends when the particle reappears in the spray zone after having subsequently traveled through the Wurster tube, the fountain, downbed, and horizontal transport regions. When one particle coating cycle ends, the next one begins. The particle cycle time is thus defined to be the time it takes for a particle to complete one cycle.

In the PEPT experiments, it was difficult to detect the tracer particle in the lower region of the fluidized bed. In the experiments, it was therefore not possible to identify a spray region and the particle cycle was defined so that it begins and ends in the Wurster tube. A spray zone is therefore not defined in the simulations either, and the particle cycle is defined in the same way as in the experiments.

It is noted that in the PEPT experiments, only one tracer particle was incorporated into each run and continuously followed for 1.5 hours. In the simulations, on the other hand, different particle trajectories are followed for 25 s. As a result, the cycle and residence time distributions were calculated from different cycles for the same tracer particle in the experiments, while they are obtained from many different particles in the DEM-CFD simulations. For the pseudo-steady particle dynamics in this fluidized bed, it is thus assumed that the information obtained (such as the particle velocity or the cycle time distribution) by following one particle for a long time is equivalent to the information obtained by following all particles for a short period of time (Li L. et al., 2015b).

4. Results and discussion

In Table 4, the mean cycle time and mean residence times of particles in different regions are summarized for the batch size of 200 g; the corresponding experimental data is also provided. Figs. 5 and 6 show the cycle time distributions and mean residence times of particles in different regions, respectively, including values both measured in the PEPT experiment and calculated from the DEM-CFD simulations for different drag models. For all selected drag models, the general feature of the particle residence times in different regions is captured. As can be seen in Table 4, the fraction of the cycle time that particles spend in different regions corresponds closely to PEPT measured values for the different drag models.

Table 4 Summary of particle cycle and residence time distributions for the batch size of 200 g

		Residence time in				Fraction of cycle time in
	Cycle time (s)	the Wurster tube (s)	the fountain region (s)	the downbed region (s)	the horizontal transport region (s)	the Wurster tube (–)	the fountain region (–)	the downbed region (–)	the horizontal transport region (–)
PEPT	8.5	1.33	0.27	0.30	6.61	0.19	0.05	0.04	0.72
Gidaspow	7.4	1.27	0.27	0.25	5.61	0.19	0.05	0.04	0.73
Min	7.7	1.50	0.27	0.30	5.68	0.21	0.04	0.04	0.71
HKL	3.6	0.46	0.30	0.14	2.72	0.14	0.09	0.04	0.72
Beetstra	4.9	0.73	0.28	0.25	3.67	0.16	0.07	0.05	0.71
Tang	7.1	1.27	0.27	0.22	5.36	0.19	0.05	0.04	0.73

Fig. 5

The cycle time distribution of particle cycles shorter than 25 s for the base case.

Fig. 6

The residence times of particles in different regions for the base case (the secondary y-axis is used for the residence time in the horizontal transport region due to the higher values).

However, the cycle time and the residence time in the different regions, particularly in the Wurster tube and in the horizontal transport region, exhibit distinct differences for the different drag models. It can be seen that both the mean cycle time in Table 4 and the CTD in Fig. 5 calculated with the Gidaspow and Tang drag models agree quantitatively with the experimental data. For the Min drag model, there is a greater fraction of particles with a longer cycle time, which can be attributed to a somewhat slower acceleration in the Wurster tube and a longer residence time in the Wurster tube (see in Table 4). Interestingly, the CTDs predicted using the HKL and Beetstra drag models peak at a shorter cycle time and are narrower than the others. These two effects result in a much shorter mean cycle time for these two drag models.

In order to understand the cause of these differences, the vertical particle velocity at selected heights above the bottom of the fluidized bed, i.e. 45 mm, 75 mm and 105 mm, is plotted in Fig. 7. At 45 mm, which is the lowest location where reliable PEPT data could be obtained, the particle velocity is close to the PEPT measured velocity in the Wurster tube for all drag models studied, while the HKL drag model shows better agreement with the PEPT data in the downbed region. At 75 mm, all the drag models give very similar results. At 105 mm, only the particle velocity in the Wurster tube is shown, and it can be seen that all the drag models overpredict the particle velocity to a great extent in a central region corresponding to a radial location smaller than approximately 5 mm. This latter difference is probably due to the fact that in the experiments, there is turbulence in this central region induced by the atomization air flow, while the simulations employ a laminar model to predict the air flow.

Fig. 7

The time-averaged vertical particle velocity at 45 mm, 75 mm and 105 mm, for the base case.

Nevertheless, since this central region contains only a few particles, it is not expected to affect the general particle motion in the fluidized bed and the predicted CTD and RTD. It can also be seen that the HKL and Beetstra drag models give a distinctly higher velocity than the other drag models in this central region. Outside of this central region, the particle velocity for all of these drag models shows better agreement with the PEPT measured velocity. But a larger overprediction of the particle velocity compared to the PEPT measurement—though smaller than that in the central region—can be found for the HKL drag model than for the other drag models. This is similar to the results presented by Link J.M. et al. (2005) for a pseudo-two-dimensional spout-fluid bed, where it was concluded that the HKL drag model has the best predictive capabilities.

Fig. 8 shows the particle volume fraction examined at two heights above the bottom of the fluidized bed, corresponding to the inlet (15 mm) and outlet (165 mm) of the Wurster tube, respectively. It can be seen that the particle volume fraction is lower for the HKL and Beetstra drag models at the inlet of the tube and vice versa at the outlet of the tube. Since the Beetstra, Tang and HKL drag models give a larger interphase momentum exchange coefficient in the dense downbed region (cf. Fig. 3), there is a larger pressure drop in the downbed region for these models. As a result of the larger pressure drop in the downbed region, a certain amount of fluidization air flow is diverted from this region into the Wurster tube. This is supported by a lower air velocity at 15 mm in the downbed region for these models, as shown in Fig. 9. Since the particles are accelerated by the air flow through the Wurster tube, the particle volume fraction at the outlet of the tube decreases compared to that at the inlet for each drag model. For the HKL and Beetstra drag models, the larger interphase momentum exchange coefficient in the dilute Wurster tube (cf. Fig. 3) results in a higher particle velocity in the tube compared to other drag models. Correspondingly, a larger mass flow rate through the Wurster tube is expected and a higher particle volume fraction at the outlet of the tube is obtained for the HKL and Beetstra models than for the other drag models.

Fig. 8

The time-averaged particle volume fraction at 15 mm (the tube inlet) and 165 mm (the tube outlet) for the base case.

Fig. 9

The time-averaged air velocity in the downbed region at 15 mm (the tube inlet) for the base case.

The higher particle volume fraction at the inlet of the Wurster tube for the Tang, Gidaspow and Min drag models can be attributed to a lower air flow rate in the Wurster tube and a lower interphase momentum transfer coefficient. The resulting slower acceleration of particles in the Wurster tube leads to particles having a greater risk of recirculating in the Wurster tube and moving back towards the horizontal transport region from the lower edge of the Wurster tube. This type of behavior was also observed in the PEPT experiments (Li L. et al., 2015a) and is sketched in Fig. 10. Thus, it is of interest to examine the number of particle recirculations in the Wurster tube, as shown in Fig. 11. It is significant in Fig. 11 that for the HKL and Beetstra drag models, the fraction of ideal cycles, i.e. the fraction of particle cycles without any recirculation, is approximately 63 % and 43 %, which is much higher than the PEPT measured value of 12 %. The reduced number of particle recirculations for the HKL and Beetstra drag models can thus decrease the residence time in the Wurster tube and prevent a very long residence time in the horizontal transport region. This again emphasizes how important the detailed particle motion in the Wurster tube is to the operation of this type of fluidized bed.

Fig. 10

An illustration showing particle recirculation. The particle moves back towards the horizontal transport region from the lower edge of the Wurster tube.

Fig. 11

The fraction distribution of particle recirculations in the Wurster tube for the base case.

In order to explore a scenario where the dense region is likely to have a greater effect, simulations of larger batch sizes of 400 g and 600 g were performed for the HKL and Tang drag models. The results obtained for the mean particle cycle times are presented in Fig. 12 for different batch sizes. In agreement with the recent PEPT experiments, all selected drag models predict a shorter cycle time for a larger batch size. Notably, the Tang drag model shows the closest agreement for the increased batch size of 400 g. The results also indicate that, as the batch size increases, the deviation from the experiments is reduced for the HKL drag model, since the effect due to particle recirculations in the Wurster tube is reduced.

Fig. 12

The mean particle cycle times for different batch sizes.

5. Conclusions

In this study, DEM-CFD simulations of particle motion in a Wurster fluidized bed have been performed in order to study the effect of drag models on particle cycle and residence times. A comparison between the calculated particle cycle and residence time distributions from the simulations and recent PEPT experimental data was made. It was found that the Gidaspow, Min and Tang drag models predict the particle cycle time, the CTD and the RTDs qualitatively and quantitatively. The Tang drag model also demonstrates better agreement with experimental data for a larger batch size. The HKL and Beetstra drag models slightly overestimate the particle velocity in the Wurster tube. This overestimation of the particle velocity in the Wurster tube results in a much smaller number of recirculations, and therefore a much shorter residence time in the horizontal transport region and a much shorter cycle time. For the larger batch sizes, where the effect of particle recirculations is less pronounced and the dense region becomes more dominant, the agreement between the experimental data and the predictions of the HKL drag model improves. The present results thus suggest that the drag model must be selected carefully for the purpose of studying residence time distributions and, hence, coating processes in fluidized beds.

Acknowledgements

This work was undertaken within the EU seventh framework program, the PowTech Marie Curie Initial Training Network (Project no. EU FP7-PEOPLE-2010-ITN-264722). Dr. Andy Ingram at the University of Birmingham is gratefully acknowledged for his contribution to the PEPT experiments. The authors would also like to thank the high-performance computing (HPC) team at AstraZeneca R&D for their kind support during computations.

Nomenclature

DEM

Discrete Element Method

CFD

Computational Fluid Dynamics

TFM

Two Fluid Model

PEPT

Positron Emission Particle Tracking

CTD

Cycle Time Distribution

RTD

Residence Time Distribution

LBM

Lattice Boltzmann Method

IBM

Immersed Boundary Method

C_D

drag coefficient (−)

d_p

particle diameter (m)

${\overrightarrow{F}}_{\mathrm{c,i}}$

contact force (N)

$\overrightarrow{g}$

gravity vector (m/s²)

m_i

particle mass (kg)

I_i

moment of inertia (kg·m²)

P

pressure (N/m²)

Re_p

particle Reynolds number (−)

${\overrightarrow{T}}_i$

torque (N·m)

t

time (s)

${\overrightarrow{u}}_g$

gas velocity (m/s)

${\overrightarrow{v}}_{\text{i}}$

particle velocity (m/s)

V_i

particle volume (m³)

${\overrightarrow{\omega }}_{\text{i}}$

rotational velocity of the particle (rad/s)

ε_g

gas volume fraction (−)

ε_s

particle volume fraction (−)

β

interphase momentum transfer coefficient (kg/ (m³·s))

μ_g

dynamic viscosity (Pa·s)

ρ_g

air density (kg/m³)

Author’s short biography

Liang Li

Liang Li obtained a B.Eng. from Sichuan University and a M.Sc. in fluid mechanics from Chongqing University, China. Since Sept. 2011, he has received a Marie Curie fellowship within EU 7^th framework program and has worked as a Ph.D. candidate, under the supervision of Professor Anders Rasmuson, in chemical engineering at the Chalmers University of Technology, Sweden. His current research is focused on developing an advanced modeling tool for fluidized bed coating. In order to understand the link between process conditions and film formation at the level of a single pellet, the discrete element method (DEM), together with computational fluid dynamics (CFD), is selected.

Johan Remmelgas

Johan Remmelgas is an associate principal scientist in pharmaceutical development at AstraZeneca R&D, Mölndal, Sweden. He is currently working on developing and applying mechanistic models to predict product and process performance. He holds degrees in chemical engineering from Stanford University and the University of California, Santa Barbara, where he worked on developing and evaluating constitutive equations for polymer solutions. Before joining AstraZeneca, Johan worked on the development of heavy-duty diesel engines at Volvo Powertrain, Gothenburg, Sweden. He is the author or co-author of patents for applications ranging from diesel engines to inhalation devices.

Berend van Wachem

Berend van Wachem is a professor of fluid mechanics at the Department of Mechanical Engineering of the Imperial College London. He has been active in developing, employing and validating advanced computational and modelling methods to predict the behaviour of complex flows. The ultimate aim of his research is understanding the underlying physics at the various length and time scales involved. He has presented comprehensive reviews on computational techniques for the modelling of multiphase flows, including droplets, interfaces and particles in fluids. He is co-editor and co-author of two handbooks on multiphase flow computations. He has also written highly cited papers on gas-liquid interface modelling. He is also the author of computational codes, such as www.multiflow.org.

Christian von Corswant

Dr. Christian von Corswant studied to be a chemical engineer at the Chalmers University of Technology and wrote his thesis at the Lund University of Technology. The title of the thesis is “Lecithin-based microemulsions for pharmaceutical use”. He has worked at AstraZeneca R&D for more than 20 years in different positions within the pharmaceutical development unit and since 2008, has a specialist role as principal scientist for advanced drug delivery systems with a special focus on orally controlled release formulations. Christian has authored/co-authored about 20 publications to date.

Staffan Folestad

Dr Staffan Folestad is a senior principal scientist in AstraZeneca with global responsibilities in manufacturing sciences and pharmaceutical innovation. On a part-time basis, he was adjunct prof at Uppsala University (from 1999), and is currently affiliated to the Chalmers University of Technology. His research focus is in advanced in-situ process sensors and in-silico process modeling for powder processes. He is the author of > 85 scientific papers, 5 book chapters, and > 20 patents and is a member of the Royal Swedish Academy of Engineering Sciences, on the board of the Chalmers Center for Chemical Process Engineering, and initiator of the pan-European science conference series EuPAT. He received the “New Safe Medicines Faster Award 2007” from the European Federation of Pharmaceutical Sciences.

Mats O. Johansson

Mats Johansson is a senior scientist in the pharmaceutical development unit at Astra-Zeneca R&D, Mölndal, Sweden. Mats has a background in chemical engineering. Mats has played a major role in establishing fluidized bed technology for orally controlled release products at AstraZeneca over the past decades. He is currently leading a network on developing and further optimizing fluidized bed equipment and processes for modified release coatings for high-performance pharmaceutical products. He is the author or co-author of several patents for pharmaceutical applications.

Anders Rasmuson

Anders Rasmuson has been a professor in chemical engineering at the Chalmers University of Technology in Gothenburg since 1990. He has published approximately 150 papers in international scientific journals with nearly 2000 citations in the areas of particulate processes, multi-phase flows and transport phenomena, mixing and separation processes. He is a member of the Royal Swedish Academy of Engineering Sciences, a member of the EFCE working party of mixing, and a member of the board of the Chalmers Centre of Chemical Process Engineering. He has recently participated in the EU Marie Curie project POWTECH as an examiner/supervisor for three PhD candidates.

References

•  Beetstra  R.,  van der Hoef  M.A.,  Kuipers  J.A.M., Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres, AIChE Journal, 53 (2007) 489–501.
•  Benyahia  S.,  Syamlal  M.,  O’brien  T.J., Extension of Hill–Koch–Ladd drag correlation over all ranges of Reynolds number and solids volume fraction, Powder Technology, 162 (2006) 166–174.
•  Bharadwaj  R.,  Smith  C.,  Hancock  B.C., The coefficient of restitution of some pharmaceutical tablets/compacts, International journal of pharmaceutics, 402 (2010) 50–56.
•  Cheng  X.X.,  Turton  R., The prediction of variability occurring in fluidized bed coating equipment. I. The measurement of particle circulation rates in a bottom-spray fluidized bed coater, Pharmaceutical Development and Technology, 5 (2000) 311–322.
•  Christensen  F.N.,  Bertelsen  P., Qualitative description of the Wurster-based fluid-bed coating process, Drug Development and Industrial Pharmacy, 23 (1997) 451–463.
•  Crowe  C.T., Multiphase flow handbook, CRC, Boca Raton, FL, 2006.
•  Crowe  C.T.,  Schwarzkopf  J.D.,  Sommerfeld  M.,  Tsuji  Y., Multiphase Flows with Droplets and Particles, Second Edition, CRC Press Inc, GB, 2012.
•  Cundall  P.A.,  Strack  O.D., A discrete numerical model for granular assemblies, Geotechnique, 29 (1979) 47–65.
•  Darabi  P.,  Pougatch  K.,  Salcudean  M.,  Grecov  D., DEM investigations of fluidized beds in the presence of liquid coating, Powder Technology, 214 (2011) 365–374.
•  Darelius  A.,  Lennartsson  E.,  Rasmuson  A.,  Niklasson Björn  I.,  Folestad  S., Measurement of the velocity field and frictional properties of wet masses in a high shear mixer, Chemical Engineering Science, 62 (2007) 2366–2374.
•  Deen  N.G.,  Annaland  M.V.,  van der Hoef  M.A.,  Kuipers  J.A.M., Review of discrete particle modeling of fluidized beds, Chemical Engineering Science, 62 (2007) 28–44.
•  Du  W.,  Bao  X.,  Xu  J.,  Wei  W., Computational fluid dynamics (CFD) modeling of spouted bed: Assessment of drag coefficient correlations, Chemical Engineering Science, 61 (2006) 1401–1420.
•  Ergun  S., Fluid flow through packed columns, Chemical Engineering Progress, 48 (1952) 89–94.
•  Forster  R.,  Seville  J.,  Parker  D.,  Ding  Y., Tracking Single Particles in Process Equipment or Probing Processes Using Positrons, KONA Powder and Particle Journal, 18 (2000) 139–148.
•  Fries  L.,  Antonyuk  S.,  Heinrich  S.,  Dopfer  D.,  Palzer  S., Collision dynamics in fluidised bed granulators: A DEM-CFD study, Chemical Engineering Science, 86 (2013) 108–123.
•  Fries  L.,  Antonyuk  S.,  Heinrich  S.,  Palzer  S., DEM–CFD modeling of a fluidized bed spray granulator, Chemical Engineering Science, 66 (2011) 2340–2355.
•  Gidaspow  D., Multiphase flow and fluidization: continuum and kinetic theory description, Academic, Boston, 1994.
•  Hertz  H., Ueber die Berührung fester elastischer Körper, Journal für die reine und angewandte Mathematik, 92 (1882) 156–171.
•  Hill  R.J.,  Koch  D.L.,  Ladd  A.J.C., The first effects of fluid inertia on flows in ordered and random arrays of spheres, Journal of Fluid Mechanics, 448 (2001a) 213–241.
•  Hill  R.J.,  Koch  D.L.,  Ladd  A.J.C., Moderate-Reynolds-number flows in ordered and random arrays of spheres, Journal of Fluid Mechanics, 448 (2001b) 243–278.
•  Jajcevic  D.,  Siegmann  E.,  Radeke  C.,  Khinast  J.G., Large-scale CFD–DEM simulations of fluidized granular systems, Chemical Engineering Science, 98 (2013) 298–310.
•  Kafui  K.D.,  Thornton  C.,  Adams  M.J., Discrete particle-continuum fluid modelling of gas–solid fluidised beds, Chemical Engineering Science, 57 (2002) 2395–2410.
•  Karlsson  S.,  Rasmuson  A.,  van Wachem  B.,  Björn  I.N., CFD modeling of the Wurster bed coater, AIChE Journal, 55 (2009) 2578–2590.
•  Ketterhagen  W.R.,  Am Ende  M.T.,  Hancock  B.C., Process modeling in the pharmaceutical industry using the discrete element method, Journal of Pharmaceutical Sciences, 98 (2009) 442–470.
•  Kriebitzsch  S.H.L.,  van der Hoef  M.A.,  Kuipers  J.A.M., Drag force in discrete particle models—Continuum scale or single particle scale?, AIChE Journal, 59 (2013) 316–324.
•  Li  L.,  Rasmuson  A.,  Ingram  A.,  Johansson  M.,  Remmelgas  J.,  Von Corswant  C.,  Folestad  S., PEPT study of particle cycle and residence time distributions in a Wurster fluid bed, AIChE Journal, 61 (2015a) 756–768.
•  Li  L.,  Remmelgas  J.,  van Wachem  B.G.M.,  von Corswant  C.,  Johansson  M.,  Folestad  S.,  Rasmuson  A., Residence time distributions of different size particles in the spray zone of a Wurster fluid bed studied using DEM-CFD, Powder Technology, 280 (2015b) 124–134.
•  Link  J.M.,  Cuypers  L.A.,  Deen  N.G.,  Kuipers  J.A.M., Flow regimes in a spout–fluid bed: A combined experimental and simulation study, Chemical Engineering Science, 60 (2005) 3425–3442.
•  Link  J.M.,  Deen  N.G.,  Kuipers  J.A.M.,  Fan  X.,  Ingram  A.,  Parker  D.J.,  Wood  J.,  Seville  J.P.K., PEPT and discrete particle simulation study of spout-fluid bed regimes, AIChE Journal, 54 (2008) 1189–1202.
•  Mann  V.,  Crosby  E., Cycle time distribution measurements in spouted beds, Canadian Journal of Chemical Engineering, 53 (1975) 579–581.
•  Parker  D.J.,  Allen  D.A.,  Benton  D.M.,  Fowles  P.,  Mcneil  P.A.,  Tan  M.,  Beynon  T.D., Developments in particle tracking using the Birmingham Positron Camera, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 392 (1997) 421–426.
•  Parker  D.J.,  Broadbent  C.J.,  Fowles  P.,  Hawkesworth  M.R.,  Mcneil  P., Positron emission particle tracking—a technique for studying flow within engineering equipment, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 326 (1993) 592–607.
•  Parker  D.J.,  Forster  R.N.,  Fowles  P.,  Takhar  P.S., Positron emission particle tracking using the new Birmingham positron camera, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 477 (2002) 540–545.
•  Roberts  R.,  Rowe  R.,  York  P., The Poisson’s ratio of microcrystalline cellulose, International journal of pharmaceutics, 105 (1994) 177–180.
•  San José  M.J.,  Alvarez  S.,  Peñas  F.J.,  García  I., Cycle time in draft tube conical spouted bed dryer for sludge from paper industry, Chemical Engineering Science, 100 (2013) 413–420.
•  Shelukar  S.,  Ho  J.,  Zega  J.,  Roland  E.,  Yeh  N.,  Quiram  D.,  Nole  A.,  Katdare  A.,  Reynolds  S., Identification and characterization of factors controlling tablet coating uniformity in a Wurster coating process, Powder Technology, 110 (2000) 29–36.
•  Ström  D.,  Karlsson  S.,  Folestad  S.,  Niklasson Björn  I.,  Laurell  T.,  Nilsson  J.,  Rasmuson  A., A new device for coating single particles under controlled conditions, Chemical Engineering Science, 60 (2005) 4647–4653.
•  Tang  Y.,  Peters  E.A.J.F.,  Kuipers  J.A.M.,  Kriebitzsch  S.H.L.,  van der Hoef  M.A., A new drag correlation from fully resolved simulations of flow past monodisperse static arrays of spheres, AIChE Journal, 61 (2015) 688–698.
•  Teunou  E.,  Poncelet  D., Batch and continuous fluid bed coating—review and state of the art, Journal of Food Engineering, 53 (2002) 325–340.
•  Thornton  C.,  Yin  K.K., Impact of elastic spheres with and without adhesion, Powder Technology, 65 (1991) 153–166.
•  Tsuji  Y., Discrete particle simulation of gas-solid flows (From dilute to dense flows), KONA Powder and Particle Journal, 11 (1993) 57–68.
•  Tsuji  Y.,  Kawaguchi  T.,  Tanaka  T., Discrete particle simulation of two-dimensional fluidized bed, Powder Technology, 77 (1993) 79–87.
•  Turton  R., The application of modeling techniques to film-coating processes, Drug Development and Industrial Pharmacy, 36 (2010) 143–151.
•  Van der Hoef  M.A.,  Beetstra  R.,  Kuipers  J.A.M., Lattice-Boltzmann simulations of low-Reynolds-number flow past mono- and bidisperse arrays of spheres: results for the permeability and drag force, Journal of Fluid Mechanics, 528 (2005) 233–254.
•  Van der Hoef  M.A.,  van Sint Annaland  M.,  Deen  N.G.,  Kuipers  J.A.M., Numerical simulation of dense gas-solid fluidized beds: A multiscale modeling strategy, 40 (2008) 47–70.
•  Wen  C.,  Yu  Y., Mechanics of Fluidization, Fluidization, Chemical Engineering Progress Symposium Series 62, (1966) 100–111.
•  Wurster  D.E., Air-suspension technique of coating drug particles. A preliminary report, Journal of the American Pharmaceutical Association, 48 (1959) 451–454.
•  Yang  S.,  Luo  K.,  Fan  J.,  Cen  K., Particle-scale investigation of the solid dispersion and residence properties in a 3-D spout-fluid bed, AIChE Journal, 60 (2014) 2788–2804.
•  Zastawny  M.,  Mallouppas  G.,  Zhao  F.,  van Wachem  B., Derivation of drag and lift force and torque coefficients for non-spherical particles in flows, International Journal of Multiphase Flow, 39 (2012) 227–239.
•  Zhu  H.P.,  Zhou  Z.Y.,  Yang  R.Y.,  Yu  A.B., Discrete particle simulation of particulate systems: A review of major applications and findings, Chemical Engineering Science, 63 (2008) 5728–5770.